Fifty population samples of apterous viviparae of black cherry aphids from nine European countries and Turkey were collected in 2004–2013, mostly from Prunus cerasus, Prunus 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.

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 PCR-amplified 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 (2 mM each), 4–8 µl of 25 mM MgCl₂ and 1.25 U of AmpliTaq Gold 360 polymerase (5 U/µl) and ddH₂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 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.

In addition to the sequences from 50 samples of Myzus cerasi, COI sequences of Myzus borealis from subgenus Myzus sensu stricto (the same subgenus as Myzus cerasi) and Myzus 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 Myzus cerasi and Myzus borealis using TCS v1.21 (Clement et al. 2000).

Samples representing different clades in the molecular tree and haplotype network were used for stepwise discriminant analysis followed by canonical analysis: 10 samples from sour cherry, Prunus cerasus, and 10 samples from sweet cherry, Prunus avium (shown in bold in Table 1).

Based on earlier taxonomic work (Heinze 1961, Dahl 1968), 19 metric (in mm) characters were studied: A2L – length of antennal segment 2; A3L – length of antennal segment 3; A4L – length of antennal segment 4; A5L – length of antennal segment 5; A6BL – length of basal part of antennal segment 6; A6TPL – length of terminal process of antennal segment 6; BL – body length (excluding cauda); Bwant3 – basal width of antennal segment 3; CL – length of cauda; CW – basal width of cauda; DT3L – length of the second segment of hind tarsus; F3L – length of hind femur; FF – depth of the frontal furrow; SL – length of siphunculus; T3L – length of hind tibia; URL – length of ultimate rostral segment; URW – basal width of ultimate rostral segment; VBSLmax – maximal length of the ventral body hairs; VBSLmin – minimal length of the ventral body hairs.

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).

Fifty partial COI sequences of Myzus cerasi and one of Myzus 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 Myzus borealis, two for samples from Prunus cerasus and two for samples from Prunus avium (Table 2). Aphids collected from Prunus mahaleb, Prunus maackii and Prunus serrulata had the same haplotype (No. 3) as the majority of samples from Prunus avium. COI haplotypes detected among samples from Prunus cerasus (No. 1 and 2) and the remaining Prunus species (No. 3 and 4) differed in the following nucleotide positions of 616 bp alignment: 300 (between No. 1–2 and No. 3–4), 321 (between No. 3 and No. 4) and 390 (between No. 1 and No. 2). The range of the intraspecific pairwise sample divergences (K2P model) was 0.0 – 0.5% (average 0.2%), 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 Myzus 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 Prunus avium, plus aphids collected from Prunus mahaleb, Prunus maackii and Prunus serrulata. The other clade contains all samples from Prunus cerasus and also includes the sample of Myzus borealis collected from Galium rubioides.

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 Prunus avium (n = 79) and 98.6% from Prunus cerasus (n = 73). The post hoc classification of the remaining thirty samples gave 92.37% correct specimen identification of Myzus cerasi cerasi (n = 118) and 93.64% of Myzus 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 Myzus 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.