The secondary contact zone of phylogenetic lineages of the Philaenus spumarius (Hemiptera: Aphrophoridae): an example of incomplete allopatric speciation.

Previous studies on the phylogeography of the meadow spittlebug Philaenus spumarius (L.) (Hemiptera: Aphrophoridae) suggest the existence of a contact zone of its main phylogenetic lineages along mountain chains in Europe and western Asia. This study presents a detailed examination of the population genetics of P. spumarius within the Carpathian Mountains. The main objective was to determine whether the populations inhabiting that area consist of individuals belonging to different genetic units and whether the observed pattern could be an example of secondary contact zone which formed after incomplete allopatric speciation. Specimens from six transects across the Carpathian arc were examined. The mitochondrial phylogeography of the meadow spittlebug in the examined area clearly shows that individuals from both main clades meet and mix there. Representatives of all three main EF1-α clades were also found. The present distribution of the main clades with a zone of overlap along the mountain ranges may suggest that these phylogenetic lineages form a young hybrid zone. Moreover, a limited number of individuals were shown to possess heteroplasmic mitochondrial DNA, which gives additional support to intraspecific hybridization. P. spumarius could be used in future work as an excellent model species in investigating population genetics, intraspecific hybridization, and speciation in progress.

The glacial and postglacial history of temperate species in Europe has been intensively studied leading to formation of paradigms of species range contraction to the southern refugial areas during glaciation and population expansion in warmer periods (basically during Holocene) ( Taberlet et al. 1998 ; Hewitt 1999 , 2011 ). This postglacial expansion of many species was disturbed by latitudinal distribution of Mediterranean and Black Seas and main mountain ranges (the Pyrenees, the Alps, the Carpathians, and the Caucasus) ( Provan and Bennett 2008 , Stewart et al. 2010 ). Moreover during unfavorable conditions during glaciations, species were divided into isolated populations which accumulated genetic differences as a result of genetic drift or natural selection. The level of such differences was partially correlated with the duration of population separation. This may have led to allopatric speciation and the formation of new taxa ( Gavrilets et al. 2000 ). When the climatic and habitat conditions became suitable for particular species, these formerly isolated populations could expand and meet in some areas. There are three basic scenarios of such events: 1) mixing of individuals from different populations without any genetic, ecological, or behavioral limitations (lack of any isolation); 2) hybridization and producing progeny with reduced fertility or viability (postzygotic barriers); and 3) no crossbreeding of individuals due to incompatibility in morphology, ecology or behavior of the newly arisen species (prezygotic barriers).

Recent studies on the phylogeny and population genetics of spittlebugs of the genus Philaenus show that Philaenus spumarius (L.) (Hemiptera: Aphrophoridae) is an excellent subject for such research ( Seabra et al. 2010 , Maryańska-Nadachowska et al. 2012 ). Currently, this genus is believed to consist of 9 or 10 species ( Drosopoulos 2003 , Tishechkin 2013 ). Most of them are distributed in the Mediterranean area, and only P. spumarius is widespread and occurs naturally throughout the entire temperate and warm Holarctic region ( Drosopoulos 2003 , Drosopoulos et al. 2010 ). The results of phylogeographic studies on the mitochondrial DNA (mtDNA) of P. spumarius showed that it is divided into two highly distinct clades: northeastern (NE) (north-central Europe and Asia) and south-western (SW) (western Europe and the Mediterranean), which meet along European mountain ranges ( Maryańska-Nadachowska et al. 2012 ). Similarly, analysis of a nuclear marker (elongation factor 1 alpha gene, EF1-α) suggests that there are three main clades: NE (Eurasian), south-eastern (east Mediterranean-Caucasian), and SW (Italo-Iberian), which probably overlap along European mountain ranges ( Maryańska-Nadachowska et al. 2012 ). P. spumarius inhabits different kinds of meadows (in lowlands, mountain valleys, and pastures) and its distribution is mainly limited by forests. In mountainous areas, it is distributed in meadows below and above the timberline and the main breaks in its range are fall along highly forested slopes.

It is probable that P. spumarius mtDNA and nuclear clades indicate separate evolutionary units. However, it is not certain if and how members of these lineages interact where they come into contact. The present distribution of the main clades suggest existence of a narrow zone of overlap along mountain ranges and recent study show no example of penetration of one clade into the range of other ( Maryańska-Nadachowska et al. 2012 ). This may suggest that these phylogenetic lineages form a young hybrid zone. As this species is highly mobile and divergence of its phylogenetic lineages happened recently (possibly in the end of Pleistocene, Maryańska-Nadachowska et al. 2012 ), the interesting question is if this allopatric speciation completed or not and if in contact zone exist individuals of mixed genotypes.

