All samples were collected into and cryopreserved in 95% or 100% ethanol. DNA from one individual per sample was isolated for molecular studies and three to five individual aphids per collection were mounted on microscope slides for morphological examination. Preserved aphid colonies were examined prior to preparation to ensure that they did not consist of multiple species. Voucher specimens for each sample were identified by G.X. Qiao based on morphological diagnostic features using standard literature-based keys (esp. Blackman and Eastop 1994, Pintera 1987, Wieczorek 2010) and by a comparison with previously identified specimens in the National Zoological Museum of China, Beijing. To avoid mutual influence and to ensure the independence of the different research methods, the morphological identification and molecular research were performed independently. All samples and voucher specimens were deposited in the National Zoological Museum of China, Institute of Zoology, Chinese Academy of Sciences, Beijing, China. Details of the sequenced taxa and voucher information are listed in Suppl. material 1.

Three aphid genes were targeted: mitochondrial cytochrome oxidase c subunit I (COI), cytochrome oxidase c subunit II (COII), and cytochrome b (Cytb), and one aphid endosymbiont Buchnera gene gluconate-6-phosphate dehydrogenase (gnd) (Kim and Lee 2008, Wang and Qiao 2009, Zhang et al. 2011, Chen et al. 2012).

Total genomic DNA was extracted from single aphid. Individual aphids were selected from the ethanol-preserved candidates with a destructive DNA extraction procedure. Plump adults are the ideal experimental material, but they must be examined under a microscope (Leica DM 2500) to eliminate parasitized individuals. Total DNA was extracted by following the Quick-Start protocol of DNeasy Blood & Tissue Kit (QIAGEN, Dusseldorf, Germany) with a single individual. The DNA solution was then stored at -20 °C for subsequent molecular experiments.

The polymerase chain reaction (PCR) mixture for the amplification of COI, COII, Cytb, and gnd genes comprised 22 μl of double distilled water (ddH₂O), 3 μl of 10 ×EasyTaq Buffer (+ Mg²⁺) (TransGen Biotech, Beijing, China), 2.4 μl of 2.5 mM/800 μl dNTPs (TransGen Biotech), 0.6 μl of 10 pmol/μl forward and reverse primers, 0.4 μl of 5 U/μl EasyTaq DNA Polymerase (TransGen Biotech), and 1 μl of DNA solution for a total volume of 30 μl.

The PCR conditions differed according to the gene and the specific primers, especially the annealing temperature, which was the most critical factor influencing product quality. The detailed primer information is shown in Suppl. material 2. The thermal setup of primer LepF/LepR (Foottit et al. 2008) or LCO1490/HCO2198 (Folmer et al. 1994) for COI gene fragment was: a 5-minute initial denaturation at 95 °C followed by 35 cycles of 30-second denaturation at 95 °C, 30 seconds of annealing at 50 °C, a 1-minute extension at 72 °C, and a 10-minute final extension at 72 °C. The protocol for primer mt2993+ (Stern 1994)/A3772 (Normark 1996) for tRNA/COII molecular marker was as follows: a 5-minute initial denaturation at 95 °C followed by 35 cycles of 1-minute denaturation at 95 °C, 1 min at 42 °C , 1 min at 72 °C, and a 7-minute final extension at 72 °C. The parameters of primer CP1/CP2 (Harry et al. 1998) for Cytb amplification was simplified as: 94 °C with 5 min, and 40 cycles of 94 °C with 50s, 48 °C with 1 min, 72 °C with 1.5 min, and 72 °C with 10 min. The setup of primer BamHI/ApaI (Clark et al. 1999) for Buchnera gnd gene was predigested as: 94 °C with 5 min, and 30 cycles of 94 °C with 1 min, 55 °C with 30s, 72 °C with 1 min, and 72 °C of 10 min final extension.

The amplification products were detected by 1.5% agarose gel electrophoresis (AGE), and then purified using EasyPure Quick Gel Extraction Kit (TransGen Biotech). The eligible products were then sent to TsingKe Biological Technology, Beijing, China or BGI, Shenzhen, China for sequencing, which was required to be bidirectional.

