Temporal genetic variation and dispersal patterns of Aedes albopictus (Diptera: Culicidae) in three different climatic regions of China

and Compared its insecticide its genetic which corresponds to present study examines how Ae. varies among different in and its dispersal The genetic variation and population structure of 17 Ae. albopictus populations collected from three climatic regions of China were investigated with 11 microsatellite loci and the mitochondrial COI gene.

in Ae. albopictus population studies on a global scale, and they continue to be a popular choice of genetic marker [15,[28][29][30]. Simultaneously and independently from Ae. albopictus microsatellites, the mitochondrial COI gene has frequently been used for mosquito barcoding and invasive species monitoring due to its conservation and high accuracy when distinguishing sequence variation at the interspecies level [31][32][33][34][35].
To evaluate how Ae. albopictus populations vary with environmental factors and decipher the potential dispersal patterns of Ae. albopictus among different climatic regions, a total of 17 Ae. albopictus populations were collected and examined via 11 microsatellite loci and the COI gene. The results provide guidelines for the future control of Ae. albopictus in China.

Materials And Methods
Eco-climatic characteristics of different climatic regions All mosquitoes were sampled from the three main climatic regions of China, including tropical, subtropical and temperate regions [36,37]. The annual accumulated temperature for the tropical region is above 8000℃, the lowest monthly average temperature is not less than 15℃, and the annual precipitation is approximately 1500~2000 mm; the annual accumulated temperature for the subtropical region is between 4500~8000℃, the lowest monthly average temperature is 0~15℃, and the annual precipitation is above 1000 mm; and the annual accumulated temperature for the temperate region is between 3000 and 4500℃, the coldest monthly average temperature is -28~0℃, and the annual precipitation is approximately 700 mm.

Mosquito sampling and DNA isolation
In the present study, ca. 600 Ae. albopictus larvae were sampled from seventeen geographically diverse sites across three environmentally distinct regions of China from June to August 2018; all the sampling site information is described in Fig. 1 and Additional le 1: Table S1. To avoid inbreeding interference, each pool of female mosquitoes for a given locality was collected from at least ve wild breeding places within 500 meters. For each of the samples, the larvae were reared separately and emerged at 25℃ ± 1℃ at 75% ± 5% relative humidity (RH) under a 14 h light/10 h dark (LD) photoperiod, and all female adult mosquitoes were identi ed under a microscope before DNA isolation [38,39]. The DNA isolation was conducted with the Qiagen DNeasy Blood and Tissue Kit (No. 69504, Qiagen, Germany) according to the manufacturer's protocol, and all the DNA samples were stored at -80°C.
Microsatellite isolation, processing, and COI gene ampli cation Referring to Chambers et al. [31], microsatellite markers were isolated from the whole genome of Ae. albopictus by magnetic-bead enrichment and PCR screening, and all markers were found to be highly polymorphic via denatured polyacrylamide gel electrophoresis (D-PAGE). A set of 11 microsatellite loci was employed for genotyping the 17 Ae. albopictus populations. Detailed microsatellite primer information is listed in Table 1. All PCRs were performed on a T100 Thermal Cycler (Bio-Rad, USA) in a 50 μl reaction system containing 10 ng of DNA, 0.25 U of PrimeSTAR HS DNA Polymerase (10 pmol/µl, TaKaRa, Japan), 6 μM of dNTPs(2.5 mM each, TaKaRa, Japan) and ca. 5 μM of each primer, with the program set to 35 cycles of 95°C for 30 sec, 57°C for 30 sec and 72°C for 1 min followed by a nal elongation at 72°C for 10 min. All products were then checked with 2% agarose gel electrophoresis under UV light and run on a 3730XL DNA Genetic Analyzer (Applied Biosystems, USA).
As described by Bonacum et al. [32], COI sequence polymorphism was investigated among 497 individuals (22-30 individuals per locality). Brie y, DNA ampli cation of a 550-bp fragment of COI was performed on a T100 Thermal Cycler (Bio-Rad, USA) with the following primer set: 5'-GGAGGATTTGGAAATTGATTAGTTC-3' (F-COI) and 5'-CCCGGTAAAATTAAAATATAAACTTC-3' (R-COI) in a 50 µl reaction mix containing 10 µl PCR Reaction Buffer (TaKaRa, Japan), 4 µl of dNTPs (2.5 mM each, TaKaRa, Japan), 1 µl of primers (10 pmol/µl, TaKaRa, Japan), and 0.5 µl PrimeSTAR HS DNA Polymerase (TaKaRa, Japan). The PCR ampli cation program was set as follows: pre-denaturation at 94°C for 3 min, followed by 35 cycles of denaturation at 94°C for 30 sec, annealing at 54°C for 45 sec, and elongation at 72°C for 1 min, with a nal elongation at 72°C for 10 min. All the PCR products were detected and separated by 2% agarose gel electrophoresis. The target fragments for COI were then cut from the gel under UV light and puri ed with the GenElute™ PCR Clean-Up Kit (NA1020, Sigma-Aldrich, USA). Each puri ed PCR product was then cloned into the pCR™2.1 Vector with a TA Cloning™ Kit (K202040, Invitrogen, USA) and selected by bacteria liquid PCR with the T7 promoter primers. Finally, at least 20 positive clones for each PCR product were sequenced on both strands using an ABI 3730XL automatic sequencer (Applied Biosystems, USA).

Population structure analyses, phylogenetic genotyping and migration analyses
The microsatellite data were processed with GeneMapper v.4.0 (Applied Biosystems, USA). All markers were tested for polymorphism by determining their polymorphism information content (PIC) values via PIC-Calc 0.6 [40], and the null allele frequency of each locus was assessed with Micro-Checker version 2.2.3 [41]. Allele frequency indices, including na and ne, were assessed with Cervus version 3.0.7 [42]. The Ho, He, and FIS values of all populations were calculated with Arlequin version 3.5.2.2 [43]. Departure from Hardy-Weinberg equilibrium (HWE) and heterozygosity de ciency were assessed via Bottleneck version 1.2.02 [44]. Additionally, STRUCTURE version 2.3.4 [45] and the ΔK method of Evanno et al. [46] were employed to calculate the optimal K value; the parameters were set as follows: (1) K ranged from 1 to 20 with ten iterations for each K value; (2) an admixture model was chosen with a Markov chain Monte Carlo algorithm run for 100,000 iterations and 1,000,000 repetitions. The optimal K value was calculated via Structure Harvester: http://taylor0.biology.ucla.edu/structureHarvester/ and depicted with Distribut version 1.1. The R packages "Adegenet 2.1.3" [47] and "Genpop 4.7.5" [48] were used to illustrate the population structure and its relationship with genetic distance, respectively; meanwhile, the genetic variation and F ST value for each population were evaluated with the AMOVA test via Arlequin version 3.5.2.2 [43].
The haplotypes of all the populations were screened by DnaSP version 6.0 [49], and the genetic relationships among all the haplotypes were displayed as a TCS network constructed by Network 10.0.0.0 [50,51]. Moreover, BEAST version 1.8.4 [52]and the R package "pheatmap" [53] were used to build a phylogenetic tree of all the haplotypes and examine the distribution of all the haplotypes among the different climatic regions, with the best model selected by JModelTest 2.1.10 [54]. The mismatch distribution and Bayesian skyline plot analyses were conducted with Arlequin version 3.5.2.2 and BEAST version 1.8.4, respectively, to investigate the current dispersal incidences of each population.
Concurrently, the possible migration routes were rebuilt via the R package "divMigrate" with the number of bootstrap replicates set to 3, the alpha value set to 0.05, the Nm method chosen for the migration statistic, and the lter threshold value set to 0.25 [55]. Finally, the relationships between the genetic indices and environmental factors were described via principal component analysis (PCA) and multiple factor analysis (MFA) with the R package "FactoMineR" [56].

Microsatellite marker isolation and assessment
In the present study, a total of 44 pairs of microsatellite markers were isolated from the whole genome of Ae. albopictus, 11 pairs of which were highly polymorphic and were therefore chosen for microsatellite genotyping analysis ( Table 1). The allele number of each locus ranged from 10 to 33, with a mean of 17.545 alleles per locus. The PIC values of each locus ranged from 0.334 to 0.925, with a mean value of 0.713. According to the de nition of PIC values by Allah et al. [57], nearly all the markers selected were highly informative (PIC value>0.5); the exception was BW-P18, which had a value of 0.357 and was considered to be an informative marker. The Micro-Checker results suggested that null alleles were present at all loci, with frequencies ranging from 0.064 to 0.157 (Additional le 8: Table S5). The linkage disequilibrium (LD) test showed that 302 pairs of loci of a total of 1870 (16.15%) across all locations were in signi cant LD after Bonferroni correction, while no consistency was found among them (Additional File 5: Figure S1).

Genetic diversity and variation
According to the results of the microsatellite data analysis, the observed number of alleles (n a ) in each Ae. albopictus population was very high, and the mean n a value in each temperature region ranged from 6.909 to 8.091 without a signi cant difference. In contrast, the effective number of alleles (n e ) ranged from 3.501 to 4.525, and the n e value increased from the temperate region (3.876) to the tropical region (4.144). The mean value of observed heterozygosity (Ho) for all the climatic regions was ca. 0.557, which was signi cantly lower than the expected heterozygosity (He) (ca. 0.684). The F IS value of each temperature region ranged from 0.266 to 0.359, and the Ae. albopictus populations signi cantly departed from HWE, except in two locations, SHJD and KZXZ (subtropical region, Table 2). Heterozygosity tests of all 17 Ae. albopictus populations based on the stepwise mutation model (SMM) revealed that nearly all the populations from the temperate and subtropical regions displayed signi cant population expansion with p-values < 0.05, while no signi cance was observed among all the tests for the Ae. albopictus populations from the tropical region after Bonferroni correction (Additional le 2: Table S2).
In all, a total of 25 COI haplotypes were observed among 497 Ae. albopictus individuals (GenBank ID: MN651301-MN651325) based on an analysis of the COI sequences, and the haplotype indices changed dramatically across Ae. albopictus populations from different climatic regions (Additional le 3: Table  S3). The haplotype diversity (Hd) ranged from 0.074 (ZGND, temperate region) to 0.750 (JYJB, tropical region), with the nucleotide diversity (π) ranging from 0.014 (BJLG, temperate region) to 0.190 (JYJB, tropical region). The average number of nucleotide differences (k) ranged from 0.067 (BJLG, temperate region) to 0.929 (JYJB, tropical region), and the number of polymorphic sites ranged from 1 to 5 across all the populations. Overall, the tropical Ae. albopictus populations showed the highest diversity, with mean values of Hd, π, and k reaching 0.657, 0.156, and 0.763, respectively, whereas the temperate populations showed the lowest diversity.

Population structure and differentiation based on microsatellite analysis
In the present study, according to the microsatellite analysis results, all the Ae. albopictus populations were adequately allocated to two clades with signi cant genetic differences, and the best K value, as determined via the ΔK method of Evanno et al., was also equal to two (Fig. 2a). Combined with the STRUCTURE bar plot analysis, the Bayesian clustering analysis showed that all the Ae. albopictus populations were adequately allocated to two clades with certain locations from the subtropical and temperate regions genetically isolated from the other locations. (Fig. 2b&c). Moreover, a total of 86.4% of the variation was explained by 50 PCs in the discriminant analysis of principal components (DAPC) analysis; the results revealed two genetically isolated groups, and there was no clear relationship between the Ae. albopictus population structure and its distribution across climatic regions (Fig.2d).
The AMOVA results are presented in Table 3. Of the total genetic variation partitioned, 31.40% could be attributed to differences among individuals within populations and 63.04% to differences within individuals (F IS =0.33253, F IT =0.36962, and all p<0.00001). Meanwhile, the pairwise F ST values between populations ranged from 0.008 (ZGND, temperate region and NNXD, subtropical region) to 0.141 (SHJD, subtropical region and JKCH, tropical region). Genetic differentiation was signi cant between all the sampled populations after Bonferroni correction (p < 0.05) except for ve pairs of F ST values among the subtropical populations NNXZ, NJTH, NJDX, KZXZ, and temperate populations ZGND, BHBG, QDDX, and HBSD (Additional le 4: Table S4). In contrast to the high individual variation, a slightly signi cant positive correlation was observed among only the temperate populations via isolation by distance (IBD) analysis (R 2 = 0.6614, p = 0.048), while no such evidence was observed in the other two regions (Additional le 6: Figure S2).

Haplotype network and phylogenetic analysis based on COI sequences
Three major haplotype clades, distributed across the three central climatic regions and closely related to each other, were reconstructed via the TCS network with 497 COI sequences. Haplotype (H1) of Clade was the most frequent haplotype and was distributed from the tropical to the temperate region with an increasing trend. Nearly all the other haplotypes derived from H1 with one or two mutations, with Clade mainly distributed in the tropical region and Clade mainly distributed in the southern subtropical region (Fig. 3b).
A phylogenetic tree of all 25 haplotype sequences, combined with a heatmap analysis, demonstrated that the 25 haplotypes were divided into three major well-supported clades. As expected, Clade was separated from Clade and Clade with 100% bootstrap support and included haplotypes H9, H10, H11, H12, H19, and H20, which were observed in only the tropical region. In comparison, Clade diverged from Clade with lower bootstrap support (70.74%) and included ve haplotypes (H2~H3 and H5~H7) that were distributed in only the southern subtropical region. Clade III was an admixture group containing all the remaining 14 haplotypes, two of which were observed in the tropical region and four in the temperature region (Fig. 3a).

Migration and correlation analyses between genetic indices and environmental factors
As illustrated in Fig. 4, all ve genetic indices and two environmental factors (i.e., temperature and rainfall) contributed equivalently to the rst axis of the PCA, explaining up to 95.2% of the variation; the exception was the environmental factor latitude, which contributed more to the second axis, with a proportion of 17.5%. Combining these results with hierarchical clustering via MFA, all 17 Ae. albopictus populations were clustered into three groups, of which cluster II and cluster III were closely related to each other. When environmental factors were regarded as the major in uencing factors, the molecular diversity indices of the populations in the tropical region were signi cantly higher than those in the other regions. These results also demonstrated that environmental factors, especially temperature and rainfall, may be the leading causes of differences in genetic diversity among different climatic regions.
Migration patterns were assessed using divMigrate networks representing all 17 Ae. albopictus populations. As expected, Ae. albopictus was observed to have migrated frequently between the tropical and temperate areas in both directions (Fig. 1). A total of four major migration trends were observed among the different climatic regions with high gene ow (Nm>0.4). Two routes were found to have originated from the tropical and subtropical areas and reached Hebei and Beijing in the temperate area, while another two arrived in Guangdong and Guangxi from the southern subtropical area. Compared with the south to north routes, the latter migration routes were substantially larger. Meanwhile, the Tajima's D and Fu's Fs values were negative, with no statistically signi cant p-values except in population JKCH (tropical region, Additional le 3: Table S3). Based on the COI sequences, the mismatch analysis results showed that the Harpending raggedness indices for all three haplotype clades were relatively low (ranging from 0.1054 to 0.1812, p>0.05), and unimodal mismatch distributions were observed among the different Ae. albopictus populations (Additional le 7: Figure S3).

Discussion
Due to their high mutation rates, codominant expression, and universal distribution throughout the eukaryotic genome, microsatellite loci have been widely used to evaluate the genetic variation and population structure of mosquito, especially for those without fully annotated genomes [29,31]. Although different sets of microsatellite loci have been identi ed by previous studies [28,29,58], the absence of PCR ampli cation products of mosquito samples from certain temperature zones has resulted in insu cient data. Therefore, in this study, we developed new primers for our samples to ensure the integrity of the experimental data. However, the potential impact of null alleles, such as reducing population genetic diversity and increasing genetic differentiation among populations [59][60][61], must be taken into consideration when microsatellite markers are used, because null alleles are more frequent in arthropods than in other species [62][63][64]. In the present study, even though nearly all the loci were considered highly informative (PIC > 0.5), null allele was still observed at all loci with frequencies ranging from 0.02 to 0.158, with an average of 0.078; however, frequency values smaller than 0.2 are considered to have no signi cant effect on the accuracy of data analysis by many studies [59,65]. Mitochondrial genes are also good markers for studying population diversity and are widely used in many studies. To improve the accuracy of the ndings in this study, the mitochondrial COI gene was also employed, and the genetic variation of all the individuals was investigated using both COI and microsatellite data.
Consistent with the studies of Zhong et al. [66] and Wei et al. [67], which were based on only mitochondrial and microsatellite diversity, respectively, the tropical Ae. albopictus populations showed higher diversity than the populations of the other two regions in China. An environment more conducive to survival and continuous dispersion may be the best explanation for this phenomenon. According to studies by Wei et al. and others [67][68][69][70], the hot and humid climate of Southern China is more suitable for the breeding and development of Ae. albopictus than a temperate climate, which has relatively low temperatures and dry air. Notably, during the winter, freezing temperatures bring about diapause in Ae. albopictus for the whole winter period, resulting in lower allele richness and population diversity in temperate Ae. albopictus populations [62][63][64]. Meanwhile, Southern China shares borders with many countries in Southeast Asia, where Ae. albopictus originates geographically, including Laos, the Philippines, and Myanmar [20]. Frequent border trades and resident travel among these areas result in the continuous dispersion of Ae. albopictus, enriching the diversity of the local population; this result was con rmed by the results of the bottleneck analysis, which showed that all the populations from temperate and subtropical regions had undergone bottleneck effects, while the populations of the tropical region had not. As suggested by Takezaki et al [71], an He value between 0.5 and 0.8 indicates high genetic polymorphism in the population. In the present study, the He values of all the sampling sites ranged from 0.606~0.733, which indicates high polymorphism. Simultaneously, compared with those in the tropical region, the mean Ho values of the temperate and subtropical regions were signi cantly lower than the mean He values, and nearly all the populations signi cantly departed from HWE. This phenomenon may be closely related to the application of mosquito control programs in these areas [72], according to recent insecticide-resistance studies of Ae. albopictus [37,[73][74][75][76], as compared with temperate regions, tropical and subtropical regions have denser populations and more developed agricultural industries; thus, the overuse of insecticides within these areas will tend to cause insecticide resistance in Ae. albopictus within these regions, resulting in differences in genetic variability among different regions. The continuous dispersion of Ae. albopictus from Southeast Asian countries to the tropical region of China also augments the local mosquito population.
Corresponding to the results of Kotsakiozi and Manni et al. [29,77], which showed no genetic differentiation within the native Asian range, our population differentiation analysis results showed that even though nearly all the pairwise F ST values between populations were signi cant, the values were not high, which indicates potential anthropogenic exchange between the populations. A study by Schmidt et al. [78] of the genetic structure of Ae. albopictus collected from 12 sample sites in Guangzhou also found that human transportation networks, particularly shipping terminals, in uenced the genetic structure of Ae. albopictus populations. With the development of China, the rapid expansion of high-speed rail, aircraft routes, and highways between different temperature zones has not only facilitated trade among different cities but also accelerated the spread of mosquitoes. Accordingly, only two genetically separated clades were observed among these populations, and molecular variations within populations and individuals were found to contribute to the genetic differentiation between these populations. As a global invasive species, Ae. albopictus has developed several strategies to cope with a broad range of temperatures and adapt to local thermal conditions [9,68], which helps this mosquito disperse and colonize different locations successfully. Here, the migration analysis results displayed four major migration trends with high gene ow (Nm>0.4) among the different climatic regions, corresponding with China's main trade and tourism routes. All these results revealed that human-aided dispersion might be the main reason for the similarity of populations across the different regions [22,[67][68][69][70][71][72][73][74][75][76][77][78]. This hypothesis was also veri ed via IBD analysis, which indicated that long distances were not signi cantly associated with individual genetic variation.
Corresponding to the microsatellite analysis results, a total of 25 haplotypes were detected among the 497 COI sequences in the present study, while nearly 44% (11 out of 25) of these haplotypes were observed in only the tropical and subtropical regions. The p-value of the Mantle test for the tropical region was signi cant. The same result for Chinese Ae. albopictus populations was also observed by both Guo et al. and Zhang et al. [79,80]. Guo et al. identi ed three predominant haplotypes in Ae. albopictus populations: H01 (21.5%) from Guangdong Province (subtropical region), H19 (22.1%) from Yunnan and Hainan provinces (tropical region), and H30 (10%) from Hainan Province, while Zhang et al. detected 25 haplotypes, including 10 shared haplotypes and 15 private haplotypes. This means that continuous dispersal from neighboring Southeast Asian countries may be the reason for the higher diversity of the Ae. albopictus populations in the tropical and subtropical regions and that the origins of the Ae. albopictus populations of these two regions may be different. Furthermore, the mismatch analysis results revealed recent population dispersal among the different climatic regions, which is the main reason for the universal distribution of the remaining haplotypes. As an important and rapidly adapting insect, Ae. albopictus has biological characteristics that are strongly in uenced by environmental conditions [9,10]. In tropical and subtropical regions, the climate is hot and wet, making it more suitable for the survival of Ae. albopictus; when this mosquito migrates from a tropical region to a temperate region, it needs time to adapt to the new environment, especially to survive during the winter period, which causes diapause in Ae. albopictus [4]. According to a study by Roiz et al. about the effects of environmental factors on Ae. albopictus, the activity of host-seeking females was positively affected by temperature and rainfall [81]. The PCA in the present study examined genetic indices and environmental factors (such as temperature, rainfall and altitude) and revealed that temperature and rainfall may be the leading causes of differences in genetic diversity among Ae. albopictus populations from different climatic regions.
Overall, strong dispersal patterns were observed in Southern China, and frequent dispersals have occurred from the tropical and subtropical regions to the temperate region. Monitoring the population structure and potential migration routes of Ae. albopictus in different climatic regions will help in the selection of appropriate mosquito control strategies, and since dengue fever continues to break out in Southern China every year, monitoring the potential migration routes of Ae. albopictus can also predict the potential transmission routes of this virus and help prevent its further spread.

Conclusion
The present study systematically evaluated the genetic variation, population structure, and haplotypebased phylogenetic relationships of Ae. albopictus populations across different climatic regions of China. All 17 Ae. albopictus populations were genetically structured in accordance with their locality or region of origin, and three major haplotype clusters were observed via phylogenetic analysis of COI. These results suggest different evolutionary histories in changing environments, especially differences in temperature and rainfall. Meanwhile, four major migration trends were observed among the different climatic regions with high gene ow; these trends contribute to the similarity of the Ae. albopictus populations.    Three climatic regions were marked with different colors, four major destinations were marked with stars of different colors, and migration routes were marked with different types of arrows.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.