Strong population genetic structure of an invasive species, Rhynchophorus ferrugineus (Olivier), in southern China

Abstract The red palm weevil (RPW), Rhynchophorus ferrugineus (Olivier), was initially reported in China in the 1990s and is now considered one of the most successful invasive pests of palm plants in the country. A total of 14 microsatellite loci and one mitochondrial cytochrome oxidase subunit Ι (cox I) gene fragment were used to investigate the genetic characteristics and structure of R. ferrugineus in southern China. High levels of genetic differentiation among populations and significant correlations between genetic and geographical distances indicated an important role of geographical distance in the distribution of the RPW in southern China. High gene flow between Fujian and Taiwan province populations illustrated the increased effects of frequent anthropogenic activities on gene flow between them. Genetic similarity (i.e., haplotype similarity) indicated that RPW individuals from Taiwan and Fujian invaded from a different source than those from Hainan. To some extent, the genetic structure of the RPW in southern China correlated well with the geographic origins of this pest. We propose that geographical distance, anthropogenic activities, and the biological attributes of this pest are responsible for the distribution pattern of the RPW in southern China. The phylogenetic analysis suggests that the most likely native sources of the RPW in southern China are India, the Philippines, and Vietnam.


| INTRODUCTION
Invasive species colonizing new environments not only have an economic impact but also often have a negative effect on biodiversity (Sakai et al., 2001). Genetic variation and population genetic structure are closely related to the evolutionary potential of a species and its resistance to a hostile environment (Grant, 1991;Hughes & Dorn, 2002;Nei, 1978). Characterizing the current population genetic structure of an invasive species can facilitate the identification of the source of the infestation, possible invasion paths, and genetic consequences of invasion (Bock et al., 2015). Genetic studies of invading pest species are crucial to understanding how organisms are able to successfully adapt to and colonize new environments.
The exact time that the RPW invaded China is unknown. It was not until the 1990s that scientists began to pay attention to this pest. In southern China, ornamental palm plants are widely used as roadside shade trees, especially in some tourist cities. The first detection of the RPW in Guangdong Province (Zhongshan City) was reported in 1997 (Liu, Zhao, Xu, Chen, & Huang, 2009). This pest was later detected in Hainan (first detected in Wenchang City) in 1998 and rapidly spread across the whole island (Liu, Peng, & Fu, 2002;Wu & Yu, 2001). Due to palm plant transportation, the RPW has been widely dispersed throughout southern China over the last thirty years (He, Yu, Chen, Zhang, & Chen, 2013;Hu, 2006;Li & Yao, 2006;Liu, 2007;Xu, 2005;Zhang, 2007). Considering the serious economic impact and the ecological destruction caused by this weevil, a better understanding of its genetic diversity and of the relationships among different geographical populations is necessary for the development of effective pest management strategies. Various factors, such as geographical barriers, climate conditions, and historical processes, as well as the dispersal abilities of the species, would influence the population genetic structure of an organism (Sun, Lian, Navajas, & Hong, 2012). RPW adults rely on flight for their short-distance dispersal but can also be carried for long distances by human activities (Qin et al., 2009), which may result in a complex population structure. Furthermore, great distances or geographical barriers are major factors that can prevent gene flow between populations, leading to genetic divergence (Campbell, Mrazek, & Karlin, 1999). There are two straits located in southern China that are likely to influence the distribution of this pest: Taiwan Strait and Qiongzhou Strait. Therefore, special geographical locations and frequent urban landscape construction increase the possibility of RPW dispersal over a wider range in southern China. We hypothesize that large geographical distances combined with frequent human activities may interact to affect the structure of RPW population genetics.
Examining how dispersal or movement patterns govern the genetic structure of this pest can help to better understand its ability to colonize rapidly.
The genetic variation of the RPW has been studied in most infested areas, including the Middle East, the Mediterranean Basin, and South Asia (Abulyazid, Kamel, Sharawi, & El-Bermawi, 2002;El-Mergawy, Nasr, Abdallah, & Silvain, 2011;Rugman-Jones, Hoddle, Hoddle, & Stouthamer, 2013). High genetic similarity was found between the RPW populations from the Kingdom of Saudi Arabia (KSA) and those from Indonesia, while no similarity was detected between the Egyptian and KSA or Indonesian populations (Abulyazid et al., 2002). Genetic variation studies of the RPW in the Middle East and the Mediterranean Basin using the cox I gene revealed that RPW populations can be subdivided into different subpopulations under the influence of genetic drift favored by founder events following three different invasive routes over 30 years (El-Mergawy et al., 2011).
Anthropogenic activities were considered one of the important factors affecting the distribution of the RPW (Qin et al., 2009). Rugman-Jones et al. (2013) investigated the large-scale genetic structure of the RPW over the distribution areas. Although there have been many studies about the RPW conducted around the world, fewer studies have focused on the genetics of the RPW in China, which is adjacent to the native range of this pest. In this study, samples were obtained from a total of eight infested locations in southern China, and the samples were analyzed to describe the genetic structure of the RPW and to determine genetic differentiation among populations using microsatellite loci and cox I sequences. Elucidating the population genetic structures of this pest will enhance our understanding of its invasion pathways and subsequently aid in developing an effective control strategy for the RPW in southern China.

