Genetic diversity and population structure of Aedes aegypti after massive vector control for dengue fever prevention in Yunnan border areas

Dengue fever is a mosquito-borne disease caused by the dengue virus. Aedes aegypti (Ae. Aegypti) is considered the primary vector of Dengue virus transmission in Yunnan Province, China. With increased urbanization, Ae. aegypti populations have significantly increased over the last 20 years. Despite all the efforts that were made for controlling the virus transmission, especially on border areas between Yunnan and Laos, Vietnam, and Myanmar (dengue-endemic areas), the epidemic has not yet been eradicated. Thus, further understanding of the genetic diversity, population structure, and invasive strategies of Ae. aegypti populations in the border areas was vital to uncover the vector invasion and distribution dynamic, and essential for controlling the infection. In this study, we analyzed genetic diversity and population structure of eight adult Ae. Aegypti populations collected along the border areas of Yunnan Province in 2017 and 2018. Nine nuclear microsatellite loci and mitochondrial DNA (mtDNA) sequences were used to achieve a better understanding of the genetic diversity and population structure. One hundred and fourteen alleles were found in total. The polymorphic information content value, together with the expected heterozygosity (He) and observed heterozygosity (Ho) values showed high genetic diversity in all mosquito populations. The clustering analysis based on Bayesian algorithm, the UPGMA and DAPC analysis revealed that all the eight Ae. aegypti populations can be divided into three genetic groups. Based on the mtDNA results, all Ae. aegypti individuals were divided into 11 haplotypes. The Ae. aegypti populations in the border areas of Yunnan Province presented with high genetic diversity, which might be ascribed to the continuous incursion of Ae. aegypti.


Results
Microsatellite genetic diversity. The polymorphic information content (PIC) value is one of the indicators used to measure allele richness of genes. A total of 114 alleles were found for the nine genetic markers; all alleles were sequenced. As shown in Table 1, the microsatellite locus SQM 6 had the highest number (20) of alleles, while marker SQM 1 had the lowest number (7) of alleles. The PIC values were high, ranging from 0.392 to 0.886, and the average value was 0.672, which indicated that sites were highly polymorphic and can reflect the genetic characteristics of all Ae. aegypti populations 24 . The microsatellite locus SQM 6 had the highest PIC value (0.886), whereas the locus SQM 1 the lowest PIC value (0.392) ( Table 1).
Normally, the genetic diversity of the mosquito population is positively related to the expected heterozygosity (He) and observed heterozygosity (Ho) value of the population. He and Ho value of all Ae. aegypti populations ranged from 0.385 to 0.605, which suggests that the genetic diversity of all mosquito populations is relatively high (Fig. 1). Except for CY, the Ho value in other regions is lower than the He value, which indicates that there may be inbreeding within the species. In addition, CY may have a large amount of individual external supplements or may experience bottlenecks.
Scientific RepoRtS | (2020) 10:12731 | https://doi.org/10.1038/s41598-020-69668-7 www.nature.com/scientificreports/ Nearly all the F IS values in all Ae. aegypti populations, except population YJ, were positive, ranging from 0.05882 to 0.24657 (Table 2), which further indicated that these populations contained different degrees of inbreeding and Heterozygote deficiency that also may be the reason for the deviation of all populations from HWE.
Microsatellite genetic structure. The F IT value over all populations was 0.309, which indicated that there were significant population differences among the individuals screened. The AMOVA results indicated that the largest proportion of genetic variation in Ae. aegypti population existed in individuals and individuals within populations, accounting for 69.14% and 15.69% of the variation, respectively (Table 3). Although the sampling site is an important factor (P < 0.0001), the ratio was relatively low ( Table 3). The Bottleneck effect analysis of all populations revealed that nearly all populations, except population ML and RL, are mutation-drift equilibrium (Table 4).   www.nature.com/scientificreports/ In order to analyze the specific genetic structure of all eight Ae. aegypti populations, based on the Bayesian algorithm, a clustering analysis was carried out (Fig. 2). Combining with The UPGMA and DAPC analyses, revealed that all eight Ae. aegypti populations can be divided into three genetic groups of the populations from Lincang city, Xishuangbanna prefecture and Dehong prefecture were genetically correlated, except for RL population, which was highly related to the population Xishuangbanna prefecture (Figs. 3, 4).
The IBD analysis displayed that the genetic distance of all eight Ae. aegypti populations were positively related to geographic distance, which meant the geographical isolation was the primary cause of genetic diversity of Ae. aegypti (Fig. 5). While the Ae. aegypti can only move hundreds of meters around their larval habitats, which suggests that the transmission of Ae. aegypti in Yunnan Province does not depend on its own activities, but on other factors, such as human activities.
The pairwise F ST values of Ae. aegypti ranged from 0.061 to 0.220 (Table 5), showing significant genetic differences. All P values were significant (P < 0.05) after Bonferroni corrections were applied.
Haplotype networks and diversity. Based on the mitochondrial COI, ND4, and ND5 areas of Ae.
aegypti, all Ae. aegypti individuals from nine populations were divided into 11 haplotypes, among which H1, 2, 3, 7, and 8 were the main haplotypes. The distribution of these haplotypes in a specific population is shown in Fig. 6. All the new sequences generated in this study are available from GenBank (Table 6). Among the eleven haplotypes, H1, 2, 3, 7, and 8 were the main haplotypes in all populations. The haplotype H1 had the most individuals, mostly from MD, MH, CY, LC, and RL, while H3 and H8 were mainly distributed at Xishuangbanna prefecture and H2 and H7 were only distributed at YJ and JH, respectively (Fig. 6). The vast majority of individuals have shared haplotypes, and sampling sites had individuals that belonged to haplotypes. Neutral test and mismatch analysis based on the mitochondrial genes ND4 and ND5 of Aedes aegypti showed that the distribution of nucleotide mismatches in this population had a single peak structure, indicating that the population of Ae. aegypti had experienced at least one significant population expansion (Fig. 7).

