Genomic changes in the biological control agent Cryptolaemus montrouzieri associated with introduction

Abstract Biological control is the main purpose of intentionally introducing non‐native invertebrate species. The evolutionary changes that occur in the populations of the introduced biological control agents may determine the agent's efficiency and the environmental safety. Here, to explore the pattern and extent of potential genomic changes in the worldwide introduced predatory ladybird beetle Cryptolaemus montrouzieri, we used a reduced‐representation sequencing method to analyze the genome‐wide differentiation of the samples from two native and five introduced locations. Our analyses based on a total of 53,032 single nucleotide polymorphism loci showed that beetles from the introduced locations in Asia and Europe exhibited significant reductions in genetic diversity and high differentiation compared with the samples from the native Australian range. Each introduced population belonged to a unique genetic cluster, while the beetles from two native locations were much more similar. These genomic patterns were also detected when the dataset was pruned for genomic outlier loci (52,318 SNPs remaining), suggesting that random genetic drift was the main force shaping the genetic diversity and population structure of this biological control agent. Our results provide a genome‐wide characterization of polymorphisms in a biological control agent and reveal genomic differences that were influenced by the introduction history. These differences might complicate assessments of the efficiency of biological control and the invasion potential of this species but also indicate the feasibility of selective breeding.

. Similar to other biological introductions or invasions, unexpected evolutionary changes can occur in populations of biological control agents, which might decrease the efficiency of pest control or increase the risk to local environments (Fauvergue, Vercken, Malausa, & Hufbauer, 2012;Guillemaud, Ciosi, Lombaert, & Estoup, 2011). Recently, there have been more reports of undesired side effects on local biodiversity caused by introduced biological control agents through rapid expansion, nontarget attack, and competition, and the invasiveness of these agents is usually related to population genetic changes (e.g., Harmonia axyridis in North and South America and Europe, Roy et al., 2016 andCactoblastis cactorum, Zimmermann, Bloem, &Klein, 2004 and Coccinella septempunctata, Losey et al., 2012 in North America). Thus, understanding the population genetics of biological control agents can aid in the exploration of their patterns of potential evolutionary changes and may be beneficial for biological control programs.
The mealybug destroyer, Cryptolaemus montrouzieri Mulsant (Coleoptera, Coccinellidae), is native to Australia and is a predator specialist of mealybug (Ślipiński, 2007). Over the last century, this predatory ladybird has been introduced to at least 64 countries or regions for classical or augmentative biological control purposes (Kairo, Paraiso, Gautam, & Peterkin, 2013). There is currently only a single available report of the side effects of such introductions on local environments, which showed interference with other biological control agents (Annecke, Karny, & Burger, 1969). However,

C. montrouzieri populations with different introduction histories can
show changes in several life history traits, such as development time, oviposition number, and performance under adverse conditions (Li, Zou, De Clercq, & Pang, 2018). We also detected range expansion in introduced populations in southern China (mentioned by Li et al. 2017). Furthermore, host range tests in a laboratory population revealed nontarget attack abilities (Maes, Grégoire, & De Clercq, 2014). These potential risks, including range expansion and nontarget attack by ladybird populations, might result from evolutionary changes that occurred during and after their introduction.
The evolution of introduced biological control agents is usually explained by processes such as population founding or admixture events, rather than being attributed to phenotypic plasticity or selection (Andraca-Gómez et al., 2014;Kajita, O'Neill, Zheng, Obrycki, & Weisrock, 2012;Lombaert et al., 2011;Retamal et al., 2016). Furthermore, high genetic differentiation is often observed among introduced populations (Fauvergue et al., 2012). We previously conducted several genetic analyses of native and introduced populations of C. montrouzieri; using 12 simple sequence repeat (SSR) markers, we found high and significant genetic differentiation between native and introduced populations (Li et al., 2017).
Additionally, we detected signals of positive selection by scanning for nonsynonymous mutations across mitochondrial protein-coding genes (Li, Liang et al., 2016). Improvements in next-generation sequencing (NGS) and bioinformatic tools have spurred the development of genome-wide genetic markers for studying the ecology and evolution of nonmodel organisms (Davey et al., 2011). Recently, this genomic methodology has been applied for the detection of genetic variation in field and mass-reared populations of biological control agents and to elucidate the genetic basis of target traits (Lommen, Jong, & Pannebakker, 2017;van de Zande et al., 2014). In this work, we aim to explore the pattern and extent of potential evolutionary changes in worldwide introduced C. montrouzieri. We conducted reduced-representation sequencing of samples from two native and five introduced locations of C. montrouzieri and performed a genome-wide survey of intra-specific polymorphism. All of the non-native beetle populations analyzed here result from artificial transfers, and their introduction histories at the studied locations are well documented.
Analyses of genomic diversity and structural differences between the populations partially revealed the evolutionary forces functioning during human-mediated species introduction. Finally, the contributions of these genomic findings to biological control applications are discussed.

| Sampling and DNA extraction
A total of 128 individuals of C. montrouzieri from two native and five introduced locations were analyzed in this study (Figure 1 and   Supporting Information Table S1). Individuals from the native lo- For each individual, total genomic DNA was extracted using a CTAB-based procedure (Milligan, 1988). DNA quality and quantity were determined using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, USA) and via electrophoretic separation in an agarose gel. Only samples with a DNA concentration >20 ng/μl and slight degradation were used.