The returned forward and reverse chromatograms were loaded and then assembled and edited by SeqMan in DNAStar software (DNASTAR, Madison, Wisconsin, USA). The nucleotide sequences were first examined in NCBI by Basic Local Alignment Search Tool (BLAST) (Altschul et al. 1990) to test their affiliations. Concurrently, for the encoding gene fragments, we translated the assembled contigs into amino acids by MEGA6 (Tamura et al. 2013) to examine whether the sequences were correct and accurate. Multiple alignments were accomplished by MAFFT (Katoh and Standley 2013), and the sequences were then adjusted and trimmed manually in MEGA6. It is noteworthy that the sequences amplified with primer mt2993+/A3772 covered the COII gene fragment as well as a tRNA, which was then removed for subsequent analysis.

In addition to sequences from 425 samples, we downloaded 245 COI and 1 COII sequence from NCBI. Here, we defined the datasets as COI-670 (including the whole research group and NCBI sequences), COII-376 (including 375 internal sequences and 1 NCBI sequence), Cytb-413 (newly gotten for this study), gnd-396 (newly obtained sequences), and COI-338, COII-338, Cytb-338, gnd-338, which contained only the specimens that acquired all 4 gene sequences.

A neighbor-joining (NJ) (Saitou and Nei 1987) tree was constructed by MEGA6 based on the aligned sequences. One thousand bootstrap replications were calculated to assess the credibility of the NJ analysis. The Kimura 2-parameter (K2P) model of base substitution (Kimura 1980) was selected in pairwise distances calculation, and for the more accurate comparison between sequences, the pairwise deletion pattern was selected for gaps/missing data treatment. After bootstrap consensus trees with bootstrap values at each node were obtained, we computed the condensed tree with a 50% cut-off value for the consensus tree. When analyzing the COI-670 tree, we chose a threshold of 2% (Foottit et al. 2009b) for a cluster standard, which has been well used in aphids. With regard to the COII-376, COII-338, Cytb-413, Cytb-338, gnd-396 and gnd-338 NJ trees, we calculated only the cluster topologies.

The Automatic Barcode Gap Discovery (ABGD) (Puillandre et al. 2012a) approach is a model-based method for delimiting species. Based on the existence of a barcoding gap (namely the intraspecific divergences are smaller than interspecific divergences) and a prior intraspecific divergence (p), the ABGD procedure first sorts the dataset into a hypothetical species, and then computes recursively with the previous groups to obtain a result optimized until there are no better partitions. We ran the ABGD with a graphic web version (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html). First, we calculated the distance values among samples by using MEGA6 with p-distance, Jukes-Cantor (JC69) model, and K2P model separately, and the result data were saved as CSV format file. We then chose 0.055, which had been suggested for Aphididae (Foottit et al. 2009b), as the prior intraspecific divergence for COI-670 and COI-338 datasets, and we used a p value of 0.1, which was the default and shown to be sufficient for analysis, for the other data sets. The other parameters were maintained by default for all analyses.

