Mitochondrial COI and morphological evidence for host specificity of the black cherry aphids Myzus cerasi (Fabricius, 1775) collected from different cherry tree species in Europe (Hemiptera, Aphididae)

Abstract Partial sequences of the mitochondrial COI gene of forty eight European and two Turkish population samples of Myzus cerasi from different winter hosts (Prunus spp.) were subjected to phylogenetic analyses. The analysed M. cerasi samples emerged as paraphyletic relative to a Myzus borealis sample used as an out-group, and formed two major clades in neighbor joining, maximum parsimony, maximum likelihood and Bayesian inference trees, corresponding to subspecies living specifically on Prunus avium and P. cerasus. Multivariate discriminant analysis (method of canonical variates) was applied to find out if morphological variation of samples correlated with mitochondrial COI and host plant information. Mean scores on the first two canonical variables clustered samples fully in accordance with their COI haplotypes and host plants confirming the existence of two morphologically similar winter host - specific subspecies of M. cerasi in Europe. No single morphological character enabled satisfactory discrimination between apterous viviparous females of the two subspecies. A three-character linear discriminant function enabled 92.37% correct identification of apterous viviparous females of M. cerasi cerasi (n = 118) and 93.64% of M. cerasi pruniavium (n = 110). A key for the morphological identification of the two subspecies is presented and their taxonomic status is discussed.


