Population Genetics of the Black Citrus Aphid Aphis aurantii (Hemiptera, Aphididae) in China

The black citrus aphid, Aphis aurantii Boyer de Fonscolombe, 1841, is one of the most destructive pests in commercial tea plantations and gardens in China. In this study, we investigated the population genetic structure of A. aurantii based on the concatenated sequences of two mitochondrial genes, cytochrome c oxidase I (cox1) and cytochrome b (cytb). A total of 166 haplotypes were identified from 177 individuals collected at 11 locations in China. The whole Chinese A. aurantii population showed a low nucleotide diversity (0.00968) and a high population diversity (haplotype diversity; 0.9991). The haplotypes of the 11 local populations were widely distributed in the neighbor-joining phylogenetic tree and haplotype network diagram, whereas no apparent lineages were detected. Gene flow analysis showed gene exchanges among local populations. The pairwise Fst values revealed a certain amount of genetic difference among local populations. Analysis of molecular variance (AMOVA) reflected genetic differences both within and among populations. The isolation by distance (IBD) analysis revealed a high positive correlation between the geographic distance and genetic distance of the different populations. Neutral test and mismatch distribution suggested that A. aurantii may have experienced recent population expansion events.


INTRODUCTION
The tea aphid (Aphis aurantii Boyer de Fonscolombe, 1841), also known as the "black citrus aphid, " is widely distributed in tropical and subtropical regions, including the Mediterranean region, England, Africa, India, Southeast Asia, Australia, South America, Central America, and North America (Tahori and Hazan, 1970;Carver, 1978). A. aurantii is an extremely polyphagous species with a very wide host range, inhabiting plants from over 190 genera of 80 families. Many of these hosts are economically important plants, such as citrus, coffee, tea, cacao, avocado, loquat, litchi, mango, fig, Camellia spp., Cinchona spp., Annona spp., Macadamia spp., Piper spp., and Artocarpus spp. (Carver, 1978). A. aurantii shows a preference for members of the families Rutaceae, Rosaceae, Apocynaceae, and Rubiaceae. In China, A. aurantii frequently occurs throughout the tea regions and causes direct feeding damage to fresh leaves and tender shoots, seriously reducing the output and lowering the quality of commercial teas. The excreted honeydew of A. aurantii can also incur mildew on the leaves, which reduces photosynthesis of the tea leaves.
In Chinese tea gardens, the main strategy for controlling these aphids is using chemical pesticides. However, this method has resulted in three main issues: residue, resistance, and resurgence. Very few genetic studies have been conducted for A. aurantii. The partial mitochondrial cytochrome c oxidase I (cox1) genes of some A. aurantii specimens were sequenced from Guangxi Province in southern China (Wang and Qiao, 2009), and apparent intraspecific genetic diversity was found within A. aurantii. The mitochondrial genome of A. aurantii from Guangxi Province was sequenced (Wang et al., 2019) and the phylogeny of Aphidoidea was rebuilt, revealing the close relationship between A. aurantii and Aphis craccivora Koch, 1854. A recent study sequenced the mitochondrial genome of A. aurantii from Sichuan Province and constructed the phylogeny of Aphididae (Pu et al., 2020), supporting the close relationship of A. aurantii and other aphids in Aphis. The transcriptome of A. aurantii has also been sequenced and analyzed (Hong et al., 2020), providing a basic transcriptomic dataset to identify simple sequence repeats (SSRs). However, these molecular studies did not involve population genetic analysis for the Chinese A. aurantii. To better understand the population genetics of this important pest and provide new insights for controlling this pest, we sampled 11 sites (Figure 1) and examined the population genetic structure of Chinese A. aurantii based on the concatenated sequences of two mitochondrial genes: cox1 and cytb.

Sampling and DNA Extraction
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. R Tissue DNA Kit (Omega, Norcross, GA, United States) following the manufacturer's protocol and preserved at −20 • C.

PCR Reaction and Sequencing
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).

Analyses of Genetic Variation
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 Analyses
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).

DNA Sequencing
The 220 sampled individuals of A. aurantii yielded 196 highquality 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.

Genetic Variation and Diversity
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 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.

Neutrality Test
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).

Mismatch Distribution Analysis
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).

Phylogenetic Relationships and Genetic Structure
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.

Genetic Variations
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)    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.

Demographic History
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.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession numbers can be found below: NCBI (accession: MW265715-MW265910 and MW289592-MW289792).

ETHICS STATEMENT
All experiments and procedures for this study complied with the current animal ethics guidelines and did not involve any protected animals.

AUTHOR CONTRIBUTIONS
H-LL and D-QP conceived and designed the experiments. CL and X-LW collected and identified the samples. H-LL, Z-TC, and K-JX performed the experiments, analyzed the data, and wrote the manuscript. Z-TC prepared the figures. All authors reviewed the manuscript and contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We are grateful to the editor and reviewers for their helpful comments.