Discussion
The information on invasion and spread of mosquito vectors is essential for understanding vector-borne disease outbreak, transmission dynamics among human populations and implementing effective mosquito control programs 25 . These are all important factors influencing the mosquito population dynamics, genetic structure patterns, and pathogen transfer through vector populations 26 . Ae. aegypti is the most important epidemic vector that can cause DF and dengue hemorrhagic fever (DHF) in human and is mainly distributed in southeast China. The most suitable habitats of Ae. aegypti include Hainan, Guangdong, Guangxi, the western and southern border areas of Yunnan, and parts of the southern Guizhou region 27 . Yet, due to climate changes and increased urbanization, a significant northward shift occurred in the northern Chinese region over recent years 28 .
Ae. aegypti is an invasive species and potential vector of disease agents in China, which has a significant impact on public health. Ae. aegypti-associated infection was first reported in Yunnan Province (Jiegao Port, near Ruili City, Dehong prefecture) in 2002 29 . In 2009, the Ae. aegypti was detected for the first time in Guanlei Port, Mengla city, Xishuangbanna prefecture 30 , and later on (2014) in Lincang, Mengding county 31 . The distribution range and abundance of Ae. aegypti species have significantly increased, and were established in at least eight cities in Yunnan Province. Therefore, the monitoring of Ae. aegypti species is essential for preventing and controlling vector-borne infectious diseases. In our study, all the Ae. Aegypti samples were collected from eight sampling places in three prefectures of Yunnan Province (Xishuangbanna prefecture, Lincang city, and Dehong prefecture). The DF cases in Yunnan Province mainly originated from these prefectures.
Our population genetics analyses of the populations from Yunnan border area that were based on two types of genetic markers (Microsatellites and mtDNA) revealed the genetic structure and the population distribution within this region. The PIC value, He and Ho value were important parameters for measuring the genetic Table 3. Hierarchical analysis (AMOVA) of the genetic variation in the Ae. aegypti samples. www.nature.com/scientificreports/ diversity of a population; the higher the value, the more complex the population structure is. Our results revealed that Ae. aegypti species in the Yunnan border region had a great allelic variation. The Ae. aegypti mosquitoes may easily transmit the virus to humans and usually find shelter in indoor habitats. Their flight range is limited, www.nature.com/scientificreports/ which means they can only move hundreds of meters around their larval habitats 25,32 . The relatively high genetic diversity of all mosquito populations is most likely caused by invasion events and human activities 33 . The results of IBD analysis support this conclusion that the dispersal of Ae. aegypti species is aided by human activities and transportation in Yunnan Province. Bayesian algorithm-based population analysis showed that all Ae. aegypti populations could be divided into three genetic groups. The first group had four populations in the Xishuangbanna prefecture (JH, MH, and ML) and RL city, which might be related to the close tourism and commercial trade exchanges between these two regions. The second group represented two populations from Lincang City. The third group was composed of LC and YJ. Contrary to other six populations, YJ population is closely related with LC population, which is significantly different from the UPGMA result. The differences may come from the sample size and range.
Except for Yingjiang, the inbreeding coefficient (F IS ) values of eight other Ae. aegypti populations were positive. Combined with the UPGMA and DAPC analysis, the results indicated that there may have a recent invasion and colonization of the Ae. aegypti in YJ city. Due to the limited flight range, this phenomenon is common for Ae. aegypti population on a small spatial scale 33 . In 2016, the Xishuangbanna Prefecture established a "Spring  www.nature.com/scientificreports/ Patriotic Health Movement" with the scope to provide the integrated control for infectious disease vectors. This may explain the bottleneck effect observed in ML and RL. MtDNA markers have been widely used to evaluate the genetic diversity of Ae. aegypti populations 34,35 . In our study, the degree of polymorphism found in the COI and ND4 sequences were relatively high (eight populations were divided into eleven haplotypes). The H1, which was the dominant haplotype, was found in five places. The analysis of all mosquito samples from two localities in Lincang City (MD and CY) showed only one haplotype (H1) for each gene. Dehong Prefecture is close to Myanmar border, and the intensive personnel activities have led to a large number of invasion events. The abundant waters and commercial activities in Xishuangbanna prefecture have also contributed to many invasion events. This idea was supported by the high levels of polymorphism detected in Xishuangbanna and Dehong prefectures (six haplotypes and seven haplotypes, respectively), which may be the main entry points of Ae. aegypti in Yunnan Province. The H2 haplotype was only distributed at YJ independent from other regions. Combined with the negative F IS value of population YJ, the Ae. aegypti species likely invaded Yunnan Province from this region over recent years.
Our research shows that in Dehong and Xishuangbanna prefectures, Ae. aegypti population invade these areas because of the continuous tourist and business activities. Inspection and quarantine need to be strengthened at the border ports and further investigation and research on mosquito vectors should be carried out. The government needs to designate effective prevention and control measures, strengthen environmental governance in the border areas and implement mosquito control measures.  www.nature.com/scientificreports/ conclusion The nuclear microsatellite markers and mtDNA sequences (COI, ND4, and ND5) were used to uncover the population genetics of the Ae. aegypti in the border area of Yunnan Province. Although several attempts have been made by the government of Yunnan Province to control the mosquito vectors, the Ae. aegypti populations in this region showed high genetic diversity and genetic structure due to the continuous invasion, and increased urbanization. Our research confirms that, over recent years, a significant Ae. aegypti invasive event occurred in YJ City; and that the Xishuangbanna and Dehong prefectures were important areas for the Ae. aegypti invasion. In summary, our results suggest that the control of Ae. aegypti in Yunnan Province is still a demanding task that needs to be taken seriously. Thus, monitoring of suspected cases of DF and the vectors should be enhanced.