There are many examples of hybridization and abutment zones in European mountains. They have been observed in the Carpathians between the newts Lissotriton vulgaris and Lissotritonmontandoni ( Babik et al. 2003 , Gherghel et al. 2012 ), the butterflies Pieris napi and Pierisbryoniae ( Varga 2008 ), and the fire-bellied toads Bombina bombina and Bombina variegate ( Szymura and Barton 1986 , Hofman 2007 , Fijarczyk et al. 2011 ). In the Alps, examples include the Valais shrew Sorex antinorii and the common shrew Sorex araneus ( Brünner et al. 2002 ), and also P. napi and P. bryoniae ( Porter 1997 ). In the Pyrenees, there is evidence for a hybrid zone of the meadow grasshopper subspecies Chorthippus parallelus parallelus and Chorthippus parallelus paralleluserythropus ( Butlin and Hewitt 1985 , Buño et al. 1994 ). These examples, as well as many others, suggest that mountain ranges in Europe are suture zones where distinct phylogenetic lineages of many taxa meet.

The main objective of this study was to investigate whether individuals belonging to different phylogenetic lineages are present in the Carpathian Mountains and if this mountain arc constituted a barrier to gene flow among P. spumarius population belonging to these clades.

Moreover, we examined the degree of diversification of P. spumarius individuals across the entire Carpathian range and the adjacent areas to verify if in this area exist populations consisting on genetically distinct individuals or are present individuals showing mixed genotypes, what may help in understanding if allopatric speciation in P. spumarius is complete or not. Specifically, we tested following hypotheses: 1) in the Carpathians exist contact zone between distinct phylogenetic lineages of P. spumarius ; 2) the Carpathians constitute barrier for dispersal of spittlebugs belonging to distinct phylogenetic and geographic lineages; 3) individuals of P. spumarius in the Carpathians show genotypes specified for particular phylogenetic lineages (= in the Carpathians are absent individuals of mixed genotypes), and consequently; and 4) P. spumarius finished its allopartic speciation and therefore should be divided into distinct taxons after examination of other characters (morphological, biological, and ecological).

Materials and Methods

Study Area

Six transects across the Carpathian arc were designed to cover all parts of these mountains ( Fig. 1 ). Each transect consists of three sampling plots: one in the outer part of the Carpathian arc, a second in the higher part of the mountains (in the center of the mountain belt), and a third in the inner part of the Carpathian arc. Beginning from the north, these transects run across: 1) the western Carpathians (WC); 2) the border between the western and eastern Carpathians (WEC); 3) the eastern Carpathians (EC); 4) the border between the eastern and southern Carpathians (ESC); and 5) the southern Carpathians (SC). One of the transects (6), although designated in a similar manner, is located entirely within the Carpathian Basin (in the inner part of the arc) in the Apuseni Mountains (AM) (Romania). Detailed information about sampling plots is presented in Table 1 .

Fig. 1.

Localization of P. spumarius sampling sites in the Carpathians. The color of circles mark transect of examined populations (west Carpathians, dark blue; west/east Carpathians, light blue; east Carpathians, green; east/south Carpathians, yellow; south Carpathians, orange; AM, red). Map of the Carpathian after and OpenStreetMap ( www.openstreetmap.org ).

Table 1.

Analyzed samples of P. spumarius with description of the sampling sites, symbols and numbers (applied to CytB) of studied specimens (updating with heteroplasmatic individuals)

Locality symbol	Country	Locality	Mountain range	Number of specimens	Transect
PS5	Slovakia	Vélke pole	Vtáčnik	5	WC
PS7	Slovakia	Malatina	Chočské vrchy	5	WC
PS17	Poland	Kamienica River valley	Gorce Mst.	5	WC
PS32	Poland	Jaśliska	Beskid Niski	5	WEC
PS35	Slovakia	Olka	Nizke Beskydy	5	WEC
PS38	Hungary	Regec	Zemplenyhegység	5	WEC
PS42	Romania	NE from Baia Mare	Gutâi Mst	6	EC
PS45	Romania	Tihuta Pass	Bârgău Mst	7	EC
PS47	Romania	Petru Voda Pass,	Stânişoarei Mst	6	EC
PS53	Romania	Predeal	Bucegi Mst	8	ESC
PS54	Romania	Sinaia	Bucegi Mst	5	ESC
PS55	Romania	Sinca Veche	Făgăraşului Mst	5	ESC
PS59	Romania	Voineasa	Lotrului Mst	12	SC
PS61	Romania	Crasna	Parângului Mst	5	SC
PS65	Romania	Băiţa	Metaliferi Mst	5	SC
PS67	Romania	Buceş	Bihorului Mst	5	AM
PS68	Romania	Marişel	Gilăului Mst	9	AM
PS69	Romania	Şuncuiuş	Pădurea Craiului Mts.	5	AM

Sampling

