Aim: A growing number of studies suggest that adaptation of invasive species play key roles in their successful establishment in novel environments. However, adaptation of invasive species to climatic conditions remains poorly characterized. This study aimed to understand the population genetic structure produced by the cotton mealybug Phenacoccus solenopsis invasion and to identify preliminary signals of selection during its range expansion.

Location: China.

Methods: We examined genetic structure of 11 population across China using SNPs, microsatellites and a segment of mitochondrial cox1 gene. ADMIXTURE, STRUCTURE and DAPC were used to infer population genetic structure; the dispersal routes were reconstructed by the DIYABC; SNPs potentially related to climate adaptation were identified by using four population differentiation methods and three environmental association methods.

Results: Strong genetic differentiation was found among populations with F_ST values ranging from 0.097 to 0.640 based on SNPs. Populations located at the northern-expansion edge exhibited the highest genetic differentiation and the lowest genetic diversity. Demographic analyses indicated that all populations were introduced from a single source population with small effective size and low recent gene flow. RDA analysis showed that climatic variables explained a higher proportion of genetic variance (43%) compared to population structure variables (15%). The top climatic variables associated with genetic differentiation were precipitation of the mean temperature of warmest quarter, mean temperature of driest quarter and isothermality. Genes related to climate candidate SNPs were mainly enriched to pathways of development, energy, and xenobiotic metabolisms.

Main conclusions: We found that extremely rapid and strong population genetic differentiation among populations appears to have developed after introduction in the cotton mealybug. Our study points to rapid neutral evolution and suggests possible climatic adaptation despite low genetic diversity in this invasive species.

Sample collection and DNA extraction

Female adults of P. solenopsis were collected from 11 locations across its distribution areas in China in 2017 (Table S1, Fig. 1a). Genomic DNA was extracted from the whole individual of P. solenopsis using DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) with an additional step of RNase treatment. Voucher DNA of specimens was stored in the Integrated Pest Management Laboratory of Beijing Academy of Agriculture and Forestry Sciences at -80°C.

SNP genotyping

We used the double-digest restriction site-associated DNA (ddRAD) method to develop genome-wide SNP markers (Peterson, Weber, Kay, Fisher, & Hoekstra, 2012). To select the endonuclease enzymes for digestion, we simulated digestion using the draft genome of P. solenopsis (Ma et al., 2019) as a reference through the DDsilico program (Rašić, Filipović, Weeks, & Hoffmann, 2014). Libraries were constructed for each individual as described elsewhere (Peterson et al., 2012). Twenty individuals from each population were used for genotyping. In brief, 120 ng of genomic DNA in a 50 µL volume was digested with NlaIII and Acil restriction enzymes (New England Biolabs, Beverly MA, USA) for 3 hours at 37 °C. Restricted DNA was purified with 75 µL (1.5X) SpeedBeads (GE, USA) and each individual was ligated to a unique pair of modified Illumina P1 (5 bp) and P2 adapters (4 bp) overnight at 16°C, followed by a heat-deactivation step under the conditions of 65 °C for 10 min, with 22 cycles at 20 °C of 1 min. These products with unique adapter of one population were pooled and cleaned with beads. We then selected cleaned DNA fragments with a size of 420-540 bp by BluePippin on a 2% gel cassette (Sage Sciences, USA). The selected pooled DNA was amplified in 20 µL reactions containing 5 µL DNA, 10 µL of PCR Phusion 2X master mix, 0.8 µL of forward primer (50 µM) and reverse primer (50 µM) as well as 3.4 µL H₂O using the conditions: 30 sec at 98°C, followed by 12 cycles at 98°C for 10 sec, 65°C for 30 sec, 72°C for 45 sec, a post-cycle incubation at 72°C for 5 min and holding at 4 °C. Final library was purified by 16 µL 0.8X speedbeads. Quality of ddRAD libraries was evaluated by Qubit 3.0 and Agilent Bioanalyses 2100, and sequenced by Illumina Hiseq 4000 platform to generate 150-bp paired-end reads.

We used Stacks version 2.0 for SNP calling (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013). FastQC version 1.1.5 (Andrews, 2013) was used to estimate the content of GC and quality of raw sequences. Raw reads with sufficiently high sequencing quality and correct barcode were retained. After quality control, ddRAD barcodes were removed by the process_radtags program (Catchen et al., 2013). Sequences were aligned to draft genome of P. solenopsis using Bowtie version 2.2.9 (Langmead & Salzberg, 2012) with default parameters, and exported the highest aligned core of SAM files. After removal of lower mapping rate individuals (< 80%), 210 samples were kept and processed for SNP identification by the ref_map.pl program (Catchen et al., 2013) with a minimum depth of five and a maximum likelihood framework.

SNPs meeting the following criteria were exported using populations: (i) SNPs that were included in at least 80% samples of a population; (ii) SNPs with a minor allele frequency (MAF) higher than 0.05; and (iii) SNPs with observed heterozygosity > 0.8. The final data sets were filtered using R package vcfR (Knaus & Grünwald, 2017) by retaining SNPs with coverage at least 10 X, a missing rate in one sample less than 5%, and samples with a missing rate of SNPs less than 20%.

Microsatellite genotyping and mitochondrial DNA sequencing

As a complement to SNPs, we genotyped 21 microsatellite loci from 24 individuals of P. solenopsis per population following the method described in Ma et al. (2019) (Table S2). Briefly, a universal PC tail (5′CAGGACCAGGCTACCGTG3′) (Blacket, Robin, Good, Lee, & Miller, 2012) was used to identify forward primer candidates for amplification. All PCR products were analyzed using an ABI 3730xl DNA Analyzer (Applied Biosystems, Foster, CA, USA) with the GeneScan 500 LIZ size standard (Applied Biosystems). Microsatellites were genotyped in GENEMAPPER version 4.0 (Applied Biosystems, USA). We sequenced a segment of the mitochondrial cox1 gene to validate the morphological identification of the specimens and level of genetic diversity on the mitochondrial genome, using the same primers, sequencing methods and analysis methods as mentioned previously (Ma et al., 2019).