Skip to main content
Download PDF

Global patterns of genomic and phenotypic variation in the invasive harlequin ladybird

BMC Biology volume 21, Article number: 141 (2023) Cite this article

Abstract

Background

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.

Results

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.

Conclusions

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.

Background

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 [3] and have become an important element of global change [4]. 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 [30]. 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 [39] 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 [30].

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.

Results

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).

Fig. 1
figure 1

Harmonia axyridis mtCOI haplotype phylogeny and distribution across its native and non-native regions. a Haplotype network of H. axyridis populations across the world, inferred using median-joining algorithm and colored by geography. b Phylogenetic tree of H. axyridis mtCOI haplotypes using neighbor-joining method, with the grey circles representing the bootstrap value, the polymorphic nucleotide positions of 22 haplotypes were shown in Table S1. c Global distribution of H. axyridis mtCOI haplotypes, with the country names in red representing a total of 566 mtCOI gene sequences downloaded from NCBI database including non-native regions of Germany (42 sequences, 658 bp), Canada (331 sequences, 563 bp), USA (61 sequences, 566 bp), Italy (20 sequences, 566 bp), and the UK (34 sequences, 566 bp), as well as native regions of Russia (58 sequences, 566 bp) and Japan (20 sequences, 566 bp)

Full size image

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).

Fig. 2
figure 2

Global geographic location and genetic population structure of Harmonia axyridis sample populations divided into four major groups (super-populations) represented by colors of purple for the west China (WCN), red for the east China (ECN), blue for the Europe (EUR), and green for the America (AME). a Global distribution of all sampled populations, see Additional file 1: Table S10 for detailed information about ladybird samples collected across the world. b A biplot of the first two principle components obtained from the principal component analysis of all sampled populations clustered into four major groups. c Genetic diversity (π) and differentiation (FST) across major groups, with the radius of a circle proportional to the genetic diversity and the length of connecting dashed line to the FST value between two groups. d Comparison of individual heterozygosity rates between the major groups, see additional file 1: Table S4 for detailed information about the heterozygosity rate. e A phylogenetic tree of all sampled populations, with Agrilus planipennis as the outgroup and WCN group as the closest to the outgroup. f Map of the structure (K = 2 and 3) matching the phylogenetic tree, with different colors representing the components from the respective major groups in each sample population

Full size image

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).

Fig. 3
figure 3

Population structure and geographic expansion of the Harmonia axyridis populations within China. a A biplot of the first two components obtained from the principal component analysis of geographic sample populations clustered into two major groups (super-populations) represented by colors of purple for the west China (WCN) and red for the east China (ECN). b A phylogenetic tree of H. axyridis populations reconstructed with Agrilus planipennis as the outgroup and with the WCN group as the closest to it. c Population structure of sample populations (K = 2, 3, and 4) matching the phylogenetic tree with colored (same as above) triangles representing respective components in each population from the two major groups. d Splits and migrations among sample populations with shades of color representing the weight of migration. e A schematic diagram showing putative origin and dispersal (with arrows representing the direction) of the H. axyridis populations represented by triangles with colors of purple for the west China (WCN) and red for the east China (ECN) super-populations

Full size image

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).

Fig. 4
figure 4

Demographic history and genome-wide patterns of differentiation of major groups (super-populations) of Harmonia axyridis. a The most likely evolutionary demographic history of global invasion of H. axyridis as estimated from FSC26 analyses under the 3-population model, suggesting a recent invasion (~ 14 ybp) of the Europe super-population from an Asian source, with significant yet low continued gene flow from the America super-population. Additionally, the estimated population size of the America super-population is considerably smaller compared to the European and two Asian super-populations, the detailed information was shown in additional file 1: Table S5-S7. b The most likely evolutionary demographic history of global invasion of H. axyridis as estimated from FSC26 analyses under the 4-population model, indicating recent (~ 13 ybp) hybrid origins of the European super-population as a result of admixture between the Eastern China and the America super-populations, which themselves diverged from their common ancestor ~ 24 ybp. The Western China and the Eastern China super-populations were estimated to have diverged around 43 ybp. The America, Europe, and Eastern China super-populations have undergone significant increase in effective population sizes after recovery from bottlenecks, while the Western China has decreased in size since divergence, the detailed information was shown in additional file 2: Table S8-S9. c Historical variation in effective population sizes estimated from the whole genome resequencing data sourced from Chen et al., 2021 (native population from China), and Boyes et al., 2021 (invasive, non-native population from the UK, Europe). All estimates were scaled by a generation time of 1 generation per year, and a mutation rate of 3.5e − 9 per site per generation

Full size image

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 [43] and another its non-native range in Europe [44], 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.

Table 1 Top ten Weir and Cockerham’s Fst outlier loci identified across all sampled populations using Arlequin, and their homology to known functional genes from GenBank
Full size table

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).