Individuals from 18 populations of P. spumarius were collected in 2011 ( Table 1 ). The spittlebugs were caught in a sweep net, immediately preserved in 99% ethanol, and finally stored at −20°C. External samples were single specimens from populations located outside the Carpathians, which were used previously in a phylogeographic study of this species ( Maryańska-Nadachowska et al. 2010 , 2012 ). In total, 108 specimens of P. spumarius from the Carpathian basin were analyzed. In addition, 23 specimens of P. spumarius and Philaenustesselatus from other localities across Europe were used for cytochrome B (CytB) analyses and 15 specimens were used for EF1-α analyses (samples from Maryańska-Nadachowska et al. 2012 ; Table 1 ). Moreover, single specimens of Philaenus italosignus and Philaenussignatus were used for the construction of an EF1-α tree (data from Maryańska-Nadachowska et al. 2012 ). All voucher specimens are preserved at the Institute of Systematics and Evolution of Animals, Polish Academy of Sciences.

Laboratory Techniques

Total genomic DNA was extracted with a NucleoSpin Tissue kit (Macherey-Nagel, Düren, Germany), using five individuals from each population. Two DNA fragments were amplified. We managed to amplify and sequence CytB from 108 individuals; however, 18 of them showed to include possible heteroplazmic mtDNA (see following chapters) and therefore, only 90 specimens were used in population genetic analyses. On the other hand, due to presence of multiple indels in EF1-α gene of most individuals we were able to obtain sequences only from single representatives from each population (18 in total). DNA isolation, amplification, purification, and sequencing were performed in the same way as in Maryańska-Nadachowska et al. (2012) . In some populations, additional specimens were used to fulfill the criterion of five specimens per population, as some individuals gave CytB chromatograms with double peaks, and these individuals were not considered in analyses. New haplotypes were deposited in GenBank ( http://www.ncbi.nlm.nih.gov/genbank/ ) under KC793208–KC793240 for CytB and KC793241–KC793258 for EF1-α.

Sequence Analyses

DNA sequences were edited using the BioEdit Sequence Alignment Editor 5.0.9. ( Hall 1999 ) and aligned using ClustalX 1.8. ( Thompson et al. 1997 ). Appropriate nucleotide substitution models were determined using MrModeltest v2. ( Nylander 2004 ) in conjunction with PAUP* ( Swofford 2002 ).

For EF1-α, Bayesian phylogenetic analyses were conducted, with two independent runs of four Metropolis-coupled Monte Carlo Markov chains (three of them “heated”) were conducted for 3 × 106 generations and sampled every 100 generations ( Huelsenbeck and Ronquist 2001 , Huelsenbeck et al. 2001 ). Convergence of Bayesian analyses was assessed using Tracer v. 1.5.0 ( Rambaut and Drummond 2009 ); ∼7,500–8,000 of the sampled trees were discarded as “burn-in”, while the remaining ones were used to reconstruct the majority rule consensus tree. The phylogenetic tree was visualized with TreeView 1.6.6 ( Page 1996 ).

CytB haplotypes were identified and standard genetic indices, such as haplotype diversity (h), nucleotide diversity (π), and number of segregating sites (S), were computed for each species and each population using the DnaSP v5 software ( Librado and Rozas 2009 ). Pairwise genetic distances for both markers were calculated using MEGA v5 ( Tamura et al. 2011 ).

Relationships between the Carpathian P. spumarius populations and samples from the other parts of the species range were described by building haplotype networks for CytB sequences using the statistical parsimony method ( Templeton et al. 1992 ) and TCS 1.21 software ( Clement et al. 2000 ). In addition, in order to determine whether the Carpathian populations are genetically distinct, we applied analysis of molecular variance (AMOVA) as implemented in ARLEQUIN 3.1. ( Excoffier et al. 2005 ), grouping populations in several ways. First, differentiation among transects was tested. Second, the distribution of genetic variation between inner and outer populations was examined (central populations were excluded, populations from the AM were assigned as inner). Third, populations were grouped according to their localization within the Carpathian arc (western: WC+WEC, eastern: EC+ESC, and southern-central: SC+AM). A Mantel test ( Mantel 1967 ) was performed using ARLEQUIN 3.1 (with pairwise Fixation index F ST values and straight-line geographic distances in kilometers) to check if the genetic structure of the populations fits an isolation by-distance (IBD) model ( Slatkin 1993 ).

Heteroplasmy Testing

Double peaks in CytB chromatograms were found in 18 P. spumarius individuals and a series of analyses were conducted to verify that this represented heteroplasmy and to exclude the possibility of: 1) DNA contamination and 2) the presence of nuclear pseudogenes. 1) Additional PCRs and sequencing were performed using templates from the individuals as well from other individuals and additional blank samples (without DNA). 2) Sequencing with both forward and reverse primers was done for the individuals. 3) The presence of stop codons was verified using MEGA v5 ( Tamura et al. 2011 ). 4) Specific primers were designed to amplify only CytB from both main mitochondrial clades (in previous work termed “south-west” and “north-east”). The CytB_PS_NE_F 5′ GGG CGA GGA ATA TAT TAT GGA TC 3′ and CytB_PS_NE_R 5′ GAT TTG CTG GAA TGA AAT TAT C 3′ primers were used for the NE clade and the CytB_PS_SW_F 5′ GGA CGA GGA ATA TAT TAT GGG TC 3′ and CytB_PS_SW_R 5′ GAT TTG CTG GGA TAA AAT TGT C 3′ primers for the SW clade. However, it was not possible to find CytB fragments with more than two or three polymorphic base pairs within a stretch of 20–25 bp, so these specific primers could probably also amplify the other CytB variant. To prevent this, a higher annealing temperature was used. 5) PCRs and sequencing of another mitochondrial gene—cytochrome oxidase I (COI) on templates from individuals were performed using C1-J-2195 and TL2-N-3014 primers ( Simon et al. 1994 ), which were also applied for this species by Seabra et al. (2010) and proved to be of clearly mitochondrial origin. 6) Alignments of CytB in individuals from both main mitochondrial clades were visually analyzed and nucleotide sites with double peaks were compared with the localization of polymorphic sites differing between the mitochondrial clades.