Materials and methods
Mosquito sampling and DNA isolation. All the adult Ae. aegypti samples were collected from following eight locations along the border area of Yunnan Province between May 2017 and September 2018 (Fig. 8,  Table 7). Each collection site covered an area of approximately 500 m in diameter. According to the Surveillance Methods for Vector Density-Mosquito (GB/T 23797-2009), a hand-held aspirator was used to collect the adult mosquitoes (intercepted before biting). All the samples were identified through the analysis of morphological characteristics in the wild field 36 and preserved in 100% ethanol at 4 °C for the isolation of genomic DNA 37 .
According to the standard DNA extraction procedure, genomic DNA was isolated from individual mosquito sample with the TaKaRa Mini-BEST Universal Genomic DNA Extraction Kit (Takara, Dalian, China); the quality and quantity of extracted DNA were analyzed using NANOdrop1000, after which samples were stored at − 20 °C until further analysis.  www.nature.com/scientificreports/ PCR amplification and microsatellite genotyping. Nine microsatellite polymorphic loci were screened from 58 loci, which were described in previous studies by denaturing polyacrylamide gel electrophoresis 38,39 . The primer sequences and information are summarized in Table 8 (Table 8) for 45 s) and 72 °C for 45 s. The final extension was performed at 72 °C for 10 min. All PCR amplification products were verified by electrophoresis of 3 μL on a 1.5% agarose gel. The formamide was mixed with LIZ 500-labeled size standard using a ratio of 100:1, and 15 μL mixture was added into the sample plate. The PCR amplification products were diluted at 1:10, 1 μL was added into the reaction, and then run on an ABI3730XL (Applied Biosystems, Foster City, USA) capillary sequencer. All microsatellite alleles were evaluated using GeneMapper software (Applied Biosystems) 14 .

