Global patterns of genomic and phenotypic variation in the invasive harlequin ladybird
BMC Biology volume 21, Article number: 141 (2023)
The harlequin ladybird Harmonia axyridis (Coleoptera: Coccinellidae), native to Asia, has been introduced to other major continents where it has caused serious negative impacts on local biodiversity. Though notable advances to understand its invasion success have been made during the past decade, especially with then newer molecular tools, the conclusions reached remain to be confirmed with more advanced genomic analyses and especially using more samples from larger geographical regions across the native range. Furthermore, although H. axyridis is one of the best studied invasive insect species with respect to life history traits (often comparing invasive and native populations), the traits responsible for its colonization success in non-native areas warrant more research.
Our analyses of genome-wide nuclear population structure indicated that an eastern Chinese population could be the source of all non-native populations and revealed several putatively adaptive candidate genomic loci involved in body color variation, visual perception, and hemolymph synthesis. Our estimates of evolutionary history indicate (1) asymmetric migration with varying population sizes across its native and non-native range, (2) a recent admixture between eastern Chinese and American populations in Europe, (3) signatures of a large progressive, historical bottleneck in the common ancestors of both populations and smaller effective sizes of the non-native population, and (4) the southwest origin and subsequent dispersal routes within its native range in China. In addition, we found that while two mitochondrial haplotypes-Hap1 and Hap2 were dominant in the native range, Hap1 was the only dominant haplotype in the non-native range. Our laboratory observations in both China and USA found statistical yet slight differences between Hap1 and Hap2 in some of life history traits.
Our study on H. axyridis provides new insights into its invasion processes into other major continents from its native Asian range, reconstructs a geographic range evolution across its native region China, and tentatively suggests that its invasiveness may differ between mitochondrial haplotypes.
Invasive species have caused severe, mostly negative impacts on biodiversity and ecosystems globally and led to immeasurable losses to economies and ecosystems [1, 2]. With climate change and globalization, biological invasions are continuing apace  and have become an important element of global change . Thus, management and prevention of invasive species have become major long-term challenges [5,6,7].
The field of invasion biology has long been seeking an understanding of traits that make a species invasive and thus predicable. These traits encompass a suite of interacting phenotypes, such as growth, reproduction, dispersal, and defense against predators, collectively termed the “invasion syndrome” [8,9,10]. An understanding of evolutionary processes that promote invasion success is essential to illustrate the concept of an invasion syndrome, and thus to develop sound, long-term approaches to preventing future invasions, and to manage extant ones [11, 12]. High-throughput sequencing is now widely appreciated as an important tool for monitoring, managing, and mitigating the impact of invasive species [13,14,15]. Population genomics can be used to obtain a detailed knowledge of the invasion history, including assessing source populations, routes of spread, number of independent introductions, the effects of genetic bottlenecks and admixture on the establishment success, adaptive potential, and further spread [16,17,18]. Genome-wide rapid rates of gene loss and gain via duplications in detoxification, gustatory, and chemosensory receptor gene families are also common in some invasive insects [19,20,21,22]. Novel allelic combinations due to extensive hybridization have led to increased hybrid fitness and enhanced environmental tolerance in invasive species [23,24,25]. Similarly, novel mutations, selective sweeps, and rapid fixation in new environments resulting in insecticide and parasite resistance have also been reported [26, 27]. In addition, serial bottleneck effects and reduced genomic diversity, by accelerating the pace of genetic drift and local adaptations to novel environments, may lead to successful invasion [12, 28, 29].
A powerful biological model for studying rapid evolution of phenotypic traits associated with invasions is provided by the invasive harlequin ladybird Harmonia axyridis (Coleoptera: Coccinellidae). It is widely distributed in its native eastern Asian range and generally appreciated as an effective natural enemy suppressing populations of aphids and other agricultural pests . Therefore, it was intentionally introduced as a biological control agent to North America and therefrom to Europe and other continents [9, 31, 32]. Its rapid establishment and wide geographical dispersal cause serious concerns for the threat it poses to biodiversity as a generalist predator [33, 34]. The past two decades have seen considerable research effort invested to understand its global invasion, with a special focus on the traits that promote its invasion success [35, 36]. Numerous phenotypic traits, especially life history traits, have been suggested as causal factors, including larger body size, increased physiological defenses against pathogens and predators, opportunistic predation, shorter development time and pre-oviposition period, higher reproductive potential, and higher temperature tolerance [10, 37, 38]. Analyses of mitochondrial  and 18 microsatellite nuclear loci [40,41,42] have suggested (1) complex population structure and heterogeneity between western and eastern North American populations, (2) multiple introductions and bottlenecks prior to establishment in the western and eastern USA, (3) subsequent colonization and establishment of South American and European populations from an eastern north American source, and (4) genetically admixed populations within Europe composed of genes originating from eastern North America, a European biocontrol population (which itself was founded from individuals from China), and western North America. These conclusions remain to be confirmed using genome-scale data which may provide new insights into the invasion success of H. axyridis. Furthermore, though constant efforts have been put into understanding H. axyridis life history traits that contribute to its invasion success, few studies have directly compared life history traits of genotypes within and between their native and non-native ranges .
Here, we present the most comprehensive global range-wide analysis of the nuclear and mitochondrial population genomics of H. axyridis. Specifically, we describe (1) global movement patterns using mitochondrial and genome-wide nuclear loci, (2) genome-wide patterns of selection signals to novel non-native environments, (3) genomic diversity, differentiation, and evolutionary demographic history, and (4) fitness differentials between native (eastern and western China) versus non-native (America) mitochondrial haplotypes. Our insights from this combination of genomic and ecological investigations improve our understanding of the invasion success of H. axyridis worldwide.
Haplotype identification by mtCOI sequencing
We amplified mtCOI gene fragments from 1025 H. axyridis individuals sampled from native (China) and non-native (North and South America, Europe) regions, obtaining a 620-bp fragment multiple sequence alignment without insertions, deletions, or stop codons. The sequencing yielded 22 haplotypes with at least two identical sequences (Additional file 1: Table S1).
Higher mitochondrial haplotypic diversity was shown in H. axyridis populations in the native range in China (Hd = 0.458, calculated after combining data of all the samples; π = 0.00204) and the non-native range in South America (Hd = 0.689; π = 0.00235), while lower diversity was found in North American (Hd = 0.370; π = 0.00183) and European populations (Hd = 0.147; π = 0.00048) (Additional file 1: Table S2). Moderate genetic differentiation was detected between native and non-native regions, based on among-site FST values (Additional file 1: Table S3). Haplotype network analysis indicated that mitochondrial haplotypes had a starlike pattern with Hap1 being ancestral and dominant and Hap2 subordinate (Fig. 1a). This pattern was congruent with the phylogenetic tree, in which all 22 haplotypes were clustered into two independent clades (Fig. 1b). The analysis, in combination with all the mtCOI sequences that have been deposited in the NCBI database, indicates that Hap1 is dominant, and Hap2 absent in the invaded regions of H. axyridis (Fig. 1c).
Genomic diversity and differentiation
Sequencing 2b-RAD libraries produced 1230.71 million raw reads from 159 individuals with an average of 7.74 million reads with restriction sites per individual (Fig. 2a). After quality filtering and variant calling using the H. axyridis reference genome, we identified 7824 polymorphic SNPs for the further analysis. A principal component analysis (PCA) revealed four putative super-populations corresponding to sample populations from west China (WCN), east China (ECN), Americas (AME), and Europe (EUR). The two principal components accounted for 36.69% of the total variance with the first (PC1) differentiating ECN, AME, and EUR and the second (PC2) distinguishing WCN from the others (Fig. 2b).
To assess the genetic differentiation among the putative super-populations, we measured genetic distances with Weir and Cockerham’s FST between the native (WCN and ECN) and non-native super-populations (AME and EUR) (Fig. 2c). The FST between WCN and AME was 0.085, and between WCN and EUR 0.083, both being higher than FST between ECN and AME (0.035) or between ECN and EUR (0.026); the FST between the WCN and ECN was 0.068. The genome-wide nucleotide diversity (π) in both native and non-native ranges indicated little genomic diversity within WCN (π = 5.78 × 10−5), ECN (π = 5.29 × 10−5), AME (π = 5.84 × 10−5), and EUR (π = 5.37 × 10−5) (Fig. 2c). The heterozygosity estimates within non-native ranges (AME: Ho = 0.175; EUR: Ho = 0.165) were higher than that within the native ranges (WCN: Ho = 0.143; ECN: Ho = 0.159) (Fig. 2d; Additional files 1: Table S4).
To further establish the phylogenetic relationship among all 27 geographic populations sampled across the world, we constructed a neighbor-joining phylogenetic tree, using Agrilus planipennis (Coleoptera: Buprestidae) as an outgroup. All sample populations clustered into four independent clades corresponding to WCN, ECN, AME and EUR super-populations as deduced by the PCA analysis (Fig. 2e). These super-populations were also supported by population structure analysis using ADMIXTURE (Fig. 2f). When using the minimum cross validation (CV) error value (K = 2), we found that the H. axyridis populations from China and America exhibited significant population structure, while those from Europe had shared ancestry with the Chinese and American populations. When using K of 3, WCN, ECN, AME, and EUR super-populations could be readily distinguished (Fig. 2f; Additional file 2: Fig. S1a).
Population origin and expansion within China
The PCA of the 7824 polymorphic SNPs within H. axyridis populations sampled from its native range China (ECN and WCN super-populations) indicated that the first two principal components accounted for 35.18% of the total variance, with PC1 reflecting the separation of XZLZ in Tibet from the other Chinese populations (Fig. 3a). This pattern was in line with the phylogenetic relationship and population structure estimates (at K = 2) (Fig. 3b, c; Additional file 2: Fig. S1b). Two populations (YNQJ and SCLS) in southwestern China were admixed with XZLZ as well as with other Chinese populations (Fig. 3c). At K = 3, populations in northeastern China (e.g., JLCC and LNSY populations) were identified as their own structured cluster. At K = 4, the XJWL, GSLZ, and SXYL populations in northwestern China could be readily distinguished (Fig. 3c). Moreover, the analysis of splits and migration rates among geographical populations in China showed that H. axyridis originated in southwestern China and spread from there to the other regions (Fig. 3d, e).
Patterns of contemporary and historical global migration
Reconstruction of demographic history using FSC26 indicated that the best fitted 3-population model suggested asymmetric migration under topology t1, with a most recent divergence of the Europe super-population from the Asia (~ 14 ybp) and an earlier divergence of the America from the common ancestral source population in Asia (~ 45 ybp, Fig. 4a; Additional file 1: Table S5, S6). The America super-population had a lower effective estimated population size than the Asia and Europe super-populations, suggesting a series of bottleneck events leading to reduced genomic diversity in the non-native populations. The estimation of contemporary migration under the island model using BA3-SNPs indicated that a majority of recent migrations occurred within sampled main ranges (Americas, Europe, and China) and not between them (Additional file 1: Table S7).
The estimation of evolutionary demographic history under the 4-population topology with recent hybrid origins of the Europe super-population provided strong support for the model of isolation with migration and exponential change in population sizes post divergence (Fig. 4b; Additional file 1: Table S8, S9). The 4-population model estimates were also largely in line with those under the 3-population model, in that (a) the America and Europe super-populations are estimated to be smaller than the east and west China super-populations, (b) contemporary migration rates between the four super-populations were low (congruent with the BA3-SNPs estimates), (c) the hybrid origins of the Europe super-population were estimated to be recent, of the order of 13 ybp, while the divergence of the America super-population from the east China super-population was estimated to have occurred around 24 ybp, and (d) the divergence between the east and west China super-populations was estimated at 43 ybp. Additionally, all non-native (introduced) populations were determined to have had significant population size expansions after the initial introduction bottleneck, whereas the native, west China super-population was noted to have undergone a ten-fold bottleneck from its ancestral size.
Historical variation in effective population sizes
Our PSMC analyses of historical variation in effective population sizes using high quality H. axyridis whole genomes, one being representative of its native range in China  and another its non-native range in Europe , indicated (1) signatures of a large progressive, historical bottleneck in the common ancestors of both populations, (2) commensurating in similar bottlenecked populations around 10,000 ybp, (3) smaller effective sizes of the non-native population through time (~ 50,000; 100,000 ybp) compared to the native population (~ 75,000; 100,000 ybp), and (4) more “recent” coalescence of the European H. axyridis genome compared to the Chinese H. axyridis genome (Fig. 4c).
Analyses of signatures of natural selection
Analyses of outliers using (a) genome-wide distribution of Weir & Cockerham’s FST across four global superpopulations, (b) estimates of FST outliers by Arlequin, and (c) OutFLANK analyses across all sample populations indicated that a single outlier locus (scaffold 1381, site 494,868) was significantly differentiated across all sample populations and was mapped back to the Histone-lysine-N-methyltransferase SETMAR locus in the H. axyridis genome (Additional file 2: Fig. S2a, b).
Several of the top ten outlier loci identified by our genome-wide scan of FST and Arlequin also mapped back to the neighborhood of the highly variable pannier locus (Additional file 2: Fig. S2c), with a couple of loci homologous to the carbamoyl-phosphate synthetase gene (CADII), which is involved in hemolymph production in coccinellid beetles. Other loci that were identified to be among the top 95% percentile of outliers by Arlequin included the procathepsin L-like locus on the X-chromosome, ATP-dependent RNA helicase spindle-E, and piggyBac transposable-element derived protein (Table 1). We acknowledge that outliers on the X-chromosome could potentially be false-positives, owing to us not controlling for sex of genotyped individuals.
Life history trait variation among H. axyridis haplotypes
Under controlled laboratory conditions, we determined if Hap1, the only mtDNA haplotype we detected in all invaded regions, had any phenotypic traits that might be associated with the invasion success of H. axyridis. The fecundity between Hap1 and Hap2 H. axyridis from China was not significantly different (P = 0.11) (Fig. 5a), while the Hap1 completed larval development by half a day faster (P = 0.0002) (Fig. 5b) and formed heavier pupae by 3.9 mg for male (P = 0.01) and 4.1 mg for female (P = 0.001) than the Hap2 (Fig. 5c, d). H. axyridis Hap1 ladybirds from Princeton, KY, USA, had significantly higher fecundity than that from Nanjing (P < 0.001) or Yangling (P < 0.001) in China, while the Hap1 from Lexington, KY, USA, was not different from that from either Nanjing (P = 0.79) or Yangling (P = 0.06) (Fig. 5e, f). Male and female H. axyridis Hap1 pupae from the USA were heavier (P < 0.001) than their counterparts from China (Fig. 5g, h).
Our comprehensive analyses of global patterns of H. axyridis invasion (1) established a timeline and sequence of invasion events indicating pervasive spread of a single mitochondrial haplotype (Hap1) from China; (2) characterized genome-wide signatures of selection in non-native and native populations, particularly at loci linked to the color morph pannier locus; (3) discovered in laboratory experiments made in both China and USA that H. axyridis mtDNA haplotype (Hap1), dominant in the putative original populations in the native range and extant only in the non-native range, had higher lifetime fecundity and improved performance in other life history traits as compared to another main haplotype (Hap2) derived from Hap1 in the native range; (4) delineated the geographical origin of H. axyridis and dispersal routes across mainland China in its native range; and (5) established historical variation in demographic history of establishment and spread of H. axyridis.
Harmonia axyridis is generally known to be native to continental, temperate and subtropical parts of east and central Asia [45, 46], yet its geographic history within this region remains largely unknown. Our analyses of divergence and migration rates among Chinese populations suggest that it might originate in southwestern China, thence spread eastward and then both north- and westward across China. The southwest part of China is widely acclaimed for its Hengduan Mountains region (HDM) as one of the earth’s 34 biodiversity hotspots . The HDM is part of the Tibeto-Himalayan region which includes a series of parallel, mostly north to south oriented high mountains separated by deep valleys with elevations ranging from 1000 m on valley floors to > 6000 m on mountain peaks [48, 49]. The HDM is not only a natural historical “museum” that has preserved plant diversity since the Cenozoic era but also a “cradle” where many new species were formed . Such highly diverse habitats and flora are assumed to be coevolved with a rich phytophagous insect fauna. Rich prey species concurred with diverse habitats may facilitate speciation in predatory lady beetles in this region. As supporting evidence, the genus Harmonia is most speciose in and most of its species are endemic to the southwestern part of China . We thus assume that H. axyridis emerged in the HDM region or nearby and spread from there to other regions in its current native range.
Previous studies that analyzed patterns of microsatellite variation and estimated migratory history suggest the establishment of North American populations from Eastern and Western China in 1988–1991, followed by subsequent importation of the American populations into South America and Africa, and a separate admixture event between the Eastern Chinese populations and American populations to establish the European populations . Our study using genome-wide SNP variation describes the evolutionary history of global populations, suggesting a recent divergence and putative hybrid admixed origins of the modern European populations, composed of genomes originating from eastern North American and eastern Chinese populations, with little continued contemporary gene flow between them and with significant bottlenecks in all invaded populations. It needs to be noted that a few populations from central Europe were clustered with Eastern and Western China groups. Our supposition is that the ladybird might have dispersed from the originally arrival places (France and Belgium) where the introduction of biological control populations was first made. Another supposition is that there was a gene exchange between ladybird populations from which our samples were collected and those introduced populations. Comparative analyses of historical effective population sizes using two whole genomes [43, 44] also indicate a large historical bottleneck event in the pre-establishment past of both native and non-native H. axyridis, along the timelines of the last glacial maximum (100,000 ybp–10,000 ybp). Similar declines with glacial recession have been reported in other insects, such as the scarce heath, Coenonympha hero , and the darkling beetle, Nyctelia confusa . However, we found that the heterozygosity values (nuclear) were higher in non-native than in native populations. A similar pattern was identified in a genetic analysis using microsatellite markers, which demonstrated a trend toward an increase in heterozygosity in some of invasive H. axyridis populations, closest to the North European core of invasion . Similarly, we discovered a high mtCOI haplotype diversity based on a single population collected from South America with relatively small sample size, where H. axyridis was introduced ca. 20 year ago from North America. Perhaps, the single population is insufficient to represent the invasive population in the South America and/or it is likely that there are many invasive populations with different origins in South America.
The FST outlier loci indicated by several lines of evidence (OutFLANK and genome-wide windowed Weir and Cockerham’s FST estimation) showed that the histone-lysine-N-methyltransferase (SETMAR) locus has putatively undergone selection pressure across both the native and the non-native regions. This locus, which has been implicated in dsDNA break repair and control of replication, is involved in reproductive regulation, epigenetic maternal imprinting, and polymorphism regulation in other insects, such as the rice brown planthopper Nilaparvata lugens  and the leafcutter ant Acromyrmex echiniator . Interestingly, this locus is found to be associated with range expansion and local adaptation across a variety of novel environments in the damselfly, Ischnura elegans . Another DNA methyltransferase (DMAP1) complex is also directly associated with female fecundity in H. axyridis , wherein gene silencing was strongly implicated in ovarian degeneration and reduced fecundity.
However, an important discovery revealed by our analyses of homology of the scaffolds containing these significantly differentiated variants as aligned with the annotated contigs clearly indicates genomic synteny with intronic regions linked and in the immediate vicinity of the pannier (pnr) locus . This locus has been previously reported to be closely associated with color morph variation in H. axyridis [58, 59]. This suggests that the non-native and native H. axyridis populations are undergoing strong selection pressure in the intronic regions linked to the pnr locus, with nearby allelic variants hitchhiking to fixation within the three global super-populations. Yet, recombination between pannier alleles may be reduced by a highly divergent sequence (~ 170 kb) in the cis-regulatory regions of pannier, leading to highly variable discrete color forms in natural populations . Interestingly, a previous study also suggests the exclusive presence of non-melanic morphs of H. axyridis in its invaded American range . This implies the linked selection effects in non-native populations leading to fixation of genomic loci around the pnr locus. Also, a previous study have shown that Other loci that were putative FST outliers included a locus linked to the procathepsin L-like protein, which is a digestive proteinase in Coleoptera  and the piggyBac locus which shows incredible diversity among insect orders . We would however encourage additional whole genomic analyses to test these linked selection effects further, in that FST outlier analyses such as those we utilized are potentially confounded by demography and selection [63,64,65,66].
We found the sole occurrence of the mitochondrial haplotype Hap1 across all globally non-native H. axyridis populations, while both Hap1 and hap2 were dominant in native populations from China. Our subsequent laboratory observations in both USA and China showed that Hap1 ladybirds were more fecund, with a higher pupal mass and a shorter larval duration than Hap2 ones. Studies of insects such as fruit flies [67, 68] and planthoppers  suggest that mitochondrial DNA variation can be associated with life history traits. So, we assume that Hap1 of H. axyridis may be more invasive than the Hap2, leading to its invasion success in North America by life history advantages. Yet, an alternative possibility is that Hap2 has been eliminated by strong negative selection in the non-native range due to lacking the standing genetic variation for adaptive potential in novel environments. Regardless, our genomic outlier analyses indicate the presence of selection pressure and linked selection across potential causal loci to differential fitness phenotypes, such as in color polymorphism and reproductive traits. The absence of Hap2 in the non-native range may be also potentially due to serial bottleneck effects and stochasticity in its establishment. Additionally, extra caution must be exercised about the life history difference performed by H. axyridis ladybirds between Hap1 and Hap2 due to uncontrollable differences in some factors (e.g., experimental conditions, diets) between the two observation sites.
As discussed in a recent review, research attention is urgently needed to focus on the characteristics of and interactions among H. axyridis and co-existent species of ladybirds in their native Asian habitats . Our current study has documented significant genetic and phenotypic differences among native and non-native populations of H. axyridis. Non-native populations of H. axyridis are currently spreading from Europe into western Asia . Therefore, hybridization between native and non-native populations will likely occur, but may be complex due to a number of differences between the populations . Based on our comparisons of the life history traits of native and non-native populations of H. axyridis, field research is required to examine the potential effects of the non-native population on the biodiversity of ladybirds in Asia and their role in the biological control of agricultural pests.
Our analyses of genome-wide nuclear population structure of H. axyridis samples from across its native and non-native ranges indicated that an eastern Chinese population could be the source of all non-native populations, revealed several putatively adaptive candidate genomic loci, and reconstructed a geographic range evolution within its native range in China. These findings provide new insights into its invasion processes into other major continents from its native Asian range. Our comparative observations of life history traits of the two main mitochondrial haplotypes in laboratories of both China and North America tentatively suggests that H. axyridis invasiveness may differ between haplotypes.
For mtCOI sequencing, 1025 H. axyridis adults were collected from 29 geographical localities across mainland China, from nine sites across North and South America and from four across Europe (Additional file 1: Table S10). Adults were collected using sweep-nets from June through October 2017–2019 mainly from various habitats including grassland, vegetable fields, garden shrubs, orchards, and soybean fields. All samples were stored in vials containing 95% ethanol at – 20 °C/ − 80 °C prior to DNA extraction.
Our previous analyses conclude that six ladybird individuals for each locality can yield a reliable estimate of population genomic parameters such as genetic diversity and estimates of population genetic differentiation with IIB-digest restriction site-associated DNA (2b-RAD) sequencing and library construction . So, we applied a sample size of n = 6 individuals (regardless of their morph variation and habitats) taken at random from the collection with > 6 individuals from each locality at 26 localities, while applying a sample of n = 5 individuals from a locality with < 6 ones for the other 3 localities (Additional file 1: Table S10). These locations were chosen to maximize the coverage of H. axyridis populations across its native and non-native ranges, at altitudes varying from a few hundreds to ca. 4000 m in Tibet in China.
The colonies of H. axyridis mtCOI haplotype for life history observations were established in insectaries in both China and USA. In 2018, two colonies of both Hap1 and Hap2 were established at an insectary at Nanjing Agricultural University, Weigang, Nanjing, from collections in two localities: one in Yangling District (latitude 34° 27′, longitude 108.08°), Xi’an city in the Shanxi province, western China; another one in Pukou district (latitude 32° 05′, longitude 118.67°), of Nanjing city in Jiangsu province, eastern China. In 2020, two Hap1 colonies were established at an insectary at University of Kentucky, Lexington, from collections in Lexington and Princeton (separated by 338 km), respectively, in Kentucky, USA. About 50 pairs of adults from each locality were individually transferred to Petri dishes (9 cm in diameter) for mating and oviposition. After laying eggs, these adults were reserved for subsequent mtCOI haplotype identification; a random sample of the 1st instar larvae hatched from those eggs were also taken for mtCOI haplotype identification. Ladybird larvae were subsequently grouped by haplotype and reared in cages at the temperature of 25 ± 1 °C, RH 65 ± 5%, and 16:8 light to dark photoperiod. The haplotype cohorts were maintained for three generations prior to the life history observation, with provision of the bean aphid Megoura japonica (Matsumura) in China and the pea aphid Acyrthosiphon pisum (Harris) in the USA, that were reared on the same host plant Vicia faba.
A total of 1025 unsexed adults from China and Europe were individually subjected to total genomic DNA extraction using the EasyPure Genomic DNA Kit (Transgen Biotech, China) by following the manufacturer’s protocol, while those from the non-native locations in the Americas were processed using the Dneasy Blood & Tissue Kit (Qiagen, Germany). The extracted genomic DNA from individual beetles was stored at – 20 °C until mtCOI fragment amplification and haplotype identification as well as 2b-RAD sequencing.
Mitochondrial haplotype identification
The mtCOI gene was amplified using the LCO1490 (5′-GGTCAACAAATCATAAA GATATTGG-3′) and HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) primers . PCRs were performed in 25 μL solutions containing 22 μL 1.1 × buffer (Tsingke, Biotech), 1 μL (10 μM) of each primer, and 1 μL of template DNA using the following thermocycling program: initial denaturation at 98 °C for 2 min, followed by 35 cycles of denaturation at 98 °C for 10 s, primer annealing at 50 °C for 10 s, and extension at 72 °C for 1 min, and a final extension at 72 °C for 2 min. Amplified products from samples from China, Europe, and South America were sequenced bi-directionally using an ABI 3730 DNA Analyzer at Genscript Biotech (Nanjing, China), and those from North America were sequenced using Sanger sequencing method at Functional Biosciences, Inc. (Wisconsin, USA).
All mtCOI sequences were aligned using ClustalW in MEGA 6.0 . Gaps and indels that were identified across all alignments were manually checked and 620 bp of mtCOI sequence was obtained from each sample. Assignment and measurements of mitochondrial haplotypes, including haplotype diversity (Hd), nucleotide diversity (π), and the average number of nucleotide differences (K), were performed using DnaSP 5.0 . The confirmation of haplotypes was done following the technique devised by De Barro and Ahmed . The neighbor joining phylogenetic tree using Kimura-2-Parameter algorithm with bootstrap values indicated on the branch, and the mtCOI sequence (NCBI Accession: KU918664) of the variegated ladybird Hippodamia variegata (Goeze) [previously Adonia variegata (Goeze)] as an outgroup. Pairwise differentiation (FST) between sampling locations was estimated with the Arlequin v.220.127.116.11 program  using the Tamura–Nei estimator . The haplotype network of the mtCOI genes was inferred using the median-joining algorithm in the software Network v.18.104.22.168 (Fluxus Technology Ltd., England) .
Life history observations
Observations of life history traits were made in the laboratory on two mitochondrial haplotypes, Hap1 and Hap2, of H. axyridis in both China (Hap1 and Hap2, at Nanjing Agricultural University, Nanjing, Jiangsu province) and USA (Hap1, at University of Kentucky, Lexington, KY). To control for the potential effects of confounding factors, the environments for rearing H. axyridis cohorts in the two insectaries were identical (temperature 25 ± 1 °C, RH 65 ± 5%, and photoperiod 16:8 light to dark). In addition, due to a lack of same prey species available for rearing the ladybird in different insectaries, two different but related aphid species were employed as diets: M. japonica in China and A. pisum in the USA. These two aphid species are in the same tribe of aphids (Macrosiphini), both feed on plants in the family Fabaceae, and both aphid species were reared on the same host plant species the broad bean V. faba for the experiments in China and the USA. Moreover, all the observations and insect-rearing procedures in both China and USA were conducted by the same individual (the first author H. L.). Neonate ladybird larvae were individually transferred into vials where aphids were supplied ad libitum as prey. They were reared at the temperature of 25 ± 1 °C, RH 65 ± 5%, and 16:8 light to dark photoperiod. Larval development was observed daily until pupation, and pupae were individually weighed (Mettler Toledo, model AL204-IC, to an accuracy of 0.01 mg). Emerged adults were paired in Petri dishes where aphids were provided ad libitum as food. After oviposition started, eggs were counted each day for 20 days, during which female ladybirds lay most of their eggs; thus, the total number of eggs laid during this period was used as a measure of fecundity.
A generalized linear model was fit to fecundity and larval duration (count variable), using quasi-Poisson error distribution when data were overdispersed , and a general linear model was fit to pupal body mass (continuous variable), as a function of independent variables of the haplotype (Hap1 and Hap2) and the site (China and USA). All models were visually checked for assumptions of normal distribution of residuals and homoscedasticity. After detecting a significant effect, the Tukey comparison procedure was used to make pairwise contrasts of least square means estimated in the model using the R package emmeans . Statistical significance was determined at α = 0.05 (two-tailed). All statistical data analyses were done using the R statistical software version 3.3.1 .
2b-RAD library construction, sequencing, and SNP calling
Approximately 200 ng of whole genomic DNA from 159 representative individuals were digested using Type IIB BsaXI restriction endonuclease, followed by the addition of five unique adapter sequences with T4 DNA ligases. The adapter ligated restriction products were PCR amplified, barcoded, and pooled prior to paired-end sequencing on an Illumina NovaSeq at OE Biotech (Qingdao, China) . Raw reads were filtered by quality (PHRED Q score > = 30 for at least 75% of each the length of each read, < 8% missing data or N bases, and the reads that did not contain BsaXI restriction sites), and all paired-end raw reads were merged using the Pear (v.0.9.6) program and assembled using the ustacks (v.1.3.4) program with the H. axyridis reference genome (http://v2.insect-genome.com/Organism/418)  using the SOAP aligner (-r0 –M4 –v2) [83, 84], allowing for a maximum mismatch of 2 nucleotide bases, and filtering for unique and optimal autosomal and X chromosome alignments. Picard Tools (http://broadinstitute.github.io/picard/; v1.118) was used to estimate, sort, remove PCR duplicates, and build BAM indices. Ambiguously mapped reads were removed from further analyses. Alignment files were then converted to BAM files using SAMtools v0.1.18 . Variant calling was performed for all samples using the Genome Analysis Toolkit (GATK, version 3.6–0-g89b7209) by filtering out sites based on the following criteria : (a) minor allele frequency (MAF) of < 0.05, (b) more than 2 alleles, and (c) > 20% of missing data. This resulted in a total of 7824 SNP sites across all samples. The SnpEff (v 5.0) program was used to annotate all SNPs against the annotated reference genome GFF3 file . The RAD-seq data analysis procedure are shown in additional file 2: Fig. S3.
Genetic diversity and population structure analysis
The phylogenetic tree of all H. axyridis geographical populations was constructed based on the neighbor-joining algorithm with 1000 bootstrap replicates using the Treebest (v1.9.2) program . A principal component analysis (PCA) was also conducted to evaluate the genetic structure of the populations using Plink v1.90) . Additionally, population structure and admixture proportions were estimated using ADMIXTURE v.1.30 with K = 2–15 subpopulation clusters , terminating at convergence at an objective function delta of < 0.0001, and using the Quasi-Newton acceleration with 3 secant conditions. The best fitting number of superpopulations was determined using the K value with the lowest cross-validation error. Heterozygosity (Ho), nucleotide diversity (π), genetic differentiation (Fst) across the four global superpopulations (ECN, WCN, AME, and EUR) as estimated by our ADMIXTURE analyses (at K = 3) were calculated using VCFtools (v0.1.16) in 10-kb nonoverlapping sliding windows along each chromosome . We utilized a conservative 10-kb window, based on the Coleopteran recombination rate of ~ 3 cM/Mb , based on a ~ 400 Mb genome, which indicates approximately one recombination event every 13-kb of the genome.
For the Chinese populations, we conducted the population genetic analysis including PCA, phylogenetic tree, and population structure analysis using identical SNP dataset (7824 SNPs). Notably, to further infer population-level splits and mixtures for H. axyridis populations from China, we filtered SNPs obtained by GATK with the following criteria: the max missing rate was 90%, collating a total of 4429 high-quality genotyped SNPs. Then, we used the Treemix 1.13 to investigate the admixture between populations, with migration edges ranging from 1 to 20 .
Estimation of contemporary and historical migration events
Prior to population genomic analyses, further filtering was applied using BCFtools v.1.11 to remove (a) the sites that were not in Hardy–Weinberg equilibrium (P < 10−4) for analyses of population structure, genetic diversity, and estimates of contemporary and historic migration and (b) the sites in linkage disequilibrium (R2 > 0.6 in 1000 sites) that were filtered for FST outlier analyses. All population genomics analyses were performed with the PPP pipeline  using the BCFtools v.1.11 (https://samtools.github.io/bcftools/bcftools.html).
Estimates of contemporary migration rates were estimated using the Bayesian MCMC method of the BayesAss3-SNPs (BA3-SNPs) [95, 96] by using a burn-in of 107 steps, followed by sampling 5 × 107 steps (sampling every 1000th iteration). Convergence and mixing of the MCMC runs were assessed using the Tracer v.1.7.1 program , and posterior density estimates of contemporary pairwise migration rates and their 95% confidence intervals were obtained.
Estimation of evolutionary demographic history
Demographic modeling of the coalescent evolutionary history was performed using the Fastsimcoal26 package (FSC26) . We combined individual beetles based on their sampling locations into three super-populations; this was done in congruence with the most optimal population structure estimated by ADMIXTURE, at K = 2, with the European population comprising of admixed Asian-American beetles: Asia, Europe, and America, and generated two-dimensional folded (minor allele) site frequency spectra (SFS) from the complete 2b-RAD dataset. With three super-populations, we considered four possible population tree topologies, shown in the NEWICK format: (1) t1: (((Asia, Europe), MRCA1), America), MRCA2; (2) t2: (((Asia, America), MRCA1), Europe), MRCA2; (3) t3: (((Europe, America), MRCA1), Asia), MRCA2; and (4) t4: (Asia, Europe, America), MRCA. We considered four possible models under each topology: (1) no migration between sampled populations (isolation with no migration), (2) asymmetric migration between sampled populations (isolation with migration, IM), (3) asymmetric migration between sampled populations (IM) with exponential growth of each population since divergence, and (4) no migration between sampled populations with exponential growth of each population since divergence (Additional file 2: Fig. S4). Additionally, based on the most likely topology and demographic history estimates from the three super-population FSC26 runs, we constructed four more scenarios: (1) no migration between populations, (2) asymmetric migration between populations, (3) no migration with exponential growth, and (4) asymmetric migration between populations and exponential growth of a four-population topology (Additional file 2: Fig. S5), such that the Asian group was split into ECN and WCN groups, with the European group formed as a result of hybridization between the American and ECN groups, which fits the K = 3 super-population structure estimates from the ADMIXTURE.
Subsequently, we then performed a series of analyses in the FSC26 to obtain estimates of demographic history: (1) 100 replicate runs of 200,000 coalescent simulations in each cycle, with 40 optimization cycles, discarding monomorphic sites, and sites with < 10 SNP’s in the SFS for estimating demographic parameters; (2) under each such model (4 three population topologies × 4 models × 100 runs each = 1600 total runs, 1 four-population topology × 4 models × 100 runs each = 400 total runs), computing AIC from the “best” likelihood scenario; (3) picking the “best” most likely model based on the AIC estimate, and simulating a total of 100 parametric bootstrap datasets of 200,000 non-recombining 100 bp DNA segments each; and (4) re-running estimation of demographic parameters, similar to step (1) from all 100 of these parametric bootstrap datasets to obtain confidence intervals around the estimates. All these estimates were scaled using an estimated mutation rate of 3.5 × 10−9 per site per generation for Drosophila melanogaster . Goodness of fit of the observed 2D site frequency spectra to that expected under the best fitting three-population and four-population demographic models were visually assessed to ensure congruence (Additional file 2: Fig. S6, S7).
Analyses of signatures of natural selection
We used the outlier detection F-statistics implemented in Arlequin v.22.214.171.124  by grouping individuals according to their global superpopulations (AME, EUR, WCN and ECN), and performing 20,000 simulations of 100 demes per group, with 10 groups of 100 demes to establish the neutral expectations of F-statistics under a hierarchical island model. Subsequently, FST outlier loci in our empirical distribution of 7824 SNPs were identified as loci in the top 1 percentile tail of the estimated distribution after the FDR (Benjamini and Hochberg) correction for multiple tests. Additionally, after thinning the dataset for LD using BCFtools, we performed analyses of FST outlier loci using OutFLANK . Weir and Cockerham FST estimates were computed across all sampled populations (i.e., without any grouping into global superpopulations), and uncorrected FST estimates were used to ensure that loci were not deviating from the corrected FST estimates. The OutFLANK function was then used to compute deviating loci from the expected neutral FST distribution, and outliers were estimated at a P-value cutoff of 0.05 using the FDR (Benjamini and Hochberg) correction for multiple testing. Outlier loci were mapped back to the raw reads from our 2b-RAD sequencing, and NCBI BLAST-n was used with publicly available Coleoptera: Coccinellidae genomes to obtain highly similar “hits.”
Analyses of historical effective population size variation in native versus non-native ranges
We inferred historical variation in effective population sizes of H. axyridis in its native range China versus in its non-native range Europe using two published resequenced whole genomes (Accession ID: SRR1348220, Beijing Academy of Agriculture and Forestry Sciences; Accession ID: ERR6054990, Wytham Woods, Oxfordshire, UK) [43, 44]. Briefly, we utilized the consensus reference genome to construct whole genome alignments using Bwa-mem2 v.2.0 [43, 101] and then sorted and constructed a consensus sequence using Samtools v.1.16.1 . The consensus sequence was then converted into the PSMC FASTA format using the fq2psmcfa tool, filtering for sites with a PHRED Q score greater than 30. We then utilized PSMC v.0.6.5-r67  with a generation time of 1 year, and a per locus mutation rate of 3.5 × 10−9 per site per generation, an upper limit of 5 for TMRCA, using 34 free atomic intervals . Results of historical variation in effective population sizes were then visualized and compared between the native (China) and non-native (Europe) genomes using the psmc_plot tool.
Availability of data and materials
All data generated or analyzed during this study are included in this published article, its supplementary information files, and publicly available repositories. All raw data of 2b-RAD sequencing generated during this study have been deposited at the National Center for Biotechnology Information (NCBI) under the accession numbers of PRJNA850112 . All sequences of 22 mtCOI haplotypes have been deposited at GenBank with the accession number: MT826903-MT826924. All scripts, files, and code utilized in the population genomic analyses have been deposited in the Zenodo database (https://doi.org/10.5281/zenodo.7796464) .
Principal component analysis
Tibetan population from western China
Eastern Chinese population
First principal component
Second principal component
Carbamoyl-phosphate synthetase gene
Hengduan Mountains region
- pnr :
IIB-digest restriction site-associated DNA
- Hd :
- π :
- K :
Average number of nucleotide differences
- H o :
- F ST :
Mitochondrial haplotypes 1
Mitochondrial haplotypes 2
Genome Analysis Toolkit
Site frequency spectra
Bellard C, Cassey P, Blackburn TM. Alien species as a driver of recent extinctions. Biol Lett. 2016;12:20150623.
Simberloff D, Martin JL, Genovesi P, Maris V, Wardle DA, Aronson J, et al. Impacts of biological invasions: what’s what and the way forward. Trends Ecol Evol. 2013;28(1):58–66.
Seebens H, Blackburn TM, Dyer EE, Genovesi P, Hulme PE, Jeschke JM, et al. No saturation in the accumulation of alien species worldwide. Nat Commun. 2017;8:14435.
Diagne C, Leroy B, Vaissière AC, Gozlan RE, Roiz D, Jarić I, et al. High and rising economic costs of biological invasions worldwide. Nature. 2021;592(7855):571–6.
Bock DG, Caseys C, Cousens RD, Hahn MA, Heredia SM, Hübner S, et al. What we still don’t know about invasion genetics. Mol Ecol. 2015;24(9):2277–97.
Fournier A, Penone C, Pennino MG, Courchamp F. Predicting future invaders and future invasions. P Natl Acad Sci USA. 2019;116(16):7905–10.
Sol D, Maspons J, Vall-Llosera M, Bartomeus I, García-Peña GE, Piñol J, et al. Unraveling the life history of successful invaders. Science. 2012;337(6094):580–3.
Tayeh A, Hufbauer RA, Estoup A, Ravigné V, Frachon L, Facon B. Biological invasion and biological control select for different life histories. Nat Commun. 2015;6:7268–73.
Roy HE, Brown PMJ, Adriaens T, Berkvens N, Borges I, Clusella-Trullas S, et al. The harlequin ladybird, Harmonia axyridis: global perspectives on invasion historyand ecology. Biol Invasions. 2016;18(4):997–1044.
Foucaud J, Hufbauer RA, Ravigné V, Olazcuaga L, Loiseau A, Ausset A, et al. How do invasion syndromes evolve? An experimental evolution approach using the ladybird Harmonia axyridis. Peer Community J. 2021;1: e33.
Colautti RI, Lau JA. Contemporary evolution during invasion: evidence for differentiation, natural selection, and local adaptation. Mol Ecol. 2015;24(9):1999–2017.
Estoup A, Ravigné V, Hufbauer R, Gautier M, Facon B. Is there a genetic paradox of biological invasion? Ann Rev Ecology Evol Syst. 2016;47(1):51–72.
Hamelin RC, Roe AD. Genomic biosurveillance of forest invasive alien enemies: a story written in code. Evol Appl. 2020;13(1):95–115.
North HL, McGaughran A, Jiggins CD. Insights into invasive species from whole-genome resequencing. Mol Ecol. 2021;30(23):6289–308.
Olazcuaga L, Loiseau A, Parrinello H, Paris M, Fraimout A, Guedot C, et al. A whole-genome scan for association with invasion success in the fruit fly Drosophila suzukii using contrasts of allele frequencies corrected for population structure. Mol Biol Evol. 2020;37(8):2369–85.
Sherpa S, Després L. The evolutionary dynamics of biological invasions: a multi-approach perspective. Evol Appl. 2021;14(6):1463–84.
Schrader L, Kim JW, Ence D, Zimin A, Klein A, Wyschetzki K, et al. Transposable element islands facilitate adaptation to novel environments in an invasive species. Nat Commun. 2014;5:5495.
Schoville SD, Chen YH, Andersson MN, Benoit JB, Bhandari A, Bowsher JH, et al. A model species for agricultural pest genomics: the genome of the Colorado potato beetle, Leptinotarsa decemlineata (Coleoptera: Chrysomelidae). Sci Rep. 2018;8(1):1931.
Wan FH, Yin CL, Tang R, Chen MH, Wu Q, Huang C, et al. A chromosome-level genome assembly of Cydia pomonella provides insights into chemical ecology and insecticide resistance. Nat Commun. 2019;10(1):4237.
Sethuraman A, Tovar A, Welch W, Dettmers R, Arce C, Skaggs T, et al. Genome of the parasitoid wasp Dinocampus coccinellae reveals extensive duplications, accelerated evolution, and independent origins of thelytokous parthenogeny and solitary behavior. G3. 2022;12(3):jkac001.
Pearce SL, Clarke DF, East PD, Elfekih S, Gordon KHJ, Jermiin LS, et al. Genomic innovations, transcriptional plasticity and gene loss underlying the evolution and divergence of two highly polyphagous and invasive Helicoverpa pest species. BMC Biol. 2017;15:63.
McKenna DD, Scully ED, Pauchet Y, Hoover K, Kirsch R, Geib SM, et al. Genome of the Asian longhorned beetle (Anoplophora glabripennis), a globally significant invasive species, reveals key functional and evolutionary innovations at the beetle–plant interface. Genome Biol. 2016;17(1):227.
Dlugosch KM, Parker IM. Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions. Mol Ecol. 2008;17(1):431–49.
Winklert DE, Chapin KJ, Francois O, Garmon JD, Gaut BS, Huxman TE. Multiple introductions and population structure during the rapid expansion of the invasive Sahara mustard (Brassica tournefortii). Ecol Evol. 2019;9(14):7928–41.
van Boheemen LA, Hodgins KA. Rapid repeatable phenotypic and genomic adaptation following multiple introductions. Mol Ecol. 2020;29(21):4102–17.
Guan F, Zhang JP, Shen HW, Wang XL, Padovan A, Walsh TK, et al. Whole-genome sequencing to detect mutations associated with resistance to insecticides and Bt proteins in Spodoptera frugiperda. Insect Sci. 2021;28(3):627–38.
Hill T, Unckless R. Adaptation, ancestral variation and gene flow in a “Sky Island” Drosophila species. Mol Ecol. 2020;30(1):83–9.
Wu NN, Zhang SF, Li XW, Cao YH, Liu XX, Wang QH, et al. Fall webworm genomes yield insights into rapid adaptation of invasive species. Nat Ecol Evol. 2019;3(1):105–15.
Fraimout A, Debat V, Fellous S, Hufbauer RA, Foucaud J, Pudlo P, et al. Deciphering the routes of invasion of Drosophila suzukii by means of ABC random forest. Mol Biol Evol. 2017;34(4):980–96.
Li HR, Li BP, Lövei GL, Kring TJ, Obrycki JJ. Interactions among native and non-native predatory Coccinellidae influence biological control and biodiversity. Ann Entomol Soc Am. 2021;114(2):119–36.
Koch RL. The multicolored Asian lady beetle, Harmonia axyridis: a review of its biology, uses in biological control, and non-target impacts. J Insect Sci. 2003;3(1):1–16.
Koch RL, Costamagna AC. Reaping benefits froman invasive species: role of Harmonia axyridis in natural bio-logical control of Aphis glycines in North America. Biocontrol. 2017;62(3):331–40.
Brown PMJ, Thomas CE, Lombaert E, Jeffries DL, Estoup A, Lawson HL. The global spreadof Harmonia axyridis (Coleoptera: Coccinellidae): distribu-tion, dispersal and routes of invasion. Biocontrol. 2011;56:623–41.
Evans EW, Soares AO, Yasuda H. Invasions by ladybugs, ladybirds, and other predatory beetles. Biocontrol. 2011;56(4):597–611.
Facon B, Hufbauer RA, Tayeh A, Loiseau A, Lombaert E, Vitalis R, Guillemaud T, Lundgren J, Estoup A. Inbreeding depression is purged in the invasive insect Harmonia axyridis. Curr Biol. 2011;21(5):424–7.
Rondoni GI, Borges J, Collatz E, Conti A, Costamagna F, Dumont E, et al. Exotic ladybirds for biological control of herbivorous insects—a review. Entomol Exp Appl. 2021;169(1):6–27.
Raak-van den Berg CL, Hemerik L, van der Werf W, de Jong PW, van Lenteren JC. Life history of the harlequin ladybird, Harmonia axyridis: a global meta-analysis. BioControl. 2017;62(1):283–96.
Lopatina EB, Reznik SY, Ovchinnikov AN, Ovchinnikova AA, Bezman-Moseyko OS, Gritsenko EV. Phenotypic plasticity of the thermal reaction norms for development in the multicolored Asian lady beetle, Harmonia axyridis (Pallas) (Coleoptera, Coccinellidae). Entomological Rev. 2020;100(6):727–44.
Blekhman A, Goryacheva I, Schepetov D, Zakharov I. Variability of the mitochondrial COI gene in native and invasive populations of Harmonia axyridis Pall. comparative analysis. PLoS One. 2020;15(4).
Lombaert E, Guillemaud T, Cornuet JM, Malausa T, Facon B, Estoup A. Bridgehead effect in the worldwide invasion of the bio-control harlequin ladybird. PLoS ONE. 2010;5(3): e9743.
Lombaert E, Guillemaud T, Thomas CE, Handley LJ, Lawson Handley LJ, Li J, Wang S, et al. Inferring the origin of populations introduced from a genetically structured native range by approximate Bayesian computation: case study of the invasive ladybird Harmonia axyridis. Mol Ecol. 2011;20(22):4654–70.
Lombaert E, Guillemaud T, Lundgren J, Koch R, Facon B, Grez A, et al. Complementarity of statistical treatments to reconstruct worldwide routes of invasion: the caseof the Asian ladybird Harmonia axyridis. Mol Ecol. 2014;23(24):5979–97.
Chen MY, Mei Y, Chen X, Chen X, Xiao D, He K, et al. A chromosome-level assembly of the harlequin ladybird Harmonia axyridis as a genomic resource to study beetle and invasion biology. Mol Ecol Res. 2021;21(4):1318–32.
Boyes D, Crowley LM, University of Oxford and Wytham Woods Genome Acquisition Lab, et al. The genome sequence of the harlequin ladybird, Harmonia axyridis (Pallas, 1773) [version 1; peer review: 1 approved with reservations]. Wellcome Open Res. 2021;6:300.
Dobzhansky T. Geographical variation in ladybeetles. Am Nat. 1933;67:97–126.
Kuznetsov VN. Lady beetles of the Russian Far East. Gainesville: The Sandhill Crane Press; 1997.
Sun H, Zhang J, Deng T, Boufford DE. Origins and evolution of plant diversity in the Hengduan Mountains. China Plant Divers. 2017;39(4):161–6.
Favre A, Päckert M, Pauls SU, Jähnig SC, Uhl D, Michalak I, et al. The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. Biol Rev Camb Philos Soc. 2015;90(1):236–53.
Muellner-Riehl AN, Schnitzler J, Kissling WD, Mosbrugger V, Rijsdijk KF, Seijmonsbergen AC, et al. Origins of global mountain plant biodiversity: testing the ‘mountain-geobiodiversity hypothesis.’ J Biogeogr. 2019;46(12):2826–38.
Ren S, Wang X, Pang H, Peng Z, Zeng T. Colored pictorial handbook of ladybird beetles in China. Beijing: Science Press; 2009. p. 194–200.
Sherpa S, Kebaïli C, Rioux D, Guéguen M, Renaud J, Després L. Population decline at distribution margins: assessing extinction risk in the last glacial relictual but still functional metapopulation of a European butterfly. Divers Distrib. 2022;28(2):271–90.
Zúñiga-Reinoso Á, Jerez V, Avaria-Llautureo J, Hernández CE. Consequences of the last glacial maximum on Nyctelia confusa (Coleoptera: Tenebrionidae) in Patagonia. Biol J Linn Soc. 2016;117(4):705–15.
Goryacheva II, Schepetov DM, Blekhman AV, Zakharov IA. On the genetic structure of Harmonia axyridis (Coleoptera, Coccinellidae) populations in native and invasive ranges, view from a position of molecular genetics. Russ J Genet. 2022;58(9):1118–28.
Liu X, Chen G, He J, Wan G, Shen D, Xia A, et al. Transcriptomic analysis reveals the inhibition of reproduction in rice brown planthopper, Nilaparvata lugens, after silencing the gene of MagR (IscA1). Insect Mol Biol. 2021;30(19):253–63.
Howe J, Schiøtt M, Li Q, Wang Z, Zhang G, Boomsma JJ. A novel method for using RNA-seq data to identify imprinted genes in social Hymenoptera with multiply mated queens. J Evol Biol. 2020;33(12):1770–82.
Dudaniec RY, Yong CJ, Lancaster LT, Svensson EI, Hansson B. Signatures of local adaptation along environmental gradients in a range-expanding damselfly (Ischnura elegans). Mol Ecol. 2018;27(11):2576–93.
Gegner J, Gegner T, Vogel H, Vilcinskas A. Silencing of the DNA methyltransferase 1 associated protein 1 (DMAP1) gene in the invasive ladybird Harmonia axyridis implies a role of the DNA methyltransferase 1-DMAP1 complex in female fecundity. Insect Mol Biol. 2020;29(2):148–59.
Gautier M, Yamaguchi J, Foucaud J, Loiseau A, Ausset A, Facon B, et al. The genomic basis of color pattern polymorphism in the harlequin ladybird. Curr Biol. 2018;28(20):3296–302.
Ando T, Niimi T. Development and evolution of color patterns in ladybird beetles: a case study in Harmonia axyridis. Dev Growth Differ. 2019;61(1):73–84.
Honek A, Brown PM, Martinkova Z, Skuhrovec J, Brabec M, Burgio G, et al. Factors determining variation in colour morph frequencies in invasive Harmonia axyridis populations. Biol Invasions. 2020;22(6):2049–62.
Beton D, Guzzo CR, Ribeiro AF, Farah CS, Terra WR. The 3D structure and function of digestive cathepsin L-like proteinases of Tenebrio molitor larval midgut. Insect Biochem Mol Biol. 2012;42(9):655–64.
Wang JJ, Du YZ, Wang SZ, Brown SJ, Park Y. Large diversity of the piggyBac-like elements in the genome of Tribolium castaneum. Insect Biochem Mol Biol. 2008;38(4):490–8.
Cruickshank TE, Hahn MW. Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. Mol Ecol. 2014;23(13):3133–57.
Hey J, Chung Y, Sethuraman A. On the occurrence of false positives in tests of migration under an isolation-with-migration model. Mol Ecol. 2015;24(20):5078–83.
Felsenstein J. Accuracy of coalescent likelihood estimates: do we need more sites, more sequences, or more loci? Mol Biol Evol. 2006;23(3):691–700.
Ahrens CW, Rymer PD, Stow A, Bragg J, Dillon S, Umbers KDL, et al. The search for loci under selection: trends, biases and progress. Mol Ecol. 2018;27(6):1342–56.
Ballard JW, Melvin RG, Katewa SD, Maas K. Mitochondrial DNA variation is associated with measurable differences in life–history traits and mitochondrial metabolism in Drosophila simulans. Evolution. 2007;61(7):1735–47.
Ballard JWO, Pichaud N. Mitochondrial DNA: more than an evolutionary bystander. Fun Ecol. 2014;28(1):218–31.
Sun JT, Duan XZ, Hoffmann AA, Liu Y, Garvin MR, Chen L, et al. Mitochondrial variation in small brown planthoppers linked to multiple traits and probably reflecting a complex evolutionary trajectory. Mol Ecol. 2019;28(14):3306–23.
Belyakova NA, Ovchinnikov AN, Bezman-Moseyko OS, Reznik SY. Comparative study of the phenotypic structure and photoperiodic responses of female Asian ladybirds Harmonia axyridis (Pallas) (Coleoptera, Coccinellidae) from Moscow, Belgorod, and Sochi. Entomol Rev. 2021;101(6):733–42.
Li HR, Qu WM, Obrycki JJ, Meng L, Chu D, Zhou X, et al. Optimizing sample size for population genomic study in a global invasive lady beetle, Harmonia axyridis. Insects. 2020;11(5):290.
Zakharov IA, Goryacheva I, Suvorov A. Mitochondrial DNA polymorphism in invasive and native populations of Harmonia axyridis. Eur J Environ Sci. 2011;1(1):15–8.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30(12):2725–9.
Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25(11):1451–2.
De Barro PJ, Ahmed MZ. Genetic networking of the Bemisia tabaci cryptic species complex reveals pattern of biological invasions. PLoS ONE. 2011;6(10): e25579.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Res. 2010;10(3):564–7.
Tamura K, Nei M. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol Biol Evol. 1993;10(3):512–26.
Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16(1):37–48.
Faraway JJ. Extending the linear model with R. 2nd ed. FL: CRC Press; 2016.
Lenth R, Singmann H, Love J, Buerkner P, Herve M. emmeans: estimated marginal means, aka least-squares means. R package version 1.5.0. 2020. https://CRAN.R-Project.Org/Package=emmeans. https://cran.r-project.org/web/packages=emmeans.
R Core Team. R: a language and environment for statistical computing. R foundation for statistical computing, Vienna. 2018; https://www.R-project.org.
Wang S, Meyer E, Mckay JK, Matz MV. 2b-RAD: a simple and flexible method for genome-wide genotyping. Nat Methods. 2012;9(8):808.
Li R, Li Y, Kristiansen K, Wang J. SOAP: short oligonucleotide alignment program. Bioinformatics. 2008;24(5):713–4.
Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol. 2013;22(11):3124–40.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303.
Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012;6(2):80–92.
Ruan J, Li H, Chen Z, Coghlan A, Coin LJ, Guo Y, et al. TreeFam: 2008 Update. Nucleic Acids Res. 2008;36:735–40.
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75.
Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19(9):1655–64.
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo M, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8.
Wilfert L, Gadau J, Schmid-Hempel P. Variation in genomic recombination rates among animal taxa and the case of social insects. Heredity. 2007;98(4):189–97.
Pickrell JK, Pritchard JK. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 2012;8(11): e1002967.
Webb A, Knoblauch J, Sabankar N, Kallur AS, Hey J, Sethuraman A. The pop-gen pipeline platform: a software platform for population genomic analyses. Mol Biol Evol. 2021;38(8):3478–85.
Wilson GA, Rannala B. Bayesian inference of recent migration rates using multilocus genotypes. Genetics. 2003;163(3):1177–91.
Mussmann SM, Douglas MR, Chafin TK, Douglas ME. BA3-SNPs: contemporary migration reconfigured in BayesAss for next-generation sequence data. Methods Ecol Evol. 2019;10:1808–13.
Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67(5:):901–4.
Excoffier L, Dupanloup I, Huerta-Sánchez E, Sousa VC, Foll M. Robust demographic inference from genomic and SNP data. PLoS Genet. 2013;9(10): e1003905.
Keightley PD, Trivedi U, Thomson M, Oliver F, Kumar S, Blaxter ML. Analysis of the genome sequences of three Drosophila melanogaster spontaneous mutation accumulation lines. Genome Res. 2009;19(7):1195–201.
Whitlock MC, Lotterhos KE. Reliable detection of loci responsible for local adaptation: inference of a null model through trimming the distribution of F (FST). Am Nat. 2015;186:24–36.
Vasimuddin M, Misra S, Li H, Aluru S. Efficient Architecture-Aware Acceleration of BWA-MEM for Multicore Systems. Proceedings of the 2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS). p. 314–24.
Li H, Durbin R. Inference of human population history from individual whole-genome sequences. Nature. 2011;475(7357):493–6.
Nadachowska-Brzyska K, Burri R, Smeds L, Ellegren H. PSMC analysis of effective population sizes in molecular ecology and its application to black-and-white Ficedula flycatchers. Mol Ecol. 2016;25(5):1058–72.
2b-RAD squencing of Harmonia axyridis populations. NCBI. 2023; https://www.ncbi.nlm.nih.gov/bioproject/ PRJNA951539 (Last accessed April 5, 2023).
Sethuraman A. arunsethuraman/harmoniapopgen: Final (Final). Zenodo; 2023. https://doi.org/10.5281/zenodo.7796464.
We appreciate the help with sample collection by Shaokui Wen, Hehe Cao, Haitian Song, Huanhuan Zhang, and Yuyang Shen. We also thank the reviewers and editors for constructive feedback on versions of the manuscript.
This work was funded by the National Key R&D Program of China #2017YFE0104900 to LM, the China Scholarship Council (CSC) grant to HL, NSF CAREER Award #2042516 to AS, NSF-ABI: 1564659 to AS and co-PI Jody Hey (Temple University), USDA-REEU: 2017–06423 to PI Vourlitis (CSUSM) and co-PI AS, NSF-REU: 1852189 to PI Betsy Read (CSUSM) and co-PI AS, the USDA-HSI:2022-77040-38529 to PI Sethuraman and co-PIs Vourlitis and Jancovich and the University of Kentucky Bobby C. Pass Research Professorship to JJO. YV, HV, NG, KK, and RC were supported by the CSUSM Summer Scholars program. This research includes calculations carried out on HPC resources supported in part by the NSF MRI: 1625061 and by the US Army Research Laboratory under contract number W911NF-16–2-0189. Computational analyses were also performed on the HPC cluster at San Diego State University, which was supported by startup funds to AS.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. Polymorphic nucleotide position of mtCOI gene defining the 22 haplotypes identified in native and non-native Harmonia axyridis populations. Table S2. Genetic diversity indices of native and non-native Harmonia axyridis populations based on 620-bp mtCOI fragment. Table S3. Pairwise FSTand migration rateof Harmonia axyridis populations from native and non-native regions. Table S4. Individual estimates the number of observed and expected homozygote sites, observed heterozygosity, and a method of moments estimate of the individual inbreeding coefficient. Table S5. Estimates of Maximum likelihood and Akaike’s information criterionacross different 3-population evolutionary demographic models. Table S6. Parameter estimates of genetic diversity, migration rate, divergence time, and population size change under the most likely 3-population topologyunder the IM model with population size changes across three major groups of Harmonia axyridis populations worldwide. Table S7. Estimates of contemporary migration between three major groups of Harmonia axyridis populations as estimated withBA3-SNPs. Table S8. Estimates of Maximum likelihood and Akaike’s information criterionacross different 4-population evolutionary demographic models. Table S9. Parameter estimates of genetic diversity, migration rate, divergence time, and population size change under the most likely 4-population topologyunder the IM model with population size changes across four major groups of Harmonia axyridis populations worldwide. Table S10. Collection information of Harmonia axyridis samples and their haplotypes from native and non-native ranges.
Fig. S1. Cross validation errors computed with ADMIXTURE v.1.3.0 by varying the number of superpopulations between K = 1 and 15. Fig. S2. Differentiation measured as Fstversus all chromosomal locations from OutFLANK analyses. Fig. S3. The workflow analysis for population genomics of the harlequin ladybird, Harmonia axyridis in the present study. Fig. S4. Demographic models of global invasion history of H. axyridis from an ancestral Asian source population tested in FSC26 between 3 superpopulations. Fig. S5. Demographic models of global invasion history of H. axyridis from an ancestral Eastern China source population tested in FSC26 between 4 superpopulations. Fig. S6. Goodness of fit estimates of the observed and estimated two dimensional site frequency spectra from our best fitting four-population model comprising asymmetric migration and exponential population size change. Fig. S7. Goodness of fit estimates of the observed and estimated two dimensional site frequency spectra from our best fitting three-population modelcomprising asymmetric migration and exponential population size change.
About this article
Cite this article
Li, H., Peng, Y., Wang, Y. et al. Global patterns of genomic and phenotypic variation in the invasive harlequin ladybird. BMC Biol 21, 141 (2023). https://doi.org/10.1186/s12915-023-01638-7