Herbivore range expansion triggers adaptation in a subsequently-associated third trophic level species and shared microbial symbionts

Invasive species may change the life history strategies, distribution, genetic configuration and trophic interactions of native species. The diamondback moth, Plutella xylostella L., is an invasive herbivore attacking cultivated and wild brassica plants worldwide. Here we present phylogeographic analyses of P. xylostella and one of its major parasitoids, Cotesia vestalis, using mitochondrial markers, revealing the genetic diversity and evolutionary history of these two species. We find evidence that C. vestalis originated in Southwest China, then adapted to P. xylostella as a new host by ecological sorting as P. xylostella expanded its geographic range into this region. Associated with the expansion of P. xylostella, Wolbachia symbionts were introduced into local populations of the parasitoid through horizontal transfer from its newly associated host. Insights into the evolutionary history and phylogeographic system of the herbivore and its parasitoid provide an important basis for better understanding the impacts of biological invasion on genetic configuration of local species.

www.nature.com/scientificreports www.nature.com/scientificreports/ its most important biocontrol agents 10,11,21 occurring in 38 countries 10 . Whilst C. vestalis has been introduced to Australia, North America, and the Caribbean in over 20 classical biological control programs 10,22 , there are no records of it being introduced to Japan, Vietnam, Malaysia (Cameron Highlands) and China [22][23][24][25] . Yet, C. vestalis is reported to be among the most predominant parasitoids of P. xylostella across parts of East Asia [25][26][27] , most likely associated to its tolerance to high temperatures 28 and insecticides 25,26 , and suggesting it may be native to this region.
The evolutionary success and diversification of insects has been aided, to some extent, by their complex associations with microorganisms 29,30 . Over the recent decades, Wolbachia symbionts have attracted intensive research effort due to their diverse behavioral effects in a broad range of insect hosts, ability to manipulate the host reproductive system, and potential role in biological control of pests 31 . Wolbachia are generally assumed to be maternally inheritable, with vertical transfer from egg cytoplasm to offspring, though recent studies have demonstrated horizontal transfer from infected to uninfected species 32 .
Invasive species may change distribution, genetic configuration and trophic interactions of other native species and trigger rapid adaptation 33,34 , which is increasingly becoming a major focus of research in evolutionary biology. In the present study, using samples of P. xylostella and C. vestalis collected from five Asian countries of China, Nepal, Thailand, Malaysia and Vietnam, we used a set of mitochondrial genes to analyze the phylogeographic relationships to (1) reveal the genetic diversity and demographic history of the two species; (2) identify the geographic origin of C. vestalis; and (3) address the adaptation of C. vestalis to the invasive P. xylostella in this region. Considering the Wolbachia-arthropod associations and the potential role in intra-specific interactions, we also investigated the impacts of horizontal transfer of Wolbachia on the genetic configuration and co-evolution of local assemblages.