| Reduced-representation genome sequencing
Approximately 500 ng of genomic DNA from each sample was processed for reduced-representation genome sequencing using the specific-locus amplified fragment sequencing (SLAF-seq) technique (Sun et al., 2013), with slight modifications. Briefly, in a pilot SLAF experiment, the restriction enzymes Hpy166II and EcoRV-HF (NEB, USA) were selected because they generate abundant and highquality nonrepetitive SLAFs that are relatively evenly distributed across the genome of the model beetle species Tribolium castaneum (Richards et al., 2008). In addition, the DNA of a plant species, Oryza sativa ssp. japonica, was used as a control to assess the normal rate of enzyme digestion. Genomic DNA from each sample was completely digested with these two enzymes, and a single-nucleotide A overhang was added to the digested fragments. The A-tailed fragments were then ligated to dual-index sequencing adapters (Kozich, Westcott, Baxter, Highlander, & Schloss, 2013). PCR was performed in a mixture of diluted restriction-ligation samples including Taq DNA polymerase (NEB), dNTPs, and adapter primers. The PCR products were purified using Agencourt AMPure XP (Beckman, UK) and subjected to electrophoresis through a 2% agarose gel. Fragments in the size range of 250-350 bp were excised and purified using a The ratio of high-quality reads exhibiting Q30 (representing a quality score of 30, indicating a 0.1% chance of an error and, thus, 99.9% confidence) and the guanine-cytosine (GC) content of the raw reads were calculated for quality control.
Only one SNP per SLAF tag was retained to remove the effects of physical linkage. Deviations from Hardy-Weinberg equilibrium (HWE, p < 0.05) for each SNP were calculated using PLINK, and those loci that significantly deviated from HWE at more than three studied locations were excluded. The filtered SNPs were mapped to the annotated transcriptome of the species (Li, Pan, De Clercq, Ślipiński, & Pang, 2016;Zhang et al., 2012) using BLASTn, with a cut-off E-value of 10e-5 and a minimum length of 70 bp. Finally, the filtered SNPs were converted into Arlequin, BayeScan, FASTA, and PED formats for subsequent analyses using PGDSpider v2.1.0.0 (Lischer & Excoffier, 2012).