Results

Nuclear DNA Pattern

The results should be interpreted with caution, because a biased sample of individuals (those lacking the indel) are analyzed. Similarly to previous phylogeographic work on P. spumarius ( Maryańska-Nadachowska et al. 2012 ), clear EF1-α chromatograms could be obtained from only ∼20% of specimens due to the presence of some indels within this gene (a variable number of repeats in short microsatellites located at several sites in the presumed intron), which prevent sequencing. Some parts of sequences showed numerous double peaks from both directions (forward and reverse). Consequently, only a single representative of each studied population from the Carpathians was genotyped. There were 41 variable sites in 14 sequences (each of ∼1,005 bp length) and 12 of them were parsimony informative (Hd = 0.954, SD = 0.039;  π = 0.01042, SD = 0.00258). In total, 14 genotypes were identified among 18 individuals. The Carpathian genotypes were compared with sequences obtained from other (non-Carpathian) populations ( Fig. 2 ). The phylogenetic EF1-α tree showed that most genotyped individuals belong to the NE EF1-α clade (E1) (with the closest genotypes found in Ukraine, Crimea, Russia, Norway, and the Alps). Exceptions were found only in five populations. The largest genetic distances between obtained EF1-α sequences reached ∼4%.

Fig. 2.

Bayesian EF1-α phylogenetic tree of P. spumarius from the Carpathians (population codes as in Table 1 with addition of I, inner population; C, central; and O, outer, the color of squares corresponds to different transects according to Figure 1 , and the color of vertical lines corresponds to distribution of EF1-α clades in Figure 4 ) and selected haplotypes found in the rest of the species range (see Fig. 1 from Maryańska-Nadachowska et al. 2012 ) with P. signatus, P. italosignus , and P. tesselatus as outgroups. Numbers above lines indicate posterior probabilities of Bayesian inference (showed only when above 0.50). Symbols: I, inner; C, center; O, outer of the Carpathians arch; N, North; S, South.

Fig. 4.

Distribution of CytB clades (left map: red circles, south-west clade; blue circles, north-east clade; violet stars, localities with detected probable heteroplasmic individuals) and EF1-α clades (right map: blue triangle, E1 clade; red triangle, E2 clade; yellow triangle, E3 clade).

The mtDNA Pattern

Indels and stop codons were not observed in the 664 bp of CytB. In total, among 90 sequences (excluding sequences containing double peaks in chromatograms) from the Carpathians we found 33 haplotypes, 48 variable sites, and 12 parsimony informative sites (Hd = 0.940, SD = 0.012; π = 0.01947, SD = 0.00401). Similarly like for EF1-α largest genetic distances between CytB clades reached ∼4%. There are 17 sites discriminating members from two main clades (17 sites are fixed between these two phylogenetic lineages). Standard genetic indices for each population as well as for each transect are presented in Table 2 . The most variable populations were found to be PS17 and PS38, and the least variable ones PS5 (monomorphic), PS45, PS54, PS59, PS65, 263 PS68, and PS67 ( Table 2 ). Phylogenetic CytB networks ( Fig. 3 ) showed that a similar number of individuals belong to the NE clade (42 individuals) and SW clade (48 individuals) ( Fig. 4 ). Haplotypes belonging to the SW CytB clade were found in five out of six transects (except EC), and those belonging to the NE CytB clade were found in all transects. Haplotypes from the SW clade were found in 11 populations. The SW clade was distributed mainly in the inner part of the Carpathian arc, in some central populations (in the WC, WEC, ESC, and SC transects), and in some populations in the outer part of these mountains (mainly in SC, but also in WC and WEC). The pattern of the NE clade distribution was the opposite of the SW clade.

Table 2.

Standard genetic indices of CytB sequences and numbers of individuals belonging to particular CytB and EF1-a clades

Transect/locality	CytB polymorphism								CytB clade		EF1-α clade
	N	V	S	H	h ± SD	π ± SD	Np.	Hp	NE	SW	E1	E2	E3
