Disentangling phylogeography, polyploid evolution and taxonomy of a woodland herb (Veronica chamaedrys group, Plantaginaceae s.l.) in southeastern Europe.

Southeastern Europe is a centre of European biodiversity, but very little is known about factors causing the observed richness. Here, we contribute to fill this gap by reconstructing the spatio-temporal diversification of the cytologically variable and taxonomically intricate complex of Veronica chamaedrys (Plantaginaceae s.l.), growing in open forests, forest edges and grasslands, with flow cytometry, molecular markers (AFLPs, plastid DNA sequences) and morphometry. Our results show that both diploid and tetraploid cytotypes are widespread, but diploids predominate on the southern Balkan Peninsula. Plastid sequences suggest a first split into three main lineages in the mid-Pleistocene and a continuous diversification during the last 0.4 my. Two of the identified plastid lineages coincide with geographically distinct AFLP clusters. Altogether, the genetic data suggest forest refugia on the southern-most Balkan Peninsula (Greece), in Bulgaria, Istria (Croatia and Slovenia) and maybe the southeastern Carpathians (Romania). Morphometric and genetic data show little congruence with current taxonomy.

Introduction

It is long and widely acknowledged that southeastern Europe, and especially the Balkan Peninsula, is a centre of European biodiversity (Turrill, 1929). One early recognised factor is the refugial character of the Balkan Peninsula with a high proportion of relic taxa, even if in many cases the claimed relic status still needs to be confirmed. Much less is, however, known about diversification processes and their spatio-temporal patterns on the Balkan Peninsula especially at lower taxonomic levels, i.e., within and among closely related species. This is partly due to the fact that in molecular studies the Balkan Peninsula is often neglected, such as in large-scale phylogeographical studies, where only a few (e.g., Taberlet, 1994; Santucci et al., 1998; Trewick et al., 2002) or no samples were included (Dumolin-Lapègue et al., 1997; Petit et al., 2002). Detailed studies focussing on the Balkan Peninsula are few and mostly deal with vertebrates (Podnar et al., 2004; Kryštufek et al., 2007; Sotiropoulos et al., 2007; Ursenbacher et al., 2008), butterflies (Schmitt et al., 2006) or mountain plants (Frajman and Oxelman, 2007; Stefanović et al., 2008; Albach et al., 2009).

Southeastern Europe is recognised as a prime refugium for temperate European forest vegetation during the cold stages of the Pleistocene, together with the Iberian and Apennine Peninsulas (Comes and Kadereit, 1998; Taberlet et al., 1998; Gömöry et al., 1999; Hewitt, 2000; Hampe et al., 2003; Petit et al., 2003; Magri et al., 2006; Médail and Diadema, 2009). On the Balkan Peninsula, Pleistocene glaciation was restricted to the high massifs (Hughes et al., 2007; Milivojević et al., 2008), but as the climate was drier and more continental than at present (Horvat et al., 1974), survival of tree species was likely restricted to small areas with favourable conditions—“refugia within refugia”—as has been hypothesised for the Iberian Peninsula (Gómez and Lunt, 2007). These forest refugia were previously assumed to have been restricted to the southern tips of the Southern European peninsulas (Horvat et al., 1974; Hewitt, 2000), but recent studies found evidence for survival of tree species significantly further north than previously assumed (Stewart and Lister, 2001).

In southeastern Europe refugia have been suggested for numerous temperate tree species. Refugia of the beech (Fagus sylvatica) were suggested in Istria and adjacent areas on the Dalmatian coast, southern Bulgaria to northwestern Greece and maybe parts of the Carpathian arc (Magri et al., 2006). Possible refugia for the hornbeam (Carpinus betulus) were proposed in Romania and northern Greece (Grivet and Petit, 2003). Caucasian and European ash (Fraxinus angustifolia and F. excelsior) survived in at least two possible refugial areas, a western one in the Dinaric Alps and an eastern one stretching from the Rhodopes to the Carpathians (Heuertz et al., 2006). The hypothesis of several unconnected forest refugia during the Pleistocene is also supported by the understory vegetation that is much more diverse, regionally differentiated, and richer in endemics than in forests in central and northern Europe (Meusel and Jäger, 1992; Willner et al., 2009). However, the hypothesis of multiple forest refugia in southeastern Europe has never been tested with a herbaceous species in a phylogeographic framework.

A major force in plant evolution and diversification is polyploidy (Ramsey and Schemske, 1998, 2002; Wendel, 2000), which may in many cases be the result of secondary contact of populations differentiated in phases of allopatry, e.g., during restriction to different refugia (e.g., Petit et al., 1999). Polyploidy is recognised as an important mode of diversification by, for instance, promoting adaptation to new ecological niches or conferring reproductive isolation, which may eventually lead to speciation (Otto and Whitton, 2000). While allopolyploids may differ conspicuously from their diploid progenitors in morphology and physiology, autopolyploids that arise from the crosses within or between populations of a single species (Ramsey and Schemske, 1998) are often more difficult to distinguish on the basis of morphology alone (Levin, 1983, 2002). Recent cytogeographical studies not only indicate a higher incidence of autopolyploidy than previously thought, but also that autopolyploids often co-exist with their diploid parental populations (Husband and Sabara, 2004; Kron et al., 2007; Kolář et al., 2009). Despite the widely recognised importance of polyploidization in plant diversification and speciation, very little is known about its contribution to the high diversity on the Balkan Peninsula.

A good system to investigate diversification patterns on the Balkan Peninsula in the contexts of putative differentiation due to isolation in refugia and of polyploidy is the Veronica chamaedrys group (Plantaginaceae s.l.). Although it is widely distributed from western Europe to western Siberia, the Caucasus and Syria (Riek, 1935; Walters and Webb, 1972) and has a rather broad ecological amplitude (Dale and Causton, 1992a,b,c,d), it is a characteristic and widespread element of southeastern European forest vegetation, growing at forest margins and open forests dominated by, e.g., oaks, hornbeam or beech, in grasslands, thickets and hedges (Walters and Webb, 1972). This perennial herb is outbreeding (Goyder, 1983) with at least central European genotypes being self-incompatible (Albach, unpublished) and has the ability of clonal growth (Boutin and Harper, 1991). The V. chamaedrys group is a member of V. subgenus Chamaedrys section Chamaedrys subsection Chamaedrys (Albach et al., 2008) and comprises V. ch. subspp. chamaedrys, chamaedryoides and micans as well as V. krumovii, V. micrantha, V. orbelica and V. vindobonensis (Albach et al., 2004), whose phylogenetic relationships have not been resolved so far (Albach, 2006). With the exception of V. micrantha, which is endemic to the northwestern and central-western Iberian Peninsula (Martínez Ortega et al., 2009), all taxa occur in southeastern Europe and V. ch. subsp. chamaedryoides, V. krumovii and V. orbelica are restricted to that area. Apart from morphological differences concerning, among others, indumentum characters, these taxa do also differ karyologically. In particular, V. ch. subsp. chamaedrys is mainly tetraploid (Fischer, 1970, 1973b; Mirek and Fischer, 1986) with only a few diploids recorded from southern Austria (Fischer, 1973a), whereas V. ch. subsp. chamaedryoides and subsp. micans, as well as V. krumovii, V. orbelica and V. vindobonensis were suggested to be exclusively diploid (Fischer, 1970, 1973b, 1974; Peev, 1972; Strid and Franzén, 1984; Mirek and Fischer, 1986). Since traditional cytotaxonomy, which is often based on a few chromosome counts only, may grossly underestimate the actual intricacy of polyploid complexes in general and of cytotype distribution patterns in particular (e.g., Suda et al., 2004), the association of ploidy level and taxonomy in the V. chamaedrys group remains to be tested, in particular if diploid and tetraploid taxa within the same genetic group are spatially segregated.

Here, we explore diversification patterns within the cytologically polymorphic V. chamaedrys group in southeastern Europe employing genetic (plastid sequences and AFLP fingerprints), ploidy level and morphometric data. Specifically, we want to assess the distribution of cytotypes in this region to test the hypothesis that polyploids are more frequent at higher latitudes, were range shifts of taxa due to climatic oscillations were more pronounced than in the South. We also want to infer mode and minimum number of polyploidization events to test (i) whether polyploids originated via autopolyploidy, as frequently observed in angiosperms (Otto, 2007), or via allopolyploidy, as possible after secondary contact of once geographically isolated diploid lineages (Petit et al., 1999) and (ii) whether polyploid cytotypes originated once or multiple times. Furthermore, we want to test whether the phylogeographical pattern of the woodland herb V. chamaedrys group agrees with those of tree species found in the same vegetation types. Finally, we want to assess whether and to which extent current taxonomy reflects genetically and/or morphometrically defined lineages.

Materials and methods

Plant material