Fig. 5
figure 5

Life history traits between Harmonia axyridis haplotypes in native and non-native ranges. a–d Life history trait distribution of the two dominant haplotypes Hap-1 and Hap-2 of H. axyridis in its native range. Legend: G = geographic site, H = haplotype, G × H = interaction between geographic site and haplotype. e–h Life history traits of Hap-1 H. axyridis in its native (NJHap-1 and SXHap-1, China) and invaded (KYHap-1 and KYWHap-1, USA) ranges

Full size image

Discussion

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 [47]. 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 [47]. 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 [50]. 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 [42]. 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 [51], and the darkling beetle, Nyctelia confusa [52]. 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 [53]. 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 [54] and the leafcutter ant Acromyrmex echiniator [55]. 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 [56]. Another DNA methyltransferase (DMAP1) complex is also directly associated with female fecundity in H. axyridis [57], 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 [58]. 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 [58]. Interestingly, a previous study also suggests the exclusive presence of non-melanic morphs of H. axyridis in its invaded American range [60]. 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 [61] and the piggyBac locus which shows incredible diversity among insect orders [62]. 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 [69] 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 [30]. 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 [70]. Therefore, hybridization between native and non-native populations will likely occur, but may be complex due to a number of differences between the populations [30]. 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.

Conclusions

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.

Methods

Sample collection

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 [71]. 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.

DNA extraction

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 [72]. 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 [73]. 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 [74]. The confirmation of haplotypes was done following the technique devised by De Barro and Ahmed [75]. 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.3.5.1.2 program [76] using the Tamura–Nei estimator [77]. The haplotype network of the mtCOI genes was inferred using the median-joining algorithm in the software Network v.4.6.1.0 (Fluxus Technology Ltd., England) [78].

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 [79], 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 [80]. Statistical significance was determined at α = 0.05 (two-tailed). All statistical data analyses were done using the R statistical software version 3.3.1 [81].

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) [82]. 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) [43] 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 [85]. 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 [86]: (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 [87]. 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 [88]. A principal component analysis (PCA) was also conducted to evaluate the genetic structure of the populations using Plink v1.90) [89]. Additionally, population structure and admixture proportions were estimated using ADMIXTURE v.1.30 with K = 2–15 subpopulation clusters [90], 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 [91]. We utilized a conservative 10-kb window, based on the Coleopteran recombination rate of ~ 3 cM/Mb [92], 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 [93].

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 [94] 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 [97], 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) [98]. 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 [99]. 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.3.5.2.2 [76] 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 [100]. 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 [88]. 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 [102] 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 [103]. 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 [104]. 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) [105].

Abbreviations

PCA:

Principal component analysis

WCN:

Tibetan population from western China

ECN:

Eastern Chinese population

AME:

Americas populations

EUR:

Europe populations

PC1:

First principal component

PC2:

Second principal component

CV:

Cross validation

FSC26:

fastsimcoal26

CADII:

Carbamoyl-phosphate synthetase gene

HDM:

Hengduan Mountains region

SETMAR:

Histone-lysine-N-methyltransferase

DMAP1:

DNA methyltransferase

pnr :

pannier

2b-RAD:

IIB-digest restriction site-associated DNA

Hd :

Haplotype diversity

π :

Nucleotide diversity

K :

Average number of nucleotide differences

H o :

Heterozygosity

F ST :

Pairwise differentiation

Hap1:

Mitochondrial haplotypes 1

Hap2:

Mitochondrial haplotypes 2

GATK:

Genome Analysis Toolkit

BA3-SNPs:

BayesAss3-SNPs

SFS:

Site frequency spectra