All Carpathians	90	50	48	33	0.940 (0.012)	0.0195 (0.004)		18	42	48	13	4	1
transect WC	15	28	28	6	0.790 (0.079)	0.0181 (0.005)		0	5	10	2	1	0
Transect WEC	15	37	36	11	0.952 (0.040)	0.0155 (0.006)		0	2	13	2	1	0
Transect EC	15	9	9	7	0.781 (0.102)	0.0042 (0.001)		4	15	0	2	1	0
Transect SEC	15	25	25	9	0.933 (0.040)	0.0161 (0.004)		3	10	5	2	1	0
Transect SC	15	25	24	7	0.838 (0.068)	0.0159 (0.004)		7	5	10	2	1	0
Transect AM	15	24	24	7	0.838 (0.068)	0.0147 (0.004)		4	5	10	2	0	1
PS5	5	0	0	1	0.000 (0.000)	0.0000 (0.000)	0	0	0	5	1	0	0
PS7	5	25	25	3	0.700 (0.218)	0.0154 (0.009)	0	0	1	4	1	0	0
PS17	5	28	28	5	1.000 (0.126)	0.0178 (0.010)	3	0	4	1	1	0	0
PS38	5	10	10	5	1.000 (0.126)	0.0078 (0.004)	3	0	0	5	0	1	0
PS35	5	5	5	4	0.900 (0.161)	0.0042 (0.002)	2	0	0	5	1	0	0
PS32	5	26	26	4	0.900 (0.161)	0.0232 (0.009)	1	0	2	3	1	0	0
PS42	5	7	7	4	0.900 (0.161)	0.0051 (0.003)	3	1	5	0	1	0	0
PS45	5	1	1	2	0.400 (0.237)	0.0006 (0.001)	1	2	5	0	0	1	0
PS47	5	5	5	3	0.700 (0.218)	0.0033 (0.002)	1	1	5	0	1	0	0
PS55	5	3	3	4	0.900 (0.161)	0.0018 (0.001)	3	0	0	5	0	1	0
PS53	5	5	5	3	0.800 (0.164)	0.0042 (0.002)	1	3	5	0	1	0	0
PS54	5	2	2	2	0.600 (0.175)	0.0018 (0.001)	0	0	5	0	1	0	0
PS65	5	6	6	2	0.600 (0.175)	0.0054 (0.002)	0	0	0	5	1	0	0
PS59	5	1	1	2	0.400 (0.237)	0.0006 (0.001)	1	7	5	0	1	0	0
PS61	5	7	7	4	0.900 (0.161)	0.0045 (0.003)	2	0	0	5	0	1	0
PS67	5	6	6	2	0.400 (0.237)	0.0036 (0.002)	0	0	0	5	1	0	0
PS68	5	1	1	2	0.400 (0.237)	0.0006 (0.000)	0	4	5	0	0	0	1
PS69	5	3	3	4	0.900 (0.161)	0.0021 (0.001)	2	0	0	5	1	0	0

N , number of individuals (without heteroplasmatic individuals); V, variable sites; S, segregating sites; H, haplotype number; h, haplotype diversity; π, nucleotide diversity; SD, standard deviation; Np, number of private haplotypes; and the distribution of CytB and EF1-a clades.

Fig. 3.

Haplotype CytB networks of P. spumarius from the Carpathians (population codes as in Table 1 ), and selected haplotypes found in rest of the species range (see Fig. 1 from Maryańska-Nadachowska et al. 2012 ). Red circles, south-west clade; blue circles, north-east clade. Numbers in square brackets show numbers of individuals harboring particular haplotypes within populations. *, haplotypes in each clade which are most closely related (they are differ by 17 substitutions).

Some Carpathian CytB haplotypes belonging to the NE clade were the same as those from Georgia and very close to those from Ukraine (mainly from the ESC, EC, WEC, and WC transects) or the same as those from Norway (SC and AM). Other Carpathian NE haplotypes were distinct from any closest by three or more mutation points. Any haplotypes from the SW clade were common with any non-Carpathian populations but some of them (from WEC, SC, and AM transects) were closest with some haplotypes from Hungary. Moreover, most individuals belonging to the SW clade from WC and WEC formed a clearly distinct cluster (0.95 posterior probability), which was not directly related to any other SW haplotypes from Eurasia.

The results of AMOVA for three methods of population grouping are presented in Table 3 . In the first method (according to transects), the greatest amount of genetic variation was found between populations within groups (∼58%), in the second method (inner and outer localization) also the greatest amount of genetic variation was revealed between populations within groups (58%), the same was in the third method (western-eastern-southern gradient) where the greatest amount of genetic variation was revealed between populations within groups (∼48%).

Table 3.

Results of AMOVA for CytB sequences grouping in few ways

AMOVA
Source of variation	A (%)	B (%)	C (%)
Among groups	14.8	13.7	26.3
Among populations	57.5	58.0	47.7
Within groups
Within populations	27.7	28.3	26.0
Total	100.0	100.0	100.0
F SC ( P )	0.674 (<0.001)	0.673 (<0.001)	0.647(<0.001)
F ST ( P )	0.723 <0.001)	0.717 (<0.001)	0.740 (<0.001)
F CT ( P )	0.148 (0.145)	0.139 (0.110)	0.263 (0.022)