In the summers of 2006 and 2007 the V. chamaedrys group was sampled in 121 sample sites (in the following referred to as “populations”). Leaf material was collected and immediately stored in silica gel. Voucher specimens are deposited at the Faculty Centre of Biodiversity, University of Vienna, Austria (herbarium WU; voucher numbers given in Table 1). Plants were determined by Manfred A. Fischer based on Fischer (1970, 1973b, 1974, 1991), Mirek and Fischer (1986) and Peev (1972, 1995) as well as on personal experience. A detailed description of taxa delimitation as used in our study (in particular subsuming V. orbelica and V. ch. subsp. chamaedrys var. eglandulosa, Table 1) is given in Appendix 1, an overview over the characters used for determination in Appendix 2.

Table 1

Population number, taxon, collection information, co-ordinates, ploidy, cpDNA haplotype and GenBank accession numbers of 129 sampled populations of the Veronica chamaedrys group from southeastern Europe. Populations that could not be assigned to a taxon are indicated with “–”. cha, V. chamaedrys subsp. chamaedrys var. chamaedrys; coi, V. chamaedrys subsp. chamaedryoides; egl, V. chamaedrys subsp. chamaedrys var. eglandulosa; kru, V. krumovii; mic, V. chamaedrys subsp. micans; kru, V. krumovii; vin, V. vindobonensis. In the column “GenBank accession no.” the first number refers to the trnH-psbA spacer sequence, the second to the rps16-trnK spacer and the third to the rpl32-trnL spacer.