Introduction
Black cherry aphid Myzus cerasi (Fabricius, 1775) is reported to be a serious pest of cherries all over the world (Barbagallo et al. 1997, Blackman and Eastop 2000, Holman 2009) and its morphology, life cycle, host specificity and potential harmfulness have therefore been the subject of intensive studies (Ross 1918, Wimshurst 1925, Pokrovskyj 1932, Vereshchagina 1966, Fomicheva 1967, Karczewska 1970, Rakauskas 1984, Gruppe 1990, Cichocka and Goszczynski 2004). Nevertheless, the species level classification of black cherry aphids has not been satisfactorily resolved. The black shiny aphids inhabiting cherry trees were originally described as a single species, Myzus cerasi (Fabricius, 1775), but European populations inhabiting sweet cherry (Prunus avium) were later separated as Myzus pruniavium Börner, 1926. Börner's species has been accepted by some (Börner 1943, 1952, Heinze 1961), but synonymised with cerasi by others (Miyazaki 1971, Eastop and Hille Ris Lambers 1976, Remaudière and Remaudière 1997, while others have treated it as a subspecies of cerasi (Shaposhnikov 1972, Favret 2014. Differences in host specificity of the two taxa have been documented in experimental transfer studies (Karczewska 1970, Dahl 1968, Vereshchagina 1966, showing M. cerasi cerasi as heteroecious, alternating between cherries (both Prunus cerasus and P. avium) and herbaceous hosts (Galium, Euphrasia, Odontites, Veronica). M. cerasi pruniavium Börner, 1926 differs from the nominative subspecies in having P. avium as the only winter host (it is incapable of living on P. cerasus), and also has a somewhat different summer host list. Enzyme electrophoretic studies indicated a reduced gene exchange between the subspecies (Gruppe 1988a), perhaps due to differences in the phloem sap of sour and sweet cherry causing divergent selection (Gruppe 1989). Morphological characters for discrimination between the two taxa have been suggested (Heinze 1961), but none of these enabled satisfactory separation between viviparous females of M. c. cerasi and M. c. pruniavium when applied to independent aphid material (Dahl 1968, Gruppe 1988b. The taxonomic status of the black cherry aphid in the western Palaearctic thus remains unclear (Blackman and Eastop 2000, Lampel and Meier 2007, Holman 2009, Osiadacz and Halaj 2010. They are possibly members of a complex of cryptic aphid species that includes M. cerasi umefoliae (Shinji, 1924), which overwinters on P. mume in Japan with Artemisia as its summer host (Miyazaki 1971). In northern Europe there is also another potentially distinct taxon that lives all year without host alternation on herbaceous plants (Dahl 1968), and to which the name M. cerasi veronicae (Walker, 1848) may be applicable.
A similar species complex is found in the mealy plum aphid, Hyalopterus spp. (Lozier et al. 2008). Separation of three Hyalopterus species was eventually justified by their distinctness at molecular level (Lozier et al. 2008), but they still remain difficult to separate by their morphological characters (Basky and Szalay-Marszό 1987, Blackman and Eastop 1994, including those used in the most recent identification keys (Lozier et al. 2008, Rakauskas et al. 2013. A similar recent case is that of birch-and oak-inhabiting aphid species of the genus Stomaphis Walker, 1870 (Depa et al. 2012).
The M. cerasi complex has not yet been the subject of detailed molecular study, although certain DNA sequences have been isolated , Valenzuela et al. 2007, Clements et al. 2000, 2000a. Preliminary data on the partial sequences of the mitochondrial COI gene support the existence of subspecies of M. cerasi (Voronova et al. 2011). The aim of this study was to clarify the taxonomic status of the host-alternating taxa in the M. cerasi complex by a combined study of partial sequences of the mitochondrial COI gene and morphological characters of the European samples collected from different species of cherries.

Material studied
Fifty population samples of apterous viviparae of black cherry aphids from nine European countries and Turkey were collected in 2004-2013, mostly from P. cerasus, P. avium, but also from three other Prunus species (Table 1). Twenty samples were used for morphology -based stepwise discriminant analysis. The remaining 30 samples were used for subsequent evaluation of the derived discrimination functions. Samples of Myzus borealis Ossiannilsson, 1959 and Myzus persicae (Sulzer, 1776) (GenBank Accession No AB506741) were used as out-groups for the phylogenetic analyses.

DNA extraction, PCR amplification and sequencing
For molecular analysis, a single aphid individual from one sampled plant was considered as a unique sample. Total genomic DNA was extracted from a single aphid using the DNeasy Blood & Tissue kit (Qiagen), which involved at least a 2 h digestion of tissue with proteinase K. Partial sequences of mitochondrial COI gene were PCRamplified using earlier published primers (Turčinavičienė et al. 2006). PCR amplification was carried out in a thermal cycler (Eppendorf ) in 50 µl volumes containing 2 µl genomic DNA, 5 µl of each primer (10 µM), 5 µl of PCR-reaction buffer, 5 µl of dNTP mix (2mM each), 4-8 µl of 25mM MgCl 2 and 1.25 U of AmpliTaq Gold 360 polymerase (5U/µl) and ddH 2 O 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" (37 cycles in total), and a final extension for 5 min (1 cycle). PCR products were subjected to electrophoresis on 2% TopVision agarose (Fermentas, Lithuania), stained with GelRed and sized against a MassRuler Low Range DNA ladder (Fermentas, Lithuania) under UV light. PCR products were purified and sequenced at Macrogen Europe (Amsterdam, the Netherlands) and Institute of Biotechnology of the Vilnius University (Vilnius, Lithuania). 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 table 1. Aphid material used in the present study. Samples used for morphology-based discriminant analysis are given in bold.
the BioEdit Sequence Alignment Editor (Hall 1999). Partial sequences were tested for stop codons and none were found. The sequence data have been submitted to the GenBank, accession numbers are given in Table 1.

Analysis of DNA sequences
In addition to the sequences from 50 samples of M. cerasi, COI sequences of M. borealis from subgenus Myzus sensu stricto (the same subgenus as M. cerasi) and M. persicae from subgenus Nectarosiphon Schouteden, 1901 were selected as out-groups for the phylogenetic analyses, which included neighbor joining (NJ), maximum parsimony (MP), maximum likelihood (ML) and Bayesian inference in phylogeny (BI). NJ, MP and ML analyses were performed using MEGA 5 (Tamura et al. 2011). For NJ analysis Kimura 2-parameter (K2P) model of base substitution was used. Bootstrap values for NJ, MP and ML trees were generated from 1000 replicates. For ML analysis Tamura 3-parameter model with Gamma distribution (T92+G) was selected by MEGA 5 model selection option (Tamura et al. 2011). Bayesian analysis was conducted in MrBayes 3.2.1 (Ronquist and Huelsenbeck 2003) using Hasegawa-Kishino-Yano model with Invariable sites and Gamma distribution (HKY+I+G), which was selected by jModeltest (Posada 2008). Four simultaneous chains, 3 heated and 1 "cold", were run for 3 000 000 generations with tree sampling every 1000 generations. The topologies obtained by NJ, MP, ML and BI were similar, so only ML tree is shown with values of NJ/MP and ML/BI bootstrap support and posterior probabilities over 50% indicated above and below branches respectively. Statistical parsimony haplotype networks were constructed for samples of M. cerasi and M. borealis using TCS v 1.21 (Clement et al. 2000).

Morphological study and discrimination analysis
Samples representing different clades in the molecular tree and haplotype network were used for stepwise discriminant analysis followed by canonical analysis: 10 sam-ples from sour cherry, P. cerasus, and 10 samples from sweet cherry, P. avium (shown in bold in Table 1).
Measurements of the slide-mounted apterous viviparous females were performed by means of interactive measurement system Micro-Image (Olympus Optical Co. GmbH). STATISTICA 8 version software (Statsoft 2007) was exploited for data analysis. Pearson's correlation coefficients were calculated to evaluate the correlation of morphometric characters with body length. Characters with strong (| r | ≥ 0.70) statistically significant (p<0.05) correlation with body length were removed from the further analysis: BL (r=1.00), F3L (r=0.83), T3L (r=0.82), A2L (r=0.7), A3L (r=0.75), A4L (r=0.71), A5L (r=0.7). The remaining 12 characters were used for the forward stepwise discriminant analysis followed by canonical analysis, with sample collection number as the grouping variable, thus excluding information about the host plant from the analysis. Mean canonical scores of the first two canonical variables were represented as bivariate scatter plots, in order to show any clustering of samples.
Morphological characters that contributed most to canonical discrimination functions were evaluated as having potential for separation of taxa. An identification key was constructed based on these discrimination functions. The key was then tested on the 30 aphid samples that had not been used in its construction (listed in normal font in Table 1).

Partial sequences of mitochondrial COI gene
Fifty partial COI sequences of M. cerasi and one of M. borealis from 11 countries were included in analysis. The alignment contained 616 bases in the final set with three variable sites, all of which appeared parsimony informative. The average base composition was A = 34.0%, C = 12.7%, G = 12.3% and T = 41.0%. The overall transition/ transversion ratio R = 1.221 for all sites.
Five COI haplotypes were detected (Fig. 1): one for M. borealis, two for samples from P. cerasus and two for samples from P. avium (Table 2). Aphids collected from P. mahaleb, P. maackii and P. serrulata had the same haplotype (No. 3) as the majority of samples from P. avium. COI haplotypes detected among samples from P. cerasus  whilst interspecific pairwise sample divergences between three species of Myzus ranged from 0.2 to 6.8% (Table 3). The maximum parsimony (MP) analysis of partial COI sequences resulted in 930 equally parsimonious trees (length = 43, CI=1.00, RI=1.00). The ML tree (T92 model) showed similar topology, as did NJ (K2P distances) and BI (HKY+I+G model) analyses. NJ, MP and ML bootstrap values over 50% together with BI posterior probabilities over 0.50 are given at respective nodes of the same tree in Fig. 2. Thus the M. cerasi samples form two major clades corresponding to two host-specific black cherry aphid taxa. One clade consists of all but two of the samples from P. avium, plus aphids collected from P. mahaleb, P. maackii and P. serrulata. The other clade contains all samples from P. cerasus and also includes the sample of M. borealis collected from Galium rubioides.

Morphology
When morphometric data of apterous viviparous females from 20 different geographical localities were subjected to discriminant analysis with sample collection number as the grouping variable, the first two canonical variates (Fig. 3) clearly separated sour cherry samples (COI haplotype No. 2) from those collected from sweet cherry (COI haplotype No. 3). Length of terminal process of antennal segment 6 (A6TPL), length of siphunculus (SL) and maximal length of the ventral body hairs (VBSLmax) appeared to be important predictors for separation of the two taxa (Table 4).
To discriminate between apterous viviparous females of host-specific black cherry aphid samples representing different clades in the haplotype network and the phylogenetic tree (Figs 1-2), the following linear discriminant function (LDF) was obtained: 3.924682×SL -5.6667×A6TPL -32.5504×VBSLmax + 1. Using this LDF, 97.37% individuals from the whole dataset were reclassified correctly into their a priori specified groups with host plant species as grouping variable, including 96.2% of apterous viviparous females from P. avium (n=79) and 98.6% from P. cerasus (n=73). The post hoc classification of the remaining thirty samples gave 92.37% correct specimen identification of M. cerasi cerasi (n=118) and 93.64% of M. cerasi pruniavium (n=110). The scatterplot of the mean LDF and body length values calculated for each of 30 samples representing different host specific subspecies of M. cerasi is shown in Fig. 4. The following key is therefore suggested for the identification of apterous viviparous females of the two subspecies of Myzus cerasi when sampled from winter hosts. M. persicae (1) 6.6 -6.8 M. borealis (1) M. persicae (1) 6.8   Table 1). Samples cluster in accordance with winter host plant and COI haplotype (haplotype number is given in parentheses, see Table 2 for other haplotypes).  Table 1) plotted against the mean body length for 30 samples of Myzus cerasi (normal font in Table 1) used to evaluate effectiveness of the eventual identification key. The icons are color-coded to match the COI haplotypes. Samples cluster in accordance with winter host plant and COI haplotype (haplotype number is given in parentheses, see Table 2 for haplotype information).