A, according to six transects across the Carpathian arch; B, into two groups: localities inner and outer of the Carpathian arch; C, into three groups according to their localization in the Carpathian arch.

The results of the Mantel test showed that there was a correlation between genetic (F ST ) and geographical (km) distances for the Carpathian populations of P. spumarius . However, this relationship was weak ( r  = 0 . 143) and on the verge of significance ( P  = 0.05).

Chromatograms of CytB with double peaks in several nucleotide positions were found in 18 individuals ( Fig. 5 ). These individuals were from the following populations: PS42 (1 individual), PS45 (2), and PS47 (1) (transect EC); PS53 (2) (transect ESC, central); PS59 (7) (transect SC, central); and PS68 (2) (transect AM, central), so exclusively from populations where other specimens showed to harbor only NE haplotypes. These heteroplasmic sites correspond to polymorphic positions between the two main mitochondrial clades ( Table 4 ). The heteroplasmic individuals always contained nucleotides typical of both mitochondrial clades.

Fig. 5.

Selected polymorphic parts of CytB chromatograms from haplotypes belonging to north-east (upper) and south-west (central) CytB clades and sequences with double nucleotide peaks (lower) which are probable examples of heteroplasmatic individuals; asterisks indicate polymorphic sites which produce double peaks in heteroplasmic sequences.

Table 4.

Selected polymorphic parts of CytB sequences from north-east and south-west CytB clades and few sequences which double nucleotide peaks which are probable examples of heteroplasmatic mitochondrial sequences

Specimen	Polymorphic parts of sequences	CytB clade
PS38C	GACGGTCAATGATCTACTACCGAGGGTTACTCATTACTATCCGGTCGTTCACCAATAGCA	NE
PS32D	.........................A .. G..............A...........C....	NE
PS38D	...........................................A...........C....	NE
PS17D	.G ..A.. G.....T ..T..T..G.. T.....A ..C..C..T..A..A..T.. T.....T.	SW
PS53D	.G.....G.....T .. T.....G .. T.....A ..C..C..T..A..A..T.. T.....T.	SW
PS47A	.G ..A..G..A..T..T..T..G.. T.....A ..C..C..T..A..A..T.. T.....T.	SW
PS42A	.R ..R.. R.....Y ..Y..T..R..K..R..M..Y..C..Y..A..A..Y.. T.....Y.
PS68D	.R.....G.....Y ..Y..T..R..K..R..M..Y..C..Y..A..A..Y.. T.....Y.
PS59F	.R.....G.....Y ..Y..T..R..K..R..M..Y..C..Y..A..A..Y.. T.....Y.
PS47D	.G ..A..G..R..T..T..T..G..K..R..A..C..C..T..A..A..T.. T.....T.

Such results were consistently repeated only in these several individuals, and not in isolates from other specimens. Blank samples did not show any PCR products, which excludes contamination of reagents or laboratory equipment. The same results were obtained by forward and reverse sequencing. No stop codons were detected in any of these 18 sequences. Some PCRs performed on DNA from individuals with double peaks using semispecific primers for the SW or NE CytB clades gave one type of sequences, but not explicitly from the NE or SW clade. However, in some samples these primers again produced sequences with double peaks, most probably because the primers were different in only two or three nucleotide positions and could amplify DNA of both mtDNA variants, despite the higher annealing temperature used for these PCRs. Amplification and sequencing of the COI mtDNA gene showed the same pattern as in CytB, but not for isolates from other specimens. And finally, alignments of “pure” haplotypes from the SW and NE clades, as well as sequences from presumably heteroplasmic individuals, clearly showed that double nucleotide peaks appeared only in those nucleotide positions which are polymorphic between these main mitochondrial clades ( Table 4 ).

Discussion

Some substantially differentiated phylogenetic lineages were found across the entire natural range of the meadow spittlebug P.spumarius ( Seabra et al. 2010 , Maryańska- Nadachowska et al. 2012 ). According to mtDNA, this species is divided into two groups: the NE clade ranging from eastern Asia to central and northern Europe and the SW clade, which is distributed in western and southern Europe and in the Middle East. The latter clade is additionally subdivided into several subclades distributed in SW and SE Europe, around the Black Sea, and in the Caucasus. A similar pattern was found for nuclear DNA—besides the NE clade (E1), there are two lineages: one in southeastern Europe and the Middle East (E2), and the other in SW and western Europe (E3). The ranges of these clades (mitochondrial and nuclear) are generally concordant which the exception of the Black Sea area, where individuals belong to the SW mtDNA clade (the Black Sea subclade) and to the NE nuclear clade. A similar population structure—SW/NE division has also been documented for other species, especially those with large Eurasian ranges extending across grassland habitats, e.g., the European grasshopper Chorthippus parallelus ( Cooper et al. 1995 ) and the buprestid Coraebus elatus ( Kajtoch et al. 2013a ). This indicates that grassland species (including steppic and those dependent on high mountain meadows) had two main refugial areas: one eastern European-Asian and the other Iberian ( Ribera and Blasco-Zumeta 1998 , Stewart et al. 2010 ). Moreover, the main P. spumarius clades (both mitochondrial and nuclear) seem to be distributed allopatrically with contact zones detected along the main mountain ranges of Europe and western Asia (the Pyrenees, the Alps, the Carpathians, and the Caucasus). It seems that populations from NE and the southern clades meet along these mountains but do not cross them. No populations apart from those inhabiting these mountains consist of individuals from two or more clades.