| Samples and DNA extraction
A total of 112 RPW adults were collected from palms in different areas of southern China during 2013 and 2014 (Table S1). Specimens were preserved in absolute ethanol and stored at −20°C. Total genomic DNA extraction and DNA quality detection were performed according to Wang, Zhang, Hou, and Tang (2015).

| Microsatellite genotyping and cox I sequencing
Individuals were genotyped with 14 microsatellite loci (Capdevielle-Dulac, 2011, personal communication, Table S2). Forward primers were labeled at the 5′ end with fluorescent dyes (FAM, HEX, or ROX; applied by JieLi Biology, Shanghai, China). The PCR was performed according to Wang et al. (2015). Amplifications were performed with an initial denaturation step of 3 min at 94°C; 40 cycles of 30 s at 94°C, 30 s at the primer-specific annealing temperature, and 1 min at 72°C; and a final elongation step of 5 min at 72°C.
We amplified and sequenced cox I to identify the genetic relationships among the different populations sampled.
Following PCR product detection by 2% agarose gel electrophoresis and sequencing (Sangon, Shanghai, China), the preliminary sequence analyses were assessed. The allele sizes of microsatellite loci were manually checked and scored based on trace data for further analysis, and mitochondrial PCR products were directly sequenced from both strands with both primers.