Discussion and conclusions
The combination of genetic distance evaluation with phylogenetic tree-building methods and multivariate analyses of morphometric data has been successfully applied to solve taxonomic problems in aphids, particularly in the genera Hyalopterus (Lozier et al. 2008), Pentalonia (Foottit et al. 2010), Aulacorthum and Neoaulacorthum (=Pseudomegoura) (Lee et al. 2011a). Based on the global data set, the average genetic divergence of COI barcode sequences between aphid species within the same genus was reported to be 5.84% (range 0.46 -11.3%), and that within species 0.05% (0.00-1.00%) , Lee et al. 2011b. Interspecific divergence of six species representing three subgenera of Myzus calculated for COI barcode sequences was reported as ranging from 5.55 to 11.3%   (Table 5). Therefore, the range of genetic divergence between the two clades of M. cerasi emerging in phylogenetic trees presented in this paper (0.0 to 0.5%) appears to be of intraspecific level. Based on the available COI data, black cherry aphids inhabiting sour and sweet cherries should therefore still be regarded as a single species. M. borealis is clearly closely related to M. cerasi and differs by only 0.2-0.5% of the COI sequences involved in the analysis. This suggests that it may also belong to the same species level taxon. More samples of M. borealis are needed to confirm this hypothesis. However, it should be noted that partial COI sequences of two biologically distinct Macrosiphum species, M. rosae (Linnaeus, 1758) and M. knautiae Holman, 1972 are very similar (Turčinavičienė and Rakauskas 2009), and low divergence levels have also been reported for Bursaphis species (Rakauskas et al. 2011) and adelgids (Žurovcová et al. 2010). It seems probable that in rapidly speciating aphid groups one may expect to find low levels of COI sequence divergence between taxa that are nevertheless functioning as distinct species.