According to mtDNA (CytB), members of the two main clades meet in the Carpathians. The NE clade was found in all parts of the Carpathians with exception of some most western and southern localities. On the other hand, SW clade was present mostly in inner part of WC and in SC, lacking from EC. The results of AMOVA showed that regardless of grouping method the highest percentage of genetic variation is consistently found within groups but not among these groups. This pattern suggests that the Carpathians were colonized from several sources (directions) and that within Carpathians there is no clear boundary between main mitochondrial lineages. Mitochondrial data support the NE source (an eastern European-Asian refugium) and also another source for SW haplotypes but it seems that these haplotypes were not derived from any of Mediterranean refugium. In the CytB network, there is a large group of individuals which form a separate subclade within the SW clade (distinct by four mutation steps from populations from southern Russia and five steps from Spanish populations). Members of this clade were found exclusively in the WC and on the border of the WEC (Poland and Slovakia), in the northern part of these mountains. This suggests that also in the WC there existed a cryptic refugium of this species. Such cryptic northern refugia within the Carpathian basin have been found for several species, mainly for water-dependent taxa such as newts ( Babik et al. 2003 , Gherghel et al. 2012 ) and fire-bellied toads ( Szymura and Barton 1986 , Hofman 2007 ), and also for grassland species such as grasshoppers ( Butlin and Hewitt 1985 , Hewitt 1993 , Buño et al. 1994 ), and meadow katydids ( Shapiro 1998 ). The refugium for a distinct mtDNA clade was not necessarily within the Carpathians—this area could have been colonized from the north, where, in the Polish highlands, there probably existed other refugia, mostly for grasslands species (e.g., Pawłowski 1999 , 2005 ; Kajtoch 2011 ; Kubisz et al. 2012 ; Kajtoch et al. 2013a , b ). A similar general pattern was observed for the nuclear EF1-α marker; however, those data should be taken with caution as only ∼20% of the studied individuals could be genotyped. The sampling was biased by the selection of only those individuals that lacked the indels. Consequently, the observed pattern could change if more samples are added or if other nuclear markers are used. EF1-α supports expansion from the north-east direction and also from south-east, probably from the Balkans. Moreover, EF1-α suggests an admixture of individuals from the SW clade (in this case from northern Italy).

Multiple sources of origin of P. spumarius populations within the Carpathians, as well as a probable cryptic refugium within this mountain range or north of the Carpathians, reveal a highly complicated history and many migration routes. The highest genetic distances among members of the main clades, both mtDNA and nuclear, are ∼4%. Assuming standard molecular clock rates for insect mtDNA protein genes to be 1.4–2.3% per million years ( Farrell 2001 , Ribera et al , Barraclough and Vogler 2002 , Borer et al. 2010 ).This distance suggests a split between the main mtDNA clades between ∼1.7 and 2.8 million years ago. This period falls within the Early Pleistocene glaciations. P. spumarius must have been divided during the Pleistocene into several isolated populations, most probably successively in consecutive glaciations—in eastern Europe/Asia, the Balkans/Asia Minor, Crimea, and the western Mediterranean. It is very likely that P. spumarius is in the process of speciation into two or three taxa. Populations from NE nuclear clade (≈E1) and SW mitochondrial clades (and probably also from nuclear clades E2 and E3, which correspond to the main mitochondrial SW clades) should be investigated in respect to their morphological, biological, ecological, or behavioral characters to verify if they deserve to be described as separate subspecies or different taxonomic units.

