Population analysis reveals genetic structure of an invasive agricultural thrips pest related to invasion of greenhouses and suitable climatic space

Abstract Biological invasions of pests into climatically unsuitable areas can be facilitated by human‐regulated environments, in which case there may be an impact on genetic structure through population processes and/or adaptation. Here, we investigated the population genetic structure of an invasive agricultural pest, Thrips palmi, in China, which has expanded its distribution range through using greenhouses. Early invaded populations showed a relatively higher level of genetic diversity than recently expanded greenhouse populations. Strong population genetic structure corresponded to a pattern of isolation by distance, with no recent gene flow and low historical gene flow among populations, reflecting limited ongoing dispersal. A genetic signature of population expansion was detected in early invaded populations and three northern populations from greenhouses, suggesting that the greenhouse environments facilitated expansion of this species. Redundancy analysis showed that the independent effects of environment and geography could explain 51.68% and 32.06% of the genetic variance, respectively. These findings point to climate‐ and greenhouse‐related spatial expansion, with the potential for adaptation by T. palmi. They emphasize the contribution of human‐regulated environments on the successes of this invasive species, a situation likely to apply to other invasive species that use greenhouse environments.

thrips (Gerson & Weintraub, 2012). Several studies have shown evidence of evolutionary adaptation of species under natural climatic conditions (Csilléry, Rodríguez-Verdugo, Rellstab, & Guillaume, 2018;Hoffmann, 2017). These factors are not necessarily independent, with many small invertebrates rapidly adapting to controlled environmental conditions when they are reared in the laboratory (Hoffmann & Ross, 2018), suggesting that species might also adapt to controlled greenhouse conditions. The rapid development of insecticide resistance in many pests under greenhouse conditions has previously been documented (Gholam & Sadeghi, 2016), illustrating the potential for evolutionary shifts in such environments.
Genetic variation across populations reflects the effects of population processes like genetic drift under limited gene flow, as well as (at loci under selection or linked to loci under selection) adaptive evolution to local environmental conditions. For invasive species in their introduced areas, population genetic structure can be formed by multiple introductions from genetically differentiated source populations, and demographic events such as bottlenecks and founder effects (Barrett, 2015;Bock et al., 2015;Cao et al., 2017;Konecny et al., 2013;Lee, 2002;Tsuchida, Kudô, & Ishiguro, 2014). Invasive species with a clear population history provide opportunities to examine climatic adaptation through geographic or temporal comparisons (Egizi, Fefferman, & Fonseca, 2015;Hoffmann, 2017;Lee, 2002;Tsuchida et al., 2014;Yamanaka, Tatsuki, & Shimada, 2008), and the effects of human-regulated environments on the genetic structure of such populations can also be investigated.
Thrips species represent small insects from the order Thysanoptera. Several species of thrips cause serious damage to crops worldwide (Mouden, Sarmiento, Klinkhamer, & Leiss, 2017;Reitz, Gao, & Lei, 2011). The cosmopolitan pest of western flower thrips, Frankliniella occidentalis (Thysanoptera: Thripidae), is a temperate species originating from southwestern USA and has invaded other areas around the world since the late 1970s (Kirk & Terry, 2003). It was introduced into China around 2003 and rapidly spread into most areas of China within 10 years. Population genetic analysis suggested multiple independent invasions of F. occidentalis into China . A lack of isolation by distance and the fact that distant populations were genetically similar suggested patterns of movement in this thrips linked to human activities . Spatial spread of F. occidentalis seems not restricted by environmental conditions, which contributes to its cosmopolitan status as a pest and its lack of population genetic structure in introduced areas.
Another important pest thrips, the melon thrips, Thrips palmi (Thysanoptera: Thripidae), was first described in tropical regions of Sumatra in Indonesia (Karny, 1925). Early reports of T. palmi were exclusively from southeastern and southern Asian countries of Thailand, India, Pakistan, Malaysia, and Philippines (Ananthakrishnan, 1955;Bhatti, 1980). The species was introduced and became established during the second half of the twentieth century across South-East Asia, South America, the Caribbean, Florida, Australia, and West Africa (Cannon, Matthews, & Collins, 2007;Kawai, 1990Kawai, , 2001. Climatic conditions are thought to have restricted the spatial expansion of T. palmi (McDosald, Bale, & Walters, 1999). In northern Japan, this species could not overwinter outdoors, while in Australia it was restricted to warmer northern regions (Layland, Upton, & Brown, 1994). However, T. palmi in temperate regions in Japan and China occurs in greenhouses, where it has become one of the most severe pests of vegetables (Kawai, 1990;Reitz et al., 2011). Spatial distribution of T. palmi expanded northward rapidly in recent years in China in greenhouses. With clear historical records and sensitivity to climatic variables, T. palmi is an ideal species to examine the effects of greenhouse environments on genetic structure.
In this study, we examined the population genetic structure of T. palmi across China and tested forces that shape this structure.
Based on its tropical origin and stepping-stone dispersal into temperate region, we hypothesized a structure formed by neutral processes such as genetic drift under geographic isolation. Based on sporadic colonization and the likelihood of strong selection by novel environmental conditions in temperate regions, we also hypothesized that greenhouse conditions should influence genetic structure of T. palmi as reflected by levels of diversity and genetic divergence. By understanding population genetic variation of T. palmi as well as evolutionary processes affecting it, factors that facilitate the spread of this invasive species and its potential to become a wider pest can be identified. This is one of the few studies that focus on population genetic structure of an agriculture insect pest in greenhouse conditions. The results shed light on understanding the impact of artificial environments on population genetic structure of pests more generally.

| Sample collection and DNA extraction
We collected specimens of T. palmi from 14 geographic populations on eggplants and cucumbers across areas of its distribution in China and one population from Japan (Table 1, Figure 1a). Specimens were collected from 20 to 30 points in a field or greenhouse, separated by a distance of about 5-10 m. One individual was used from each collection point to reduce the likelihood of genotyping siblings.
Populations from northern China were collected from greenhouses, where T. palmi is not known to overwinter in fields, while populations from southern China were collected from fields. In total, 348 female individuals from 15 populations were used for genotyping. We used female individuals for population genetic analysis due to the haplodiploid sex determination of T. palmi, in which haploid individuals develop to males while diploid individuals develop to females.
Genomic DNA was extracted for each specimen with a DNeasy Blood & Tissue Kit (Qiagen).

| Microsatellite genotyping and DNA sequencing
Twenty-six genome-wide microsatellite loci for T. palmi were used for genotyping (Table S1). PCR products were labeled by fluorescence for length determination with the method of Blacket, Robin, Good, Lee, and Miller (2012). Conditions for PCR amplification were described in previous publications (Cao, Li, et al., 2016;Song, Cao, Wang, Li, & Wei, 2017;. PCR products were analyzed using ABI 3730xl DNA Analyzer (Applied Biosystems) with the GeneScan 500 LIZ size standard (Applied Biosystems) by Tsingke Biotechnology Co. Ltd. Microsatellite loci were genotyped using GENEMAPPER 4.0 (Applied Biosystems) and checked for stuttering and large allele dropout using MICRO-CHECKER version 2.2.3 (van Oosterhout, Hutchinson, Wills, & Shipley, 2004).
To characterize mitochondrial variation and validate correct identification of the specimens, we sequenced a fragment of cox1 gene on DNA barcoding region of animals by the primer pairs FO-AF (5′ TTTCGTCTAACCATAAAGATATCGG 3′) and FO-AR (5′ TAAACTTCTGGGTGCCCAAAAAATCA 3′) modified from previously published primers (Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994) based on the mitochondrial genome sequence of Franklienella occidentalis (Yan et al., 2012). Polymerase chain reaction (PCR) was conducted using the Mastercycler pro system (Eppendorf) under standard conditions with an annealing temperature of 52°C. PCR components were added as recommended by the manufacturer of Takara LA Taq (Takara Biomedical). Amplified products were purified and sequenced on an ABI 3730xl DNA Analyzer by Tsingke Biotechnology Co. Ltd.

| Genetic diversity analysis
For microsatellite loci, null allele frequency was estimated using FreeNA with 10,000 bootstraps (Chapuis & Estoup, 2007). Then, we examined Hardy-Weinberg equilibrium (HWE) for each locus at each population and linkage disequilibrium between each pair of loci at each population by GENEPOP version 4.2.1 (Rousset, 2008).
Heterozygosity excess and deficit were tested by GENEPOP version 4.2.1. GENCLONE version 2.0 (Arnaud-Haond & Belkhir, 2007) was used to estimate the total number of alleles (A T ) and the unbiased expected heterozygosity (H ET ) (Nei, 1978). Furthermore, we compared the number of alleles (A S ) and standardized expected heterozygosity (H ES ) among samples with different sample sizes using a rarefaction method in GENCLONE. Allelic richness (A R ) and allelic richness of private alleles (P AR ) were calculated with a rarefaction approach in HP-RARE version 1.1 (Kalinowski, 2005) on a minimum sample size of 15 diploid individuals.

| Population genetic structure analysis
For microsatellite loci, we used two Bayesian cluster methods and discriminant analysis of principal components (DAPC) to investigate genetic structure across the populations.
The first Bayesian method, Bayesian analysis of population genetic structure (BAPS), was conducted using BAPS version 6.0 (Cheng, Connor, Sirén, Aanensen, & Corander, 2013). This method can infer the clustering of individuals by incorporating spatial information. The maximum number of genetically diverged clusters (K) was set to 5, 8, 10 or 14, to ensure convergence and consistency of the results. For each K, 15 repeat runs were performed.  (Pritchard, Stephens, & Donnelly, 2000). This method can estimate the ancestral gene frequencies and the admixture proportions for each individual (Lawson, van Dorp, & Falush, 2018). An admixture model with correlated allele frequencies was chosen. Thirty replicates for each K (from 1 to 10) were run with 200,000 Markov Chain Monte Carlo (MCMC) iterations after a burn-in of 100,000 iterations.
The optimal value of K was determined using the Delta (K) method (Evanno, Regnaut, & Goudet, 2005) by submitting the outputs of STRUCTURE to STRUCTURE HARVESTER WEB version 0.6.94 (Earl & Vonholdt, 2012). The membership coefficient matrices (Q-matrices) of replicated runs for each K were combined using CLUMPP version 1.1.2 (Jakobsson & Rosenberg, 2007) with the Greedy algorithm and then visualized using DISTRUCT version 1.1 (Rosenberg, 2004).
Third, discriminant analysis of principal components (DAPC) was performed using adegenet version 2.0.1 (Jombart, 2008) in the R F I G U R E 1 Sampling locations (red points in figure a) and population genetic structure of the 15 Thrips palmi populations inferred from BAPS analysis (pie charts in figure a) and STRUCTURE (b) based on 26 microsatellite markers.
In the BAPS analysis, distributions of the seven identified clusters are shown by different colors. The color proportions in the pie charts represent the frequency of each cluster in a population. The red circle shows the central and eastern group (CE), and the blue circle shows the northern group (NT) of populations collected from greenhouses (a). In the STRUCTURE analysis, the optimal K determined based on delta K was 4. Clusters of individuals when K is 4, 5, and 6 are presented. Four outlier populations of JANP, HNSY, SCPZ, and YNXS, as well as two population groups of EC and NT, are presented. Codes for collection sites are shown in Table 1 100 This method does not rely on biological models that can be a complementary analysis of model-based methods such as BAPS and STRUCTURE.
For mitochondrial DNA, haplotypes and their distribution in each population were analyzed in DnaSP version 5.10 (Librado & Rozas, 2009). Phylogenetic relationships among mitochondrial cox1 haplotypes were constructed using a split networking method implemented in SPLITSTREE version 4.13.1 (Huson & Bryant, 2006) with the neighbor-net method under a distance model of K2P after 1,000 bootstraps.

| Gene flow analysis
Recent and historical gene flow among populations was estimated using Bayesian methods implemented in BayesAss version 3.0.4 (Wilson & Rannala, 2003) and Migrate version 3.7.2 (Beerli & Felsenstein, 2001), respectively. BayesAss was modeled to estimate gene flow for the past 1-2 generations. We ran 100 million steps with different start seeds after preliminary runs for adjustment of mixing parameters for allele frequencies and inbreeding coefficients.
The trace outputs of ten longer runs were combined using Tracer

| Demographic history analysis
We estimated the variation of effective population size for each population based on microsatellite loci using coalescent algorithms implemented in Migraine version 0.5.4 (Rousset, Beeravolu, & Leblois, 2017). This software uses PAC-likelihood (product of approximate conditional likelihoods) based on quantities inherent to the importance of sampling algorithms (Cornuet & Beaumont, 2007 Then, 15 iterations with 800 points and 1,000 runs per point were conducted starting from points of last iteration. In the third round, 3 iterations with 800 points and 20,000 runs per point were conducted starting from points of the last iteration. Estimation of past demographic changes is usually challenging (Beaumont, 2010;Hey, 2010;Rousset et al., 2017). When the OnePopVarSize model was used, estimation became instable. We evaluated the confidence of OnePopVarSize-based estimation by comparing the similarity of 2Nμ values estimated by OnePop and OnePopVarSize as well as the confidence intervals of the estimates.

| Partitioning geographic and climate effects on genetic variation
We analyzed geographic/climatic effects on population genetic variation using two methods. First, isolation by distance (IBD)  bio5, max temperature of warmest month; bio8, mean temperature of wettest quarter) and two precipitation-related (bio15, precipitation seasonality; bio18, precipitation of warmest quarter) variables.
Both full (environmental and geography) and partial (environmental or geography) models of RDA were analyzed in the R package

| Gene flow and effective population size
Among the 15 populations, seven showed an expansion of effective population size, two showed a reduction (JANP and GDZD), and six showed no significant difference between current and ancestral effective population sizes (Table 5). Based on our criteria, estimates using OnePopVarSize model fitted well in at least five populations, three of which showed population expansion. Three of the five greenhouse populations showed population expansion (Table 5).

| Geographic and climate effects on genetic variance
A Mantel test indicated a significant correlation between genetic distance and geographic distance (r = .6305, p = .001).  (Figure 3a). When geographic variables were constrained in the RDA analysis, habitat, isothermality of temperature (bio3), and maximum temperature of warmest month (bio5) were highly correlated with genetic distance (Figure 3b).

| Origin of Thrips palmi
Population genetic analyses provide novel approaches for investigating dispersal and invasion routes of species compared to historical records (Boissin et al., 2012;Estoup & Guillemaud, 2010;Lombaert et al., 2014;Rollins, Woolnough, Wilton, Sinclair, & Sherwin, 2009 showing that this species originated from tropical counties of Asia and spread into pantropical regions in its early stages of range expansion (Cannon et al., 2007).

| Population processes shaping the genetic structure of Thrips palmi
We tested possible population processes that may shape the genetic

| Greenhouse-related population structure and climatic adaptation in Thrips palmi
The genetic structure of T. palmi populations in China follows geography, whereas five populations collected from greenhouses in northern China form a genetic group as inferred from the microsatellites.  Note: 95% of coverage confidence interval are provided in bracket; pGSM, parameter for GSM model of microsatellite mutation; 2Nμ, the current effective population size scaled by mutation rate; 2N anc μ, ancestral effective population size scaled by mutation rate; type of population size variation was determined by overlap of 95% of coverage confidence interval between 2Nμ and 2N anc μ estimated by OnePopVarSize model. The bolded and underlined numbers show confidence estimation from OnePopVarSize model determined based on differences between estimation of 2Nμ by two models and confidence intervals of estimated 2N anc μ. Five populations collected from greenhouses are underlined. NA, indicates the value is either extremely low or high.
F I G U R E 3 RDA analysis on genetic variance explained by the environmental effects of climate and habitat. A full RDA model was run by considering environmental and geographic effects simultaneously (a), and a partial RDA model was run by constraining geographic effects to analyze the correlation of environmental variables (b). Individuals from the same population are indicated by circles with the same color. PCNM2-5, geographic variables; bio3, isothermality; bio5, maximum temperature of warmest month; bio8, mean temperature of wettest quarter; bio15, precipitation seasonality; bio18, precipitation of warmest quarter; habitat, samples collected from greenhouse or field. Correlations of each variable are indicated by an arrow. Long arrows indicate a high correlation between the variable and genetic distance to selected loci; however with a limited number of microsatellites,  (Kawai, 1985;McDosald et al., 1999). Based on limited temperature adaptation, introduced T. palmi were successfully eradicated in England after they were introduced (MacLeod, Head, & Gaunt, 2004). In northern Japan, the species could not overwinter outdoors and survived cold winters in greenhouses. The population dynamics of T. palmi can also be influenced by humidity (Su, Chiu, & Lin, 1985). In Australia, the distribution of T. palmi is restricted to northern regions with warm temperatures in winter, and its southern distribution may be limited by prevailing aridity (Layland et al., 1994).
These limits point to the likelihood of migration and population size of the thrips being limited by environmental conditions in China.
Greenhouse conditions in northern China can help T. palmi to persist despite low temperatures outside, and they also provide humid conditions. The ecological opportunity provided by greenhouses may explain why there was a signature of populations expanding in three of five greenhouse populations. Such ecological opportunities could promote adaptive radiation by generating genetic changes in organisms (Stroud & Losos, 2016;Yoder et al., 2010). Further studies involving well-designed experiments are needed to explore this from a biological and genomic perspective.

| Distinct genetic structure between two thrips species
Thrips are tiny insects that are difficult to detect in quarantine.
Invasion and dispersal of such small insects are usually mediated by human activities such as plant transport. Both F. occidentalis and T. palmi are serious invasive pest of agriculture. Nevertheless, our studies showed distinct genetic structure in these two thrips across their distribution range of China. The F. occidentalis thrips showed population genetic differentiation unrelated to its geographic distribution in China, pointing to multiple introductions and human-mediated dispersal in sporadic directions . In contrast, in T. palmi there was a high level of genetic differentiation across its spatial distribution in China-related partly to geographic distance.
This pattern of population genetic structure is similar to that found in species with stepping-stone dispersal (Cao, Wei, Wei, Hoffmann, Wen, & Chen, 2016;Kimura & Weiss, 1964;Wei et al., 2015) rather than sporadic long-distance dispersal . The pattern is also congruent with historical records in China; T. palmi was first reported in southern China and occurred in northern areas much later (Yi & Liang, 2001;Zhang, Han, & Fu, 1985). Thrips palmi is particularly abundant on the foliage of cucumbers and eggplants, and dispersal of this thrips may have been slow because these vegetables are transported as fruit rather than foliage (Huang, 1989

| Implications for pest management
Greenhouses extend the ability to produce crops in different seasons and regions. The protected conditions in greenhouses lead to serious pest problems that are different from outdoor fields, such as outbreaks of whiteflies, mites, and thrips (Gerson & Weintraub, 2012). Our findings illustrate that they can also influence the genetic structure and variation of the pest organism through affecting population processes. These, in turn, could influence processes like the evolution of pesticide resistance that are influenced by climatic conditions (Maino, Kong, Hoffmann, Barton, & Kearney, 2016;Wimmer, Hoffmann, & Schausberger, 2008). Compared to the temperate species F. occidentalis, which can be readily dispersed across a wide range through human activities, dispersal of the tropical species into temperate areas will be restricted by environmental conditions. However, environmental restrictions can be overcome by using artificial environments and there may also be evolutionary adaptation to the relatively constant conditions present in glasshouses. Thus, evolutionary adaptation by T. palmi should be further considered both to cold conditions and to the relatively constant conditions of greenhouses given that adaptive shifts can influence predictions around pest distributions and pest abundance.

| CON CLUS IONS
In this study, we revealed population genetic structure and potential evolutionary forces affecting genetic differentiation of an invasive agricultural pest, T. palmi. We found that environmental factors and geographic isolation correlated with genetic differentiation. The ecological opportunity provided by greenhouses may contribute to recently expanded populations of T. palmi.

ACK N OWLED G EM ENTS
We

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

DATA AVA I L A B I L I T Y S TAT E M E N T
All data were achieved in Dryad under DOI: https ://doi.org/10.5061/ dryad.bp27sd7 (Cao et al., 2019).