How many species are there in the subgenus Bursaphis (Hemiptera: Sternorrhyncha: Aphididae)? CO-I evidence

Species-level problems in the Aphis (Bursaphis) complex are reconsidered based on the partial sequences of the mitochondrial cox1 gene together with morphological and ecological data. This indicates that the American species A. oenotherae is a complex of four species (A. oenotherae, A. holoenotherae, A. costalis and A. neomexicana) and the taxonomic status of the species couples A. varians – A. manitobensis and A. epilobii – A. grossulariae require further clarification. Aphis sp. (USA: California, Oregon) of Blackman & Eastop (2006, p. 415) deserves the status of a species provided there is information on its host association and life cycle. Partial cox1 sequences might be misleading when used as standard DNA barcodes of aphid species of the subgenus Bursaphis.

It is common knowledge that the entire genus Aphis (including subgenus Bursaphis) is characterised by its morphological homogeneity and for a comprehensive taxonomic treatment of the genus it is necessary use information on the life cycles and host specificity of the species (Stroyan, 1984;Heie, 1986).In the case of the subgenus Bursaphis, life cycles and host specificity have been experimentally studied in only four of twelve species in this subgenus (Table 1).Thus it is possible that there are also cryptic species, the most recently being Aphis holoenotherae Rakauskas recorded in 2007 in Europe (Rakauskas, 2007).A. holoenotherae appeared morphologically inseparable from the Nearctic A. oenotherae, although the two species have different hosts and life cycles, and different distributions.For taxa that contain complexes of cryptic species and for which there is little or no information on their ecological characters, molecular data is of special interest (Coeur d'acier et al., 2007;Foottit et al., 2008;Mitchell, 2008).The aim of this paper is to reconsider species-level problems in the Aphis (Bursaphis) complex by using morphological and ecological data, and information on partial sequences of the mitochondrial cox1 gene.
species in this subgenus are designated ill-defined due to the lack of information on their life cycles and host specificity.
When establishing aphid clones, colonies were initiated from a single parthenogenetic female and are expected to contain DNA from genetically homogenous individuals.Aphid rearing provided important information on aphid life cycles and host specificity of the aphid lineages used for DNA extraction.Detailed information on the aphid cloning methods used is already published (Rakauskas, 1993).Field samples were taken from a single colony.Aphids were either frozen at -80°C or stored in 96% ethanol.

DNA extraction, amplification and sequencing
Total genomic DNA was extracted from a single individual using the DNeasy Blood & Tissue kit (Qiagen), which involved at least a 2 h digestion of tissue with proteinase K.
The target sequence was a partial fragment from the mitochondrial gene encoding Cytochrome Oxidase subunit I (CO-I).The 620 bp region was PCR-amplified using previously published primers (Tur inavi ien et al., 2006).PCR amplification was carried out in a thermal cycler (Eppendorf) in 50 µl volumes containing 1-2 µl genomic DNA, 5 µl of each primer (10 µM), 5 µl of PCR-reaction buffer, 5 µl of dNTP mix (2 mM each), 4-8 µM of 25 mM MgCl2 and 1,25 U of Taq polymerase (recombinant) (5U/µl) and ddH2O to 50 µl.The cycling parameters were as follows: denaturizing at 95°C for 10 min (1 cycle), denaturizing at 95°C for 30˝, annealing at 49°C for 30˝ and extension at 72°C for 30˝ (32 cycles in total), and a final extension for 5 min (1 cycle).PCR products were subjected to electrophoresis on 2% Top Vision LE GQ agarose (Fermentas, Lithuania), stained with ethidium bromide and sized against a Gene Ruler 100 bp DNA ladder (Fermentas, Lithuania) under UV light.PCR products were cleaned using the QIAquick PCR Purification Kit (Qiagen).Cyclic sequencing was performed at the Institute of Biotechnology (Lithuania) using a BigDye® Terminator v3.1 Cycle Sequencing Kit and products sequenced using a 3130 × l Automated Sequencer (Applied Biosystems).The amplification primers were also used as sequencing primers.DNA sequences for each specimen were confirmed with both sense and anti-sense strands and aligned in the Bio-Edit Sequence Alignment Editor (Hall, 1999).Sequences were tested for stop codons and none were found.

DNA based species delimitation
Of the numerous methods available to delimit species boundaries (Pons et al., 2006, Vogler & Monaghan, 2007) the following were used in this study.
1.The DNA barcoding method (Hebert et al., 2003) that has been successfully used to identify Linnaean species and indicate presumably cryptic species (Footit et al., 2008).It is based on the pair wise nucleotide sequence divergences measured as distances: within-species sample divergences normally appear to be much lower than inter-specific ones.Pair wise nucleotide sequence divergences were calculated using Kimura 2-parameter (K2P) model of base substitution (Kimura, 1980), which has been evaluated as the best model for species level analysis at low distances (Hebert et al., 2003).The species Identifier tool, as implemented in TaxonDNA (Meier et al., 2006), evaluates uncorrected genetic distances and the results are similar to those obtained using the K2P model (Meier et al., 2008).The smallest inter-specific distances were used to predict cryptic species complexes in accordance with the suggestion of  (Müller, 1974;Martin, 1982;Stroyan, 1984;Cook, 1984;Heie, 1986;Holman, 1990;Remaudière, 1993;Rakauskas, 1996Rakauskas, , 1998Rakauskas, , 2003Rakauskas, , 2007Rakauskas, , 2008;;Blackman & Eastop, 2006).Life cycle characteristics confirmed by experimental studies (Patch, 1927;Rakauskas, 1993Rakauskas, , 2007) ) are given in bold.Asterisked species of Aphis (Bursaphis) are designated well-defined species in the text based on their specific morphology and at least indirect evidence on the specificity of their life cycles and host associations.Meier et al. (2008).Standardization of inter-specific sequence differences were avoided because inter-specific sequence differences appear rather variable when comparing different taxa.Pair wise nucleotide sequence divergences were calculated using MEGA version 4 (Tamura et al., 2007) and neighbour-joining (NJ) analysis (Saitou & Nei, 1987) was used to represent these distances as a tree.

