A total of 220 individuals of A. aurantii were sampled from 11 tea plantations in China in 2020 (Figure 1 and Table 1). The populations from the 11 locations were, respectively, named as P1, P2, and so on. All individuals were preserved in 99.9% ethanol until DNA extraction. Genomic DNA was extracted from each individual using the E.Z.N.A.^® Tissue DNA Kit (Omega, Norcross, GA, United States) following the manufacturer’s protocol and preserved at −20°C.

The cox1 gene was amplified with the universal primers cox1F 5′-ATTCAACCAATCATAAAGATATTGG-3′ and cox1R 5′-TAAACTTCTGGATGTCCAAAAAATCA-3′ (Hebert et al., 2004). The cytb gene was amplified with primers cytbF 5′-GATGATGAAATTTTGGATC-3′ and cytbR 5′-CTAATGCAATAACTCCTCC-3′ (Harry et al., 1998). The conditions for the PCR reactions were as follows: initial denaturation at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 30 s, annealing at 54°C for 45 s, and elongation at 72°C for 60 s, and a final elongation at 72°C for 10 min. The PCR products were separated by electrophoresis in 1.0% agarose gels and purified with the Axygen DNA Gel Extraction Kit (Axygen Biotechnology, Hangzhou, China). The purified PCR fragments were sequenced by an ABI 3730 automated sequencer (Biozeron Co., Ltd., Shanghai, China).

The gene length, conserved sites, variable sites, parsimony informative sites, singleton sites, and the average A + T contents were, respectively, calculated for the cox1 and cytb genes and their concatenated sequences by Mega v6.0 (Tamura et al., 2013). Genetic diversity parameters, including the haplotype number (H), number of polymorphic sites (S), haplotype diversity (H_d), and nucleotide diversity (Pi), were estimated using DnaSP software (Librado and Rozas, 2009). Each haplotype was named as H1, H2, and so on. Molecular variance analysis (AMOVA), pairwise F_st values, neutrality test, and mismatch distribution analysis were performed with the Arlequin v3.5 software (Excoffier and Lischer, 2010). Gene flow analysis was also conducted with the Arlequin v3.5 software to generate the N_m values and examine the gene exchanges between different populations. To test the correlation between geographic distance and genetic difference, the IBD (isolation by distance) analysis based on the Mantel test was performed using IBD v1.53 (Bohonak, 2002). Principal coordinates analysis (PCoA) was conducted from the distance matrix of the 11 populations and eight provincial groups.

Phylogenetic reconstruction using the concatenated sequences was performed with the neighbor-joining (NJ) method hosted by Mega v6.0 (Tamura et al., 2013). In the NJ analysis, the bootstrap method was used with 1,000 bootstrap replications. The Kimura two-parameter model was selected as the substitution model. Both transitions and transversions were included, and the rates among sites were regarded as uniform. Gaps/missing data were pairwise deleted. All computed trees were adjusted and visualized in FigTree v1.4.2. The relationships among the different haplotypes were inferred using a median-joining method in PopART software (Leigh and Bryant, 2015).

The 220 sampled individuals of A. aurantii yielded 196 high-quality cox1 sequences and 201 cytb sequences (GenBank accession numbers MW265715–MW265910 and MW289592–MW289792). Both genes were simultaneously obtained for 177 individuals. The remaining double-peak sequences were excluded from the analyses. For the cox1 genes, the consensus length was 702 bp for all populations (Table 2). The average A + T content for the cox1 genes was 75.5%. The cox1 sequences of the 11 populations included 623 conserved sites, 76 variable sites, 54 parsimony informative sites, and 22 singleton sites. For the cytb genes, the consensus length was 829 bp for all populations (Table 2). The average A + T content for the cox1 genes was 76.9%. The cytb sequences of the 11 populations included 666 conserved sites, 150 variable sites, 114 parsimony informative sites, and 33 singleton sites. For the concatenated sequences of the cox1 and cytb genes of the 177 individuals, the consensus length was 1,507 bp and the average A + T content was 76.2% (Table 2). There were 1,287 conserved sites, 214 variable sites, 160 parsimony informative sites, and 54 singleton sites in the concatenated sequences. These concatenated sequences were used for subsequent genetic and phylogenetic analyses.

A total of 166 haplotypes were identified for the 177 samples of A. aurantii, which had both cox1 and cytb sequences (Supplementary Material). Most of the haplotypes were unique and corresponded to their own population; only three haplotypes (H28, H34, and H63) were shared by different populations. The largest haplotype was H28, which was shared by one individual of P2 from Hunan, two individuals of P7 from Sichuan, and one individual of P10 from Anhui. The entire dataset comprised 92 polymorphic sites, with a high H_d of 0.9991 (±0.0008) and a low Pi of 0.0097 (±0.0005) (Table 3). For each population, the number of polymorphic sites (S) ranged from 37 (P7) to 88 (P11). The H_d was 1.000 for all populations. The Pi values varied from 0.0112 (±0.0008, P2) to 0.0211 (±0.0014, P9).

