Introduction
Eriophyoid mites (Acari: Eriophyoidea) have a worldwide distribution. They exist in great variety and are highly host-specific plant feeders. Nearly 80% of them have been reported on a single host species, 95% on one host genus, and 99% on one host family (Skoracka et al., Reference Skoracka, Smith, Oldfield, Cristofaro and Amrine2010). They have only a weak ability to disperse. It seems that they cannot actively seek new hosts. Their only active mode of dispersal is walking from one plant to another, where the plants touch one another. Passive dispersal may be their main form of long-distance dispersal (Michalska et al., Reference Michalska, Skoracka, Navia and Amrine2010). This high host specificity and limited ability to disperse indicates that eriophyoid mites probably become divergent and undergo speciation when they shift to host plants and encounter new geographic conditions.
The discovery of cryptic diversity did not only help researchers involved in α-taxonomy, but also provided insights in the speciation, biodiversity, phylogeography, evolutionary theory, and ecological interactions of mites (Bickford et al., Reference Bickford, Lohman, Sodhi, Ng, Meier, Winker, Ingram and Das2007). With the current DNA-based technology, data can be used to estimate the genetic distances between taxa and even delimit species boundaries, thus allowing new species to be identified (Hebert et al., Reference Hebert, Cywinska, Ball and deWaard2003). Morphology-based taxonomy of eriophyoid mites has been used for more than one hundred years. This process was challenging when cryptic speciation events were not accompanied by any obvious morphological differentiation, but these events can be detected more effectively using DNA-based methods, particularly if genetic divergence is pronounced (Carew et al., Reference Carew, Schiffer, Umina, Weeks and Hoffmann2009; Skoracka & Dabert, Reference Skoracka and Dabert2010; Skoracka et al., Reference Skoracka, Kuczyński, Santos, Dabert, Szydło, Knihinick, Truol and Navia2012; Miller et al., Reference Miller, Skoracka, Navia, Mendonca, Szydło, Schultz, Smith, Truol and Hoffmann2013). Moreover, cryptic speciation can be supported by reproductive and host-adaptive data (Skoracka, Reference Skoracka2008; Skoracka et al., Reference Skoracka, Kuczyński, Szydło and Rector2013).
However, the results of cryptic diversity analyses may conflict when methods based on different types of data are used. Although speciation is not necessarily accompanied by morphological differentiation (Calcagno et al., Reference Calcagno, Bonhomme, Thomas, Singer and Bourguet2010; Henry & Wells, Reference Henry and Wells2010) quantitative morphological methods may reveal small differences between divergent populations. These differences may have arisen as adaptations to host morphology (Skoracka et al., Reference Skoracka, Kuczynski and Magowski2002). With molecular methods, mitochondrial and nuclear DNA sequence data often lead to different phylogenetic results due to differences in their evolutionary modes and substitution rates (e.g. Shaw, Reference Shaw2002; Leaché & Mulcahy, Reference Leaché and Mulcahy2007). Nuclear DNA is transmitted biparentally and mitochondrial DNA transmitted maternally as a single non-recombining block (Moore, Reference Moore1995). However, mitochondrial DNA experiences a higher mutation rate than nuclear DNA in Acari (Navajas & Boursot, Reference Navajas and Boursot2003). The incongruence between mitochondrial and nuclear DNA data can produce different indications and arguments in population genetic studies and may confuse the interpretation of cryptic diversity in eriophyoid mites.
Here, populations of an eriophyoid mite Tetra pinnatifidae, Xue et al., (2006) collected from different host plants and localities are surveyed. Morphometric data were obtained for 32 quantitative characters and used to explore morphological differentiation within T. pinnatifidae. For molecular analysis, one mitochondrial (16S) and one nuclear (28S D2 region) marker were surveyed. The objectives of the present study were to explore: (1) cryptic diversity in T. pinnatifidae and (2) the congruence of among different data sets in detecting cryptic diversity.
Materials and methods
Sampling
Mites from one host plant species and locality were regarded as one population for analytical purposes. From September 2004 to May 2012, field surveys were conducted in 12 areas within China. A total of 15 populations of T. pinnatifidae were collected with the aid of a hand-lens (30×). All mites were vagrant on the lower surfaces of the leaves. They infested nine different host plants, all Rosaceae. Most of the host plants were fruit trees or green trees that are widely distributed in China. No obvious damage to the host trees was observed. Six other Eriophyoidea species served as outgroups. Their morphologies were significantly different from those of T. pinnatifidae. Information on the samples is shown in table 1. All mites were preserved in 70% ethanol before use.
Morphometric analysis
The morphological terminology used here is consistent with that of Lindquist (Reference Lindquist, Lindquist, Sabelis and Bruin1996). The genus-level classification follows Amrine et al. (Reference Amrine, Stasny and Flechtmann2003). Slides were mounted using modified Berlese medium (Amrine & Manson, Reference Amrine, Manson, Lindquist, Sabelis and Bruin1996). Specimens were examined with a Zeiss A2 research microscope with phase contrast. For each population, 32 characters of 11–25 female specimens of good quality were measured. Abbrevations are given in fig. 1. These characters have been used in previous morphometric analyses of eriophyoid mites (Skoracka et al., Reference Skoracka, Kuczynski and Magowski2002; Magud et al., Reference Magud, Stanisavljević and Petanović2007; Vidović et al., Reference Vidović, Stanisavljević and Petanović2010; Wang et al., Reference Wang, Kuo, Jeng and Huang2011). An analysis of variance (ANOVA) was performed for each character to determine the differences between clades inferred from the DNA data (see below). The total variation among all specimens was assessed with a principal component analysis (PCA), which is a way of reducing the dimensionality of a problem. This method calculates a linear combination of principal components that represent 100% of the variation of multiple characters, and then iteratively calculates new combinations to explain any residual variation. This procedure does not assume any a priori groupings. Both ANOVA and PCA were performed in SPSS (SPSS Inc.). After that, two axes were constructed using the combination of PC1/PC2 and PC2/PC3.
DNA extraction, PCR, and cloning
Pooled samples of 5–10 mites were processed as one DNA sample for the extraction using a DNeasy Blood and Tissue Kit (Qiagen). First, the DNA was extracted from a single mite but this failed in most samples. Hence, the protocol was modified according to previous studies (Dabert et al., Reference Dabert, Ehrnsberger and Dabert2008). DNA samples were stored at −20 °C before use. PCR primers were designed as in previous studies. For 16S, samples of populations SYB and LA were amplified using primers LR-J-122887 and LR-N-13398 (Simon et al., Reference Simon, Frati, Beckenbach, Crespi, Liu and Flook1994). For others, LR-J-122887 and WCM16 s R (Carew et al., Reference Carew, Schiffer, Umina, Weeks and Hoffmann2009) were used. For 28S, primers D2R and ND2f (Campbell et al., Reference Campbell, Steffen-Campbell and Werren1993) were used. All PCR processes were performed in 25 μl, containing 12.5 μl of Super Taq Mix (Pudi Biotech Co., Ltd., Shanghai, China), 5 μl of template DNA, and 0.5 μM of each primer. Thermocycles consisted of an initial denaturation step at 96 °C for 5 min, followed by 30 cycles of denaturation at 95 °C for 30 s, annealing at 53 °C for 16S and 52 °C for 28S for 30 s, and extension at 72 °C for 30 s. The final extension took place at 72 °C for 5 min. An aliquot of 5 μL of the PCR amplification product was run on a 1.5% agarose gel stained with ethidium bromide and visualized under UV light to ensure the correct size of the amplified fragment.
Amplified fragments were purified using a TIANgel Midi Purification Kit (Tiangen Biotech Co., Ltd., Beijing, China). Then, the distinct single-band amplicons were cloned into a pEASY-T3 Cloning Vector (pEASY-T3 Cloning Kit, TransGen Biotech Co., Ltd., Beijing, China) and one positive clone was chosen at random. Then, the inserts were sequenced in ABI 3730XL by a biological reagent company, Meiji Biotech Co., Ltd., Shanghai.
Alignment and phylogenetic analyses
Nucleotide diversity (Pi) and haplotype diversity (Hd) were estimated using DnaSP 5.0 (Rozas et al., Reference Rozas, Sánchez-DelBarrio, Messeguer and Rozas2003). Haplotypes were compared online with sequences in GenBank using BLAST. Sequences were aligned using Clustal X (Thompson et al., Reference Thompson, Gibson, Plewniak, Jeanmougin and Higgins1997). Kimura-2-parameter (K2P) (Kimura, Reference Kimura1980) distances between pairs of sequences were calculated in MEGA 5 (Kumar et al., Reference Kumar, Nei, Dudley and Tamura2008). Four different tree reconstruction methods, neighbor-joining (NJ), maximum likelihood (ML), maximum parsimony (MP), and Bayesian inference (BI), were applied on both molecular markers. NJ was run in MEGA 5 using the K2P model with pairwise deletion of gaps. MP was implemented in PAUP 4.0b10 (Swofford, Reference Swofford2002). A heuristic search procedure was repeated 1000 times with randomized taxa additions and branch-swapping using 1000 bootstrap replicates. For ML and BI, jModeltest 0.1.1 was first used to identify the best-fit model of evolution (Posada, Reference Posada2008). These were found to be GTR+G for 16S and TIM1+G for 28S according to the Akaike Information Criterion (AIC). These models were implemented in PhyML 3.0 (Guindon et al., Reference Guindon, Delsuc, Dufayard and Gascuel2009) and MrBayes 3.1.2 (Ronquist & Huelsenbeck, Reference Ronquist and Huelsenbeck2003), respectively. For ML, bootstrap values were calculated for 100 replications. For BI, four Markov chains were run for 1,000,000 generations for 16S and 2,000,000 generations for 28S, with a sampling frequency of 100. The independent runs were considered to have converged when the standard deviation of the split frequencies value dropped below 0.01. Consensus trees were then calculated after omitting the first 25% trees as burnin. Finally, trees were rooted in Diptilomiopus ligustri and details were edited in TreeGraph 2 (Stover & Muller, Reference Stover and Muller2010).
Species delimitation
Species delimitation was inferred using two methods. First, the barcoding gap was tested to establish the threshold between intra- and inter-specific differences. The use of the barcoding gap is controversial, but it is still useful in the initial delimitation of species, even if it needs to be combined with other data to decrease the probability of misidentification (Wiemers & Fiedler, Reference Wiemers and Fiedler2007). Second, Pons generalized mixed Yule coalescent was used (GMYC, Pons et al., Reference Pons, Barraclough, Gomez-Zurita, Cardoso, Duran, Hazell, Kamoun, Sumlin and Vogler2006). The GMYC analysis divides a single locus gene tree into a portion in which a Yule speciation process affects the branch lengths and a portion in which there is a shift to a coalescent branching process. The boundary between these two portions is then used to define species. For the analysis, an ultrametric tree was required. First, all identical haplotypes were removed from the ML tree to render it fully dichotomous using TreeEdit (Rambaut & Charleston, Reference Rambaut and Charleston2002). Second, the ML tree was converted to an ultrametric one in r8s (Sanderson, Reference Sanderson2002) applying penalized likelihood to search for the best smoothing parameter and fixing the age of the root node at an arbitrary value of 1.0. These trees were then used as input trees in the search for the optimal threshold, identifying independently evolving entities using the package splits 1.0–11 (https://://r-forge.r-project.org/projects/splits) in R 2.13.1 (R Development Core Team, 2011).
Intra-clade analysis
To precisely estimate relationships among closely related haplotypes, a haplotype network was estimated (Bandelt et al., Reference Bandelt, Forster and Rohl1999). Once different populations had been confirmed to exist in the same clade, their haplotypes were estimated using median-joining as implemented in Network 4.6.1.0.
Results
Morphometrics
A total of 15 populations were identified as T. pinnatifidae according to traditional morphological taxonomy. Details of morphometric data and ANOVA are provided in Supplementary Material 1. ANOVA of all single characters showed significant morphometric differences between clades except for the leg characters: b, d, g, h and j. In the PCA analysis, the first two principal components accounted for 28.545% of the total variations and the first three accounted for 36.359%. Two axes were combined to determine whether any distinct groups could be recognized. As such, the combinations of PC1/PC2 and PC2/PC3 were considered. In general, there were no clear distinction among the populations. Nevertheless, the combination PC1/PC2 (fig. 2a) showed that populations XNA and XNB were grouped quite far from the others. Then LZ, LA, and CZ were grouped and separated from the two-member group containing JYA and JYB. However, these two groups were mixed into a larger group with GH, HZ, MD, SYA, and SYB. PC2/PC3 did not provide any further separation or grouping of populations (fig. 2b).
Phylogenetic analyses
A 341–343 bp (outgroups not included) of 16S rRNA and a 581–582 bp (outgroups not included) of the 28S rRNA D2 region were amplified and sequenced (GenBank accession numbers are given in Supplementary Material 2). In total 98 sequences of 16S and 96 sequences of 28S of T. pinnatifidae were included in the analyses. BLAST showed all sequences to have considerable similarity to the sequences of species of Eriophyoidea. All the target sequences except 28S of LZ and SYB were successfully amplified. Final alignments (including outgroups) were composed of 353 bp of 16S and 620 bp of 28S for phylogenetic analyses. Pi showed that 16S was more polymorphic than 28S (table 2. Pi of 16S was higher than that of 28S, except in Clade E, which was composed of only two populations, with only 12 sequences of 16S and 13 sequences of 28S). In contrast, degrees of Hd were similar for both markers.
The trees NJ, ML, MP, and BI trees based on 16S and 28S (figs 3–4) were mainly congruent. All trees indicated the monophyletic nature of T. pinnatifidae, showing either high bootstrap values or posterior probabilities (NJ: 100/ML: 79-96/MP: 100/BI: 1.00). All the trees of 16S showed the existence of at least eight clades (Clades A–H) (100/85-95/100/1.00), which was congruent with the 28S trees except for those lacking LZ and SYB (98-100/79-99/87-100/0.99-1.00). Two more encompassing clades were distinguished, too. These two clades were not always highly supported, including clade (C(D,E) (61/48/53/0.88 in 16S and 98/68/79/0.90 in 28S) and clade (H(F,G) (−/38/46/0.74 in 16S and 46/35/37/0.84 in 28S). Moreover, the topologies of the trees of the two makers were quite different. In 16S, the clade (C(D,E) was grouped with Clade B, though with very low support values (44/15/34/-), but in 28S, the clade (C(D,E) was grouped with clade (H(F,G) (−/23/44/0.63), yet again with very low support values but in line with the result of the morphometric analyses. The topological positions of XA between trees produced using the two makers were not congruent, either.
The following populations diverged in a manner corresponding to their hosts. Populations from Prunus spp. were in the same clade (LA and CZ in Clade E; JYA in clade (C,(D,E)). Similarly, populations from Malus baccata were in the same clade (MD, JG, SYA, GH, and HZ in Clade G and JYB in clade (H,(F,G)). Exceptional cases were also observed. XA and SYB, which both infest Grataegus pinnatifida were in two very distant clades (Clade A and H). In Clade B, three populations were from three different hosts: two from Pyrus and one from Malus. LZ in Clade C was not closely related to other populations from Pyrus species.
All populations or clades from the northeast showed close relationships in clade (H(F,G) (fig. 6), and two populations from eastern areas were placed in Clade C. However, populations from northwestern areas were separated into different clades, including XA in Clade A, XNA/XNB in Clade B, JYA in Clade D, JYB in Clade F, and LZ in Clade C. Population KM, which was collected from a southwestern area, clustered with two populations from the northwest, XNA and XNB. All three were joined in Clade B. The existence of three sympatric populations, XNA/XNB, JYA/JYB, and SYA/SYB, should be emphasized here. All of these pairs were found on different hosts. They could be separated into two cases. For XNA/XNB, the two populations were joined a clade. For JYA/JYB and SYA/SYB, the members of each pair were in different clades.
Species delimitation
Average K2P distances between different clades are shown in table 3. A barcoding gap was set using a histogram of K2P distances. K2P distances were first separated into intra- and inter-clade values. There is a well-defined gap in 16S (fig. 5a), ranging from 0.04 to 0.15. There is also a much less conspicuous barcoding gap in 28S (range 0.01–0.02) (fig. 5b). The gap between in- and outgroup clades in 28S was much larger (0.05–0.21), than that found in 16S (0.30–0.32).
The results of GMYC analyses differed considerably between the two markers. In 16S, each clade in the phylogenetic trees was defined as a separate species. In 28S, these clades were all as one species.
Intra-clade analysis
There were three clades composed of different populations. Among these, Clade B and Clade E contained populations from different hosts, and Clade B included sympatric populations. In the 16S haplotype network, haplotypes in Clades B and E were completely separated, corresponding to the different hosts and different main haplotypes (H5, H11, and H18 in Clade B and H38 and H41 in Clade E). In Clade G, only one main haplotype was found. This clade contained two separate populations, GH and SYA. For 28S, the results were quite different. No independent population was separated among three clades. Clades E and G shared one main haplotype (H41 in Clade E and H26 in Clade G), and Clade B had two (H4 and H12).
Discussion
Host-associated cryptic diversity in geographic populations
Host-associated diversity has been reported in previous studies of eriophyoid mites (Carew et al., Reference Carew, Schiffer, Umina, Weeks and Hoffmann2009; Skoracka & Dabert, Reference Skoracka and Dabert2010; Skoracka et al., Reference Skoracka, Kuczyński, Santos, Dabert, Szydło, Knihinick, Truol and Navia2012; Miller et al., Reference Miller, Skoracka, Navia, Mendonca, Szydło, Schultz, Smith, Truol and Hoffmann2013). This phenomenon probably also exists in the morphospecies T. pinnatifidae. First, the diversity of parapatric populations in northeastern and northwestern China showed different host distributions. Populations from northeastern China, which infest only one host (fig. 6), were joined in clade (H(F,G). This clade included JYB, a population from northwestern China that infests the same host as the other populations. Similarly, populations from eastern China and that infested two closely related Prunus hosts, were joined in one and the same clade. Populations in northwestern China, which were collected from six kinds of host plants, were more complex. They were distributed over six clades and did not constitute a single monophyletic taxon. Second, sympatric host-associated diversity, which is considered strong evidence of natural selection exerted by host plants (Cunningham et al., Reference Cunningham, Zalucki and West1999; Peccoud et al., Reference Peccoud, Ollivier, Plantegenest and Simon2009), was observed in two sympatric populations, XNA and XNB, of Clade B. These two populations were separated in the 16S haplotype network because of their different hosts. Three exceptional cases of host-associated diversity are dealt with here. First, unlike the sympatric case given above, two sympatric cases observed in the present study, JYA/JYB and SYA/SYB, did not form clades. Second, XA and SYB, which came from the same host species, were not monophyletic. Rather, they showed considerable genetic distance and were assigned to different clades. Third, KM in Clade B showed considerable geographic isolation and infested more sorts of hosts than the other two populations in the same clade. These cases illustrate the complex of dispersal and host transition history of these mites.
Implications of results obtained using different methods
In the present study, three kinds of data were used independently to evaluate the cryptic diversity of a host-associated mite. The results were not completely congruent and they are compared in fig. 7.
The morphometric data seemed to conceal the divergence shown in the molecular phylogeny. Only two large groups were well-defined, far fewer than the number of clades defined by 16S and 28S. Moreover, the separation of these two groups was not distinct because they had the overlapping parts. Morphological characters such as the number of empodium rays and prodorsal ornamentations were re-examined, but no significant differences were found between clades or populations. Nevertheless, in previous studies of eriophyoid mites, host-associated diversity was observed using similar morphological methods (Skoracka et al., Reference Skoracka, Kuczynski and Magowski2002; Magud et al., Reference Magud, Stanisavljević and Petanović2007; Vidović et al., Reference Vidović, Stanisavljević and Petanović2010). However, Skoracka et al. (Reference Skoracka, Kuczyński, Santos, Dabert, Szydło, Knihinick, Truol and Navia2012) failed to separate the host-associated wheat curl mite (WCM, Aceria tosichella) into cryptic lineages using this method. This suggests that the divergence is too recent for any morphological differentiation to have taken place.
The results obtained from the mitochondrial and nuclear DNA data usually differ due to the differences in their evolutionary modes and substitution rates (e.g. Shaw, Reference Shaw2002; Leaché & Mulcahy, Reference Leaché and Mulcahy2007). Although the molecular phylogenies and barcoding gaps of 16S and 28S all supported the existence of eight clades, these data also showed some degree of incongruence. First, the nucleotide diversity was lower in 28S than in 16S, indicating slower accumulation of mutations at this nuclear marker than in 16S. Second, GMYC of the 16S tree defined different clades as separate species, but in the 28S tree these clades constituted only one species. GMYC of 28S may conceal cases of cryptic speciation and evidence from other host-associated species of Anthocoptini mites is consistent with this conclusion. Nucleotide diversity of 16S was higher than that of nuclear ITS in populations of Aceria guerreronis (Navia et al., Reference Navia, Moraes, Roderick and Navajas2005). The mitochondrial sequences of Abacarus hystrix and Abacarus lolii, two host-associated mites whose reproductive isolation has been demonstrated to diverge by 0.20 (COI), despite their 28S D2 regions differed by only 0.002 (Skoracka, Reference Skoracka2008; Skoracka & Dabert, Reference Skoracka and Dabert2010). Similarly, the genetic divergence of host-associated WCM showed a range of 0.002–0.021 in 28S D2, which was much lower than the 0.073–0.235 observed in mitochondrial COI (Skoracka et al., Reference Skoracka, Kuczyński, Szydło and Rector2013). Miller et al. (Reference Miller, Skoracka, Navia, Mendonca, Szydło, Schultz, Smith, Truol and Hoffmann2013) detected a mean divergence of 0.069 for 16S vs 0.019 for ITS1 and 0.011 of adenine nucleotide translocase among worldwide populations of WCM. Moreover, Skoracka et al. (Reference Skoracka, Kuczyński, Szydło and Rector2013) concluded that host-associated WCM represents different cryptic lineages based on plant bioassays. Third, 16S was found to be more informative than 28S in uncovering intra-clade haplotype diversity and as such 16S was more useful in the exploration of the host-associated diversity that leads to speciation (Peccoud et al., Reference Peccoud, Ollivier, Plantegenest and Simon2009), while 28S was not. Host-associated populations shared common haplotypes in 28S, indicating incomplete lineage sorting, potential polyphagous ability, and hybridization opportunities (Funk & Omland, Reference Funk and Omland2003). High mitochondrial divergence, as well as low nuclear and morphological divergence, indicates very recent stages of cryptic diversity in this eriophyoid mite, an understanding that might not have been reached by focusing on a single data set only.
Perspectives in systematics of host-associated eriophyoid mites
Traditional morphology-based taxonomy has an inherent problem: speciation is not always accompanied by morphological differentiation. The taxonomy of the superfamily Eriophyoidea faces such a problem, and especially cryptic speciation of these mites is often in a very recent stage. Hence, morphological characters should no longer be considered as the only criterion for taxonomy, particularly not since current molecular technology provides a powerful tool to disclose hidden, taxonomically relevant, diversity. However, the differences of results obtained from different molecular markers or algorithms should be seriously considered. Eriophyoid mites are highly host-specific and cryptic lineages of eriophyoid mite species usually have specific host plants. Hence, host association can provide substantial systematic information. In addition, cross breeding experiments and host transition bioassays are required to better understand the systematics of the Eriophyoidea.
The supplementary materials for this article can be found at http://www.journals.cambridge.org/BER
Acknowledgements
This research was funded by the National Natural Science Foundation of China (31172132). The authors would like to thank Kai-Jun Zhang, Jing-Tao Sun, Xian-Ming Yang, and Xia Rong of Nanjing Agricultural University (NJAU) for reviewing an early draft of the manuscript and for providing constructive criticism. They would also like to thank Professor Xin-Hua Li of the NJAU Department of Botany for identifying host plants.