PCR amplification of the Vgsc gene in these Cx. p. pallens specimens yielded a 521 bp PCR fragment. The sequence data for this gene were analyzed to assess kdr allele frequency distributions, kdr codon evolution, polymorphisms in the neighboring intron downstream of the kdr codon, and population genetic diversity. PCR amplification of the mitochondrial COI gene yielded a 651 bp PCR fragment. As the mitochondrial genome is maternally inherited and haploid, drift towards different haplotype frequencies will occur more rapidly within isolated populations, resulting in differentiation roughly twice that observed for nuclear markers. Accordingly, these mtDNA sequences were analyzed to assess population variations, population expansion, spatial population structure, genetic differentiation, and patterns of gene flow.
Kdr allele frequency distributions
Three alleles (wild-type L1014, L1014F, and L1014S) and eight genotypes (TCA/TCA, TTC/TTC, TTA/TTA, TTT/TTT, TTA/TCA, TTA/TTT, TTT/TCA, and TTC/TTT) were detected in this study population, corresponding to two non-synonymous codon 1014 mutations (Table 2). The kdr mutant allele frequencies were found to vary with geographic location. Both mutant alleles (L1014F and L1014S) were observed in mosquitos in the Northwest plain of Shandong Province, with the L1014F allele being predominant. Both mutant alleles were also found in the eastern hilly area of Shandong Province. However, only the L1014F mutant allele was detected in the southern mountainous area and DY at a low mutational frequency (Fig. 1).
Kdr haplotype diversity
In total, 28 polymorphic sites were identified when analyzing these samples. Genetic diversity indices and neutrality tests (Fu’s Fs and Tajima’s D) were performed using the neighboring intron downstream of the kdr allele in Cx. p. pallens, with 520 kdr intron sequences being included in these analyses. As these polymorphic sites exhibited heterozygous phenotypes, a phase analysis approach was used for haplotype predictions, leading to the detection of 31 haplotypes from 9 populations (Table 3). Overall respective haplotype diversity (Hd) and nucleotide diversity (Pi) values were 0.664 and 0.0521. The HZ populations exhibited low Hd (0.351) and Pi (0.01261), whereas these values were high in all other analyzed populations. No significant departures from neutrality were detected in any of these populations using Tajima’s D test, whereas Fu’s F test detected significant departures in the QD and ZB tests, with positive F test results indicating low levels of high-frequency populations consistent with a population size reduction and/or balancing selection (Table 3).
Genealogical analyses of kdr mutations
To gain further insight into the evolutionary relationships among the detected kdr mutations, the 31 haplotypes identified above were subjected to further analysis. Two haplotypes were predominant in the mosquito study populations: H01 (52.5%) and H02 (23.8%) (Fig. 2). The other detected haplotypes exhibited only limited geographic distributions or were unique to individual populations. Overall, 18 shared and 13 unique haplotypes were identified in this study, with the former accounting for 58.1% of all haplotypes. The H01 and H02 haplotypes were shared among all study sites, while substantial haplotype sharing was also observed between the LC and DZ populations (including H06, H11, H13, and H15). Many shared haplotypes were also observed among the QD, DY, YT, and RZ populations. In addition to these eastern coastal cities, the QD population exhibited only limited sharing of haplotypes with other regions suggesting that this population is relatively differentiated from other populations, in line with the genetic differentiation coefficient Fst. Other than the H01 and H02 haplotypes, only one other shared haplotype was present in the ZB population. However, Except for the shared haplotypes of all study sites, only one shared haplotype exists in ZB. The geographically non-adjacent YT and LC populations also shared seven haplotypes, with a negative genetic differentiation coefficient (Fst) of -0.00859 between these populations consistent with no appreciable genetic differentiation and high levels of communication between these populations. Given the limited sampling at each study site and associated time constraints, further analyses will be needed to confirm and expand upon these findings.
Network analysis results revealed a complex reticulate pattern and multiple independent mutational events resulting in the observed kdr haplotypes (Fig. 3). The H11 haplotype is largely in the center of this network as it is most closely related to the other haplotypes such that it may be a more evolutionarily ancient haplotype that can undergo mutation through one or more steps to yield the other haplotypes. Haplotypes H01-F and H01-S comprise the majority of samples in all populations and may be derived from ancestral H01-L, H08-L, H03-L, and H12-L haplotypes through one or two mutational steps. Likewise, the H02-F and H02-S haplotypes may be derived from the ancestral H02-L, H08-L, H010-L, and H118-L haplotypes through one or two mutational steps. The H05-F haplotype found in HZ stems from H01-L through one mutational step, while the H26-F haplotype from TY and LC is likely to have arisen from H11-L through three mutational steps. The H16-F and H29-F haplotypes are derived from a shared ancestor lacking the kdr haplotype. Overall, these genealogical analyses suggest that as few as one mutational step may differentiate between resistant and susceptible phenotypes among Cx. p. pallens.
Analyses of mtDNA sequence variation
Sequence polymorphism analyses of 520 COI mtDNA sequences were next conducted, leading to the identification of 8 variable sites and 9 haplotypes were detected (Table 4). The respective overall haplotype diversity (Hd) and nucleotide diversity (Pi) values were 0.306 and 0.00052, respectively. These relatively low values suggest that the individual geographic Cx. p. pallens populations may have been subject to population bottleneck effects. The DZ population exhibited the highest levels of diversity, whereas the DY and ZB populations were the least diverse. The Median-joining method was used to construct a haplotype network graph consisting of a radiative evolutionary center and a star-shaped distribution (Fig 4). In total, 2 and 7 shared and unique haplotypes were identified, respectively. Of these, the H1 haplotype was the most widely distributed primitive haplotype, radiating out to one shared haplotype and multiple unique haplotypes. This star-shaped network phylogeny of this network also supports the notion that these Cx. p. pallens populations expanded following population bottleneck events, but that these populations are still undergoing expansion. The H4 haplotype was shared among the RZ, QD, DZ, and YT populations, while other haplotypes were only present at low frequencies. H2, H3, and H9 were respectively unique to the HZ, LY, and LC populations, while H5-H8 were unique to the DZ population. The observation that these haplotypes were unique to these respective regions may be a consequence of limited sample selection or may be a result of genetic variation among these geographic environments. To better assess the Cx. p. pallens population structure and phylogenetic relationships among these populations, a Neighbor-joining (NJ) tree incorporating these 9 COI sequence-based haplotypes was constructed (Fig 5), revealing that all haplotypes clustered in a single branch.
Genetic differentiation among populations
An analysis of molecular variance (AMOVA) revealed that genetic differences in the COI gene among and within populations respectively accounted for 18.3% and 81.7% of the observed variance, with a fixation index (Fst) value of 0.18293. As such, within-populations variance is the primary source of genetic variability among the Cx. p. pallens populations of Shandong Province. The pairwise Fst and Nm value estimates further highlight genetic differentiation among these populations (Table 5). A range of factors can contribute to the genetic differentiation of particular populations, including transport or other human activities, landscape barriers, and geographical distances among these populations. Fst values for the genetic differentiation index among these different mosquito populations ranged from -0.00802 to 0.46939. The only negative Fst values were those between the DZ and YT populations. No significant differentiation was observed among 14 pairs of populations (Fst<0.05; 38.9%), while moderate differentiation (0.05<Fst<0.15), high differentiation (0.15<Fst<0.25), and very high differentiation (Fst>0.25) were respectively observed among 7 (19.4%), 8 (22.2%), and 7 (19.4%) pairs of populations. The highest Fst value was observed between the DZ-YT and QD-ZB/DY populations (Table 5). Nm values corresponding to the gene flow among these populations were between -31.43 and 134.26. Together, these findings were thus consistent with roughly half of the populations exhibiting high or very high levels of genetic differentiation and low levels of gene flow among these populations. Fixed coefficient Fst values between pairs of populations ranged between -0.00802 and 0.46939 (Table 5), and the highest levels of genetic differentiation were observed when comparing the mosquito populations from QD and either DY or ZB (FST = 0.46939; Nm = 0.28), owing to the presence of the H4 haplotype in QD and its absence in DY/ZB.
Population demography history
In an effort to infer the demographic history of Cx. p. pallens populations in Shandong Province, e Tajima’s D and Fu’s Fs tests were next conducted, with the results and simulated P-values being shown in Table 4. Negative Tajima’s D and Fu’s Fs values were obtained for the LC, DZ, LY, and HZ populations, suggesting high numbers of low-frequency mutations in these populations and that they are actively undergoing expansion. The frequency distributions for pairwise differences between sequences were also assessed to evaluate historical demographic expansions. Mismatch distributions can be used to visually assess historical population dynamics. Two models correspond to expected DNAsp values, with one assuming that the population size remains constant and the other assuming that the population grows or declines. When observations conform to one of these models, it indicates that the population is in line with the expected model-derived values and the population is constant or shrinking, whereas if the values deviate from expected values this indicates the population is expanding with a single-peak curve. The present results indicated that actual observations were consistent with hypothetical model values, thus indicating that the Cx. p. pallens populations in Shandong were in a state of genetic equilibrium.