| Detection of F ST outliers
We detected F ST outliers among the total SNPs using methods implemented in Arlequin 3.5.1.1 (Excoffier & Lischer, 2010) and BayeScan 2.1 (Foll & Gaggiotti, 2008). We scanned for outlier loci deviating from the expected null distribution of F ST in a system of hierarchically structured populations (Arlequin's hierarchical island model, AH, Excoffier, Hofer, & Foll, 2009), which was identified as the best structure in the following analyses (six-group structure: QL+BM, SY, SA, GT, SZ, and TP). In addition, we used Arlequin's nonhierarchical island model (ANH) to test for F ST outliers. In both the AH and ANH analyses, we ran 20,000 simulations, with 10 simulated groups in AH, 100 demes per group, and a maximum expected heterozygosity of 0.5. We also employed the Bayesian Thus, we drew Venn diagrams to determine the intersection of high/low F ST outliers detected by AH, ANH, and BayeScan. We also drew heatmaps to present the allele frequencies of the high F ST outliers. The F ST outliers detected by the three methods were used for further analysis.

| Genetic diversity and population structure
To detect the contribution of different loci to the diversity of this biological control agent, the total SNPs were further separated into three datasets, which included the high and low F ST outliers and the remaining loci (assumed as neutral loci). Genetic diversity and population structure were analyzed based on the total loci, the neutral loci, and the high F ST outlier loci.
We employed Arlequin to calculate genetic diversity parameters at each of the studied locations, including the percentage of polymorphic loci (Poly%), nucleotide diversity (Pi), observed (H O ), and expected heterozygosity (H E ). Comparisons of these parameters between native and introduced locations were performed using t tests in SPSS 21 (IBM SPSS Statistical, Chicago, USA).
We calculated overall and pairwise F ST values using Arlequin, with significance estimated based on 10,000 permutations. Then, population structure was studied through (a) the reconstruction of phylogenetic trees, (b) the estimation of ancestry proportions, and (c) an analysis of molecular variation (AMOVA). Neighbor-joining (NJ) trees were reconstructed in MEGA6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) using the p-distance between individuals and 1,000 bootstrap replicates. The software ADMIXTURE (Alexander, Novembre, & Lange, 2009) was employed to estimate ancestry proportions using different numbers of ancestral clusters (K). We explained the genetic structure by obtaining the lowest cross-validation error from K = 2 to 10. Then, AMOVA was performed in Arlequin with 10,000 permutations to further test for the population structure identified by clustering with ADMIXTURE.

| SNP calling
A total of 504.25 million reads were obtained from the raw sequencing data. The sequencing quality of each sample is shown in Supporting Information  Figure S1).

| Detection of F ST outliers
The three genome scanning methods detected a total of 521 high

| Genetic diversity and population structure
The genetic diversity parameters of each studied location are shown in Table 2 (all 53,032 loci) and Supporting Information Table S3 ( for K = 4-7 were quite close. For K = 5, SY and SZ shared the same ancestral cluster, while for K = 4, QL+BM and TP were also grouped together. For K = 7, SA appeared to be an admixed population of two ancestral clusters. However, a too-large K value may have led to an overfitting effect, which forced the algorithm to subdivide the populations. The percentage of variation associated with differences among groups (F CT ) was higher and significant for K = 6 structuring (26.94%, p = 0.0489, Supporting Information Table S4) than for K = 5 (1.92%, p = 0.4360, Supporting Information Table S4) and K = 4 (4.29%, p = 0.2776, Supporting Information Table S4). However, these AMOVAs had limited power since some of the groups were formed by one location only.
The high F ST outlier loci generated a population structure that was in many aspects qualitatively similar but quantitatively stronger.
The NJ tree showed that SA and SY+SZ were extremely divergent  Table S5). The ADMIXTURE analysis showed that K = 7 resulted in the lowest cross-validation error (Figure 5e), with the individuals from the same locations belonging to a unique cluster.
For K = 6, the subdivision of SA was possibly an artefactual effect.
The AMOVAs showed that the K = 5 and K = 4 structuring explained nearly none of the total variance among populations (Supporting Information Table S4).

| Genomic changes in biological control agents arising from introduction
In this study, we demonstrated very high genetic differentiation among introduced populations of the ladybird beetle C. montrouzieri.
Our results based on a large number of genome-wide SNP markers are overall consistent with the differentiation patterns in this system detected with a few SSR markers (Li et al., 2017). More generally, high genetic differentiation based on traditional genetic markers (e.g., SSRs) is commonly detected among different introduced populations of biological control agents (Lombaert et al., 2014;Retamal et al., 2016;Sethuraman, Janzen, & Obrycki, 2015). In theory, genetic differentiation among populations can be largely attributed to genetic drift or selection, and the repeated process of establishment of introduced populations from native ones or other introduced populations may contribute to the strength of these effects (Fauvergue et , 2012). Because environmental conditions (e.g., temperature and diet) differ from native conditions of C. montrouzieri, we considered that not only genetic drift but also selection might have contributed to the observed differences between populations. Our reducedrepresentation genomic sequencing approach provides a relatively powerful dataset for a better understanding of the extent and the forces contributing to the genetic changes in the introduction history of this biological control agent.

