Adaptive evolution following colonization can affect the impact of invasive species. The fall webworm (FWW) invaded China 40 years ago through a single introduction event involving a severe bottleneck and subsequently diverged into two genetic groups. The well-recorded invasion history of FWW, coupled with a clear pattern of genetic divergence, provides an opportunity to investigate whether there is any sign of adaptive evolution following the invasion. Based on genome-wide SNPs, we identified genetically separated western and eastern groups of FWW and correlated spatial variation in SNPs with geographical and climatic factors. Geographic factors explained a similar proportion of the genetic variation across all populations compared to climatic factors. However, when the two population groups were analyzed separately, environmental factors explained more of the variation than geographic factors. SNP outliers in populations of the western group had relatively stronger response to precipitation than temperature-related variables. Functional annotation of SNP outliers identified genes associated with insect cuticle protein potentially related to desiccation adaptation in the western group and genes associated with lipase biosynthesis potentially related to temperature adaptation in the eastern group. Our study suggests that invasive species may maintain evolutionary potential to adapt to heterogeneous environments despite a single invasion event. The molecular data suggest that quantitative trait comparisons across environments would be worthwhile.

Sample collection

Samples of FWW were collected from 16 locations across its distribution range in China; 14 of these have previously been used for population genetics analysis using microsatellite markers (Cao et al., 2016). The other two newly collected populations were obtained from the expansion fronts of FWW in 2017-2018. Larvae of FWW were each sampled from different silk webs at each sampling location to reduce the chances of collecting siblings. In total, 306 larvae of FWW were obtained and used for DNA extraction, library construction, and genotyping, with 13-20 individuals per population.

Library construction, SNP calling, and filtering

Genomic DNA was isolated from larvae individually using a DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). We used the ddRAD method to develop genome-wide SNPs for FWW (Peterson et al., 2012). Genomic DNA from each individual was digested by the restriction enzymes NlaIII and AciI for 3 hours at 37 °C (Aguilar et al., 1979; Li et al., 2018). Then we used 67.5 µl (1.5×) SpeedBeads (GE) to purify the digested DNA. A pair of uniquely modified Illumina P1 (5 bp) and P2 adapters (4 bp) were ligated to the digested DNA at 16 °C overnight. A heat-deactivation step was used to end the ligating reaction under conditions of 65 °C for 10 min and 22 cycles at 20 °C for 1 min. We pooled ligated products with a unique adapter into one library, followed by a purifying step using (1.5×) SpeedBeads (GE). Fragments of 420 - 540 bp were selected using BluePippin on a 2% gel cassette (Sage Sciences, Beverly, MA, USA) and then amplified using 12 PCR (polymerase chain reaction) amplification cycles. We used 64 µl 0.8× SpeedBeads to purify the amplified libraries. The quantity and quality of each library were evaluated using Qubit 3.0 and Agilent Bioanalyses 2100. The Illumina NovaSeq 6000 platform was used for sequencing to obtain 150-bp paired-end reads.

We used Stacks version 2.52 to filter the low-quality sequencing data and call SNPs (Catchen et al., 2013). Raw sequencing reads were demultiplexed and trimmed using the process_radtags. Reads for each individual were mapped to the reference genome of FWW with a size of 510.5 Mb from NCBI (Assembly: GCA_003709505.1 ASM370950v1) (Wu et al., 2018) using Bowtie version 2.3.5.1 (Langmead et al., 2012). SNPs were called using a maximum likelihood framework and filtered with populations implemented in Stacks, VCFtools version 0.1.16 (Danecek et al., 2011), and the R package vcfR (Knaus et al., 2017) based on the following criteria: (a) samples with a mapping rate less than 80% were removed; (b) SNPs with a sequencing depth higher than eight and less than 500 were removed; (c) samples and SNPs with a missing rate higher than 10% in the corresponding dataset were removed; (d) SNPs with a minor allele count lower than 10 were removed; (e) SNPs with observed heterozygosity of > 0.75 across all populations were removed; (f) SNPs with a p-value of Hardy-Weinberg equilibrium (HWE) lower than 10^-7 in all populations were removed to generate dataset of neutral SNPs. In order to reduce the influence of linkage on population structure inferences, we retained only SNPs separated by at least 1000 bp (Lowry et al., 2017).

References

Aguilar, J. D., & Riom, J. (1979). Nemoraea pellucida (Meigen), A new parasite of Hyphantria cunea (Drury) [France; fall webworm]. Bulletin De La Société Entomologique De France, 84, 204-207.

Cao, L. J., Wei, S. J., Hoffmann, A. A., Wen, J. B., & Chen, M. (2016). Rapid genetic structuring of populations of the invasive fall webworm in relation to spatial expansion and control campaigns. Diversity and Distributions, 22, 1276-1287.

Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., & Cresko, W. A. (2013). Stacks: an analysis tool set for population genomics. Molecular Ecology, 22, 3124-3140.

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., anks, E. B., Depristo, M. A., . . . Sherry, S. T. (2011). The variant call format and VCFtools. Bioinformatics, 27, 2156-2158.

Knaus, B. J., & Grünwald, N. J. (2017). VCFR: A package to manipulate and visualize variant call format data in R. Molecular Ecology Resources, 17, 44-53.

Langmead, B., & Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nature Methods, 9, 357-359.

Li, B. Y., Gao, Q., Cao, L. J., Hoffmann, A. A., Yang, Q., Zhu, J. Y., & Wei, S. J. (2018). Conserved profiles of digestion by double restriction endonucleases in insect genomes facilitate the design of ddRAD. Zoological Systematics, 43, 341-355.

Lowry, D. B., Hoban, S., Kelley, J. L., Lotterhos, K. E., Reed, L. K., Antolin, M. F., & Storfer, A. (2017). Breaking RAD: an evaluation of the utility of restriction site-associated DNA sequencing for genome scans of adaptation. Molecular Ecology Resources, 17, 142-152.

Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., & Hoekstra, H. E. (2012). Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PloS ONE, 7, e37135.

Wu, N., Zhang, S., Li, X., Cao, Y., Liu, X., Wang, Q., . . . Zhan, S. (2018). Fall webworm genomes yield insights into rapid adaptation of invasive species. Nature Ecology and Evolution, 3, 105-115.

Here we provided VCF files generated and its population map generated in this study. Three VCF files were included.

1.fww_invariant+SNP_miss20_DP3.GQ20.vcf.gz, includes SNPs and invariant sites of all populations;

2.fww.ddRAD.all.vcf.gz, includes SNPs of all populations;

3.fww.4fds.vcf.gz, includes four degenerated SNPs of all populations;

4.fww_265_popmap.txt, includes a population map of all individuals.

The three VCF files and a population map file can be opened by VCFtools and used as input files for population genetic diversity, population genetic structure, demographic inference, and outlier scanning analysis.