The AMOVA showed more variation of genetic differentiation within populations than among populations of A. aurantii (Table 4). Over 87.49% variation of genetic differentiation occurred within populations, and only 12.51% of the genetic variation was found among populations. The pairwise F_st values were highly variable, ranging from 0.0229 (between populations P3 and P4) to 0.2984 (between populations P2 and P11). High genetic differentiations were present between populations P2 and P11 (F_st > 0.25, p < 0.05), whereas no significant differentiations were present between other pairwise populations (Table 5). Most of the N_m values from the gene flow analysis were more than 1, except for the negative values between P11 and populations P2, P3, P4, P7, and P8 (Table 5). The highest N_m value was found between P3 and P4.

The IBD analysis revealed a highly significant positive correlation (r = 0.4190, p < 0.01) between the geographic distance and genetic distance of the A. aurantii populations.

The first two coordinates of the PCoA explained 78.78% of the total variation among the 11 populations (Figure 2). The first and second axes explained 52.8 and 25.98% variation, respectively. Most genotypes on the PCoA graph were separately clustered, except for populations 1 and 10, which were partially overlapped. The PCoA at the population level also revealed that P2, P4, P7, and P11 were clearly distinct from the other populations. In the PCoA at the province level, most genotypes from Sichuan were closely located with those from Guizhou, Yunnan, and Anhui; the genotypes from Hunan, Fujian, and Zhejiang were clearly distinct from those of the other provinces.

The neutrality test was conducted based on the concatenated cox1 and cytb sequences (Table 6). Tajima’s D values were positive for most populations except for P11, and all these values were not statistically significant (Table 6). Fu’s F_s values were negative for most populations except for P7, and most of the values were not statistically significant, except for that of P1. For the entire dataset of A. aurantii, Tajima’s D value was positive and not significant, while Fu’s F_s value was negative and significant (Table 6).

The mismatch distribution analyses generated single bell curves for the combined dataset of all A. aurantii populations (Figure 3). In the two curves calculated by different models, the observed data generally corresponded with the simulated mismatch distribution under the spatial expansion model rather than the sudden expansion model.

For the entire dataset of A. aurantii under the sudden expansion model, the goodness-of-fit test computed the sum of squared deviation (SSD) as 0.0040, which was not statistically significant (p = 0.2300); the Harpending’s raggedness index (HRI) was 0.0003, which was also not statistically significant (p = 1.0000). For each population, the SSD and HRI values were all not significant (Table 7).

For the entire dataset of A. aurantii under the spatial expansion model, the goodness-of-fit test computed the SSD as 0.0003, which was not significant (p = 0.8500), and the HRI was 0.0003, which was also not significant (p = 0.9700). For each population, the SSD and HRI values were all not significant (Table 7).

The NJ tree was constructed using the concatenated cox1 and cytb sequences (Figure 4). However, the bootstrap values of most clades were lower than 50%. Individuals of the 11 populations were scattered in the tree and did not cluster into apparent lineages. Each geographic population cannot be divided in the phylogenetic analysis, indicating more variation of genetic differentiation within populations than among populations of A. aurantii. The only exception is P11, most samples of which grouped together in the NJ tree.

The median-joining network was constructed for populations from eight provinces (Figure 5). The 166 haplotypes were widely distributed among the sampling locations and were connected without apparently diverged lineages. Most haplotypes were unshared by different geographic groups (each province), except for haplotypes H28, H34, and H63. H28 was the mostly shared haplotype by three geographic groups from Hunan, Sichuan, and Anhui. H34 was shared by two geographic groups from Sichuan and Anhui. H63 was shared by two geographic groups from Hunan and Sichuan.

The purpose of the current study was to determine the genetic diversity and interrelationship of the Chinese populations of A. aurantii using mitochondrial gene markers. Assessment of the genetic diversity of the species can reveal its adaptability to environmental change, which is important for pest control or biological resource protection (Schmitt and Hewitt, 2004). The A. aurantii populations studied herein were characterized by a high haplotype diversity and a low nucleotide diversity. This suggests that A. aurantii might have undergone expansion from a low effective population size, or ancient population bottleneck or founder effect (Hedgecock et al., 1989; Grant and Bowen, 1998). The star-like structures of the haplotype network diagram also provided evidence for the recent expansion of A. aurantii from a small number of ancestors (Slatkin and Hudson, 1991; Grant and Bowen, 1998).