| Microsatellite data analyses
Because the sample size collected in Sichuan Province was limited and unsuitable for microsatellite analysis, a total of seven populations were genotyped for 14 microsatellite loci. The presence of null alleles was investigated using Micro-Checker 2.2.3 (van Oosterhout, Hutchinson, Wills, & Shipley, 2004). Genetic diversity statistics, including the mean allele number in each locus (A), mean observed heterozygosity (Ho), and expected heterozygosity (He), were estimated using Excel Microsatellite Tools (Park, 2001). Deviations from Hardy-Weinberg equilibrium (HWE), F IS values, and pairwise F ST values (Weir & Cockerham, 1984) were calculated using the program GENEPOP (Raymond & Rousset, 1995) (http://genepop.curtin.edu.au/). To minimize the bias of genetic diversity statistics affected by null alleles, genotypes were corrected based on the allele frequencies by the EM algorithm using FreeNA (Dempster, Laird, & Rubin, 1977). F ST values were corrected for null alleles present at microsatellite loci by implementing the ENA algorithm within the program FreeNA (Chapuis & Estoup, 2007). Isolation by distance (IBD) using F ST /(1−F ST ) and the semi-logarithmic geographical distance matrix (km) among all populations were determined by GENEPOP. The mean allelic richness (AR) was estimated using the FSTAT program (Goudet, 2001). BOTTLENECK version 1.2.02 (Cornuet & Luikart, 1996;Piry, Luikart, & Cornuet, 1999) was used to assess deviation from the expected heterozygote excess relative to allelic diversity across 14 loci under different mutation models (two-phased model of mutation, TPM; stepwise mutation model, SMM; infinite alleles model, IAM).
Deviations of observed heterozygosity relative to that expected at drift-mutation equilibrium were evaluated with the Wilcoxon signrank test with 1,000 iterations (Luikart, Allendorf, Cornuet, & Sherwin, 1998). The population structure of RPW was inferred from microsatellite data using the Bayesian clustering algorithm of STRUCTURE 2.3.3 (Pritchard, Stephens, & Donnelly, 2000), with an initial burn-in period of 5 × 10 4 iterations followed by a run of 10 6 Markov chain Monte Carlo (MCMC) repetitions. We used the admixture ancestry model (Hubisz, Falush, Stephens, & Pritchard, 2009) and the correlated allele frequency model. A series of simulations were run from K = 2 to 7 (from K = 2 to 12 when microsatellite data of Wang et al., 2015 were included) with 20 replicates for each K value. The most likely value of K was calculated according to Evanno, Regnaut, and Goudet (2005) in STRUCTURE HARVESTER (Earl & vonHoldt, 2012). DISTRUCT 1.1 (Rosenberg, 2004) was used to graphically display the genetic structure results. The distribution of genetic variance within and among clusters was determined using ARLEQUIN 3.5 (Excoffier & Lischer, 2010) through analysis of molecular variance (AMOVA). Gene flow among populations was estimated using the formula N e m = (1−F ST )/4F ST (Wright, 1931).

| Mitochondrial cox I sequence and data analyses
DNA sequence data were derived from 112 samples obtained from eight different locations (see Table S4). All sequences were aligned by CLUSTAL-X 1.81 (Thompson, Gibson, Plewniak, & Higgins, 1997).
A median joining (MJ) network analysis was performed in NETWORK 4.6.1.1 (Bandelt, Forster, & Röhl, 1999). The maximumlikelihood method (ML; best model: Kamura 3-parameter; bootstrap values: 1,000) by MEGA 5.03 (Tamura et al., 2011) was used to construct the phylogenetic tree based on the haplotype data. AMOVA was used to test the genetic differentiation among different clusters by calculating F-statistics from haplotypes with 10,000 permutations using the ARLEQUIN 3.5 software. Mismatch distribution as well as neutrality tests were also calculated using ARLEQUIN 3.5. In addition, a Bayesian skyline plot (BSP) analysis was performed using the BEAST v 2.1.3 program (Bouckaert et al., 2014) to examine the demographic history based on coalescent theory. TRACER v 1.5 (Drummond, Rambaut, Shapiro, & Pybus, 2005) was used to show the results of the BSP analysis.

| Genetic diversity of microsatellites
The mean numbers of alleles and AR at a microsatellite locus were 3.74 and 3.32, respectively. The results showed that He and Ho ranged from 0.273 to 0.650 and from 0.250 to 0.621, respectively.
Null alleles were observed at five microsatellite loci with a frequency ranging from 0.022 to 0.283. Null alleles existed in the P1F8 locus in

| Genetic variation of cox I sequences
Eight populations were used for cox I analysis. A total of 544 bp of the cox I sequences were generated after final alignment (GenBank accession numbers KY629226-326). Thirty haplotypes with global values of haplotype diversity and nucleotide diversity and an average number of nucleotide differences were found in the eight populations (total Hd = 0.914, π = 2.313%, k = 12.584) ( Table 1). The highest levels of haplotype and nucleotide diversities were exhibited in the GDSZ population (Hd = 0.964, π = 2.089%), while the lowest levels were observed in the TWTZ population (Hd = 0.222, π = 0.041%).  Table 2).

| Population genetic structure
The pairwise F ST differences based on cox I showed significant differentiation in 24 of the 28 population pairs (Table 2).
According to the method described by Evanno et al. (2005), We combined our cox I sequences with 29 sequences obtained from GenBank (Table S3)   Pop, population label; A, number of alleles; AR, allelic richness; R Ho, observed heterozygosity calculated by the raw data, C Ho, observed heterozygosity calculated by the corrected data; R He, expected heterozygosity calculated by the raw data, C He, expected heterozygosity calculated by the corrected data; R F IS , fixation index calculated by the raw data, C F IS , fixation index calculated by the corrected data; Nh, no. haplotypes; Nph, number of private haplotypes; Hd, haplotype diversity; π, nucleotide diversity; k, average number of nucleotide differences; **indicate significant deviations from HWE at p < 0.01 to the analysis, suggesting a significant correlation between genetic and geographical distances among the different populations ( Figure 5).