Assuming that P. spumarius consists of two or more distinct evolutionary units, it would be interesting to determine if members of these lineages interbreed without any pre or postzygotic difficulties. P. spumarius can be an example of the intraspecific hybridization, which is known in, e.g., grass snakes ( Thorpe 1983 ), meadow grasshoppers ( Hewitt 1993 , Buño et al. 1994 ), and stick insects ( Demontis et al. 2010 ). As such tests were not part of the present work, only some indirect conclusions can be drawn from the distribution of genetic variants among and within populations. The current distribution of the main clades may be explained in a two alternative ways: 1) P. spumarius may have just recently met along mountains in Europe (e.g., the Carpathians) and western Asia, and is currently forming a “mixing” zone as a result of dispersion or 2) this contact zone was formed at the beginning of the Holocene and has been stable ever since. The first scenario is less probable as given that P. spumarius is highly mobile ( Thompson 1984 , Halkka and Halkka 1990 , Stewart and Lees 1996 , Drosopoulos et al. 2010 ) and its expansion probably has started in the end of glaciations, there was enough time to spread all genetic variants across the entire range of the species, what was not observed ( Maryańska-Nadachowska et al. 2012 ). Second scenario is more likely, however it assume that some factors must shape the stability of this zone and prevent populations belonging to different phylogenetic lineages from mixing. Such factor(s) could include selection against hybrids or some ecological or behavioral barriers which prevent reproduction (pre or postzygotic isolation). The presented data do not give a clear answer if spittlebugs from different phylogroups breed in nature and if hybrids are viable and fertile. Only one study reports crossing P. spumarius individuals from different populations in a laboratory ( Yurtsever 2002 ). However, the study concerns populations from New Zealand and Wales, probably both of which belong to the same SW mitochondrial clade. The existence of populations within the Carpathians that consist of individuals belonging to two mitochondrial clades (mainly in the WC) and some discrepancies in the distribution of CytB and EF1-α clades in several populations (mainly in the inner part of the Carpathians) suggest that individuals belonging to different phylogenetic lineages can cross, and “hybrids” with mixed genotypes can be found. The Carpathians are a suture zone, where many species have their hybrid zones, e.g., newts, ( Babik et al. 2003 , Gherghel et al. 2012 ), fire-bellied toads ( Szymura and Barton 1986 , Hofman 2007 ), barbels ( Brun et al. 1992 , Chenuil et al. 2004 ), and butterflies ( Varga 2008 ). P. spumarius can be added to this list of hybridizing taxa in the Carpathians and other mountains. However, this phenomenon needs to be verified, both in nature and in the laboratory.

The last clue as to the existence of P. spumarius hybrids involves examples of mitochondrial sequences of presumably heteroplasmic origin. Contamination of the probe was excluded as no blank samples led to amplified products and particular products were obtained consistently from the same isolates. The nuclear origin of mtDNA sequences was ruled out as no stop codons or indels were detected in any CytB sequences, and primers semispecific to the NE or SW clade mostly produced only single sequences (without double peaks). Moreover, the same pattern was observed in a different mtDNA gene, i.e., COI, located on the opposite side of the molecule, which strongly suggests that heterozygous mitochondrial sequences cannot be of nuclear origin. If mitochondrial pseudogenes were true, this would mean that entire mtDNA or both of these genes independently were incorporated into the nuclear genome, which is a rather implausible situation. Moreover, if part or entire mtDNA has its nuclear homologue, it should be present in all individuals, and not only in a small fraction of them. In the Carpathians, only ∼17% of the sequenced individuals possessed heterozygous CytB sequences, and no such individuals had been found elsewhere in the range of this species ( Maryańska-Nadachowska et al. 2012 ) or in COI sequences ( Seabra et al. 2010 ); however, these studies had not been so detailed in sampling. Most importantly, almost all heterozygous nucleotide sites were the same as polymorphic nucleotide sites which differentiated the NE and SW clades. This last observation strongly supports the hypothesis that heterozygous mtDNA sequences are of heteroplasmic origin. If one CytB copy was of nuclear origin, mutations should appear also in other nucleotide positions not associated with DNA polymorphism between the clades. It is most probable that some individuals possess double mtDNA haplotypes, one from the maternal ancestor, and the other as a result of paternal leakage. Such a high number of double peaks in chromatograms suggests a similar frequency (density) of both CytB variants in the isolates. This indicates that the studied individuals are not F1 hybrids (as they should have much less paternal mtDNA), but rather later-generation hybrids (in which there is a similar number of both CytB variants). This assumption is also partially supported by the fact the all heteroplasmic individuals were found in populations where other members belong exclusively to NE clade. This observation suggests that recent hybridization is unlikely. Therefore, these individuals may harbor two mitochondrial types as a result of past introgression. On the other hand, it is possible that our sampling (five specimens per population) was not sufficient and we missed some heteroplasmic individuals in other Carpathian and Eurasian populations.

Heteroplasmy is very difficult to detect and most probably can be found among hybrids which possess mtDNA variants from both parental species. As P. spumarius consists of two highly diversified mtDNA haplogroups, it presents an opportunity for observing individuals harboring both mitochondrial variants. Heteroplasmy with paternal leakage has been previously detected in other insects, among others in hemipterans: the cicadas Magicicada septendecim, Magicicada septendecula , and Magicicada cassini ( Fontaine et al. 2007 ), and in Nezara viridula ( Kavar et al. 2006 ). Heteroplasmy is not so rare a phenomenon as it was assumed previously (e.g., White et al. 2008 ); however, its detection and interpretation is not easy and may be affected by pseudomitochondrion genome detection ( Parr et al. 2006 ).

The contact zone of distinct phylogenetic lineages, P. spumarius could be used in future work as an excellent model species for the investigation of population genetics, intraspecific hybridization, and possibly speciation in progress. In addition, the detected heteroplasmic individuals are examples of a rare and interesting phenomenon. Still, more variable genetic markers (e.g., microsatellites, single nucleotide polymorphisms, genome and transcriptome sequencing) are needed to enable a comprehensive study of this species.