| Role of genetic drift and selection in genomic population changes
The studied introduced populations were artificially transferred from the native areas in Australia or introduced populations on other continents. This intentional introduction process shows some similarities to biological invasions (Blackburn, Lockwood, & Cassey, 2015) but could potentially exhibit even more extreme founder effects due to the total isolation of the introduced populations from the native populations. In this study, a significant decrease in genetic diversity was detected in the introduced populations, which is consistent with founder effects together with strong genetic drift during the early stages of introduction. Additionally, our pairwise F ST results, NJ trees, and clustering results supported the notion that all the introduced populations are deeply divergent from the native populations. Since these patterns are very similar in the datasets of all loci and only neutral loci, we suggest that genetic drift-in addition to allele frequency changes during the founding event-was the main force shaping the genetic diversity and population structure of the introduced populations of this biological control agent. Only ~100 years has passed since the introduction of C. montrouzieri for biological control was initiated (Kairo et al., 2013), but the production of more than 10 generations per year can contribute to rapid genomic changes. We also found that some alleles in the introduced populations were not present in the native locations. Apart from probably rare de novo mutations, it is possible that these alleles are too rare to have been picked up in the limited sampled individuals within a population. An alternative explanation is that some of these genomic differences go back to naturally existing differentiation in the native range of this species, rather than being caused by human activities. Our analyses of native samples from two distinct locations It would require a detailed survey of the natural range (Zemanova, Broennimann, Guisan, Knop, & Heckel, 2018;Zemanova, Knop, & Heckel, 2016) to determine the full extent of natural genetic varia-  (Beaumont & Nichols, 1996;Fischer, Foll, Heckel, & Excoffier, 2014). However, with the increase in differentiation caused by genetic drift, loci with higher F ST values can be mistaken for evidence of positive selection, since the two mechanisms produce similar signals in genomic variation (Cruickshank & Hahn, 2014;Currat et al., 2006;Vera, Díez-del-Molino, & García-Marín, 2016

| Implications for biological control
The extent of heritable phenotypic differences between the introduced or native populations of C. montrouzieri has not been characterized at a larger scale, but there is likely variation in F I G U R E 5 Ancestral clustering using ADMIXTURE based on all 53,032 loci (a, b), 52,318 neutral loci (c, d), and 521 high F ST outlier loci (e, f). The cross-validation error of different numbers of ancestral clusters (K) and the individual ancestry bar plots for K = 4-7 are shown the performance as a biological control agent. In fact, ladybirds from the SY population have a significantly longer developmental time and perform better under cold stress but worse under starvation than beetles from SA when tested under controlled laboratory conditions (Li et al., 2018). The large genetic differentiation between these and among the other populations indicates several potential issues related to biological control applications.
First, intra-specific variation might complicate the assessment of biological control efficiency and invasion potential. The evolution of biological control agents not only influences their effectiveness in pest control but also generates potential side effects on local biodiversity (Roderick, Hufbauer, & Navajas, 2012). Such changes could certainly result from phenotypic plasticity that enables ladybirds to respond to varying environments (Murren et al., 2015). However, relatively often, the successful establishment of introduced species has been associated with rapid genetic changes. These genetic changes could be caused by genetic drift, as most likely in the present study, and/or adaptive responses to novel selection pressures (Fauvergue et al., 2012). For example, in a parasitoid wasp introduced to New Zealand from South America for the control of weevil pests, changes in biotype frequency and better performance were detected after 10 years, indicating local adaptation (Phillips et al., 2008). Further studies should take the potential for intra-specific variation and different reactions within C. montrouzieri into account.
Second, the extent of the detected genetic changes might indicate the feasibility of selective breeding for improvement of traits of interest (Lommen et al., 2017). Using methods such as genome-wide association studies (GWAS), it is possible to identify and characterize loci affecting target traits for further genomicbased selection (Visscher, Brown, Mccarthy, & Yang, 2012). For example, similar population genomic techniques are being applied to parasitoid wasps, predatory bugs, and predatory mites in the Breeding Invertebrates for Next Generation BioControl Training Network (BINGO-ITN, http://www.bingo-itn.eu/en/bingo.htm).
Our genome-wide description of polymorphisms in the ladybird C. montrouzieri demonstrates ample genetic variation and thus suggests that such approaches might be promising avenues for further improvements of the performance of this biological control agent.

ACK N OWLED G EM ENTS
This work was supported by the National Natural Science

CO N FLI C T O F I NTE R E S T
None declared.