Pop. No.	Taxon	Locality (voucher number)⁎	Altitude (m a.s.l.)	Co-ordinates (E, N)	Ploidy	Haplotype	GenBank Accession No.
1	—	A: Mayerhofen (Bardy_ch1)	1650	11°54′56″, 47°03′58″	4x	h27	HM370828, HM370699, HM370570
2	egl	A: Obervellach (Bardy_ch2)	850	13°09′34″, 46°56′18″	4x	h33	HM370829, HM370700, HM370571
3	egl	A: Greifenburg (Bardy_ch3)	1525	13°11′53″, 46°46′10″	4x	h27	HM370830, HM370701, HM370572
4	egl	A: Lake Weißensee (Bardy_ch4)	1000	13°18′25″, 46°43′08″	4x	h28	HM370831, HM370702, HM370573
5	egl	A: Mt. Maltaberg (Bardy_ch5)	1250	13°30′52″, 46°57′54″	4x	h33	HM370832, HM370703, HM370574
6	vin	HR: Poreč (Bardy_ch6)	140	13°39′39″, 45°15′31″	2x	h23	HM370833, HM370704, HM370575
7	egl	HR: Pula (Bardy_ch7)	60	13°51′50″, 44°53′51″	2x	h23	HM370834, HM370705, HM370576
8o	egl	SLO: Predmeja (Frajman and Schönwetter 11129)	970	13°52′31″, 45°57′12″	4x	h52	HM370835, HM370706, HM370577
9	egl	A: St. Egyden (Bardy_ch9)	640	14°04′10″, 46°35′01″	4x	h52	HM370836, HM370707, HM370578
10	egl	HR: Lanišće (Bardy_ch10)	470	14°05′49″, 45°24′48″	2x	h26	HM370837, HM370708, HM370579
11	cha	A: Sengsengebirge (Bardy_ch11)	1400	14°12′14″, 47°48′27″	4x	h27	HM370838, HM370709, HM370580
12	egl	HR: Opatija (Bardy_ch12)	270	14°16′05″, 45°19′21″	4x	h23	HM370839, HM370710, HM370581
13	egl	SLO: Pekel pri Borovnici (Frajman and Schönwetter 11113)	345	14°22′18″, 45°53′26″	4x	h53	HM370840, HM370711, HM370582
14	—	A: Rosenau am Hengstpaß (no voucher)	800	14°23′00″, 47°42′00″	4x	h27	HM370841, HM370712, HM370583
15	—	HR: Trstenik mire (no voucher)	965	14°27′29″, 45°29′19″	4x	h56	HM370842, HM370713, HM370584
16	egl	SLO: Vrhnika pri Ložu (Frajman 12022)	590	14°30′42″, 45°41′69″	4x	h23	HM370843, HM370714, HM370585
17	egl	A: Bad Eisenkappel (Bardy_ch17)	1600	14°40′44″, 46°31′00″	4x	h35	HM370844, HM370715, HM370586
18	egl	HR: Vratnik pass (Bardy_ch18)	750	14°59′9″, 44°58′43″	4x	h24	HM370845, HM370716, HM370587
19	cha	A: Mt. Hochschwab (Bardy_ch19)	1510	15°04′37″, 47°35′19″	4x	h33	HM370846, HM370717, HM370588
20	egl	HR: Žuta Lokva (Bardy_ch20)	440	15°06′41″, 44°57′08″	4x	h25	HM370847, HM370718, HM370589
21o	egl	HR: Baške Oštarije (Frajman and Schönwetter 12063)	1150	15°10′40″, 44°30′20″	4x	h52	HM370848, HM370719, HM370590
22	egl	HR: Ogulin (Bardy_ch22)	330	15°11′12″, 45°16′21″	4x	h24	HM370849, HM370720, HM370591
23	—	HR: Plitvička jezera lakes (no voucher)	700	15°36′43″, 44°52′29″	4x	h24	HM370850, HM370721, HM370592
24	egl	HR: Zagreb (Bardy_ch24)	230	15°41′09″, 45°46′42″	4x	h55	HM370851, HM370722, HM370593
25	mic	A: Mt. Schneeberg (Bardy_ch25)	1600	15°49′58″, 47°45′14″	2x	h30	HM370852, HM370723, HM370594
26	cha	A: Mt. Schneeberg (Bardy_ch26)	1700	15°50′19″, 47°45′27″	4x	h24	HM370853, HM370724, HM370595
27o	egl	HR: Gračac (Bardy_ch27)	560	15°52′08″, 44°17′48″	4x	h24	HM370854, HM370725, HM370596
28	vin	A: Mödling (Bardy_ch28)	300	16°16′40″, 48°04′49″	2x	h38	HM370855, HM370726, HM370597
29	egl	HR: Križevci (Bardy_ch29)	130	16°27′30″, 45°58′38″	4x	h23	HM370856, HM370727, HM370598
30	vin	BiH: Mt. Troglav (Bardy_ch30)		16°36′06″, 43°56′07″	2x	h39	HM370857, HM370728, HM370599
31	vin	A: Mt. Hundsheimer Berg (Frajman and Schönwetter 11096)	340	16°56′14″, 48°08′21″	2x	h55	HM370858, HM370729, HM370600
32	vin	A: Mt. Hundsheimer Berg (Englisch_ch32)	450	16°56′20″, 48°07′57″	2x	h43	HM370859, HM370730, HM370601
33	egl	HR: Bilo Gora (Bardy_ch33)	210	17°14′19″, 45°52′36″	4x	h23	HM370860, HM370731, HM370602
34	—	H: Nemesvita (Bardy_ch34)	190	17°21′56″, 46°49′25″	4x	h33	HM370861, HM370732, HM370603
35	egl	H: Balatonboglar (Bardy_ch35)	140	17°30′02″, 46°34′16″	4x	h34	HM370862, HM370733, HM370604
36	egl	BiH: Mt. Čvrsnica (Surina and Modrić_ch36)	1530	17°37′36″, 43°39′06″	4x	h54	HM370863, HM370734, HM370605
37	vin	H: Ajka (Bardy_ch37)	420	17°42′30″, 47°02′45″	2x	h43	HM370864, HM370735, HM370606
38	vin	H: Ajka (Bardy_ch38)	420	17°42′30″, 47°02′45″	2x	h39	HM370865, HM370736, HM370607
39	vin	H: Ajka (Bardy_ch38)	420	17°42′30″, 47°02′45″	4x	h55	HM370866, HM370737, HM370608
40	egl	H: Sokorópátka (Bardy_ch40)	230	17°42′33″, 47°28′27″	4x	h55	HM370867, HM370738, HM370609
41	—	H: Pécs (Bardy_ch41)	200	18°03′33″, 46°08′32″	4x	h56	HM370868, HM370739, HM370610
42	cha	HR: Slavonski Brod (Bardy_ch42)	170	18°11′35″, 45°12′40″	4x	h31	HM370869, HM370740, HM370611
43	cha	HR: Slavonski Brod (Bardy_ch42)	170	18°11′35″, 45°12′40″	2x	h55	HM370870, HM370741, HM370612
44	—	H: Dombóvár (Bardy_ch44)	220	18°19′27″, 46°25′52″	4x	h33	HM370871, HM370742, HM370613
45	vin	HR: Vladislavci (Bardy_ch45)	90	18°34′48″, 45°26′58″	2x	h37	HM370872, HM370743, HM370614
46	egl	HR: Vladislavci (Bardy_ch45)	90	18°34′48″, 45°26′58″	4x	h33	HM370873, HM370744, HM370615
47	vin	BiH: Maglić (Bardy_ch47)	1700	18°45′10″, 43°17′52″	2x	h29	HM370874, HM370745, HM370616
48o	egl	SRB: Podromanija (Bardy_ch48)	1070	18°53′19″, 44°01′27″	4x	h6	HM370875, HM370746, HM370617
49o	egl	H: Baja (Bardy_ch49)	140	18°53′20″, 46°11′41″	4x	h23	HM370876, HM370747, HM370618
50	vin	H: Esztergom (Bardy_ch50)	110	18°57′18″, 47°47′32″	2x	h47	HM370877, HM370748, HM370619
51	—	PL: Andrychow (no voucher)	348	19°20′15″, 49°50′45″	4x	h33	HM370878, HM370749, HM370620
52	vin	SRB: Kremna (Bardy_ch52)	780	19°34′40″, 43°51′43″	2x	h52	HM370879, HM370750, HM370621
53	vin	SRB: Kremna (Bardy_ch52)	780	19°34′40″, 43°51′43″	4x	h11	HM370880, HM370751, HM370622
54	vin	H: Gyöngyö (Bardy_ch54)	390	19°57′44″, 47°50′19″	2x	h50	HM370881, HM370752, HM370623
55	—	SRB: Mt. Kapaonik (Bardy_ch55)	1770	20°50′08″, 43°18′55″	4x	h5	HM370882, HM370753, HM370624
56	vin	GR: Smolikas (Frajman and Schönwetter 11638)	1400	21°03′48″, 40°06′09″	2x	h66	HM370883, HM370754, HM370625
57o	egl	H: Békéscsaba (Bardy_ch57)	130	21°12′04″, 46°43′6″	4x	h32	HM370884, HM370755, HM370626
58	kru	MK: Nižepole—Orlove bari (Frajman and Schönwetter 11634)	1800–2150	21°12′35″, 40°57′01″	2x	h60	HM370885, HM370756, HM370627
59	kru	MK: Kopanke (Frajman and Schönwetter 11636)	1700–1890	21°12′45″, 41°01′43″	2x	h59	HM370886, HM370757, HM370628
60	coi	GR: Metsovo (Calvo et al. JC0881)a	1618	21°13′39″, 39°47′32″	2x	h67	HM370887, HM370758, HM370629
61	egl	MK: Begovo pole (Frajman and Schönwetter 11698)	1800–2200	21°24′50″, 41°43′16″	4x	h1	HM370888, HM370759, HM370630
62	coi	GR: Kato Kerasovo (Bardy_ch62)	230	21°26′07″, 38°31′06″	2x	h76	HM370889, HM370760, HM370631
63	egl	SRB: Bojnik (Bardy_ch63)	480	21°35′36″, 43°05′25″	2x	h70	HM370890, HM370761, HM370632
64	egl	MK: Pletvar (Frajman and Schönwetter 11682)	1000	21°39′18″, 41°21′58″	2x	h64	HM370891, HM370762, HM370633
65	coi	GR: Hani Ponopoulou (Bardy_ch65)	620	21°39′47″, 37°48′19″	2x	h77	HM370892, HM370763, HM370634
66	coi	GR: Vlachomandra (Bardy_ch66)	200–230	21°42′06″, 38°27′05″	2x	h74	HM370893, HM370764, HM370635
67	coi	GR: Lampia (Bardy_ch67)	860–880	21°48′45″, 37°51′47″	2x	h78	HM370894, HM370765, HM370636
68o	egl	H: Debrecen (Bardy_ch68)	130	21°52′02″, 47°45′02″	4x	h51	HM370895, HM370766, HM370637
69o	egl	SRB: Mt. Rtanj (Frajman and Schönwetter 11374)	1000–1560	21°53′30″, 43°46′30″	4x	h9	HM370896, HM370767, HM370638
70	coi	GR: Aghios Georgios (Bardy_ch70)	630–640	21°55′24″, 38°56′10″	2x	h68	HM370897, HM370768, HM370639
71	egl	SRB: Dobra (Bardy_ch71)	140	21°58′36″, 44°37′50″	4x	h45	HM370898, HM370769, HM370640
72	coi	GR: Moní Philosóphou (Bardy_ch72)	680	22°02′35″, 37°33′30″	2x	h72	HM370899, HM370770, HM370641
73	egl	SRB: Mt. Suva planina (Bardy_ch73)	940	22°05′24″, 43°13′45″	4x	h12	HM370900, HM370771, HM370642
74	coi	GR: Klitoria (Bardy_ch74)	490	22°09′04″, 37°51′47″	2x	h79	HM370901, HM370772, HM370643
75	—	GR: Kalavrita (Aedo et al. MA760495)b	1980	22°11′22″, 37°59′50″	4x	h65	HM370902, HM370773, HM370644
76	coi	GR: River Styx (Bardy_ch76)	1250	22°13′42″, 37°59′41″	2x	h76	HM370903, HM370774, HM370645
77	coi	GR: Ostrakina (Bardy_ch77)	1500	22°14′29″, 37°40′09″	4x	h62	HM370904, HM370775, HM370646
78	kru	GR: Ostrakina (Bardy_ch78)	1350–1420	22°16′11″, 37°38′00″	4x	h62	HM370905, HM370776, HM370647
79	kru	GR: Ostrakina (Bardy_ch78)	1350–1420	22°16′11″, 37°38′00″	2x	h61	HM370906, HM370777, HM370648
80	coi	GR: Poulani (Bardy_ch80)	1110–1130	22°19′36″, 38°43′17″	2x	h76	HM370907, HM370778, HM370649
81	coi	GR: Amfissa (Bardy_ch81)	870–890	22°20′39″, 38°31′54″	2x	h76	HM370908, HM370779, HM370650
82	coi	GR: Olimbos (Frajman and Schönwetter 11675)	1200–1500	22°24′08″, 40°04′53″	2x	h63	HM370909, HM370780, HM370651
83	vin	GR: Litochoro (Frajman and Schönwetter 11660)	580	22°28′44″, 40°06′54″	2x	h63	HM370910, HM370781, HM370652
84	coi	GR: Litochoro (Frajman and Schönwetter 11659)	580	22°28′44″, 40°06′54″	2x	h67	HM370911, HM370782, HM370653
85	coi	GR: Litochoro (Frajman and Schönwetter 11658)	370	22°29′28″, 40°06′49″	2x	h67	HM370912, HM370783, HM370654
86	vin	MK: Mt. Rujen (Frajman and Schönwetter 11708)	1800–2252	22°30′55″, 42°09′29″	2x	h13	HM370913, HM370784, HM370655
87	vin	MK: Mt. Rujen (Frajman and Schönwetter 11708)	1800–2252	22°30′55″, 42°09′29″	4x	h1	HM370914, HM370785, HM370656
88o	egl	BG: Vidin vis Gradec (Bardy_ch88)	50	22°33′57″, 43°59′51″	4x	h4	HM370915, HM370786, HM370657
89	coi	GR: Amfiklia (Bardy_ch89)	280	22°36′41″, 38°40′07″	2x	h73	HM370916, HM370787, HM370658
90	coi	GR: Mt. Parnassos (Bardy_ch90)	1400	22°37′11″, 38°30′18″	2x	h73	HM370917, HM370788, HM370659
91	coi	GR: Regini (Bardy_ch91)	370	22°41′07″, 38°42′43″	2x	h73	HM370918, HM370789, HM370660
92	egl	BG: Belogradčik (Bardy_ch92)	450	22°41′11″, 43°38′09″	4x	h2	HM370919, HM370790, HM370661
93o	egl	BG: Vidin (Bardy_ch93)	100	22°52′60″, 43°47′35″	4x	h1	HM370920, HM370791, HM370662
94o	egl	BG: Montana (Bardy_ch94)	325	22°53′40″, 43°27′33″	4x	h57	HM370921, HM370792, HM370663
95	vin	BG: Godec (Bardy_ch95)	650	23°03′09″, 42°59′32″	2x	h36	HM370922, HM370793, HM370664
96	vin	BG: Godec (Bardy_ch95)	650	23°03′09″, 42°59′32″	4x	h15	HM370923, HM370794, HM370665
97	vin	RO: Fildu de Jos (Bardy_ch97)	350	23°04′55″, 46°55′31″	2x	h44	HM370924, HM370795, HM370666
98o	egl	BG: Petrohan pass (Bardy_ch98)	1320	23°07′20″, 43°07′15″	4x	h20	HM370925, HM370796, HM370667
99	vin	BG: Mt. Vitosha (Bardy_ch99)	1650	23°17′00″, 42°34′00″	2x	h21	HM370926, HM370797, HM370668
100o	egl	BG: Vraca (Bardy_ch100)	250	23°20′52″, 43°16′14″	4x	h3	HM370927, HM370798, HM370669
101o	egl	BG: Mt. Vihren (Bardy_ch101)	2049	23°24′50″, 41°45′27″	4x	h58	HM370928, HM370799, HM370670
102	vin	BG: Bansko (Frajman and Schönwetter 11326)	1900	23°24′59″, 41°45′36″	2x	h13	HM370929, HM370800, HM370671
103o	egl	BG: Mt. Vihren (Sternburg_ch103)	1900–1930	23°25′00″, 41°45′00″	4x	h14	HM370930, HM370801, HM370672
104	coi	GR: Manoudi (Bardy_ch104)	110	23°28′37″, 38°46′12″	2x	h75	HM370931, HM370802, HM370673
105o	egl	BG: Vraca (Bardy_ch105)	600	23°45′39″, 43°12′31″	2x	h8	HM370932, HM370803, HM370674
106o	egl	BG: Saranci (Frajman and Schönwetter 11275)	760	23°52′56″, 42°42′47″	4x	h17	HM370933, HM370804, HM370675
107	vin	RO: Valea Lunga (Bardy_ch107)	400	24°03′08″, 46°08′31″	2x	h49	HM370934, HM370805, HM370676
108	vin	BG: Batak (Bardy_ch108)	1450	24°08′44″, 41°45′50″	2x	h22	HM370935, HM370806, HM370677
109o	egl	BG: Batak (Bardy_ch108)	1450	24°08′44″, 41°45′50″	4x	h1	HM370936, HM370807, HM370678
110o	egl	RO: Sibiu (Bardy_ch110)	550	24°12′24″, 45°48′36″	4x	h48	HM370937, HM370808, HM370679
111o	egl	RO: Sibiu (Bardy_ch111)	550	24°12′24″, 45°48′36″	4x	h48	HM370938, HM370809, HM370680
112	kru	BG: Peštera (Bardy_ch112)	800	24°16′47″, 42°00′09″	2x	h41	HM370939, HM370810, HM370681
113	coi	BG: Peštera (Bardy_ch113)	700	24°16′47″, 42°00′09″	2x	h69	HM370940, HM370811, HM370682
114	vin	BG: Trojanski pass (Bardy_ch114)	1550	24°33′57″, 42°47′28″	2x	h39	HM370941, HM370812, HM370683
115o	egl	BG: Trojanski pass (Bardy_ch115)	1350	24°43′60″, 41°40′29″	4x	h1	HM370942, HM370813, HM370684
116o	egl	BG: Loveč (Bardy_ch116)	200	24°44′38″, 43°05′09″	4x	h18	HM370943, HM370814, HM370685
117o	egl	BG: Sokelovtsi (Bardy_ch117)	1570	24°46′25″, 41°41′45″	4x	h1	HM370944, HM370815, HM370686
118o	egl	BG: Trojan (Bardy_ch118)	500	24°47′01″, 42°54′34″	4x	h17	HM370945, HM370816, HM370687
119	kru	BG: Manastir Sveta Petka (Bardy_ch119)	525	24°54′49″, 41°58′42″	2x	h42	HM370946, HM370817, HM370688
120	egl	BG: Manastir Sveta Petka (Bardy_ch120)	525	24°54′49″, 41°58′42″	4x	h1	HM370947, HM370818, HM370689
121	—	BG: Kazanlăk (no voucher)	450	25°10′33″, 42°33′39″	4x	h19	HM370948, HM370819, HM370690
122o	egl	RO: Vatra Dornei (Bardy_ch122)	875	25°13′41″, 47°19′52″	4x	h51	HM370949, HM370820, HM370691
123o	egl	BG: Šipka pass (Bardy_ch123)	1250	25°19′24″, 42°44′52″	4x	h16	HM370950, HM370821, HM370692
124o	egl	BG: Veliko Tărnova (Bardy_ch124)	350	25°27′52″, 43°00′25″	4x	h7	HM370951, HM370822, HM370693
125	—	RO: Braşov (Bardy_ch125)	900–950	25°35′34″, 45°38′02″	2x	h46	HM370952, HM370823, HM370694
126o	egl	RO: Gura Humorului (Bardy_ch126)	500	25°56′05″, 47°32′06″	4x	h40	HM370953, HM370824, HM370695
127	—	BG: Tsarevo (Albach_2007_45)	247	27°45′58″, 42°07′04″	4x	h10	HM370954, HM370825, HM370696
128	—	BG: Goritsa (Albach_2007_21)	70	27°49′13″, 42°55′01″	4x	h10	HM370955, HM370826, HM370697
129	—	BG: Ahtopol (Albach_2007_52)	50	27°57′15″, 42°04′15″	4x	h71	HM370956, HM370827, HM370698