The General Mixed Yule Coalescent (GMYC) (Fujisawa and Barraclough 2013) is a tree-based approach for the delimitation of species. We ran the GMYC method in R project (available from: https://www.r-project.org/) by using the “splits” package (available from: http://r-forge.r-project.org/projects/splits). The input tree was required to be strictly ultrametric and bifurcating, which meant there was no zero-length branch. Here, we used a maximum likelihood (ML) tree as the input. Therefore, the haplotypes were calculated and generated by DnaSP (Librado and Rozas 2009), and an ML tree was constructed by RAxML (Stamatakis 2006) with haplotype data. Due to the ultrametric and bifurcating requirement, the ML tree was constructed with r8s (Sanderson 2003). The outcome tree modified by r8s was read into the “splits” R package, and the delimiting result was obtained with relevant commands.

The Poisson Tree Process (PTP) (Zhang et al. 2013) model is another tree-based method for inferring putative species. The PTP approach is a close relative of the GMYC method, but it only needs a simple phylogenetic tree as its input without requiring it to be ultrametric and bifurcating. As an updated version of the original PTP, the bPTP method was employed simultaneously to separate hypothetical species, which added a Bayesian support value to the tree. The PTP and bPTP analyses were run on a web server (http://species.h-its.org/ptp/) and the value 500,000 was selected for MCMC generations, with the other parameters set by default. The input tree was an ML tree constructed by RAxML with GTRCAT model. However, we encountered the same problem as Schwarzfeld and Sperling (2015), namely that the bPTP analysis failed to show convergence under 500,000 generations (the upper limit of the web server). Therefore, only the PTP result is displayed and discussed below.

The 425 samples collected by the group members in recent years were carefully authenticated with mounted individuals under the microscope, and all 425 samples were identified to species. The few vouchers with uncertain species identification were sorted into featured clusters and were given the epithet “sp.”, which made them convenient for further analysis. A total of 75 morphological species were determined from 670 whole samples, and 51 were identified from the 425 mounted samples.

The COI sequences were trimmed to a length of 658 bp, which included 365 conserved sites, 293 variable sites and 258 parsimony-informative sites. The sequences had an average nucleotide composition of 38.0% T, 17.1% C, 34.4% A, and 10.5% G. The COII sequences were trimmed to a final length of 672 bp, among which 399 sites were conserved, 273 sites were variable, and 251 sites were parsimony-informative. The average T, C, A, G compositions of these sequences were 38.7%, 14.0%, 39.5%, and 7.8%, respectively. The Cytb gene was 760 bp, in which there were 420 conserved sites, 340 variable sites and 303 parsimony-informative sites. The Cytb sequences consisted of 41.4% T, 15.3% C, 34.3% A, and 9.0% G. We obtained a total length of 807 bp for the gnd gene with an average nucleotide composition of 37.8% T, 9.8% C, 39.5% A, and 12.8% G, among which there were 368 conserved sites, 439 variable sites and 417 parsimony-informative sites. Across all 4 genes, a strong T and A nucleotide composition bias existed.

From a total of 425 samples, 425 COI gene fragment sequences, 375 COII gene fragment sequences, 413 Cytb gene fragment sequences, and 396 gnd gene fragment sequences were acquired. The successive amplification efficiency of those markers in order was COI (100%) > Cytb (97%) > gnd (93%) > COII (88%).

Genetic divergences were assessed by 5 disparate metrics among and within species. For the interspecific divergences of congeneric species, we chose the average interspecific distance, which was calculated within genera that contained more than one species, and the smallest interspecific distance, which meant the minimal value of interspecific distance within genera with at least two species. When evaluating the intraspecific divergences, three variables (average intraspecific distance, mean theta, and average coalescent depth) were applied. The average intraspecific distance was the average value of the genetic distances between samples within species that had at least two individuals. The mean theta signified a modified theta, which expressed the average pairwise distance scored for species with more than one obtained representative, by dislodging improper individuals concerned with the asymmetrical acquisition of samples. The average coalescent depth, namely the average value of maximum intraspecific distance, was calculated for species in which there were no fewer than two samples.

All five interspecific and intraspecific metrics were determined within genera and species (Table 1). The results of different genes and datasets showed distinctly high interspecific divergences and low intraspecific distances. For DNA barcoding the smallest interspecific distance and average coalescent depth were the most useful and intuitive parameters. Within Chaitophorinae, the genetic divergence ranges of smallest interspecific distance and average coalescent depth of COI, COII, Cytb, and gnd were (0.0693–0.1233, 0.0060–0.0218), (0.0563–0.1110, 0.0024–0.0070), (0.0703–0.1060, 0.0052–0.0230), and (0.0807–0.1427, 0.0010–0.0090), respectively. The figures above were obtained across the whole dataset. To obtain a more reliable and comparable analysis, we calculated the results of the 338-sample datasets. The smallest interspecific distance and average coalescent depth ranges of COI-338, COII-338, Cytb-338, and gnd-338 were (0.0693–0.0995, 0.0044–0.0165), (0.0563–0.1190, 0.0024–0.0090), (0.0703–0.1040, 0.0052–0.0230), and (0.0811–0.1600, 0.0010–0.0090), respectively. The computations above showed a properly high smallest interspecific distance and a comparatively low average coalescent depth.

To observe the occurrence frequency of different genetic divergences, we drew the frequency line charts of inter- and intra-specific genetic distances based on 338 datasets (Figure 1). Each gene was signified in one chart: the top half was calculated with all 338 samples; and the bottom half was scored by eliminating two samples (Nos. 25138 and 25161) of Chaetosiphella longirostris. The overlap region was indicated by the red dotted rectangle. No obvious barcoding gap was found in these samples across COI, COII, Cytb, and gnd genes. The overlap regions of COI, COII, Cytb, and gnd in the top half were 0.000–0.031, 0.000–0.023, 0.000–0.045, 0.000–0.018; and those in the bottom half were 0.022–0.031, 0.011–0.023, 0.015–0.045, 0.006–0.018, respectively. It was clear that the overlap region was much narrower by eliminating the questioned vouchers, and the total frequency within that region was also reduced significantly (Figure 1). The data in the overlap region represented samples that the DNA barcoding would fail to identify. Therefore, a lower total frequency in that region was better. The order of the total frequencies in the overlap region was COI < COII < Cytb < gnd in the top and bottom half parts.

Figure 1. 

Frequency line charts of inter- and intra-specific genetic distances based on 338 dataset. The x-axis represents the genetic distance, and the y-axis represents the occurrence times in the whole genetic distance matrix. Each peak was a data point with corresponding genetic distance and occurrence times. The data points on the green and red line were calculated with the interspecific distances, and the points on purple and blue line were calculated with the intraspecific distances. The overlap region, which was the crossing area of inter- and intra-specific divergence, is indicated by the red dotted rectangle. Each gene was signified in one chart: the top half was calculated with all the 338 samples; and the bottom half was scored by eliminating the queried samples of Chaetosiphella longirostris.

Figure 1. 

Continue.

Given that there were no unambiguous and credible references of genetic thresholds for COII, Cytb, and gnd in aphids, the method of threshold with NJ was applied only in the COI-670 dataset with a threshold of 2% (Foottit et al. 2009b). For all the COI, COII, Cytb and gnd datasets, one distance-based method ABGD (for the analysis results, see Suppl. material 3 for details) and two tree-based approaches GMYC and (b)PTP were employed simultaneously as reference counterpoints to each other.

The morphological and molecular identification results are shown in Table 2. To suitably display the analysis results, four group names were used: accurate, of which the putative species clusters were identical with morphological identifications; split, whose component vouchers were merely part of a specific unabridged morphological species without representatives of others; lumped, which was defined as an aggregation of more than one species including all samples of those species; and partial lumped, a multi-species cluster that consisted of all the vouchers of one or several species as well as part samples of other species.

Seventy-five morphological species were identified from COI-670, which included sequences downloaded from NCBI. For COI-670, we obtained 89 putative species by the GMYC approach with a 65.17% accuracy rate, 85 species using PTP with 67.06% accuracy, 81 species using ABGD with 72.84% accuracy, and 81 species by threshold-NJ with 72.84% accuracy. The COII-376 dataset with only one sequence from NCBI contained 51 morphological species, and 53 hypothetic species were gleaned using the GMYC method with an accuracy rate of 83.02%, 48 species using PTP with 83.33% accuracy, and 50 species using ABGD with 82.00% accuracy. The Cytb-413 data contained 48 morphological species and generated 54 clusters by the GMYC method with an accuracy rate of 79.63%, 49 species using PTP with 81.63% accuracy, and 48 species using ABGD with 79.17% accuracy. There were 49 morphological species within gnd-396, and the putative species found using the GMYC approach was 46 with an accuracy rate of 86.96%, using the PTP method was 45 species with an accuracy rate of 84.44%, and using the ABGD analysis was 48 species with 87.50% accuracy. An analysis of COI-338, COII-338, Cytb-338, and gnd-338 were performed to compare the results of different genes with diverse methods under the same sample composition. There were 45 morphological species within the given 338 samples. The analysis results of various genes and accuracy rate were: COI-338 (GMYC: 47, 89.36%; PTP: 49, 85.71%; ABGD: 46, 91.30%), COII-338 (GMYC: 45, 93.33%; PTP: 45, 86.67%; ABGD: 45, 86.67%), Cytb-338 (GMYC: 50, 80.00%; PTP: 46, 80.43%; ABGD: 42, 83.33%), and gnd-338 (GMYC: 42, 92.86%; PTP: 41, 90.24%; ABGD: 44, 93.18%). The final results were all displayed in NJ trees (see Suppl. material 4–11) with bootstrap values. The tree topology should be regarded only as distance affinity and not a phylogenetic relationship.

COI has not been the only gene marker used for aphid DNA barcoding, other genes from the mitochondrial genome and from endosymbionts having been used for various aphid groups (Chen et al. 2012, Chen et al. 2013a, Liu et al. 2013, Tang et al. 2015) . Additionally, many analytical methods have been used in DNA barcoding. Here, we focused on COI, COII, Cytb, and gnd genes, and applied four methods, threshold with NJ, ABGD, GMYC, and PTP, to delimitate species of Chaitophorinae.

For different genes, the amplification efficiency was COI (100%) > Cytb (97%) > gnd (93%) > COII (88%). Within the 338-sample dataset, the difference in values between the smallest interspecific distance and average coalescent depth were unequal in different groups. For Chaitophorus, Periphyllus, and Chaetosiphella (Table 1), the difference values were COI > Cytb > COII > gnd, gnd > COI > COII > Cytb, and gnd > COI > COII = Cytb, respectively. From the overlap region and total frequency (Figure 1), COI and COII were similar with a narrower overlap region and less frequency than Cytb. Although the overlap span of gnd was sufficient, its total frequency in that region was slightly larger. Therefore, the COI and COII genes may be better markers for DNA barcoding.

The most important factor in choosing the delimitation method was the identification accuracy within different genes. Therefore, a better approach means higher identification accuracy and a greater range of application with various genes. The accuracy of GMYC, PTP, and ABGD within COI-338, COII-338, Cytb-338, and gnd-338 were ABGD (91.30%) > GMYC (89.36%) > PTP (85.71%), GMYC (93.33%) > PTP = ABGD (86.67%), ABGD (83.33%) > PTP (80.43%) > GMYC (80.00%), and ABGD (93.18%) > GMYC (92.86%) > PTP (90.24%), respectively (Table 2). In Chaitophorinae, the ABGD was a much better analytical method, and GMYC was also better than PTP in tree-based approaches. Considering that the analysis of ABGD required prior intraspecific distance, a tree-based method needed to be employed concurrently. A tree-based approach should be crosschecked against a non-tree-based approach within species delimitation studies (Fontaneto et al. 2015). As different methods may yield inconformity conclusions (Carstens et al. 2013), the accurate identification of species requires further integrative analysis (Puillandre et al. 2012b, Pante et al. 2015). Herein, a brief investigation of the morphological characteristics combined with a distance-based method of ABGD and a tree-based method of GMYC may be a very suitable pattern for species delimitation and the rapid identification of Chaitophorinae.

Chaitophorus saliniger Shinji is an important pest on willows in East Asia. Based on the topology structure and results of analysis with different methods (Figure 2A, Suppl. material 4–11), all samples of this species were divided into two clades which could be regarded as two different species. However, based on the morphological characteristics, all samples should be Chaitophorus saliniger. On examining the geographical distribution information of all samples, we found that one of the clades consisted of two samples (Nos. 17651 and 33320) of C. saliniger collected from Northeast China, Heilongjiang Province, and the voucher of another sequence (C. saliniger-GU978785.1, a downloaded sequence from NCBI) from the Korean Peninsula. The location sites of these three samples were all at a relatively high latitude. All samples in the other clade were not from the aforementioned regions. So, within C. saliniger, a population differentiation emerged among the samples from different locations. The genetic divergences between the two clades for COI, Cytb, and gnd were 0.043, 0.039, and 0.020, respectively, which could be regarded as interspecific distances. Although the morphological characteristics of all samples were similar, differentiation at the gene level seems to have occurred between northern and southern populations. Similar genetic differentiation between these populations has also been demonstrated for a nuclear gene, EF-1α (Fang et al. 2016).

Figure 2. 

The analysis results of some species from the COI-670 dataset. The analysis results based on other genes were almost identical. The NJ tree was constructed based on the Kimura 2-parameter (K2P) model with a bootstrap value over 50% displayed. The gray blocks behind the tree represent the putative species, which means that the taxa in the tree corresponding to a single block are in one putative species. The number of blocks express the number of putative species using this method. A Chaitophorus saliniger B Laingia psammae.

In a similar manner, two samples (Nos. 17613 and 19950) of Laingia psammae Theobald were divided into two independent clades (Figure 2B), a result supported at all genes and using different approaches. Sample No. 17613 was from Jilin, northeastern China, whereas sample No. 19950 was from Xinjiang, northwestern China. The genetic distances of three genes (COI, COII and Cytb) between the two samples were 0.064, 0.068, and 0.062, respectively, which reach the level of species (Wang and Qiao 2009). Therefore, differentiation between northeastern and northwestern populations in L. psammae exists.

From the topology structures and the constructed consequences of threshold with NJ, ABGD, GMYC, and PTP, we observed that population differentiation was clearly present within both C. saliniger and L. psammae. Similar findings have been reported in other aphid species (Lee et al. 2011, Wang et al. 2011). The prominent differences among populations may even be an indication of cryptic species (Bickford et al. 2007). Speciation is a long and continuous process, and cryptic species are not easily explained. Within the process, the incipient species may hold for millions of years (Avise and Walker 1998). Therefore, cryptic species need further study with more samples, combined with morphological characteristics and biological information.