References

  1. Bellard C, Cassey P, Blackburn TM. Alien species as a driver of recent extinctions. Biol Lett. 2016;12:20150623.

    Article  PubMed  PubMed Central  Google Scholar 

  2. 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.

    Article  PubMed  Google Scholar 

  3. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. 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.

    Article  CAS  PubMed  Google Scholar 

  5. 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.

    Article  PubMed  Google Scholar 

  6. Fournier A, Penone C, Pennino MG, Courchamp F. Predicting future invaders and future invasions. P Natl Acad Sci USA. 2019;116(16):7905–10.

    Article  CAS  Google Scholar 

  7. 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.

    Article  CAS  PubMed  Google Scholar 

  8. 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.

    Article  CAS  PubMed  Google Scholar 

  9. 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.

    Article  Google Scholar 

  10. 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.

    Article  Google Scholar 

  11. Colautti RI, Lau JA. Contemporary evolution during invasion: evidence for differentiation, natural selection, and local adaptation. Mol Ecol. 2015;24(9):1999–2017.

    Article  PubMed  Google Scholar 

  12. 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.

    Article  Google Scholar 

  13. Hamelin RC, Roe AD. Genomic biosurveillance of forest invasive alien enemies: a story written in code. Evol Appl. 2020;13(1):95–115.

    Article  PubMed  Google Scholar 

  14. North HL, McGaughran A, Jiggins CD. Insights into invasive species from whole-genome resequencing. Mol Ecol. 2021;30(23):6289–308.

    Article  PubMed  Google Scholar 

  15. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Sherpa S, Després L. The evolutionary dynamics of biological invasions: a multi-approach perspective. Evol Appl. 2021;14(6):1463–84.

    Article  PubMed  PubMed Central  Google Scholar 

  17. 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.

    Article  CAS  PubMed  Google Scholar 

  18. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  19. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  20. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  23. 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.

    Article  CAS  PubMed  Google Scholar 

  24. 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.

    Article  Google Scholar 

  25. van Boheemen LA, Hodgins KA. Rapid repeatable phenotypic and genomic adaptation following multiple introductions. Mol Ecol. 2020;29(21):4102–17.

    Article  PubMed  Google Scholar 

  26. 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.

    Article  CAS  PubMed  Google Scholar 

  27. Hill T, Unckless R. Adaptation, ancestral variation and gene flow in a “Sky Island” Drosophila species. Mol Ecol. 2020;30(1):83–9.

    Article  PubMed  PubMed Central  Google Scholar 

  28. 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.

    Article  PubMed  Google Scholar 

  29. 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.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  31. 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.

    Article  Google Scholar 

  32. 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.

    Article  Google Scholar 

  33. 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.

    Article  Google Scholar 

  34. Evans EW, Soares AO, Yasuda H. Invasions by ladybugs, ladybirds, and other predatory beetles. Biocontrol. 2011;56(4):597–611.

    Article  Google Scholar 

  35. 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.

    Article  CAS  PubMed  Google Scholar 

  36. 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.

    Article  Google Scholar 

  37. 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.

    Article  Google Scholar 

  38. 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.

    Article  Google Scholar 

  39. 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).

  40. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  41. 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.

    Article  CAS  PubMed  Google Scholar 

  42. 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.

    Article  PubMed  Google Scholar 

  43. 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.

    Article  CAS  Google Scholar 

  44. 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.

    Article  Google Scholar 

  45. Dobzhansky T. Geographical variation in ladybeetles. Am Nat. 1933;67:97–126.

    Article  Google Scholar 

  46. Kuznetsov VN. Lady beetles of the Russian Far East. Gainesville: The Sandhill Crane Press; 1997.

    Google Scholar 

  47. 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.

    Article  PubMed  Google Scholar 

  48. 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.

    Article  PubMed  Google Scholar 

  49. 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.

  50. 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.

    Google Scholar 

  51. 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.

    Article  Google Scholar 

  52. 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.

    Article  Google Scholar 

  53. 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.

    Article  CAS  Google Scholar 

  54. 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.

    Article  CAS  PubMed  Google Scholar 

  55. 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.

    Article  CAS  PubMed  Google Scholar 

  56. 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.

    Article  CAS  PubMed  Google Scholar 

  57. 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.

    Article  CAS  PubMed  Google Scholar 

  58. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. 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.

    Article  PubMed  Google Scholar 

  60. 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.

    Article  Google Scholar 

  61. 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.

    Article  CAS  PubMed  Google Scholar 

  62. 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.

    Article  CAS  PubMed  Google Scholar 

  63. 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.

    Article  PubMed  Google Scholar 

  64. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  65. 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.

    Article  CAS  PubMed  Google Scholar 

  66. 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.

    Article  PubMed  Google Scholar 

  67. 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.

    Article  CAS  PubMed  Google Scholar 

  68. Ballard JWO, Pichaud N. Mitochondrial DNA: more than an evolutionary bystander. Fun Ecol. 2014;28(1):218–31.

    Article  Google Scholar 

  69. 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.

    CAS  PubMed  Google Scholar 

  70. 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.

    Article  Google Scholar 

  71. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  72. 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.

    Google Scholar 

  73. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25(11):1451–2.

    Article  CAS  PubMed  Google Scholar 

  75. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  76. 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.

    Article  Google Scholar 

  77. 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.

    CAS  PubMed  Google Scholar 

  78. Bandelt HJ, Forster P, Röhl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999;16(1):37–48.

    Article  CAS  PubMed  Google Scholar 

  79. Faraway JJ. Extending the linear model with R. 2nd ed. FL: CRC Press; 2016.

    Book  Google Scholar 

  80. 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.

  81. R Core Team. R: a language and environment for statistical computing. R foundation for statistical computing, Vienna. 2018; https://www.R-project.org.

  82. 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.

    Article  CAS  PubMed  Google Scholar 

  83. Li R, Li Y, Kristiansen K, Wang J. SOAP: short oligonucleotide alignment program. Bioinformatics. 2008;24(5):713–4.

    Article  CAS  PubMed  Google Scholar 

  84. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  85. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  86. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. 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.

    Article  CAS  PubMed  Google Scholar 

  88. Ruan J, Li H, Chen Z, Coghlan A, Coin LJ, Guo Y, et al. TreeFam: 2008 Update. Nucleic Acids Res. 2008;36:735–40.

    Article  Google Scholar 

  89. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19(9):1655–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. 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.

    Article  CAS  PubMed  Google Scholar 

  93. Pickrell JK, Pritchard JK. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 2012;8(11): e1002967.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Wilson GA, Rannala B. Bayesian inference of recent migration rates using multilocus genotypes. Genetics. 2003;163(3):1177–91.

    Article  PubMed  PubMed Central  Google Scholar 

  96. 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.

    Article  Google Scholar 

  97. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  99. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. 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.

    Article  Google Scholar 

  101. 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.

  102. Li H, Durbin R. Inference of human population history from individual whole-genome sequences. Nature. 2011;475(7357):493–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  104. 2b-RAD squencing of Harmonia axyridis populations. NCBI. 2023; https://www.ncbi.nlm.nih.gov/bioproject/ PRJNA951539 (Last accessed April 5, 2023).

  105. Sethuraman A. arunsethuraman/harmoniapopgen: Final (Final). Zenodo; 2023. https://doi.org/10.5281/zenodo.7796464.