Voucher specimen deposited at the Universidad de Salamanca (SALA).

Voucher specimen deposited at the Real Jardin Botanico de Madrid (MA).

Population corresponds to V. orbelica.

A, Austria; BG, Bulgaria; BiH, Bosnia and Herzegovina; GR, Greece, H, Hungary; HR, Croatia; MK, FYR of Macedonia; RO, Romania; SLO, Slovenia; SRB, Serbia.

Flow cytometry

DNA ploidy levels were estimated for five individuals per population. Flow cytometry was conducted with silica gel dried material following the protocol of Baranyi and Greilhuber (1996) with propidium iodide staining using a CyFlow ML (Partec GmbH, Münster, Germany) equipped with a green laser (Cobolt Samba 532 nm, Cobolt AB, Solna, Sweden) and using Pisum sativum cultivar ‘Kleine Rheinländerin’ as internal standard.

Molecular methods

Total genomic DNA was extracted from silica gel dried tissue (ca. 10 mg) of one individual per population. In eight populations with mixed cytotypes, one individual per cytotype was analysed and will in the following be referred to with a separate population identifier: this concerns populations 38/39, 42/43, 45/46, 52/53, 78/79, 86/87, 95/96 and 108/109 (Fig. 1 and Table 1). Extraction followed the CTAB-protocol of Doyle and Doyle (1987) with a few modifications: after precipitation with isopropanol and subsequent centrifugation, the DNA pellet was washed in 70% ethanol, dried at 37 °C and re-suspended in TE-buffer. The quality of the extracted DNA was checked on 1% TAE-agarose gels.

Fig. 1

Sampled populations of taxa of the Veronica chamaedrys group in southeastern Europe. Details of the collected populations are given in Table 1.