Results
Genetic diversity. A total of 1,621 bp DNA was obtained from concatenation of three P. xylostella mitochondrial genes (hereafter referred as p3m). From 323 P. xylostella individuals coming from 29 sampling locations (Table 1), we found 187 polymorphic sites and 212 haplotypes, representing a high haplotype diversity with an average of 0.931. We identified 174 haplotypes represented by single individuals, with the remainders represented by multiple P. xylostella. Nucleotide diversity was overall low with an average of 0.329%, except for the samples from three locations of Quanzhou in China (FJQZ, 0.606%), Cameron Highlands in Malaysia (MLCH, 0.906%) and Katmandu in Nepal (NPKT, 0.682%) ( Table 1).
For C. vestalis, based on 1,232 bp DNA from three mitochondrial genes (hereafter referred as c3m), we found a relatively low haplotype and nucleotide diversity with an average of 0.415 and 0.172%, respectively. From 324 individuals, 43 polymorphic loci and 29 haplotypes were identified, with 19 haplotypes coming separately from single individuals and the rest from multiple individuals. We also observed higher haplotype and nucleotide diversity from two sampling locations of Guiyang (GZGY: Hd = 0.667 and θ = 0.431%) and Yuxi (YNYX: Hd = 0.769 and θ = 0.795%) in Southwest China (Table 1).
Phylogeny and haplotype network. The p3m-based phylogenetic tree revealed overall low genetic differentiation among individuals in terms of the branch length, and no isolated clusters containing individuals from specific geographic regions or populations, regardless of high genetic differentiation identified between Wolbachia PlutWB1-infected and -unfected individuals (Fig. 1a). The Cytb-based network exhibited a star-like shape with many unique haplotypes present at the terminals (Fig. 1b), indicating recent population expansion by P. xylostella. In terms of this Cytb-based network, haplotype 2 (H2) was dominant and present in most sampled populations. The frequency of identified haplotypes was randomly distributed in each of the P. xylostella populations. The c3m-based phylogenetic tree (Fig. 2a) and haplotype network (Fig. 2b) showed four separate lineages with multiple basal clusters/lineages containing the samples from Southwest China. Lineage 1 comprised the samples from part of China, Malaysia and Vietnam, Lineage 2 consisted of individuals from China, Thailand and Nepal, while lineage 3 was represented by individuals from China only. Using the CoxI gene sequences of our C. vestalis samples and additional sequences from India, Kenya, Benin, Hungary, Malaysia, New Zealand and Russia, we constructed a CoxI-based phylogeny (Fig. 2c) that does not change the overall 3cm-based topology of the basal position and paraphyly of samples from Southwest of China (Fig. 2a,b).
Four primary lineages were present on the Wsp-based Phylogeny of Wolbachia (Fig. 3). Hosts involved in lineage 1 and 2 are C. vestalis (Lineage 1) and P. xylostella (Lineage 2), respectively, while lineage 4 was identified in both of these two species (Fig. 3).
Demographic history. Neutrality tests for P. xylostella were conducted using Tajima's D and Fu's Fs statistics ( Table 1). The p3m-based Tajima's D and Fu's Fs statistics were significantly negative (Tajima's D = −2.501, P < 0.001; Fu's Fs = −5.532, P < 0.05) when all sampled populations were considered as one group. A significantly negative Tajima' D value was associated with a significantly negative Fu's Fs value, suggePhysting a recent expansion of P. xylostella populations in this region. The p3m-based mismatch distribution was unimodal when all sampled individuals of P. xylostella were considered as one group for analysis (Fig. 4a1), further supporting the concept of a recent expansion of the P. xylostella populations in East Asia. The c3m-based mismatch distribution was multimodal when all sampled individuals were considered as one group, but the three defined lineages exhibited unimodal distributions ( Fig. 4a2-a5).