AschS19
2. Character based tree building methods -maximum parsimony (MP) and maximum likelihood (ML).Different clustering methods consider monophyletic groups of specimens as distinct species.Clades were suggested to be distinct species when they received strong support (> 50%) in all phylogenetic reconstructions.MP (Swofford et al., 1996) analysis was computed using MEGA 4 (Tamura et al., 2007).Heuristic searches were carried out with ten random taxon addition replicates.Bootstrap values were generated from 1000 replicates, each with ten randomaddition heuristic searches.Maximum likelihood (ML) analysis was conducted with PHYML (Guindon & Gascuel, 2003) and the best fitting model for the evolution was determined using the Akaike information criterion, as implemented in jModeltest 0.1.1 (Posada, 2008).The general time reversible model with a proportion of invariant sites and a gamma distribution (GTR + I + G) was used.Bootstrapping of the ML analysis was based on 100 replicates.Three out-group species were used in order to establish the rooting of the Aphis clade: Protaphis carlinae (Aphidini: Aphidina), which represents a closely related genus, Schizaphis rotundiventris (Aphidini: Rhopalosiphina, GenBank Accesion N AF220511), which represents a sister sub-tribe and Acyrthosiphon pisum (Macrosiphini: Macrosiphina, GenBank Accesion N AF077776), which is a more distantly related species.
3. Parsimony networks and implemented 95% parsimony connection limit (Hart & Sunday, 2007).Samples that stuck together in a single haplotype network were suggested to be con-specific.Intra-specific phylogenetic relationships among mtDNA haplotypes were inferred based on a statistical parsimony method (Templeton et al., 1992).Statistical parsimony is more convenient when investigating closely related haplotypes with a low number of substitutions.Haplotype networks were constructed using the software TCS 1.13 (Clement et al., 2000).
The sequence data for all species were submitted to the Gen-Bank under following Accession Nos HM2 45794 to HM2 45856.

