Introduction
Ceratovacuna Zehntner (Reference Zehntner1897) and Pseudoregma Doncaster (Reference Doncaster1966) are closely related genera that belong to Cerataphidini in Hormaphidinae (Hemiptera, Aphididae). Most cerataphidini aphids alternate between their primary hosts, Styrax (Styracaceae), where they produce galls, and their secondary hosts (usually Asteraceae and bamboos) (Ghosh, Reference Ghosh1988; Stern et al., Reference Stern, Aoki and Kurosu1997a). Almost all species of Cerataphidini produce soldier aphids in galls on their primary hosts, but many Pseudoregma and Ceratovacuna species also produce horned soldiers on their secondary hosts. Behavioral, morphological, and phylogenetic evidence suggest that primary and secondary host soldiers are not homologous and that horned soldiers evolved after the primary host soldiers (Aoki & Kurosu, Reference Aoki and Kurosu1989, Reference Aoki and Kurosu2010; Fukatsu et al., Reference Fukatsu, Aoki, Kurosu and Ishikawa1994). The unusual biology of the horned soldiers suggests that they are interesting groups for investigating social behavior (Aoki, Reference Aoki, Ito, Brown and Kikkawa1987; Dixon, Reference Dixon1985, Reference Dixon1998; Stern, Reference Stern1994, Reference Stern1995, Reference Stern1998; Stern et al., Reference Stern, Aoki and Kurosu1997a, Reference Stern, Whitfield and Fosterb; Costa, Reference Costa2006). However, to better understand the origin and history of sociality in these species, their evolutionary history must be clarified first.
The genus Ceratovacuna contains 21 known species (Qiao et al., unpublished data), while Pseudoregma includes 12 species (Qiao et al., unpublished data). Using mitochondrial genes, Stern (Reference Stern1994) performed the first molecular phylogenetic study of Cerataphidini and found that Ceratovacuna and Pseudoregma each formed a monophyletic group (fig. 1a). However, using the same mitochondrial sequences but replacing P. alexanderi with P. bambusicola and P. sundanica, the monophyly of Ceratovacuna and Pseudoregma was not supported (fig. 1b) (Stern, Reference Stern1995). The phylogenetic relationships of Ceratovacuna and Pseudoregma have been discussed in several works later (Stern et al., Reference Stern, Aoki and Kurosu1997a; Stern, Reference Stern1998; Fukatsu et al., Reference Fukatsu, Shibao, Nikoh and Aoki2001). However, these studies have mainly focused on the life history and evolution of soldier aphids. While the monophyly of Pseudoregma was supported in most of these studies, the status of Ceratovacuna remained inconclusive. In addition, these studies were based on mitochondrial genes, and sampling bias existed due to specific study aims (e.g. evolution of life history and soldier aphids). Therefore, a detailed investigation of the phylogenetic relationships of these two genera is needed.
Phylogeny of Cerataphidini based on mitochondrial COI and COII (Stern, Reference Stern1994) (a) and Stern (Reference Stern1995) (b).
In addition to their uncertain phylogenetic history, the species of both Ceratovacuna and Pseudoregma are widely distributed across complex habitats in South-East Asia. Their distribution ranges include Java in the south, Japan and Korea in the north, eastern India in the west, and Taiwan and Japan in the east (Stern et al., Reference Stern, Aoki and Kurosu1997a). Considering their similar distributions, the two genera provide good models for investigating the structure of inter- and intraspecific genealogical patterns and examining how aphid species diversification correlates with historical events that have shaped the speciation of co-distributed taxa.
To test the monophyly of Ceratovacuna and Pseudoregma, we used a more extensive sampling and both mitochondrial (cytochrome c oxidase subunit I (COI); cytochrome b (Cytb)) and nuclear (elongation factor-1α (EF-1α); long-wavelength opsin (LWO)) genes to reconstruct the phylogenetic relationships of these two related genera. We also estimated divergence times and tested whether geographical barriers and host evolution have been the possible driving forces for species differentiation in both genera.
Materials and methods
Taxon sampling
We sampled 12 species of Ceratovacuna and five species of Pseudoregma. Specimens for slide-mounting were stored in 75% ethanol. Ceratovacuna sp1., Ceratovacuna sp2., Ceratovacuna sp3., Ceratovacuna sp4., and Ceratovacuna sp5. have recently been recognized as new species based on morphological and molecular data, and this information will be published separately. Samples for molecular experiments were stored in 95% ethanol. All samples were preserved at −30 °C. On the basis of current knowledge of the phylogenetic relationships within Hormaphidinae and Cerataphidini (Stern, Reference Stern1998; Ortiz-Rivas et al., Reference Ortiz-Rivas, Moya and Martínez-Torres2010; Huang et al., Reference Huang, Xiangyu, Ren, Zhang, Zhang and Qiao2012), outgroups were chosen from Aleurodaphis van der Goot (1917), Astegopteryx Karsch (1890), Cerataphis Lichtenstein (1882), Ceratoglyphina van der Goot (1917), Tuberaphis Takahashi (1933) in Cerataphidini; Hamamelistes Shimer (1867) and Hormaphis Osten-Sacken (1861) in Hormaphidini; Neohormaphis Noordam (1991) and Nipponaphis Pergande (1906) in Nipponaphidini, respectively. Multiple outgroups are indispensable to reliably test the monophyly of the two genera. All voucher specimens were deposited in the National Zoological Museum of China (NZMC) at Institute of Zoology, Chinese Academy of Sciences (Beijing) (voucher numbers see table 1). Collection information for all samples, including collection localities, host plants, and collection dates are listed in table 1.
Related information of aphid samples examined in this study.
DNA extraction, PCR and sequencing
The CTAB DNA isolation technique (Vavre et al., Reference Vavre, Girin and Boulétreau1999) was used to extract genomic DNA. A 700 bp fragment of COI was amplified using primers LepF (5′-ATTCAACCAATCATAAAGATATTGG-3′) and LepR (5′-TAAACTTCTGGATGTCCAAAAAATCA-3′) from Foottit et al. (Reference Foottit, Maw, von Dohlen and Hebert2008). About 800 bp of Cytb was amplified using primers Cp1 (5′-GATGATGAAATTTTGGATC-3′) and Cp2 (5′- CTAATGCAATAACTCCTCC-3′) from Harry et al. (Reference Harry, Solignac and Lachaise1998). The primers for amplification of EF-1α were EF3 (5′-GAACGTGAACGTGGTATCAC-3′) and EF2 (5′- ATGTGAGCAGTGTGGCAATCCAA-3′) (Palumbi, Reference Palumbi, Hillis, Moritz and Mable1996). LWO sequences were amplified using primers OPSETF1 (5′-GGYRTYACNATTTTYTTCTTRGG-3′) and OPSETR1 (5′-GANCCCCADATYGTNAATAAYGG-3′) (Benjamin Ortiz-Rivas, private communication). PCR for Cytb, EF-1α, and LWO was performed in 30 μl reaction volumes: 3 μl 10×PCR buffer, 2.4 μl dNTPs, 20 μl dd H2O, 0.4 unit Taq DNA polymerase (all from TransGen Biotech, Beijing, China) and 0.6 μl 10 μM forward and reverse primers (synthesized by Invitrogen Biotech, Shanghai, China). The PCR volumes for COI were the same as those for Cytb, except for 19.4 μl dd H2O and 0.8 μM primers. All PCR reactions included a 95 °C hot start for 5 min and a 72 °C extension for 10 min at the end. The cycling conditions for COI included 34 cycles of denaturation at 95 °C for 30 s, annealing at 52 °C for 30 s, and extension at 72 °C for 1 min. The cycling conditions for Cytb included 35 cycles of denaturation at 92 °C for 1 min, 48 °C for 1.5 min, and 72 °C for 1 min. The PCR thermal regime for EF-1α and LWO included 35 cycles of denaturation at 94 °C for 1 min, annealing at 50 °C for 1 min, and extension at 72 °C for 1.5 min. Sequencing reactions were performed using the same primers that were used for PCRs using ABI 3730 automated sequencer (ABI, USA).
Sequence alignment and phylogenetic analyses
Sequences were assembled using Seqman II (DNASTAR) and checked manually. Multiple alignments were performed with ClustalX 1.81 (Thompson et al., Reference Thompson, Gibson, Plewniak, Jeanmougin and Higgins1997). Nucleotide composition, conserved sites, variable sites, and parsimony informative sites were calculated using MEGA 4 (Tamura et al., Reference Tamura, Dudley, Nei and Kumar2007). We used the coding regions of EF-1α and LWO in phylogenetic analyses. Introns of EF-1α and LWO were removed according to the GT-AG rule and comparison with the cDNA sequences of Hormaphis similibetulae [GenBank accession number: DQ493866] and Cerataphis sp. D15 [GenBank accession number: AJ539465], respectively. A partition homogeneity test was implemented with PAUP* 4.10b (Swofford, Reference Swofford2002) and used to evaluate the appropriateness of combining different genes into one dataset. The results of the partition-homogeneity test indicated that all genes can be grouped to four combined datasets (COI+Cytb, EF-1α+LWO, COI+Cytb+EF-1α and COI+Cytb+EF-1α+LWO) (P<0.01). For Bayesian analysis, partitions were made for each gene separately using their own parameters derived from Modeltest 3.7 (Posada & Crandall, Reference Posada and Crandall1998).
Phylogenetic analyses were conducted in PHYML (Guindon & Gascuel, Reference Guindon and Gascuel2003) and MrBayes 3.1.2 (Ronquist & Huelsenbeck, Reference Ronquist, Huelsenbeck and Nielsen2005). Modeltest was used to determine the best-fit nucleotide substitution model for maximum likelihood (ML) analysis. The combined data were partitioned into mitochondrial versus nuclear data, and the best model chosen under the Akaike Information Criterion (AIC) (Posada & Buckley, Reference Posada and Buckley2004) for COI, COI+Cytb, EF-1α, EF-1α+LWO, and COI+Cytb+EF-1α+LWO was GTR+I+G. ML analysis was performed using PHYML, and nodal support was evaluated by bootstrap analysis with 1000 replicates. The values of model parameters were treated as unknown variables with uniform prior probabilities, and were estimated during the analyses. Four chains were run, starting from a random tree and proceeding for 15 million Markov chain Monte Carlo generations, with tree sampling every 100 generations. Four independent runs were conducted to verify the results. The first 2500 trees (25% of the total) were discarded as burn-in. For the remaining trees, 50% majority-rule consensus trees were generated. An examination of the log-likelihood scores, the average standard deviation of split frequencies (<0.01), and potential scale reduction factor (PSRF) for branch lengths suggested that the burn-in periods were long enough for chains to become stationary.
Divergence time estimation
The new species Unicohormaphis sorini Wegierek & Zyla (Wegierek & Zyla, Reference Wegierek and ŻYŁA2011), which belongs to Cerataphidini, was found in Baltic amber and dates back to the Eocene. Since fossils generally indicate the minimum age of a species, we assigned a conservative age of 57 mya (with normal prior distribution) as a calibration point to represent the minimum age for the most recent common ancestor (MRCA) of Cerataphidini in our study. Based on a mitochondrial tRNA/COII dataset and calibrated nucleotide substitution rate, von Dohlen et al. (Reference von Dohlen, Kurosu and Aoki2002) calculated the mean of divergence estimates between East-Asian and North American Hormaphis as 25.5 mya. We also used this age as a calibration point to verify the results of fossil calibration. Based on mitochondrial and nuclear data, the divergence times among species were estimated using BEAST (v1.6.1, Drummond & Suchard, Reference Drummond and Suchard2010) under a relaxed molecular clock and uncorrelated lognormal model, respectively. A GTR model of nucleotide substitution with gamma rate heterogeneity among sites and a proportion of invariant sites and Yule speciation process were used.
After determining optimal parameters over several preliminary runs, multiple analyses were run for 10 million generations, and 25% of the samples were removed as burn-in. The results were visualized in Tracer v1.5 (Rambaut & Drummond, Reference Rambaut and Drummond2008), where the stationarity was assessed by ensuring that the effective sample size (ESS) values were high for most parameters. The results proved repeatable over multiple independent runs. Thus, a chronogram was reconstructed based on results from a representative single run and visualized using the program FigTree v1.3.1 (Rambaut, Reference Rambaut2009).
Results
Sequence characters and phylogenetic analyses
COI sequences were obtained from 85 samples. The sequence length varied from 608 to 703 bp, and the final alignment was of 574 bp. Seventy-eight Cytb sequences were obtained with length from 666 to 814 bp. Their final alignment was 665 bp. The length of the EF-1α sequences obtained from 67 samples varied between 990 and 1168 bp, with a final alignment of 806 bp. Forty-eight LWO sequences ranging from 948 to 1198 bp were obtained, and their 636 bp alignment was used for phylogenetic analyses. All sequences were submitted to GenBank (see table 1 for accession numbers). The intraspecific pairwise distances were 0–4.5% for C. lanigera and much lower in the other species (0–1.6%) (Suppl. 1).
Phylogenetic analyses
The Bayesian and ML analyses based on COI did not support the monophyly of Ceratovacuna and Pseudoregma, when all species studied here were included; however, the monophyly of the two genera as a whole was suggested (ML/BI: 78/1.00) in the phylogeny including most related species from other genera (fig. 2). Except C. silvestrii and P. alexanderi, the taxonomic status of most species in these two genera was supported. Although with low resolution, all sampled species tended to cluster into three main groups: I, C. silvestrii+C. hoffmanni+Ceratovacuna sp1.+Ceratovacuna sp2.; II, C. nekoashi+Ceratovacuna sp3.+C. graminum+Ceratovacuna sp5.+C. panici+C. lanigera; and III, P. sundanica+C. longifila+P. panicola+P. alexanderi+P. koshunensis+P. bambucicola. Ceratovacuna sp4., was located as a sister clade of all other species. The monophyly of group I was supported (ML/BI: 93/1.00). One sample (No. 14674) of C. lanigera formed a sister group with all other samples of this species (ML/BI: 70/0.96).
Phylogenetic relationships between genus Ceratovacuna and Pseudoregma based on maximum likelihood analysis of mitochondrial gene COI. Bootstrap percentages (≥50%) from ML are shown above the branch, posterior probabilities (≥50%) of Bayesian analysis are shown below the branch.
Phylogenetic analysis based on EF-1α (fig. 3) and the combined datasets of COI+Cytb (Suppl. fig. 1), EF-1α+LWO (Suppl. fig. 2), and COI+Cytb+EF-1α+LWO (fig. 4) yielded trees with similar topology to the COI tree. The monophyly of Ceratovacuna and Pseudoregma as a whole was supported in all analyses. In the EF-1α tree, monophyly of group II (ML/BI: 59/1.00) and Pseudoregma (ML/BI: 69/1.00) was supported, respectively. However, the monophyly of C. silvestrii in group I was not resolved.
Phylogenetic relationships between genus Ceratovacuna and Pseudoregma based on Bayesian analysis of nuclear gene EF-1α. Bootstrap percentages (≥50%) from ML are shown above the branch, posterior probabilities (≥50%) of Bayesian analysis are shown below the branch.
Phylogenetic relationships between genus Ceratovacuna and Pseudoregma based on Bayesian analysis of the combined datasets of COI, Cytb, EF-1α, and LWO. Bootstrap percentages (≥50%) from ML are shown above the branch, posterior probabilities (≥50%) of Bayesian analysis are shown below the branch.
The results of combined datasets improved the resolution of species relationships. The sister relationship of C. silvestrii, C. hoffmanni with Ceratovacuna sp2. (ML/BI: 91/1.00), and C. lanigera, C. nekoashi with C. panici (ML/BI:70/1.00) were supported based on the combined data from COI+Cytb (Suppl. fig. 1). Similarly, P. panicola with P. sundanica (ML/BI: 85/0.99), C. silvestrii with Ceratovacuna sp2. (ML/BI: 52/1.00), and Ceratovacuna sp3. with Ceratovacuna sp5. (ML/BI: 97/1.00) were all clustered into sister clades in the analysis based on a combined dataset of EF-1α+LWO (Suppl. fig. 2). Bayesian analysis of the combined dataset COI+Cytb+EF-1α+LWO strongly supported the monophyly of three main groups (fig. 4).
Divergence time
The divergence times were estimated separately for mitochondrial and nuclear data, and the 95% confidence intervals of most nodes were largely congruent (fig. 5; table 2). The analysis based on mitochondrial data included more taxa and was used to evaluate possible factors that influenced both species divergence and the formation of distribution patterns. The MRCA of Ceratovacuna and Pseudoregma were approximately 48.9 mya (95% highest posterior density, HPD: 35.6–54.9 mya). Most species divergences in the two genera occurred in the late Tertiary. The earliest speciation time was 24.3 mya (HPD: 9.5–37.9 mya) for C. lanigera, and the divergence of Ceratovacuna sp4. occurred at 8.8 mya (HPD: 4.9–13.1 mya).
A chronogram of Ceratovacuna and Pseudoregma based on the mitochondrial gene COI. Median age estimates are shown behind nodes (mya). Numbers 1–12 corresponding to table 2. Black dot represents two calibration points used in divergence time analyses. Locality abbreviated in capital letters: XZ, Tibet; GZ, Guizhou; YN, Yunnan; HN, Hainan; GX, Guangxi; FJ, Fujian; GD, Guangdong; ZJ, Zhejiang; SC, Sichuan; TW, Taiwan; CQ: Chongqing.
Estimated divergence times for selected nodes.
Node numbers are corresponding to that of fig. 5.
Discussion
Phylogenetic relationships of Ceratovacuna and Pseudoregma
Phylogenetic analyses based on both mitochondrial and nuclear genes have not reconstructed the monophyly of Ceratovacuna and Pseudoregma. Some previous studies suggested the monophyly of Pseudoregma (Stern, Reference Stern1994, Reference Stern1998; Stern et al., Reference Stern, Whitfield and Foster1997b; Fukatsu et al., Reference Fukatsu, Shibao, Nikoh and Aoki2001). In the present study based on more sampling and both mitochondrial and nuclear data, monophyly of Pseudoregma was supported by phylogenetic analyses based on nuclear genes and combined datasets (figs 3 and 4), while the monophyly of Ceratovacuna was unresolved in all analyses. The monophyly of the two genera as a whole was strongly supported in all analyses. These results are also supported by morphological evidences. The pleural grooves on the pronotum of apterous viviparous females were an important morphological characteristic in distinguishing Ceratovacuna and Pseudoregma (Ghosh, Reference Ghosh1988). The species of Pseudoregma have two pleural grooves on their pronotum that are separated by a median ridge, while Ceratovacuna species have no obvious pleural grooves. As more species are identified and described, distinguishing these two genera based on the occurrence/absence of pleural grooves on pronotum may be problematic, because this morphological characteristic varies within and among species. Molecular phylogenetic analyses and morphological characters, such as the arrangement and number of wax plates as well as pleural grooves on the pronotum of apterous viviparous females, provide evidence that Ceratovacuna and Pseudoregma should be redefined. In the future, more data are needed to test this assumption.
The sampled species of these two genera tend to cluster into three main groups. Groups I and II consist of Ceratovacuna species, except for P. sundanica, which clusters into group I in the analyses based on COI+Cytb (Suppl. fig. 1). All Pseudoregma species cluster into group III. Although the individual mitochondrial genes yield low resolution in resolving the topology of deeper nodes, the analyses based on nuclear genes and combined dataset of COI+Cytb+EF-1α+LWO (fig. 4) robustly supported the monophyly of the three main groups. Morphological evidence also supported molecular phylogenetic relationships between species. For example, the relationships of C. silvestrii with Ceratovacuna sp2., and the species of group II are further supported by their wax plate characteristics. The abdominal tergites of all species in group II have only marginal wax plates, while C. silvestrii and Ceratovacuna sp2. in group I, along with all species of Pseudoregma in group III, all exhibit spinal wax plates and some have pleural wax plates. In addition, wax plates on the head dorsum of C. silvestrii and Ceratovacuna sp2. are composed of 8–16 wax cells, and both of them are rich in small round or polygonal wax plates on their abdominal tergites. On the other hand, marginal wax plates on the abdominal tergites of all species in group II are composed of fewer than 12 wax cells. Compared to Pseudoregma, the wax plates of Ceratovacuna are more strongly developed. In contrast to species of Ceratovacuna, the posterior margin of the pronotum of Pseudoregma species, where some wax plates are distributed, is swollen; while this part of other species of Ceratovacuna is flat. Thus, morphological traits supported these three main groups of Ceratovacuna and Pseudoregma.
In the phylogenetic analyses based on EF-1α (fig. 3), the monophyly of C. silvestrii and P. alexanderi samples was not supported. The monophyly of C. silvestrii was not in the EF-1α+LWO tree. However, these two species could be distinguished in the trees using COI as well as the combined datasets of COI+Cytb and COI+Cytb+EF-1α+LWO. This may be explained by the differential phylogenetic information content of these genes. Owing to their slower evolutionary rate, nuclear markers are usually applied to investigate deeper phylogenetic relationships of aphids above genus level, whereas mitochondrial genes are more suitable for resolving relationships at the intra- and interspecific levels (Zhang & Qiao, Reference Zhang and Qiao2006). Ceratovacuna sp2. and C. silvestrii are likely to have diverged recently and hence nuclear genes may not yet have sorted completely.
Divergence times and implications for species diversification
The estimated divergence times (fig. 5) indicate that the diversification of most sampled species of Ceratovacuna and Pseudoregma occurred at approximately 10 mya in the late Tertiary (table 2). Short time intervals are evident among different species. Phylogenetic analyses indicate low resolution of the relationships among the three main groups of the sampled species of Ceratovacuna and Pseudoregma. Therefore, the divergence times and phylogenetic relationships suggest that Ceratovacuna and Pseudoregma may have undergone a rapid burst of speciation events in the late Tertiary. The radiation is likely to be related to the evolution of their secondary host plants. The secondary hosts of Ceratovacuna and Pseudoregma are plants of Poaceae. Owing to climatic cooling, herbaceous angiosperms such as Poaceae and Asteraceae appeared to be dominant in the middle to late Tertiary period (Heie, Reference Heie1996; von Dohlen & Moran, Reference von Dohlen and Moran2000). In fact, the divergence times in the present study are consistent with the results of Huang et al. (Reference Huang, Xiangyu, Ren, Zhang, Zhang and Qiao2012), which suggest that Cerataphidini most likely acquired its secondary hosts during the middle to late Tertiary period.
Except for C. lanigera, intraspecific genetic distances of mitochondrial COI were 0–1.6% in the two genera. However, the distances for C. lanigera samples were 0–4.5%. The genetic distances between one sample (No. 14674) from Jiuwan Mountain with other samples varied from 4.0 to 4.5%, which were much larger than the differences between the other species. Phylogenetic analyses based on mitochondrial genes all supported sample No. 14674 forming a separate clade from all other C. lanigera samples. Sample No. 14674 has no wax plates on the dorsal side of the head and pronotum, and the wax cells of the marginal wax plates on mesonotum and abdominal tergites are loosely arranged. All other samples of C. lanigera exhibit a pair of marginal wax plates on the head dorsum and pronotum, along with a much tighter arrangement of wax cells on the marginal wax plates of the mesonotum and abdominal tergites. Current evidences of genetic distances, phylogenetic topology, and morphological characteristic indicate that there may have been one cryptic species in C. lanigera. Geographical isolation of Jiuwan Mountain is likely to be a main factor for intraspecific divergence. However, further investigation based on more samples especially from the Jiuwan Mountain region is needed to test whether this putative cryptic species is valid.
The supplementary material for this article can be found at http://www.journals.cambridge.org/BER
Acknowledgements
We thank Jing Chen for providing some sequences of outgroups, and are grateful to Ms. Fendi Yang for preparing aphid slides. The work was supported by the National Science Funds for Distinguished Young Scientists (No. 31025024), National Natural Sciences Foundation of China (Nos. 31061160186 and 31272348), National Science Fund for Fostering Talents in Basic Research (No. J1210002), and the Ministry of Science and Technology of the People's Republic of China (MOST GRANT No. 2011FY120200).