Microsatellite data analysis. Genetic diversity.
The PIC values of all nine loci were calculated with PIC-Calc 0.6 14 . The genetic diversity of all Ae. Aegypti populations were characterized by expected heterozygosity (He) and observed heterozygosity (Ho), using POPGENE version 1.32. The F IS value of each mosquito population was also calculated. The statistical significance test was performed with the exact tests available in POPGENE 40 .
Genetic structure. The genetic variation was tested by the AMOVA test with Arlequin (version 3.5.2.2) for the interpretation of genetic variability and structure among different locations, mosquito populations. The AMOVA was evaluated at four different hierarchical levels: (1) all samples (non-grouped) were analyzed as a single group to test the overall genetic differences between samples; (2) the samples in one region were analyzed as a unique group; (3) the interregional populations were analyzed as an unique group; (4) the individual sample within population was analyzed as an unique group. Based on the stepwise mutation model (SMM), the recent genetic bottleneck in each mosquito populations was calculated by the software BOTTLENECK 1.2.02. The data were analyzed with the recommended settings: an index statistic closer to 1 indicates that the population is in a stable state, while a very low value indicates that the population has experienced a genetic bottleneck in the past 41 . The sign test implemented in the software was used to test for significant heterozygosis excess. The isolation by distance (IBD) was estimated with Mantel's test in R, using the correlation between genetic distance and geographic distances by the regression of pairwise FST/(1 − FST) on the natural logarithm (Ln) of straight-line geographic distance. www.nature.com/scientificreports/ For the determination of real genetic clusters (K) within all mosquito samples, a Bayesian clustering algorithm-based software STRU CTU RE 2.2 was employed. All mosquitos were divided into different populations represented by a specific number (K = 8), under the assumption of Hardy-Weinberg equilibrium and linkage equilibrium 42 . The software parameters were set as follows: the assumed populations ranged from 1 to 8, and the calculation model was set as admixture ancestry and independent allele frequency models 100,000 burn-in steps followed by 1,000,000 MCMC replicates, and each population was calculated for 10 runs. The optimum K value was estimated with Evanno's △k method based on the second-order rate of change in the log probability of the △k among 10 runs of each assumed K 43 , and all the results were uploaded to a web-based utility Harvest for the calculation of the optimum K value (https ://taylo r0.biolo gy.ucla.edu/struc t_harve st/). Furthermore, the     : AAT CGT GAC GCG TCT TTT G  CT10(TT)CT  233-239  R: TAA CTG CAT CGA GGG AAA CC   SQM 2  F: CAA ACA ACG AAC TGC TCA CG  GA15  157-183  R: TCG CAA TTT CAA CAG GTA GG   SQM 3  F: ATT GGC GTG AGA ACA TTT TG  CAT7  156-