Discussion
The paraphyly nature of samples from Southwest China, located at the early branching nodes of both c3m-based ( Fig. 2a) and CoxI-based ( Fig. 2c) phylogeny, suggests that Southwest China is the geographical origin of C. vestalis. Such an inference was also supported by the high nucleotide polymorphism of C. vestalis populations (YNYX and GZGY) in Southwest China (Table 1). We also found that all our sampled C. vestalis individuals derived from Southwest China based on the c3m-based phylogenetic tree (Fig. 2a) and haplotype network (Fig. 2b). According to the CoxI-based global phylogeny of C. vestalis (Fig. 2c, and supported by c3m tree, Fig. 2b), C. vestalis individuals from African, Oceanian and European countries were evolutionarily closely related to those from Thailand (TLPH) and Nepal (NPKT).
C. vestalis is one of the most common parasitic wasps of P. xylostella 10 and previously considered to be native in Malaysia 24 , Japan 23 , and China 25 . It has also been suggested to originate from Europe based on the species description using Ukrainian specimens 21 . However, our Cox I-based phylogenetic analysis demonstrated that C. vestalis populations in Europe, Africa and Oceania derived from East Asia (Fig. 2b). Although no samples from the New World were included in this study, we speculate that the haplotypes from the New World may derive from the haplotypes of the Old World as C. vestalis was reported to be recently introduced into North America 22,35 and South America 21,35 as a biological control agent.
P. xylostella was recorded to colonize many regions (including East Asia) of the world in recent centuries 14,15 . A significantly negative Tajima' D value was associated with a significantly negative Fu's Fs value, suggesting a recent expansion of P. xylostella populations in East Asia, which is further supported by the unimodal p3m-based mismatch   vestalis sampled individuals, L is the length of DNA fragments, S is the segregating sites, h and Hd are the number and diversity of haplotypes, θ is nucleotide diversity, "−" denotes populations with <5 individuals or with only one haplotype and they were not used for calculation of population parameters, "|" symbolizes the separation of P. xylostella (left) and C. vestalis (right), the significance of statistic tests are indicated by "*" (P < 0.05), "**" (P < 0.01) and "***" (P < 0.001). The first 25 populations listed are from China, with the other five populations from Nepal (1), Thailand (1), Vietnam (1) and Malaysia (2). The samples within a country are listed in order from the most northerly to the most southerly latitude. Full information of the populations refers to www.nature.com/scientificreports www.nature.com/scientificreports/ distribution of all sampled individuals. Our metadata analysis revealed a regional expansion of C. vestalis populations associated with the invasion and colonization of P. xylostella in East Asia. A comparable observation has also been reported for two Diadegma parasitoids of P. xylostella in Europe 36 . The results of BEAST analysis also showed that the most recent common ancestor (TMRCA) of C. vestalis was between the TMRCAs of P. xylostella in Old World (OW) and Oceania (OC), suggesting that the population expansion of C. vestalis could be related to the regional invasion of P. xylostella into Oceania through East Asia. Based on our results, we propose that the regional distribution of C. vestalis is a case of ecological sorting 7 that add P. xylostella, which had become the dominant herbivore in brassica crops of East Asia 11,37 , to the hosts for C. vestalis 38 . Such a new trophic association may facilitate the rapid adaptation of C. vestalis, by the previously-documented means of altered the life history traits and physiological manipulation 39,40 , to P. xylostella.
The co-occurrence of plutWB1 41 in P. xylostella and C. vestalis (Lineage 4 in Fig. 3) suggested horizontal transfer of Wolbachia between herbivore and parasitoid. Our observation of Wolbachia-infected P. xylostella pupae and adults implied that the direction of this horizontal transfer was from P. xylostella to C. vestalis, given that C. vestalis would kill P. xylostella before pupation 35 . In addition, in concordance with findings of Delgado & Cook 41 , a distinct clade of five individuals (from different sampling locations) infected by plutWB1 were identified in the phylogenetic tree of P. xylostella (Fig. 1a), suggesting a long history of co-evolution between P. xylostella and this Wolbachia strain and a horizontal transfer scenario from the herbivore to its parasitoid. Another distinct clade (Lineage 3) in the wsp-based phylogenetic tree showed that the host of Wolbachia involved parasitoids, herbivores and predators, suggesting that the horizontal transfer of Wolbachia can occur across multiple trophic levels.
Li et al. 42 reported that Bemisia tabaci-associated Wolbachia can be horizontally transferred between infected and uninfected individuals via plants. This can be further supported by a recent study of flowers and wild megachilid bees, in which the plants can act as hubs for bacterial transmission between multiple organisms 43 . In the present study, using our wsp sequences and the NCBI-based wsp sequences, we found that the individuals infected by Wolbachia were involved with several species of herbivores (Lineage 3 in Fig. 3). We assumed that the presence of such a Wolbachia lineage (Lineage 3) in various herbivores might have come from food intake or potentially from horizontal transfer between the species. www.nature.com/scientificreports www.nature.com/scientificreports/ In this study, using an interactive system involving an invasive herbivore, P. xylostella, and its parasitoid, C. vestalis, we demonstrated how the invasion of an alien host herbivore (P. xylostella) could significantly affect the genetic variation of a higher trophic species (C. vestalis). Through ecological sorting, this parasitoid could have switched from an original local host to P. xylostella and, as an important agent of classical biological control, it underwent significant, human-aided, population expansion. In addition, during this expansion, the endosymbiont Wolbachia (plutWB1) in P. xylosetlla was also introduced into local species. Our work provides a comprehensive picture (Fig. 5) of how invasion by an alien species can trigger significant evolutionary changes in a newly associated parasitoid and the transmission of exotic bacteria into local species, leading to the formation of new biological interactions and the genetic configuration of local species.

Methods
Sample collection and species identification. We collected P. xylostella and C. vestalis samples from the same or nearby cabbage and broccoli fields in China (25), Nepal (1), Thailand (1), Vietnam (1) and Malaysia (2) from 2012 to 2014 ( Fig. S1; Table S1). Twenty-eight samples of P. xylostella and C. vestalis were collected from common locations, and one sample of P. xylostella was collected in Chongqing, China (CQ) and its counterpart C. vestalis sample from nearby, Luzhou in Sichuan Province, China (SCLZ) (Fig. S1; Table S1). The second and third instar larvae of P. xylostella were maintained on vegetable leaves for emergence of parasitoids. P. xylostella pupae and adults, and C. vestalis cocoons were morphologically identified and preserved in 95% ethanol. Specimens were stored at −80 °C prior to DNA extraction. A total of 323 P. xylostella and 324 C. vestalis individuals were used in this study. We used a 600 bp mitochondrial gene sequence (CoxI) ( Table S2) and DNA barcoding criteria to individually confirm the species identity of P. xylostella and C. vestalis based on BOLD 44 . DNA Extraction and sequencing. Total genomic DNA was extracted from individual insects using DNeasy Blood and Tissue Kit (Qiagen, Germany). Three mitochondrial genes, CoxI, cytochrome b (Ctyb), and NADH dehydrogenase subunit I (NadhI) were sequenced for both P. xylostella and C. vestalis (Table S2). Primers were developed for Cytb and NadhI of P. xylostella as well as Cytb and NadhI of C. vestalis using Primer Premier version 5 (Premier Biosoft International, Palo Alto, CA, USA) based on the reference mitochondrial genomes. Primers of other gene segments were referred to published references (Table S2). www.nature.com/scientificreports www.nature.com/scientificreports/ PCR was conducted using the Mastercycler pro system (Eppendorf, Germany) under the following conditions: an initial denaturation for 2 min at 94 °C, followed by 35 cycles of 10 s at 96 °C, 15 s at specific annealing temperature of each genes (Tables S2), and 1 min at 72 °C, and a subsequent final extension for 10 min at 72 °C. Amplified products were purified and bidirectionally sequenced using the ABI 3730xl DNA Analyzer by Sanboyuanzhi Biotechnology Co., Ltd. (Beijing, China). All the sequences were deposited in Genebank database (accession number from KX604356 to KX606864).
Infection of P. xylostella and C. vestalis by Wolbachia was determined using ∼600 bp products of the wsp gene amplified with specific primers (Table S2). A positive control of PCR reaction (with DNA of Wolbachia infected samples as templates) was used to test infection of Wolbachia in samples.