The AFLP procedure followed Vos et al. (1995) with the modifications described in Schönswetter et al. (2009). Initially, selective primers were screened using 23 primer combinations. The five final primer combinations for the selective PCR (fluorescent dye in brackets) were EcoRI (6-FAM)-ACT/MseI-CTA, EcoRI (VIC)-ACC/MseI-CAA, EcoRI (NED)-ACC/MseI-CTT, EcoRI (6-FAM)-ACT/MseI-CAT and EcoRI (VIC)-AGG/MseI-CTC. Five microlitres each of 6-FAM, NED and VIC labelled selective PCR products were combined and purified using Sephadex G-50 Fine (GE Healthcare Bio-Sciences, Uppsala, Sweden) applied to a Multi Screen-HV plate (Millipore, Molsheim, France). 1.2 μl of the elution product were mixed with 10 μl formamide (Applied Biosystems, Foster City, CA, USA) and 0.1 μl GeneScan 500 ROX (Applied Biosystems) and run on an ABI 3130x automated capillary sequencer (Applied Biosystems). Twenty individuals were replicated to calculate the error rate and to allow non-reproducible fragments to be excluded from the analysis. Raw AFLP data were aligned with the internal size standard using ABI Prism GENESCAN 3.7.1 (Applied Biosystems), and imported into GENOGRAPHER 1.6.0 (available at http://hordeum.oscs.montana.edu/genographer) for scoring. The error rate was calculated as the ratio of mismatches (scoring of 0 vs. 1) over matches (1 vs. 1) in AFLP profiles of replicated individuals (Bonin et al., 2004).

Three regions of the plastid genome, the trnH-psbA spacer (primers trnH-F, psbA-R; Tate and Simpson, 2003), the rps16-trnK spacer (primers rps16x2F2, trnK(UUU)x1; Shaw et al., 2007) and the rpl32-trnL spacer (primers rpl32-F, trnL (UAG); Shaw et al., 2007), were sequenced from one individual per population. PCR conditions for all three regions were 5 min at 95 °C followed by 30 cycles of 1 min at 95 °C, 1 min at 50 °C and 3 min at 65 °C, followed by 7 min at 65 °C. Reaction volumes of 24 μl included 8 μl REDTaq ReadyMix PCR Reaction Mix (Sigma–Aldrich, Vienna, Austria), 2 μl of each primer (10 μM), 8 μl of H2O and 4 μl of 1:10 diluted template DNA of unknown concentration. The PCR products were cleaned with Exonuclease I and Calf Intestine Alkaline Phosphatase (Fermentas, St. Leon-Rot, Germany) according to the manufacturer’s instructions. All reactions were carried out on a GeneAmp 9700 thermocycler (Applied Biosystems). BigDye Terminator chemistry (Applied Biosystems) was used according to the manufacturer’s instructions for cycle sequencing following electrophoresis on a 3130xl Genetic Analyzer capillary sequencer (Applied Biosystems). Sequences were edited with SEQMAN II 5.05 (DNAStar, Madison, WI, USA) and manually aligned using BIOEDIT 7.0.4.1 (Hall, 1999). All sequences were deposited in GenBank (Table 1). Sequences from the three plastid regions were concatenated based on the assumption that the plastid forms a single linkage group.

Data analyses

Genetic structure of AFLP data was inferred using the population mixture analysis implemented in BAPS 5.2 (Bayesian Analysis of Populations Structure; (Corander et al., 2003); http://www.abo.fi./fak/mnf/mate/jc/software/baps.html). This program, which can handle dominant markers like AFLPs under the module “clustering with linked loci” (Corander and Tang, 2007), treats both the frequencies of the markers and the number of genetically divergent groups as random variables. Stochastic optimisation is used to infer the mode of the posterior distribution. As our data set included diploid as well as tetraploids individuals, the following strategy was adopted. A mixture analysis of the diploid individuals only was conducted with the maximal number of groups (K) set to 2–10. Each run was replicated 10 times and the results were averaged according to the resultant likelihood scores. Results of the mixture analysis of diploid individuals (excluding two individuals which could not be unambiguously assigned to a single gene pool, see Section 3) were used to define “training” clusters for subsequent assignment of tetraploid individuals (“admixture based on pre-defined populations”). Admixture coefficients were estimated using 500 iterations, and the significance of these coefficients was estimated by employing the simulation strategy described by Corander and Marttinen (2006) using 50 reference individuals and 10 iterations each. An Unweighted Pair Group Method with Arithmetic mean (UPGMA) tree was inferred based on Kullback–Leibler distances (Kullback and Leibler, 1951) among clusters as implemented in BAPS 5.2. A Neighbour-Joining (NJ) analysis based on a matrix of Nei–Li distances (Nei and Li, 1979) including 2000 pseudo-replicates was conducted with TREECON 1.3b (Van de Peer and De Wachter, 1997). Using the program SPLITSTREE4 version 4.6 (Huson and Bryant, 2006), a NeighbourNet was constructed based on the same distance matrix. A Principal Co-ordinate Analysis (PCoA) based on a matrix of Jaccard distances among individuals was conducted using NTSYS-PC 2.0 (Rohlf, 1998).

Phylogenetic analysis was conducted using the approach implemented in BEAST 1.4.8 (Drummond and Rambaut, 2007), as this allows taking into account the genealogical uncertainty due to the stochastic nature of the coalescence process. As inversions introduce substitutional changes which actually are due to a structural mutation, these have been inverted prior to all analyses. The best-fit substitution model was determined using the Akaike Information Criterion (AIC) as implemented in MODELTEST 3.6 (Posada and Crandall, 1998). As the set of models until the cumulative Akaike weight exceeded 0.95 included often non-nested models with at least two substitution rates, we finally used a GTR+Γ model subsuming the proportion of invariable sites in the gamma distribution and using Jeffrey’s priors for the substitution model parameters.

Since a strict clock model was rejected (Bayes factors <−13; calculated with Tracer 1.4 available from http://tree.bio.ed.ac.uk/software/tracer/), rate evolution was modelled in a relaxed clock framework using a lognormal distribution (Drummond et al., 2006) with uniform distributions for mean and standard deviation of 0–100 and 0–10, respectively. Due to the lack of external calibrations, we used a strong prior on the substitution rate, derived from previously published substitution rates for plastid regions (Yamane et al., 2003; Smith et al., 2008), and modelled it with a normal distribution with a mean of 4 × 10−3 substitutions per site per million years and a wide standard deviation of half the mean. After initial analyses, the root was constrained to a maximum age of 10 my. We used the Bayesian Skyline Plot (Drummond et al., 2005) as the most general demographic model, as it also allows fluctuations in population size to be detected. Using different group intervals (m = 3, 5, 10) gave very similar results (absolute Bayes factors <2.4), and only those with m = 3 are shown. Stationarity of the Markov chain, which was run for 3 × 107 generations with sampling every 1000th generation, was determined using TRACER 1.4. The first 10% of sampled generations was discarded as burn-in, after which all effective sample size (ESS) values were greater than 290. A second run was conducted to confirm convergence of the Markov chain on the stationary distribution. All parameter estimates were based on these two runs combined (54,000 sampling points). A chronogram was constructed from the majority rule consensus tree calculated using PAUP 4.0b10 (Swofford, 2001) with node heights being the median values of the age estimates determined with FIGTREE 1.2.3 (available from http://tree.bio.ed.ac.uk/).

A statistical parsimony haplotype network was constructed using TCS 1.21 (Clement et al., 2000). For this analysis, insertions/deletions longer than one base pair as well as inversions were re-coded as single step mutations, and then sequence gaps were treated as a fifth character state. Mononucleotide repeats of varying length were excluded, since these are prone to homoplasy at larger geographic scales (Ingvarsson et al., 2003).

In order to compare the within-group differentiation of the lineages identified by BEAST, we calculated π, the mean number of pairwise differences (Tajima, 1983) and its variance (Tajima, 1993) with Arlequin 3.11 (Excoffier et al., 2005). Population expansions of the four cpDNA lineages were tested using mismatch distribution with a unimodal distribution indicating population expansion (Rogers and Harpending, 1992). Agreement between the observed and the expected distribution under a sudden-expansion model was tested in ARLEQUIN 3.11 (Excoffier et al., 2005) via the sum of squared differences, which, if significant at p ⩽ 0.05, indicates deviation from the expansion model (Schneider and Excoffier, 1999). Significance was assessed via a non-parametric bootstrapping procedure with 10,000 replicates.

Morphometry

In 98 of the 129 sampled populations one or two flowering or fruiting individuals, whose phenological stage allowed scoring of more than half of the morphometric characters, were available for morphometric analysis. Five sampled locations (38/39, 42/43, 52/53, 78/79 and 86/87) included two cytotypes that could not be distinguished in the morphometric analysis and were subsequently treated as polymorphic in ploidy level. One qualitative and nineteen quantitative characters (including one ratio), including those deemed diagnostic for intraspecific or specific taxa in the V. chamaedrys group, were scored. Unless stated differently, length characters were measured in mm, and densities in number per mm2; leaf characters were scored on the leaf pair subtending the lowermost inflorescences while the stem indumentum was determined from the internodium below the lowermost inflorescences. Qualitative characters were scored once, whereas quantitative characters other than leaf dimensions were measured five times and averaged. Of quantitative characters that were highly correlated (p < 0.001) only the one with the fewest missing characters was retained for the final analysis, thus resulting in twelve quantitative characters (Table 2; data matrix in Appendix 3).

Table 2

Morphological characters or ratios of characters employed in a morphometric analysis of 98 populations of the Veronica chamaedrys group from southeastern Europe.

Abbreviation	Morphological character or ratio
KG	Calyx hairs glandular (1) or not (0)
KHD	Density of hairs on the calyx lobes
KHL	Length of hairs on the calyx lobes
LLW	Length/maximum width of lamina
LP	Length of the petiole
LS	Length of the style
LT	Number of teeth on one half of the lamina
LTB	Vertical distance from first tooth to the base of the lamina (negative if lamina base is rounded to cuneate, positive if base is cordate)
LTLW	Length/width of lamina tooth at maximum width line
LUBD	Density of hairs on the upper side of the lamina next to the apex of the lamina
LUBL	Length of hairs on the upper side of the lamina next to the apex of the lamina
SRD	Density of hairs between the two opposite lines of hairs on the stem
SRL	Length of hairs on the two opposite lines of hairs on the stem

The seven excluded characters were (1) the length of the stem that was correlated with the number of teeth on the lamina (LT), (2) the length of the petiole one leaf pair above and (3) one leaf pair below the leaf pair subtending the lowermost inflorescence, both correlated with the length of the petiole subtending the lowermost inflorescence (LP), (4) the length and (5) density of hairs on the upper side of the lamina next to the base of the lamina and (6) the length and (7) density of hairs on the lower side of the lamina next to the apex of the lamina, which were correlated with the length (LUBL) and density (LUBD) of hairs on the upper side of the lamina next to the apex of the lamina, respectively.

Of the remaining twelve quantitative characters (the qualitative character was treated separately) a matrix of pairwise Gower distances (S15; Legendre and Legendre, 1998) was calculated with R Package 3.0 (Legendre and Vaudor, 1991) and served as basis for a Principal Co-ordinate Analysis (PCoA) using NTSYS-PC 2.0 (modules DCENTER and EIGEN; Rohlf, 1998). To check for correspondence of the morphological characters with the first two principle co-ordinates, a test for association between paired samples using Spearman’s rho was carried out.

Additionally, a linear discriminant analysis of the twelve quantitative characters was carried out using R 2.7.2 (R Development Core Team, 2008). As pre-defined groups we used the Eastern, Western and Southern Groups as inferred from AFLP analyses; individuals pertaining to the Krumovii Group (see Section 3) were excluded. Missing values were replaced by the average values of the respective group (Eastern Group, 6%; Western Group, 7%; Southern Group, 1%). In order to ascertain characters differentiating the Southern Group from the remaining samples, we conducted a second discriminant analysis with these two groups only. The results of the coefficients that contributed most to the linear discriminator were virtually identical with the coefficients of the first discriminator resulting from the analysis with three groups and are thus not shown.

To test for significant morphological differentiation between the Southern Group and the remaining samples (see Section 3), t-tests were undertaken and subsequently boxplots produced for all characters studied. The latter analyses were conducted with R 2.7.2 (R Development Core Team, 2008).

Results

Flow cytometry analyses yielded histograms with mean CVs of G1 peaks of the sample and internal standard of 9.1% and 4.2%, respectively. DNA ploidy levels (Suda et al., 2006) inferred from measured genome sizes revealed that DNA diploids (for simplicity in the following referred to as “diploids”) and DNA tetraploids (referred to as “tetraploids”) were present throughout most of the study area (Fig. 2). However, tetraploids were more frequent in the North of the distribution area, resulting in a significant association of cytotypes with latitude (Spearman’s rho = 0.337, p < 0.001). In contrast, no altitudinal separation could be detected (rho = −0.029, p = 0.756).

Fig. 2

Distribution of ploidy level (diploid, thin outline; tetraploid, thick outline) and of three genetic clusters derived from BAPS analysis of AFLP markers in 129 populations of the Veronica chamaedrys group from southeastern Europe. The insert shows a UPGMA tree of Kullback–Leibler distances among clusters. Population identifiers are given in Fig. 1.

AFLPs

We scored 468 AFLP fragments ranging from 101 to 498 base pairs. The error rate amounted to 1.6% and was thus well within the range deemed acceptable by Bonin et al. (2004). Mixture analysis including only diploid individuals identified three clusters (data not shown). Subsequent admixture analysis assigned most individuals to a single cluster; only two individuals (from populations 58 to 59) were admixed and were excluded from the “training clusters”. The subsequent admixture analysis including diploid as well as tetraploid individuals (Fig. 2) identified three clusters, hereinafter called Southern, Western and Eastern Cluster, each containing both diploid and tetraploid individuals. Allele distributions of the Western and the Eastern Cluster were more similar to each other than they were to the Southern Cluster (insert in Fig. 2). Admixture was encountered among all clusters and mostly concerned tetraploid individuals, albeit not all tetraploids were admixed, while only three diploid individuals were admixed. The three clusters were geographically somewhat separated: The Western Cluster was most prominent in the northwest, the Eastern Cluster in the east and the Southern Cluster in the south, but especially Western and Eastern Clusters were widely geographically overlapping (Fig. 2). In the NeighbourNet (Fig. 3), two groups (Southern and Western Group) corresponding to two clusters resolved by the admixture analysis were separated along strongly weighted splits. The remaining populations fell into the weakly differentiated Eastern and Krumovii Group, the latter being positioned between the Southern and the Eastern Group. All groups comprised more than one taxonomic entity, and V. vindobonensis was found in all four groups. In the Neighbour-Joining analysis (not shown) both the Southern and the Western Group received high bootstrap support (100 BS and 96 BS, respectively). Neither the Krumovii nor the Eastern Group had bootstrap support >50% in agreement with the lack of strongly weighted splits.

Fig. 3

NeighbourNet diagram of AFLP data constructed for 129 populations of the Veronica chamaedrys group from southeastern Europe. Splits with weight <0.5 have been omitted to aid legibility. Numbers along branches are bootstrap values above 50% derived from a Neighbour-Joining analysis and are given for major branches only. Bold branches mark tetraploid individuals.

cpDNA sequencing

The trnH-psbA sequences were 295–304 bp, the rps16-trnK sequences 701–709 bp and the rpl32-trnL sequences 843–867 bp long; the length of the concatenated alignment was 1978 bp. Concatenating the sequences yielded 79 haplotypes (Table 1). Phylogenetic analysis using BEAST (Fig. 4) revealed three major clades, from here on referred to as Southern, Widespread and Eastern + Western Lineage, and haplotype h71 found in a tetraploid individual. The oldest differentiation among the lineages (age given as median and its 95% highest posterior density interval) occurred 0.81 (0.22–3.01) mya. The Southern and the Widespread Lineage (posterior probabilities of 0.82 and 0.97, respectively) diversified nearly simultaneously 0.40 (0.10–1.55) and 0.37 (0.10–1.36) mya, respectively. The Southern Lineage (h59–h70, h72–h79) was distributed on the southern Balkan Peninsula (Greece and northerly adjacent areas; Fig. 5). The Widespread Lineage (h39–h58) extended from the southeastern Alps (Austria, Slovenia) to the Bulgarian Stara Planina and was the sole lineage in the Carpathians (Romania). A deep phylogenetic split was evident in the Eastern + Western Lineage (posterior probability 0.84). The Eastern and the Western Lineages, which separated 0.47 (0.14–1.72) mya, started to diversify at 0.33 (0.09–1.20) and 0.29 (0.07–1.05) mya, respectively. The Western Lineage (h23–h38) extended from the southeastern Alps (Austria, Slovenia) to the northern Dinaric Mountains (Croatia, Bosnia and Herzegovina) and the eastern Pannonian Plain (Hungary), the Eastern Lineage (h1–h22) from western Serbia to Bulgaria. Both lineages geographically overlapped with the Widespread Lineage. Rapid population expansion for all lineages was inferred with the Bayesian Skyline Analysis (Fig. 4) to have occurred about 0.05 mya.

Fig. 4

Plastid DNA haplotype diversity encountered in 129 populations of the Veronica chamaedrys group from southeastern Europe. (Left) Simplified majority rule consensus tree from relaxed clock Bayesian analysis with the software BEAST. Node heights correspond to median ages (see text for details); nodes with age estimates younger than 0.04 my are not shown and are thus not distinguishable from unresolved lineages. Numbers along branches are Bayesian posterior probabilities; identical haplotypes or unresolved polytomies are collapsed as triangles, their vertical extension being proportional to the number of individuals. The bar depicts the ploidy level of a lineage (light grey, diploid; dark grey, tetraploid). The insert at the bottom shows the Bayesian Skyline Plot (median and 95% high posterior density limits). (Right) Statistical parsimony network of plastid DNA haplotypes. Haplotypes that were not sampled are shown with small open circles. Outline shading corresponds to the three main AFLP groups presented in Fig. 3; individuals of the Krumovii Group are marked with an asterisk. Di- and tetraploids have thin and thick outlines, respectively. The four main lineages identified by Bayesian analysis (at the left) are indicated.

Fig. 5

Geographical distribution of plastid DNA haplotypes encountered in 129 populations of the Veronica chamaedrys group from southeastern Europe. In the insert, the most frequently sampled haplotypes (>2 occurrences) are listed, and their position is indicated in the statistical parsimony network onto which the four main lineages identified by BEAST are overlaid (green, Eastern Lineage; blue, Western Lineage; grey, Widespread Lineage; yellow, Southern Lineage). The map shows the geographic distribution of the eleven most frequent haplotypes (big dots) as well as the rarer ones (small dots, coloured as the overlay in the insert).

In contrast to the clearly divergent Southern Lineage, differentiation among the Widespread, Eastern and Western Lineages identified by BEAST was only weakly reflected in the parsimony network (Fig. 4). Relatively few haplotypes were frequent and many were found in a single population only. The most frequent haplotypes were found in three separate regions (Fig. 5): (i) central to southern Greece (h67, h73 and h76), (ii) Bulgaria and westerly adjacent areas (h1), and (iii) northern Croatia to eastern Austria and southern Poland (h23, h24, h27, h33 and h55); the only exception of this pattern were the haplotypes h39 and h52 from the Widespread Lineage, which occurred in the West as well as in the East. Most haplotypes corresponding to the Western Lineage were closely related and maximally three mutational steps apart. Haplotype lineages were only partly congruent with AFLP-derived clusters (BAPS) or groups (NeighbourNet). The Southern Lineage formed a concise group in the network and comprised all individuals from the Southern Group as well as two haplotypes corresponding to the Krumovii Group (h59 and h60) and haplotype h71, which belongs to the Eastern Group. With the exception of haplotypes h36–h38 sampled in populations pertaining to the Eastern Group, the Western Lineage comprised mainly haplotypes of the Western Group, albeit not all. The Eastern Lineage included exclusively haplotypes which occurred in the Eastern Group, whereas the Widespread Lineage contained mostly haplotypes of the Eastern Group plus four sampled in the Western Group (h52, h53, h55 and h56). Haplotypes sampled in individuals from the Krumovii Group were spread over the network and were partly internal to those corresponding to the Southern Lineage and partly related to those of the Eastern as well as the Widespread Lineages (Fig. 4).

The mean number of pairwise differences (π) calculated for sets of haplotypes of the Western, Eastern, Widespread and Southern Lineages amounted to 2.13 ± 1.21 (mean ± variance), 4.06 ± 2.08, 4.85 ± 2.45 and 6.34 ± 3.09, respectively. The mismatch distributions showed a unimodal distribution only for the Eastern and Western Lineages, the distribution of the Widespread Lineage was close to unimodal, and that of the Southern Lineage was multimodal; only the Eastern distribution was significant (p = 0.02).

The PCoA of the morphometric data (Fig. 6) showed a weak segregation of di- and tetraploids along the first factor (16.8%; second factor: 10.3%), but failed to reflect the taxonomic grouping. There was, however, a weak clinal separation mainly along the first factor with V. chamaedrys subsp. chamaedryoides and subsp. chamaedrys var. eglandulosa being most distant and subsp. chamaedrys var. chamaedrys as well as V. vindobonensis and V. krumovii being intermediate. The first factor was positively correlated with different length characters (e.g., length of hairs on different organs (LUBL; SRL; KHL)), the distance between the first tooth and the lamina base (LTB), the length of the style (LS) and the number of teeth on the lamina (LT), but negatively correlated with the densities of hairs (KHD; LUBD; SRD), the length-to-width ratio of the tooth at the broadest part of the lamina (LTLW) or the length of the petiole (LP). The second factor showed positive correlation with the length of the petiole (LP) and with the length-to-width ratio of the lamina (LLW), but was negatively correlated with the density of hairs on different organs (KHD; LUBD), with the length-to-width ratio of the tooth at the broadest part of the lamina (LTLW), the length of the style (LS), as well as the number of teeth on the lamina (LT) (Appendix 4). Differentiation along the first axis was also detected when the only qualitative character, i.e., presence of glandular hairs on the calyx (KG), was plotted (plot not shown); importantly, all individuals of subsp. chamaedryoides had a glandular hairy calyx. When displaying the genetic groups obtained from the NeighbourNet of the AFLP data onto the morphometric plot, individuals of the Southern, Western and Krumovii Groups were separated but the Eastern Group blurred the separation (Fig. 6).

Fig. 6

Principal Co-ordinate Analysis of a matrix of pairwise Gower distances based on 12 morphological characters scored for 187 individuals from 98 populations of the Veronica chamaedrys group from southeastern Europe. Symbols of the upper graph correspond to taxa as used in Fig. 1; diploid individuals are white, tetraploids are shaded in grey. Sampling sites with both cytotypes are marked both grey and white. The graphs at the bottom show the correspondence with the four AFLP groups as indicated in Fig. 3.

The linear discriminant analysis (Appendix 5) showed a similar pattern as the PCoA, but separation of the groups was more pronounced. The proportion of among-group variance explained by the first linear discriminator was 67%, that of the second linear discriminator 33%. The first discriminator separated the Southern Group from the Western and Eastern Group, which showed considerable overlap along the second discriminator. The coefficients that contributed most to the first discriminator were the length of hairs on the two opposite lines of hairs (SRL, 2.00) and the ratio length to maximum width of lamina (LLW, 1.35) as well as the ratio length to width of lamina tooth at maximum width (LTLW, −1.02). The coefficients that contributed most to the second discriminator were the length of hairs on the calyx lobes (KHL, 2.41), the length of hairs on the upper side of the lamina next to the apex of the lamina (LUBL, −1.19), the length of the petiole (LP, 0.84) and the ratio length to maximum width of lamina (LLW, 0.69).

Significant differences between the Southern Group and the remaining samples (excluding the Krumovii Group, for details, see Section 4) overlapped only partly with the strength of the coefficients contributing to the first discriminator (details in Appendix 6). Altogether, individuals of the Southern Group had denser but shorter hairs and the teeth at the margin of the lamina were longer and narrower than in the Western and the Eastern Groups.

Discussion

Distribution of cytotypes within the V. chamaedrys group in southeastern Europe

Di- and tetraploid cytotypes within the V. chamaedrys group occurred all over southeastern Europe, as has been assumed previously (Mirek and Fischer, 1986). There was, however, a significant association with latitude in the study area as diploids were predominant in the south, whereas in the north both di- and tetraploids occurred (Fig. 2). Previous studies revealed similar patterns of cytotype distribution across Europe. Autotetraploids of Rorippa amphibia (Luttikhuizen et al., 2007), for example, occur in northern Europe, whereas diploids grow in Central and Western Europe. A similar pattern was also found in Plantago media in which diploids had small fragmented distributions in southern and southeastern Europe but tetraploids were more widespread and continuously distributed over the continent (Van Dijk and Bakx-Schotman, 1997).

The cytotypes of the V. chamaedrys group did neither show altitudinal segregation (Table 1), as observed in Lotus corniculatus s.l. (Gauthier et al., 1998) or Taraxacum sect. Ruderalia (Calame and Felber, 2000), nor was there obvious habitat differentiation between the cytotypes, as has been shown for Anthoxanthum alpinum (Felber-Girard et al., 1996), Claytonia virginica (Lewis and Suda, 1976), Dactylis glomerata (Lumaret et al., 1987) or Senecio carniolicus (Hülber et al., 2009). Our results are in concert with studies on Solidago altissima, where autopolyploids co-exist with diploids in a broad zone of overlap (Halverson et al., 2008), or Ranunculus adoneus with no ecological differentiation between di- and tetraploids (Baack and Stanton, 2005).

Tetraploid individuals of the V. chamaedrys group were mainly of autopolyploid origin and evolved several times independently within each genetic group. This is supported by the unambiguous assignment of tetraploid individuals to clusters including also diploids (Fig. 2), the NeighbourNet diagram uniting di- and tetraploid individuals with heavily weighted splits (Fig. 3), as well as the co-occurrence of di- and tetraploids in plastid lineages (Figs. 4 and 5). Our results are thus adding to an ample body of evidence that multiple origins of polyploids are rather the rule than the exception (e.g., Doyle et al., 1990; Soltis et al., 2003; Albach, 2007). Tetraploids mainly evolved within various northern genetic lineages, and only a few originated within the southern refugium (Figs. 2–4). Although the exact causes remain elusive, we speculate that the more stable conditions on the southern Balkans (Tzedakis et al., 2002) may have prevented the successful establishment of polyploids. More pronounced climatic oscillations on the northern Balkan Peninsula and adjacent areas in the north, in contrast, triggered massive range shifts of forest vegetation that also concerned the V. chamaedrys group and were probably a driving force allowing tetraploid establishment.

Whereas our data strongly support prevalence of autopolyploidy in the V. chamaedrys group in southeastern Europe, they are inconclusive with respect to the role of allopolyploidy. The admixed state of many tetraploid individuals in the Bayesian clustering analysis (Fig. 2) may indicate that they originated by crossings between divergent diploid lineages. Alternatively, these crosses may as well have occurred at the tetraploid level, in line with the close geographical proximity of non-admixed tetraploid populations. In any event it needs to be pointed out that the NeighbourNet, which is a method of choice for uncovering reticulate relationships (Huson and Bryant, 2006), only identified two tetraploid individuals (populations 23 and 36; Fig. 3) as being intermediate between the Western and Eastern Groups.

Spatio-temporal evolution of the V. chamaedrys group

Slowly mutating, maternally inherited (Zhang et al., 2003) plastid DNA sequence data (Figs. 4 and 5) and rapidly evolving, nearly entirely nuclear-derived (Bussell et al., 2005) AFLP marker (Figs. 2 and 3) do not yield congruent patterns of spatio-temporal evolution of the V. chamaedrys group, most probably because both marker systems trace differentiation processes taking place at different time horizons. Integrating information obtained from both marker systems, however, allows reconstructing a detailed scenario from the earliest diversification in the mid-Pleistocene over the relatively constant diversification during the last 0.4 my (Fig. 4) to the allopatric differentiation accompanying isolation in refugia during the Last Glacial Maximum and postglacial secondary contact leading to admixture of gene pools (Fig. 2).

The oldest traceable phylogeographic pattern in the V. chamaedrys group is the simultaneous differentiation of three plastid DNA lineages in the mid-Pleistocene (Fig. 4). A fourth lineage is constituted by haplotype h71 sampled in population 129 from the Bulgarian Black Sea coast. Since this sample is not differentiated in the AFLP data set (Figs. 2 and 3) and our sampling from the southeastern-most Balkan Peninsula is very scarce, we here only point out the possibility of a divergent plastid lineage probably distributed along the southern coast of the Black Sea. Reconstructing the palaeo-distribution of the three main lineages is straightforward for only one of them. The Southern Lineage is restricted to the south of the Balkan Peninsula (Fig. 5) except for two isolated occurrences in southern Serbia (population 63) and central Bulgaria (population 113). In contrast to the other lineages, the multimodal mismatch distribution indicates distributional stasis and lack of population expansion. The distinctness of the genetic entity from the southern Balkans including its two northern outliers is corroborated by AFLP data (Figs. 2 and 3), and the nearly perfect congruence between the data sets with contrasting inheritance strongly suggests that hybridisation with other entities was at least rare. It remains to be tested whether the observed integrity has only historical reasons or if also ecological or crossing reasons are involved. The mountain ranges of northern Greece were glaciated (e.g., Hughes et al., 2007) and probably acted as strong barriers, contributing to the genetic isolation of the Southern Lineage/Group. Accordingly, the southern-most Balkan Peninsula is well-known as one of the centres of endemicity in Europe (Tan et al., 2001).

The distributions of the Widespread and Eastern + Western Lineages (Fig. 4) overlap throughout their entire distribution area from the Alps to the Rhodope Mountains and both lineages sometimes co-occur in immediate vicinity (e.g., populations 52 and 53; Fig. 5). Our data do not allow discriminating if both lineages either co-occurred in sympatry (“ancient polymorphism”) after their split or if they regained secondary contact after an initial allopatric phase. The Carpathians, however, were obviously never reached by the Western + Eastern Lineage and were colonised by the Widespread Lineage in the course of the first diversification event at about 0.4 my (Fig. 4), as evidenced by the presence of micro-lineages that have their most internal haplotypes in the Carpathians (h44–h46 and h48–h51; Figs. 4 and 5).

The Eastern + Western Lineage was disrupted into the vicariant Western and Eastern Lineages at about 0.47 my (Fig. 4). The two entities occur from the Alps and the Western Hungarian Plains to the northern Dinaric Mountains and from the Stara Planina to the Rhodope Mountains, respectively (Fig. 5). Diversification within the two lineages started roughly simultaneously at about 0.29 and 0.33 my ago (Fig. 4). Longitudinal vicariance is also seen in the AFLP data (Figs. 2 and 3). The Eastern Group is distributed from the Hungarian Plain and the Carpathians to the Stara Planina and the Rhodope Mountains. It comprises the entire Eastern Lineage and most haplotypes of the Widespread Lineage (Fig. 5). In contrast to the other main AFLP groups, the Eastern Group lacks strongly weighted splits and bootstrap support (Fig. 3). It is unclear if the lack of divergence is a consequence of the relatively recent contact of two formerly differentiated genetic entities corresponding to the Widespread and Eastern plastid Lineages.

The northwestern genetic entities encountered in the V. chamaedrys group, the Western Lineage and the Western Group, are roughly congruent and occur from the Alps and the Hungarian Plain to the Dinaric Mountains. The Western Lineage possesses a set of closely related haplotypes some of which were sampled frequently (Figs. 4 and 5), and thus exhibits the lowest value of pairwise differences, amounting to only one third of that of the Southern Lineage. All this suggests either a comparatively young age or a stronger bottleneck within the contemporary distribution range. The Western AFLP group comprises mostly tetraploid individuals (Fig. 2); diploids were scattered throughout the range but were most frequent in Istria.

Allopatric differentiation during the Last Glacial Maximum presumably shaped the pattern seen in the AFLP data (Figs. 2 and 3). It is mainly governed by three refugia in the southern, northwestern and eastern parts of southeastern Europe whose existence was already previously proposed for associated tree species based on macrofossil charcoal, pollen evidence and phylogeographical studies. (1) The southern Balkan Peninsula (Pindus mountain range and central Greece) provided buffered climatic oscillations (Tzedakis et al., 2002) and refugia of tree and shrub species such as beech and hornbeam were hypothesised (Grivet and Petit, 2003; Magri et al., 2006), overlapping with the Southern Cluster/Group. (2) On the eastern Balkan Peninsula (Eastwood, 2004) and the Carpathians (Willis and van Andel, 2004), refugia for beech, hornbeam and ash were proposed (Grivet and Petit, 2003; Heuertz et al., 2006; Magri et al., 2006) within the distribution area of the Eastern Cluster/Group. (3) Refugia for beech and ash (Heuertz et al., 2006; Magri et al., 2006) have been suggested on the northwestern Balkan Peninsula and adjacent areas, namely Istria, the Pannonian Basin and the valley of River Danube in Austria (Willis and van Andel, 2004), roughly congruent with the extent of the Western Group/Cluster. Generally, the forest refugia were either separated by glaciated mountain massifs as the climatic snow line was about 1000 m lower than today (e.g., Bognar and Prugovečki, 1997; Reuther et al., 2007), or by Artemisia-steppes prevailing in the lowlands (Willis, 1994).

Triggered by most probably postglacial climatic amelioration, the V. chamaedrys group underwent massive range expansions as suggested by the Bayesian Skyline Plot (Fig. 4) for ca. 0.05 mya. As substitution rates derived from phylogenetic studies with usually deeper time coverage (older than 1–2 my) will often be gross underestimates of more recent substitution rates (Ho et al., 2005) the obtained age estimate will be biased towards older ages; thus the rapid expansion may have occurred after the Last Glacial Maximum. Subsequently, as suggested by the admixture analysis (Fig. 2) but not by the NeighbourNet diagram (Fig. 3), secondary contact between Western and Eastern Groups partially broke down previously existing genetic differentiation by hybridisation involving frequent independent polyploidization (Fig. 2).

Taxonomic considerations

Our data show that current taxonomy does not reflect genetic lineages (Figs. 2–5). Furthermore, there is no tight association between taxa and ploidy level, as previously suggested, but diploids and polyploids are found in all recognised taxa as well as in all major genetic lineages (Fig. 4). Morphometric analysis does not recover current taxonomic entities, either (Fig. 6). Therefore, our data suggest changes to the currently used taxonomic framework of the V. chamaedrys group. The only entity that is congruently identified by different data types and thus merits taxonomic recognition is V. ch. subsp. chamaedryoides from the southern Balkan Peninsula. It is mainly (but not exclusively, see Fig. 3) diploid and its genetic distinctness was confirmed by both the AFLP and the cpDNA data sets (Figs. 2, 3, 5). Although subsp. chamaedryoides is only weakly separated in the PCoA plot of the morphometric data (Fig. 6), the Southern Group is relatively well separated from Eastern and Western Groups in the discriminant analysis (Appendix 5). Using only slightly overlapping characters (Appendix 6), subsp. chamaedryoides can thus be separated from the other taxa. Specifically, indumentum characters as well as the shape of the teeth and the lamina can be used to distinguish subsp. chamaedryoides (Appendix 6). The identifiability of this taxon is evidenced by the fact that about 80% of the voucher specimens have been correctly identified in the initial phase of our study.

The two widespread taxa V. ch. subsp. chamaedrys and V. vindobonensis both comprise diploid and tetraploid cytotypes (Fig. 3) and could neither be differentiated by molecular data (Figs. 2–5) nor by PCoA of morphometric data (Fig. 6)—even in the discriminant analysis Eastern and Western Groups overlap strongly (Appendix 5). Additionally, in contrast to the northwestern distribution edge of V. vindobonensis in eastern Austria (Fischer et al., 2008), the vegetation of the sampling sites did not reveal obvious habitat differentiation between these entities and both entities were found from wet to dry habitats, in open forests or meadows (K. Bardy, personal observations). Although the AFLP data (Fig. 3) suggested two genetic entities, these can neither be paralleled with morphometric data nor are they congruent with a priori determinations of voucher specimens (Figs. 3 and 6). Although we cannot entirely exclude the possibility that morphological differential characters might still be found, we favour to treat these two lineages as informal phylogeographical groups within V. ch. subsp. chamaedrys, which thus becomes a morphologically and cytologically variable taxon comprising diploids as well as mainly autotetraploid derivatives.

Although V. orbelica was described as diploid (Peev, 1972), we found only tetraploid individuals morphologically resembling this taxon. As continuous variation of characters connects V. orbelica (Peev, 1972) and V. ch. subsp. chamaedrys var. eglandulosa (Mirek and Fischer, 1986; see Appendix 1), and it is not genetically divergent, either, V. orbelica should be included in subsp. chamaedrys var. eglandulosa. Similarly, genetic indistinctness and the lack of morphometric differentiation of V. ch. subsp. micans suggest merging this taxon with subsp. chamaedrys.

Veronica krumovii morphologically resembles V. ch. subsp. chamaedrys (Fig. 6) but differs by glandular hairs of varying length scattered over the whole plant. According to the AFLP data (Fig. 3), the majority of populations appears to be intermediate between the Southern and the Eastern Group, and may actually be hybrids between subsp. chamaedrys and chamaedryoides as circumscribed here. Plants morphologically belonging to V. krumovii also occur disjunctly in southern Greece (Fig. 1), suggesting that the erratically occurring, but conspicuous diagnostic character is under simple genetic control and consequently of no taxonomic value.