| Gene flow and demographic history
The gene flow values suggested the presence of variable levels of gene flow between populations (Table S4) (Table   S5).
Together with Tajima's D (2.099, p = .933) and Fu's Fs (−0.563, p = .604) tests, mismatch distribution analysis showed a multimodal mismatch graph when all populations were considered as a single group (Fig. S2a), indicating no evidence of recent demographic expansion. These results were confirmed by a Bayesian skyline plot (BSP) analysis based on mitochondrial data (Fig. S2b).

| Genetic diversity of RPW populations in southern China
By genotyping different populations from the major RPW outbreak centers in China using microsatellite markers, we determined that genetic variation was generally low within the southern China populations. These results were consistent with those previously found for the RPW in Fujian Province (Wang et al., 2015). The average number    (Yang, Sun, Xue, Li, & Hong, 2012).
Heterozygote excess was found in the GXNN population with a negative F IS value, which is consistent with the result of the heterozygote excess test using BOTTLENECK software (IAM, p = .00003 and TPM, p = .00131). We therefore deduced that individuals in the GXNN population had experienced a bottleneck effect and did not meet the criteria for HWE.
High haplotype diversity (Hd > 0.5) was present in most populations, especially in individuals from the two southernmost provinces, Hainan (HNWC1 and HNWC2) and Guangdong (GDSZ), indicating F I G U R E 5 Scatter plots of genetic isolation by geographical distance and genetic distance among Rhynchophorus ferrugineus populations using microsatellite and mitochondrial sequences that these populations have experienced more propagule pressure.
Multiple introduction or repeated introduction events will enhance diversity among invasive populations. The local climate seems to play a role in the expansion process of invasive species. Located in the tropical zone, with relatively high temperatures and rainfall amounts and a wide diversity of potential hosts, Hainan is considered an optimal habitat for the RPW. The high level of ecosystem invasibility of Hainan Island promotes colonization by invasive populations. The RPW has infested many coconut palms in Hainan and has had a serious economic impact on the island. Shenzhen City has the most ports in China and is also an important border city in Guangdong Province; the (international) trade of palms from native or infected areas to new areas across this province increased the possibility of infection by invasive species. Based on the palm transportation records in China, frequent human-aided dispersal events occurred in these two provinces (Li & Yao, 2006;Xu, 2005;Zhang, 2007). We can therefore deduce that individuals from these two provinces have been subjected to disproportionate invasion pressures, which enhance population diversity and consequently have resulted in an excess number of haplotypes during population colonization. Hainan provinces, gene flow between the populations was high, indicating that geographical barriers were not a main factor that prevented the distribution of the RPW in these regions. Furthermore, the Hercynian construction and the Maritime Silk Road construction may considerably increase human-aided transportation, which will enhance immigration between these regions. According to these results, it is possible to hypothesize that based on a broad perspective, geographical distances have historically played a very important role in the large-scale genetic structure of the RPW in southern China, although anthropogenic activities cannot be ignored.

| Population genetic structure and expanding events
In summary, after testing the genetic characteristics of the RPW in southern China, we found strong genetic structure, which indicated a significant connection between the geographical distribution and anthropogenic activities. More than one haplotype detected in one population revealed that multiple introduction events from other infected areas (including native or invasive regions) to China played a major role in the distribution of haplotypes in southern China. From palm plants and offshoot transportation, we can deduce that the RPW in China originated from different populations and followed various paths. The origin sources of this pest in southern China are likely from India, the Philippines, and Vietnam. To some extent, the construction of city gardens, which promoted palm plant transportation in China, widely influenced the distribution of the RPW. On the other hand, the successful invasion of the RPW in southern China is probably due to its biological attributes together with the existence of numerous suitable habitats and climates across southern China (Qin et al., 2009). The cryptic behavior, egg deposition inside plant tissue, and polyphagous nature (feeding on 40 palm species, even sugarcane [Nirula, 1956;Lever, 1969;Esteban-Duran et al., 1998;OJEU, 2008;Anonymous, 2013]), make detection difficult and facilitate the invasion of the RPW to new environments. Additionally, because of its high fecundity (total lifetime fecundities range from 58 to 531 progeny per female) (Wattanapongsiri, 1966), high population growth potential