Genetic analysis.
Sequences for each of the gene fragments were aligned using MEGA5.2 45 . All mitochondrial sequences for each of the individuals of both insect species were aligned independently using MAFFT-7.037 46 . Conservative regions selected by Gblock-0.91b 47 were used for gene concatenation, which was performed by Sequence-Matrix-1.7.8 with default parameters 48 . Parameters of genetic diversity, haplotype diversity (Hd) and nucleotide diversity (θ), were calculated using the DnaSPv5 49 . Populations with <5 individuals (3 populations of P. xylostella and 12 of C. vestalis) were not included in the calculation of the parameters related to the genetic diversity (Table 1). www.nature.com/scientificreports www.nature.com/scientificreports/ Phylogenetic and network analysis. Phylogeographic analysis can be used to explore the evolutionary history of a species [50][51][52] or to examine the temporal and spatial effects on co-evolutionary relationships of closely related species 53,54 . This type of study can help reveal the impacts of biological expansions on local communities over wide spatial scales. Using the sequences of three concatenated mitochondrial genes, phylogenetic relationships were constructed for P. xylostella (here after p3m) and C. vestalis (here after c3m). We selected Cotesia flapvis, which is from the Cotesia genus, as the outgroup for C. vestalis phylogenetic tree construction. The CoxI  www.nature.com/scientificreports www.nature.com/scientificreports/ sequences of C. vestalis and C. flapvis (outgroup) from NCBI (http://www.ncbi.nlm.nih. gov/) were also downloaded for construction of a global phylogenetic tree of C. vestalis. The phylogeny of wsp gene was developed using the NCBI-based wsp sequences with the best hit (with <3 gaps and > = 99% identity) when BLAST conducted using the wsp sequences from this study, plus additional sequences of the wsp gene of Wolbachia parasitizing P. xylostella and C. vestalis.
Phylogenetic inferences were performed using the neighbor-joining (NJ) and maximum likelihood (ML) methods by PAUP*4.0b10 55 . The software MrModeltest version 2.3 56 was used to select the best-fit nucleotide substitution model. The General Time Reversible model was used with invariable sites and a gamma-shaped distribution of rates across sites (GTR + I + G) based on the Akaike Information Criterion (AIC).
The network analysis was conducted for mitochondrial genes of P. xylostella and C. vestalis using median-joining algorithm implemented in the software Network, version 4.6.1.3 57 . We constructed the haplotype networks of both species for individual genes as well as the concatenated mitochondrial gene sequences. Haplotype type and frequency for each population were also recorded.
Demographic analysis. We calculated the Tajima's D and Fu's Fs for each of the populations (≥5 individuals) of the two species based on the concatenated mitochondrial genes using Dnasp V5 49 . Analyses of mismatch distributions were also performed for both species. For C. vestalis, mismatch distribution was analyzed for not only a collection of all samples, but also three major clusters based on the phylogenetic tree of three concatenated mitochondrial genes.
BEAST 58 was also used to calculate the coalescent time of lineages in P. xylostella and C. vestalis, based on CoxI sequences, which has a reported mutation rate 53 . For P. xylostella, more samples from the Old World and Oceania 36 were included to increase precision in coalescent time inference. As the mutation rates varied among insect lineages, we used lognormal relaxed clock while estimating the evolutionary timescales. The chain of Markov Chain Monte Carlo (MCMC) was set to 50 million with Log parameters in every 5000.

Data Availability
Data generated during the study available in Genbank with the primary accession number from KX604356 to KX606864.