Different locations in China are represented by different local tea products. Such long-term monocultures of plant hosts could lead to the low nucleotide diversities of each tea aphid population. The ability of A. aurantii to adapt to chemical insecticides is poor, so it is susceptible to regional extinction events. The periodically excessive use of chemical insecticides in tea gardens might cause the total extinction of a regional population, which could lead to a further reduction in the genetic diversity of A. aurantii.

The AMOVA suggested the existence of more variation of genetic differentiation within populations than among populations of A. aurantii. The pairwise F_st results (0.02292–0.29841) suggested that low to moderate degrees of genetic differentiations are present among the A. aurantii populations (Wright, 1990). Meanwhile, most of the N_m values generated by gene flow analysis were higher than 1, further indicating the existence of frequent gene flow events among the A. aurantii populations (Slatkin and Maddison, 1989). Both the F_st and N_m analyses supported the genetically close relationship between P3 of Guangdong and P4 of Fujian. The PCoA also indicated the existence of genetic differentiations in most populations of A. aurantii.

The IBD analysis revealed a highly significant positive correlation (r = 0.4190, p < 0.01) between the geographic distance and genetic distance of the A. aurantii populations. Geographic isolation was supported as the main factor that reduced the gene flow events among the different populations of A. aurantii. Despite the restriction by geographic distance, the strong dispersal ability of A. aurantii resulted in a certain degree of gene exchange between populations, which was supported by the AMOVA and gene flow analysis. The current cosmopolitan distribution of A. aurantii (Piron et al., 2019) supported its possession of a strong dispersal ability, which could have promoted the gene exchange between different populations. The ability of inhabiting a wide range of plants also provided intermediate hosts for the dispersal of A. aurantii. The continuous transportation of countless tea products or other plants between regions could have also accelerated the process of gene exchange.

In both the NJ tree and haplotype network diagram, no strong genetic structure was detected. Each population was widely distributed on the trees and networks, comprising different haplotypes and not forming distinct lineages. This may be explained by the long period of migration-caused gene exchange between the populations of A. aurantii, which reduced the degree of genetic differentiation and made it difficult to establish independently evolved lineages.

For the entire dataset of Chinese A. aurantii, the neutrality test showed insignificantly positive Tajima’s D values and significantly negative Fu’s F_s values, suggesting the probable existence of population expansion events in history (Fu, 1997). For each of the 11 sampled populations, most of the D values were insignificantly positive and most of the F_s values were insignificantly negative, deviating from neutrality. For the entire dataset and for each sampled population, both the tests of sum of squared deviation and Harpending’s raggedness index were not significant, which suggested that the hypothesis of population expansion cannot be rejected (Harpending, 1994). The mismatch distribution of A. aurantii was consistent with the results of the neutrality test, being consistent with the simulated mismatch distributions of the expanded populations (Harpending et al., 1993). The obtained low Harpending’s raggedness index (<0.04) was also regarded as a significant evidence of ancient expansion events (Harpending, 1994), which further supported the hypothesis of population expansion in Chinese A. aurantii populations.

The significant negative values of the neutrality test and the good fitting degree of the mismatch distribution (demographic expansion and spatial expansion models) were the signatures of the population expansion events and supported the hypothesis that Chinese A. aurantii populations have undergone a rapid population expansion during their history (Tajima, 1989; Rogers, 1995; Fu, 1997). The result of the demographic history was consistent with that of haplotype and nucleotide diversities, which also indicated that A. aurantii underwent a rapid population expansion.

According to the above multiple analyses, Chinese A. aurantii populations might have derived from a low effective ancestral population size and then expanded along with their suitable habitats, such as tea gardens or tea factories. The continuous increase of habitats and the transregional transportation of host plants provided favorable conditions for the fast expansion of the tea aphids. The dispersal and spread of A. aurantii into different areas generated different geographic populations, and each population gradually accumulated genetic differences due to geographic isolation. Usually, the long geographical distance could reduce the gene exchange and cause different genetic structures within populations (Taguchi et al., 2015; Zhang et al., 2019). However, the excessive use of chemical insecticides continuously caused regional extinction and loss of the genetic diversity of Chinese A. aurantii populations. Although geographic isolation was considered the main factor that influenced the genetic diversity of Chinese tea aphids, the gradually increased cultivation areas and frequent transportation apparently promoted the gene exchange between the different A. aurantii populations and led to fewer genetic differences. From the aspect of pest control, reduction of population genetic diversity by chemicals or gene exchanges is beneficial to the control of tea aphids.