Download references

Acknowledgements

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.

Funding

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.

Author information

Authors and Affiliations

  1. Department of Entomology, College of Plant Protection, Nanjing Agricultural University, Nanjing, People’s Republic of China

    Hongran Li, Yansong Wang, Xiaohan Shu, Ling Meng & Baoping Li

  2. Shenzhen Branch, Guangdong Laboratory of Lingnan Modern Agriculture, Genome Analysis Laboratory of the Ministry of Agriculture and Rural Affairs, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen, People’s Republic of China

    Hongran Li & Yan Peng

  3. Department of Biological Sciences, California State University, San Marcos, CA, USA

    Bryce Summerhays, Yumary Vasquez, Hannah Vansant, Christy Grenier, Nicolette Gonzalez, Khyati Kansagra, Ryan Cartmill & Arun Sethuraman

  4. Department of Life and Environmental Sciences, University of California, Merced, CA, USA

    Yumary Vasquez

  5. Empresa Brasileira de Pesquisa Agropecuária (Embrapa), Brasilia, DF, Brazil

    Edison Ryoiti Sujii

  6. Department of Entomology, University of Kentucky, Lexington, KY, USA

    Xuguo Zhou & John J. Obrycki

  7. Department of Agroecology, Flakkebjerg Research Centre, Aarhus University, Aarhus, Denmark

    Gábor L. Lövei

  8. ELKH-DE Anthropocene Ecology Research Group, University of Debrecen, Debrecen, Hungary

    Gábor L. Lövei

  9. Department of Zoology & Ecology, Hungarian University of Agriculture & Life Sciences, Godollo, Hungary

    Gábor L. Lövei

  10. Department of Biology, San Diego State University, San Diego, CA, USA

    Arun Sethuraman

Authors
  1. Hongran Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Yan Peng
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yansong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Bryce Summerhays
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xiaohan Shu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Yumary Vasquez
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Hannah Vansant
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Christy Grenier
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Nicolette Gonzalez
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Khyati Kansagra
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Ryan Cartmill
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Edison Ryoiti Sujii
    View author publications

    You can also search for this author in PubMed Google Scholar

  13. Ling Meng
    View author publications

    You can also search for this author in PubMed Google Scholar

  14. Xuguo Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  15. Gábor L. Lövei
    View author publications

    You can also search for this author in PubMed Google Scholar

  16. John J. Obrycki
    View author publications

    You can also search for this author in PubMed Google Scholar

  17. Arun Sethuraman
    View author publications

    You can also search for this author in PubMed Google Scholar

  18. Baoping Li
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

The study was conceived by B.L. and designed by B.L., J.O., and H.L. Insect sample collections were made by H.L., Y.W., X.S., and L.M. The data analyses were mainly done by H.L., Y.P., X.S., and. A.S. with contributions from B.S., Y.V., H.V., C.G., N.G., K.K., R.C., and E.S. The manuscript was written by H.L., and A.S., and finalized by B.L. with contributions from G.L., J.O., and X.Z. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Arun Sethuraman or Baoping Li.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

 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.

Additional file 2:

 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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12915-023-01638-7

Keywords

BMC Biology

ISSN: 1741-7007

Contact us