RESULTS AND DISCUSSION
The alignment of partial CO-I sequences contained 623 sites, of which 125 were variable and 91 appeared parsimony informative.The sequences were heavily biased towards A and T nucleotides.The average base composition was A = 34.9%,C = 13.9%,G = 12.3% and T = 38.9%.The overall transition/transversion ratio was R = 4.113 for all sites.
The maximum parsimony (MP) analysis of partial CO-I sequences resulted in 220 equally parsimonious trees (length = 297, CI = 0.58, RI = 0.87).The strict consensus MP tree is shown in Fig. 1.A similar tree was produced by a ML analysis using model GTR + I + G. Topologies inferred in strict consensus MP and ML trees with the best likelihood scores appeared congruent and produced the same strongly supported clades.Therefore, both MP and ML bootstrap values are given at the respective nodes in Fig. 1.The neighbour joining (NJ) tree based on Kimura 2-parameter distances (Fig. 2) revealed a somewhat different topology.Noticeably, nodes, that are absent in the NJ tree had low supports in the MP and ML trees.CO-I partial sequences indicate the subgenus Bursaphis is a monophyletic group, although ML/MP bootstrap support values (62/78 respectively) are not very strong.Earlier studies, based on traditional (Stroyan, 1984;Heie, 1986;Remaudière, 1993) and molecular (Coeur d'acier et al., 2007;Foottit et al., 2008) data have also suggested the monophyly of Aphis (Bursaphis).The smallest pair wise inter-specific divergence was between A. holoenotherae and A. oenotherae s. str., followed by the couple A. varians and A. manitobensis (Table 4).The pair wise within-species sample distances for welldefined species of the subgenus Aphis (Bursaphis) are given in Table 5.Their minimal within-species sequence divergence ranges from 0.00 to 0.16% (A.holoenotherae).The maximum within-species sample distance in well-defined species reaches 0.97% in A. varians.Of the ill-defined species, the A. oenotherae-like samples are the most numerous (8 samples altogether, see Table 2).The pair wise within-group sample distances for A. oenotherae-like samples collected on different hosts are given in Table 6, together with the Aphis sp.samples collected from similar host plants.Table 6 indicates the distances between A. oenotherae-like samples collected on similar host plants are comparable with inter-specific distances between well-defined species.For the present, A. oenotherae is reported to be a holocyclic species alternating between currants and species of herbaceous plants of the family Onagraceae (Remaudière, 1993;Blackman & Eastop, 2006).This is based on the experimental study of the life cycle of A. sanborni (Patch, 1927), which was subsequently synonymised with A. neomexicana by Robinson & Rojanavongse (1976).Patch (1927) reports A. sanborni as a holocyclic species alternating between Ribes and Epilobium spp.(E.adenocaulon, E. lineare, E. coloratum).The life-cycle of American populations of A. oenotherae, synonymised with A. neomexicana by Remaudière (1993), has never been checked.In the MP/ML and NJ trees (Figs 1-2), A. oenotherae-like samples occur at different nodes.Four of them are located in clade A, together with samples of A. costalis and new unnamed species mentioned as Aphis sp.(New Zealand) in the key of Blackman & Eastop (2006) (Fig. 1).These samples emerge in the same nest of the haplotype parsimony network, which indicates all them belong to the same species.Of the three remaining A. oenotherae-like samples, one (Aneo14) collected from currants in Arizona, might be a representative of a separate species (Fig. 1).In our opinion, this justifies the restitution of A. neomexicana (W.P. Cockerell & T.D. A. Cockerell, 1901) as a valid species.Finally, two A. oenotherae-like samples emerging close to the A. holoenotherae node in clade D2 (Aoen19 and Aoen71, Figs 1-2) might represent one more separate species, A. oenotherae Oestlund, 1887 in a narrow sense.Thus, based on the partial CO-I sequences, the A. oenotherae-like complex (together with A. costalis) might include four cryptic species: A. oenotherae, A. holoenotherae, A. costalis and A. neomexicana.This supports the doubts concerning the synonymy of A. neomexicana and A. oenotherae (Blackman & Eastop, 2000, 2006;Rakauskas, 2000).
The remaining samples tend to form a sister group to clade A and include European, East Palaearctic and North American material (Fig. 1).European species form two well supported separate clades consisting of A. grossulariae, A. schneideri and A. epilobii samples (clade B, 68/90% bootstrap support for ML/MP respectively), and A. epilobiaria samples (clade C, 93/97%).Generally, CO-I partial sequences correspond well with the available information on the ecological and morphological features of the European species of the subgenus Bursaphis (Stroyan, 1984;Heie, 1986).However, the European A. epilobii samples remain poorly resolved in the MP/ML and NJ trees (Figs 1-2).Noticeably, these samples emerge together with those of A. grossulariae and A. schneideri in the same nest of the haplotype parsimony network, which suggests they all belong to the same species.A. epilobii is reported as a holocyclic species, oligophagous on Epilobium herbs, mostly E. montanum (Stroyan, 1984;Heie, 1986;Holman, 2009).Some of the samples collected on Epilobium spp.for this study were difficult to identify using the morphological characters available in the recent widely accepted key (Blackman & Eastop, 2006), because of a mixture of or intermediate characters.For example, the sample of Aphis (Bursaphis) (Aepi43) collected in Kyiv (Ukraine) from Epilobium sp. had the key morphological characters closer to A. epilobii, but the body colour (light green) of A. grossulariae.This sample appeared inside the well supported A. grossulariae node (Figs 1-2).Sample (Aepi90) collected in D kštai (Lithuania) had the morphology of A. epilobii, but appeared in the A. grossulariae node.Sample (Aepi91) from Rabi (Czech Republic) also had morphological characters of A. epilobii, although appeared closer to the A. schneideri node.This indicates the need for a rigorous study of host specificity and life cycles, and a subsequent morphological and molecular analysis of A. epilobii in Europe.
The present data, together with that published by Tur inavi ien et al. (2006), confirm that the European species A. grossulariae and A. schneideri are very similar.The A. grossulariae and A. schneideri samples used in this study  form separate clades in MP/ML and NJ trees, but they appear in the same nest of the haplotype parsimony network.This supports an earlier hypothesis that A. grossulariae and A. schneideri hybridize naturally (Rakauskas, 1999(Rakauskas, , 2003)).Experimental crosses between these species resulted in hybrid clones that differ in their morphology, host specificity and life cycle characteristics, namely, some hybrid clones have the morphology of A. schneideri but the life cycle of A. grossulariae (Rakauskas, 1999).Subsequently, hybrid morphotypes were found in field samples (Rakauskas, 2003).In this study, a sample (Asch75) collected from black currants at Salaspils (Latvia) had the morphology of A. schneideri, but appeared in the A. grossulariae clade (Figs 1-2).This might be due to the hybrid origin of the sample and certain maternal effects.Noticeably, A. epilobii samples are located inside the A. grossulariae + A. schneideri clade.A. epilobii sample Aepi90 groups together with the above mentioned Asch75 sample of A. schneideri in the A. grossulariae cluster (Figs 1-2).This suggests that some of the A. epilobii samples were also of hybrid origin (A.grossulariae × A. schneideri).Life-cycle and host specificity data (Stroyan, 1984;Rakauskas, 1993Rakauskas, , 1998)), however, indicate that A. grossulariae and A. schneideri are species.
Due to the lack of material, it is still unclear whether A. popovi is a separate species.A. popovi was described from aphid material collected on currants at Yakutia (Russian Siberia).It is morphologically similar to A. schneideri (Rakauskas, 1996), but its life cycle and host specificity are unknown.
Several North American (A.manitobensis, A. varians, A. oenotherae s. str.) and Palaearctic (A. holoenotherae) samples tend to form a separate clade in the MP and ML trees, but it has low bootstrap support (D, 58/42% for ML/MP, respectively, Fig. 1) and is not apparent in the NJ tree (Fig. 2).In this group, the A. varians and A. manitobensis samples nest together in the haplotype parsimony network and form clearly separated clades in MP/ML and NJ trees (Figs 1-2, D1, 54/97 and 98%, respectively).Foottit et al. (2008) report that A. manitobensis samples differ from all A. varians haplotypes, although the distance between A. manitobensis and A. varians samples is about the same as that between the two A. varians clusters.When these authors were asked to enlarge on this (Foottit & Maw, 2009;pers. comm.)they replied that "given the relative uniformity of the transcontinental A. varians cluster, the observed difference between that group and the few A. manitobensis samples available may prove to be constant despite the small distance".This is in accord with the results of this study.The case of A. manitobensis and A. varians is very similar to that of A. grossulariae and A. schneideri (see above), with the ranges of pair wise inter-specific sample divergences comparable for both couples of species (0.48-1.30 and 0.81-1.63,respectively; Table 4).In addition the morphological and ecological evidence indicates that the four above mentioned species are good species (Remaudière, 1993;Remaudière & Remaudière, 1997).In contrast to the case of A. grossulariae-A.schneideri, proper evidence of the principal ecological niche characteristics (host specificity and life cycles) is lacking for A. manitobensis.Life cycle and host specificity of A. varians is holocyclic, with the species host alternating between currants and Epilobium herbs (Patch, 1927), and only indirect information indicating that A. manitobensis is monoecious on currants (Robinson & Rojanavongse, 1976).Therefore, it is necessary to rear both species and carry out morphological and molecular analyses of the lineages in order to resolve this problem.
A highly supported clade D2 (100/99 and 99% for ML/MP and NJ respectively, Figs 1-2) contains both Palaearctic and American samples.Palaearctic samples of A. holoenotherae form a uniform group, nearly all iden-tical in their sequences; with only one differing in one nucleotide.The A. holoenotherae node appears in sistergroup relationships with two A. oenotherae s. str.samples from Wyoming and Washington in the ML/MP and NJ trees (Figs 1-2).The range of pair wise inter-specific sample divergences between A. holoenotherae and A. oenotherae s. str. is the lowest among the couples of welldefined species analysed (0.48-0.81,Table 4).American samples differ from the rest of the A. holoenotherae node in three (Aoen71) and four (Aoen19) substitutions and therefore nest together in the haplotype parsimony network.Moreover, the partial CO-I sequences of the Canadian sample, Aoen97, are identical with most of the Palaearctic A. holoenotherae samples.It was not possible to morphologically identify this sample, because there were only immature nymphs in the subsample that R. Foottit kindly supplied.In reply to our request for further information R. Foottit & E. Maw (2010, pers. comm.)stated that "there is only a single adult aptera in their collection … Terminal process lengths (left/right) are 0.216/0.221mm.Measurements (left/right) for 3 specimens of instar 4 alates: 0.187/0.187,0.187/0.171,0.210/0.200".The key of Rakauskas (2008) places apterae with these characteristics in A. oenotherae, but the lengths of nymphs are more similar to those of A. holoenotherae.The case of A. oenotherae s. str.and A. holoenotherae is comparable to that of A. manitobensis and A. varians (see above). A. holoenotherae is holocyclic and monoecious on Oenothera spp. in Lithuania and Poland (Rakauskas, 2007), but there is no experimentally based information on the life cycle and host specificity of A. oenotherae s. str. in North America.In addition, both of these presumed species seem to have separate distributions (Palaearctic versus Nearctic for A. holoenotherae and A. oenotherae s. str.respectively).In order to resolve this problem it is necessary to determine the life cycle and host specificity of A. oenotherae s. str.and subject it to morphological and molecular analyses.
Sample Aepi23, collected from Epilobium in Washington State, which ran to Aphis sp.California in the key of Blackman & Eastop (2006) is not included in any grouping in MP and ML cladograms (Fig. 1), but forms a separate node together with the A. oenotherae-like sample Aneo14 in the NJ tree (Fig. 2).This tends to confirm the presupposition of Blackman & Eastop (2006) that there is a new species inhabiting Epilobium in western USA.In contrast, sample Aphis83 collected from Epilobium in Washington State (Moses Lake), which ran to the Aphis sp.(New Zealand) in the same key appears to be a member of clade A (Figs 1-2), which consists of A. costalis and A. oenotherae-like samples.The sequences of sample Aphis83 are identical with those of A. oenotherae-like sample Aoen16 from New Mexico.
In conclusion, the exact number of species in the subgenus Aphis (Bursaphis) still needs further clarification.This is due to the lack of reliable information on the host specificity and life cycles of eight of the twelve species in this subgenus.Based on the results of the present study, cox1 gene partial sequences indicate that the following Aphis (Bursaphis) species need to be revised.The A. oenotherae-like complex seems to include four or more species.Rigorous life cycle and host specificity studies on the species in this complex, together with morphological and molecular analysis of the same material, might lead to the resurrection of A. neomexicana.The same also applies to the A. varians -A.manitobensis and A. epilobii -A.grossulariae species couples.Aphis sp.California in the key of Blackman & Eastop (2006), which corresponds morphologically to our sample Aepi23 collected from Epilobium in Washington State, awaits a proper definition based on a combined ecological, morphological and molecular study.Three out of twelve currently recognized species in the subgenus Aphis (Bursaphis) (A.solitaria, A. fluvialis, A. popovi) were not included in the present study due to the lack of material.A. popovi might be of special interest due to its similarity to A. schneideri (see Rakauskas, 1996 for details).
The data presented here support the statement of Coeur d 'acier et al. (2007) that "mitochondrial DNA does not allow the differentiation of species that are difficult to identify morphologically", but clearly indicates the existence of complexes of cryptic species and the need to synonymize morphologically distinct "species" (see above).Distance and tree-building methods of species delimitation appeared congruent when dealing with closely related species of the subgenus Bursaphis.Designation of any threshold for distance based methods can be accepted only in cases when ecological data (life cycle and host specificity) for the specimens are also available, because even the smallest inter-specific distances might correctly reflect different species (e.g. the case of A. grossulariae -A.schneideri, see above).Once a species is assigned to a monophyletic group, only those that are well supported clades deserve the nomination of distinct species (Meier et al., 2006).Congruence of topologies produced by different clustering methods is deemed to be important evidence when delimiting species, despite the small inter-specific pair wise nucleotide sequence divergences.Recently, mitochondrial DNA analysis was used to indicate the possible existence of cryptic species in the genus Toxoptera (Wang & Qiao, 2009).We have reservations about the claim that partial sequence of cox1 gene (658-bp fragment from the 5' region) can be "adopted as the standard DNA barcode region for animal life" (Foottit et al., 2008).Apart from the problem of whether it can be used for all animal life, the use of partial CO-I sequences as the standard barcode for aphids is limited.A. manitobensis and A. varians are good species (see above), so partial CO-I sequences might appear misleading when applied as standard barcodes in this case.Aphid species Macrosiphum rosae (Linnaeus, 1758) and Macrosiphum knautiae Holman, 1972 differ in their ecology and morphology, but their cox1 gene partial sequences are identical [see Tur inavi ien & Rakauskas (2009) for details].It is likely that molecular data has the same limitations as any other character used in classification, including the most popular morphological ones.Therefore, ecological specificity is the most important feature when trying to classify living things, followed by morphological and molecular analyses [see Rakauskas (2009) for wider discussion].The experience of the authors leads them to support the opinion that DNA barcoding is not a replacement for morphology-(and ecology!-) based taxonomy (Coeur d'acier et al., 2007;Mitchell, 2008;Žurovcová et al., 2010;Tan et al., 2010).

CONCLUSIONS
Summarising, the results of the present study, based on an analysis of partial sequences of the mitochondrial cox1 gene: 1. Partial CO-I sequences might be misleading when used as standard DNA barcodes for aphid species of the subgenus Bursaphis.In addition to DNA studies, it is necessary to undertake ecological, morphological and molecular studies in order to determine the number of species in this subgenus.
2. A. oenotherae is likely to be a complex of cryptic species.In addition to A. oenotherae; A. holoenotherae, A. costalis and A. neomexicana should be validated.
3. Aphis sp.(USA: California, Oregon) of Blackman & Eastop (2006, p. 415) deserves the status of a good species once there is information on its host association and life cycle.

Fig. 1 .
Fig. 1.Phylogeny of the genus Aphis based on a region of mtCO-I.The maximum parsimony strict consensus topology is shown.Bootstrap values at branches are ML (left) and MP (right).Only those exceeding 50% are indicated.(-) indicates nodes that collapsed in ML consensus due to low support (bootstrap values less than 50%).Sample abbreviations as in Tables 2-3.-new, the yet undescribed species Aphis sp.New Zealand in Blackman & Eastop (2006: p. 415); -new, the yet undescribed species Aphis sp.California in Blackman & Eastop (2006: p. 415); -A.oenotherae-like samples collected on Epilobium or Oenothera; N -A.oenotherae-like samples collected on Ribes; -A.costalis.

TABLE 1 .
Species of the subgenus Bursaphis McVicar Baker in the World

TABLE 3 .
Clonal material of the genus Aphis L. studied.Morph abbreviations the same as in Table2; fx -fundatrix.

TABLE 6 .
Range of pair wise within-group sample divergences (Kimura 2-parameter model) for A. oenotherae-like and Aphis sp.samples collected from different host plants (number of samples in parentheses).