Annals of Botany 118: 1071–1088, 2016 doi:10.1093/aob/mcw142, available online at www.aob.oxfordjournals.org
On the origins of Balkan endemics: the complex evolutionary history of the Cyanus napulifer group (Asteraceae) n ova1 and Jaromır Kucera1 Katarına Olsavska1,*, Marek Slovak1, Karol Marhold1,2, Eliska Stub 1
Institute of Botany, Slovak Academy of Sciences, D ubravsk a cesta 9, SK-84523 Bratislava, Slovakia and 2Department of Botany, Faculty of Science, Charles University, Ben atsk a 2, CZ-128 01 Praha 2, Czech Republic *For correspondence. E-mail [email protected] Received: 20 January 2016 Returned for revision: 2 March 2016 Accepted: 3 June 2016 Published electronically: 21 July 2016
Background and Aims The Balkan Peninsula is one of the most important centres of plant diversity in Europe. Here we aim to fill the gap in the current knowledge of the evolutionary processes and factors modelling this astonishing biological richness by applying multiple approaches to the Cyanus napulifer group. Methods To reconstruct the mode of diversification within the C. napulifer group and to uncover its relationships with potential relatives with x ¼ 10 from Europe and Northern Africa, we examined variation in genetic markers (amplified fragment length polymorphisms [AFLPs]; 460 individuals), relative DNA content (40 ,6-diamidino-2phenylindole [DAPI] flow cytometry, 330 individuals) and morphology (multivariate morphometrics, 40 morphological characters, 710 individuals). To elucidate its evolutionary history, we analysed chloroplast DNA (cpDNA) sequences of the genus Cyanus deposited in the GenBank database. Key Results The AFLPs revealed a suite of closely related entities with variable levels of differentiation. The C. napulifer group formed a genetically well-defined unit. Samples outside the group formed strongly diversified and mostly species-specific genetic lineages with no further geographical patterns, often characterized also by a different DNA content. AFLP analysis of the C. napulifer group revealed extensive radiation and split it into nine allopatric (sub)lineages with varying degrees of congruence among genetic, DNA-content and morphological patterns. Genetic admixture was usually detected in contact zones between genetic lineages. Plastid data indicated extensive maintenance of ancestral variation across Cyanus perennials. Conclusion The C. napulifer group is an example of a rapidly and recently diversified plant group whose genetic lineages have evolved in spatio-temporal isolation on the topographically complex Balkan Peninsula. Adaptive radiation, accompanied in some cases by long-term isolation and hybridization, has contributed to the formation of this species complex and its mosaic pattern. Key words: Cyanus napulifer group, Centaurea; Compositae, Balkan Peninsula, AFLP, cpDNA, flow cytometry, multivariate morphometrics, endemics, allopatric speciation, hybridization, incomplete lineage sorting.
INTRODUCTION The Balkan Peninsula is one of the most important European centres of plant diversity. Compared with any other region of Europe, it harbours the richest flora, including a high percentage of endemics (e.g. Eastwood, 2004; Krystufek and Reed, 2004; Stevanovic et al., 2007; Hewitt, 2011). Its astonishing diversity was likely generated by a combination of factors and processes. One of the most important is undoubtedly the refugial character of the Balkan Peninsula. The paleo-environmental conditions in this region allowed the survival of relic taxa, including some of tertiary origin (Tzedakis, 2004; Peev and Delcheva, 2007). Furthermore, the territory of today’s Balkan Peninsula served as a crossroad for Asian lineages and taxa during east–west colonization of Europe since the Early Oligocene (32 million years ago; Hewitt, 2011; Manafzadeh et al., 2014). The heterogeneous climate and complex physical geography of the Balkan Peninsula provided a wide range of conditions for different species and vegetation types to evolve in. The patchy, topographically varied landscape of the Balkans, with its mosaic of
mountain ranges, islands and peninsulas that acted as migration barriers for most organisms, significantly contributed to the formation of Balkan biodiversity (Stevanovic, 1996; Reed et al., 2004; Meshinev, 2007; Peev and Delcheva, 2007; Velikov and Stoyanova, 2007; Radford and Ode´, 2009). The Balkan flora is nevertheless still largely unexplored and calls for further investigation (but see Perny et al., 2005; Frajman and Oxelman, 2007; Stefanovic et al., 2008; Frajman and Schneeweiss, 2009; Scho¨nswetter and Schneeweiss, 2009; Bardy et al., 2010; Kolarcik et al., 2010; Kucera et al., 2010; Surina et al., 2011; Kuzmanovic et al., 2013; Lakusic et al., 2013; Ronikier and Zalewska-Gałosz 2014). The Balkan Peninsula accommodates almost 27% (18 out of 67 species) of perennial taxa of Cyanus sect. Protocyanus (Dobrocz.) Olsavska (Olsavska et al., 2013) and is considered its evolutionary centre (Hellwig, 2004). A large proportion of Cyanus perennials from the Balkans belong to the taxonomically extremely complex group of C. napulifer (Olsavska et al., 2013), which is the focus of the present study. The C. napulifer group, as defined by Olsavska et al. (2013), comprises five
C The Author 2016. Published by Oxford University Press on behalf of the Annals of Botany Company. V
All rights reserved. For Permissions, please email: [email protected]
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
1072
Balkan endemics – C. napulifer (Rochel) Sojak (s.str.), C. orbelicus (Velen.) Sojak, C. velenovskyi (Adamovic) Wagenitz & Greuter, C. nissanus (Petrovic) Sojak and C. tuberosus (Vis.) Sojak; one species occurring in the Balkan Peninsula and Turkey [C. thirkei (Sch. Bip.) Holub]; and one recently described Turkish species [C. eflanensis Kaya & Bancheva (Greuter, 2006–2009; Kaya and Bancheva, 2009)]. Although their mutual phylogenetic relationships remain unresolved, previous genetic studies indicate their close evolutionary relationships. However, C. thirkei seems to occupy a separate position, and no information is available on C. eflanensis (Borsic et al., 2011; Lo¨ser, 2012). Morphologically, the C. napulifer group is clearly characterized by a suite of common features, mainly by rhizomatous and/or tuberous roots and narrow stem leaves. All the taxa are diploids and share the base chromosome number of x ¼ 10 (no data are available for C. eflanensis), save for a few triploids (2n3x30) found in C. nissanus and C. thirkei (Olsavska et al., 2013). Two other Balkan species, C. diospolitanus Bancheva & Stoyanov and C. pseudoaxillaris (Stef. & T. Georgiev) Holub, previously included in Cyanus sect. Napuliferi (Stef. & T.Georgiev) Bancheva & Raimondo (Bancheva and Raimondo, 2003; Bancheva and Stoyanov, 2009), were proved to be genetically, karyologically and morphologically unrelated to the C. napulifer group (Borsic et al., 2011; Olsavska et al., 2013) and so were not included in the present study. The C. napulifer group represents a promising model for the study of homoploid speciation promoted mainly by geographical isolation and adaptation to new ecological niches. However, previous efforts to infer the evolutionary relationships of Cyanus perennials and the mechanisms of their diversification largely failed because of the insufficient or conflicting phylogenetic signals inferred from genetic sequence markers (Borsic et al., 2011; Lo¨ser, 2012). To uncover the evolutionary relationships among members of this group and the factors responsible for their diversification, we applied a combination of approaches that has proved to be effective in similar studies (e.g. Scho¨nswetter et al., 2009; Bardy et al., 2010; Kucera et al., 2010; Olsavska et al., 2011). Specifically, we employed fast-evolving multilocus amplified fragment length polymorphism (AFLP) markers, chloroplast DNA (cpDNA) data, flowcytometric data and multivariate morphometrics. The aims of our current study were: (1) to reveal which populations/taxa constitute the C. napulifer group and to uncover their relationships with other Cyanus taxa with x ¼ 10 occurring in Europe and Northern Africa; (2) to infer the overall genetic, karyological and morphological patterns of Balkan members of the C. napulifer group and to examine their congruence; and (3) to identify the evolutionary processes and factors that triggered speciation and divergence within the C. napulifer group and to infer biogeographical implications. MATERIALS AND METHODS Sampling design
The sampling strategy was focused on members of the Cyanus napulifer group occurring predominantly in the Balkan Peninsula. A total of 38 population samples of C. napulifer (s.str.), C. orbelicus, C. velenovskyi, C. nissanus, C. thirkei and
C. tuberosus were collected from natural habitats in 2009–13 (May/August) throughout the Balkan Peninsula [Bulgaria (BU), Croatia (HR), Greece (GR), Former Yugoslav Republic of Macedonia (MK), Serbia (RS)] (Fig. 1, Table 1). To elucidate the evolutionary relationships of the C. napulifer group, material from 30 populations of other morphologically similar European and North African species of the section Protocyanus, with the base chromosome number of x ¼ 10 , was collected for AFLP and DNA content analyses (Table 1, Fig. 1): C. pindicola (Griseb.) Sojak, C. epirotus (Halacsy) Holub, C. lingulatus (Lag.) Holub, C. graminifolius (Lam.) Olsavska, C. fuscomarginatus (K.Koch) Greuter, and C. pichleri (Olsavska et al., 2013). Plant material
Intact fresh leaves were collected from three to eight plants per locality and dried in silica gel for AFLP analyses. Roots/rhizomes from up to 14 plants were transferred to an experimental garden in Banska Bystrica, Slovakia (48 450 08900 N, 19 090 29000 E, 390m a.s.l.) for flow cytometry estimation. Depending on the abundance of flowering individuals, 10–52 plants of the C. napulifer group (not including C. thirkei) were collected at original localities for morphometric analyses. The type of bedrock and biotope were recorded during the field research. Each population of the species under investigation had at least several hundreds of individuals, and the taxa under study are neither endangered nor protected by law. Voucher specimens were deposited in the herbarium of the Institute of Botany of the Slovak Academy of Sciences (SAV). DNA extraction and AFLP generation
Total genomic DNA was extracted using the DNeasy Plant Mini Kit (Qiagen, Du¨sseldorf, Germany). The AFLP procedure (Vos et al., 1995) followed a modified protocol described by Mered’a et al. (2008). To obtain highquality AFLP profiles, we tested 38 selective primer combinations using four individuals from different, geographically distant taxa. Final analyses were performed using four primer combinations labelled with different fluorescent dyes and producing clear, unambiguous and polymorphic profiles: EcoRIACA-(6-FAM)/MseI-CTG, EcoRI-AAG-(VIC)/MseI-CAG, EcoRI-AAC-(NED)/MseI-CTG and EcoRI-ACA-(PET)/MseICAC (Applied Biosystems). Amplification products were pooled and fragment analysis was performed on an ABI 3100 Avant capillary sequencer of the BITCET Consortium, Department of Molecular Biology, Comenius University, Bratislava. The internal size standard GeneScan –500 LIZ (Applied Biosystems) was employed for size calibration of the AFLP profiles obtained. AFLP trace files were read and analysed using DAX software (Van Mierlo Software Consultancy, Netherlands). Only unambiguously scorable markers in the size range of 50–500 base pairs (bp) were scored manually as present (1) or absent (0). The reproducibility of AFLP profiles was determined using 35 replicates (76% of the final dataset). This dataset comprised randomly chosen re-extracted samples of all taxa under study, which were analysed and scored independently. The error rate (Bonin
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group Cyans napulifer group C. napulifer + orbelicus BU
C. orbelicus BU/RS C. velenovskyi GR C. velenovskyi RS C. nissanus RS C. orbelicus Vich C. tuberosus GR/MK C. tuberosus MK C. tuberosus HR
1073
other Cyanus taxa: C. pindicola C. epirotus C. thirkei C. pichleri C. fuscomarginatus C. lingulatus C. graminifolius
AL - Albania, BA - Bosnia and Herzegovina, BG - Bulgaria, HR - Croatia, GR - Greece, ME - Montenegro, MK - Former Yugoslav Republic of Macedonia, RS - Serbia
FIG. 1. Distribution map of Cyanus samples used in molecular, DNA-content and morphometric analyses. Samples of C. napulifer group are assigned to genetic lineages (coloured solid lines) and sublineages (coloured dashed lines) according to the results of AFLP data analyses. Symbols with thick lines are used for populations of presumably hybrid origin.
et al., 2004) was calculated as the ratio of mismatches (1 versus 0) to matches (1 versus 1) between replicated samples.
AFLP data analyses
To reveal the genetic pattern of the taxa under study, we assembled four AFLP datasets. The first dataset, AFLP dataset 1, comprised all analysed taxa (460 individuals); the second one, AFLP-dataset 2 (264 individuals), comprised only accessions of the C. napulifer group excluding C. thirkei, which early on turned out to be genetically clearly distinct from the rest of the group (see the Results section). AFLP dataset 3 (233 individuals) comprised accessions of the C. napulifer group excluding four genetically admixed populations (TRI 151, TRI 153, TRI 265 and TRI 273; see the Results section). For the purposes of multispecies coalescent analyses, we created AFLP dataset 4 (120 individuals), composed of two randomly selected individuals for each population analysed except four genetically admixed ones (TRI 151, TRI 153, TRI 265 and TRI 273; see also below in this section). The overall genetic pattern of relationships among populations and individuals was investigated using the following approaches based on AFLP datasets 1–3: (1) neighbour-joining tree (NJ; Saitou and Nei, 1987) based on Sørensen’s ( ¼ Dice) similarities transformed into a distance matrix with bootstrap
support for clades and subclades computed using 5000 replicates; (2) split network (neighbour-net, NN) analysis based on Sørensen’s similarities transformed into a distance matrix (Huson and Bryant, 2006); and (3) principal coordinates analysis (PCoA; Krzanowski, 1990) based on Jaccard’s similarity matrix. The NJ tree analysis and PCoA were performed using Famd 1108 beta software (Schlu¨ter and Harris, 2006), while NN was conducted using SplitsTree 4 (Huson and Bryant, 2006). To search for genetically homogeneous groups within the C. napulifer group encompassing only diploid taxa (AFLP dataset 2), Bayesian model-based clustering as implemented in Structure 223 (Falush et al., 2007) was employed. The settings for the analysis were the same as those used in Slovak et al. (2012a, b). The R script Structure-sum-2009 (part of AFLPdat; Ehrich, 2006) was used to calculate the statistics necessary for determining the correct K value (Evanno et al., 2005). Graphical outputs of analyses were generated using Clumpp v11.1 (Jakobsson and Rosenberg, 2007) and Distruct (Rosenberg, 2004). To infer hierarchical relationships between genetic lineages of the C. napulifer group and related species with x ¼ 10, we computed a species trees in a multispecies coalescent framework using SNAPP v2.2.0 (Bryant et al., 2012). A priori ‘species’ delimitation was done with respect to the results of NN, NJ, PCoA and structure analyses, with the exceptions of
1074
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
TABLE 1. Locality details, including geographical coordinates, altitude, the date of collection and the collector(s) of plants of Cyanus species included in the molecular (AFLP), DNA content flow cytometric and morphometric analyses Population (voucher) code*
Locality and collection details†
Cyanus napulifer group Cyanus napulifer þ orbelicus BU (central Bulgarian) genetic lineage TRI 214 Bulgaria; Stara Planina Mts, Korduna peak; 42 450 1900 N, 24 040 3600 E; 1605 m; (A) 2.viii.2011; (B) 8.vii.2013, MP & KO TRI 215 Bulgaria; Stara Planina Mts, 2 km SE from Vezhen cottage; 42 450 47500 N, 24 240 10800 E; 1869 m; 2.viii.2011, KO TRI 216 Bulgaria; Stara Planina Mts, Troyanski prokhod, Beclemento; 42 460 3700 N, 24 360 5200 E; 1601 m; 3.viii.2011, KO TRI 217 Bulgaria; Stara Planina Mts, Troyanski prokhod, Mt Viloto; 42 460 4700 N, 24 390 5600 E; 1603 m; 3.viii.2011, KO TRI 300 Bulgaria; Stara Planina Mts, Buzludzha; 42 440 07000 N, 25 230 51800 E; 1400 m; 9.vii.2013, MP & KO TRI 222 Bulgaria; Rila Mts, Belmeken dam; 42 090 5100 N, 23 460 5700 E; 1963 m; 9.viii.2011, KO TRI 153 Bulgaria; Rila Mts; (A) 25 km NE from Granchar Hut; 42 080 0800 N, 23 360 5400 E; 2294 m; 19.vi.2010, KO; (B) 12 km N from Granchar Hut; 42 070 5100 N, 23 350 4500 E; 2277 m; 20.vi.2010, KO Cyanus orbelicus BU/RS (Bulgarian–Serbian) genetic sublineage TRI 296 Serbia; Kopaonik Mts Karaman; 43 170 34000 N, 20 490 47700 E; 1910; 6.vii.2013, MP & KO TRI 260 Serbia; Mt Besna Kobila; (A) below the peak; 42 320 57900 N, 22 130 56200 E; 1662 m; (B) N of the peak; 42 330 29900 N, 22 140 09400 E; 1559 m; 15.vi.2012, CL & KO TRI 298 Bulgaria; Ossogovska planina; (A) Kulin Kamak; 42 110 41400 N, 22 350 54300 E; (B) Gramadite; 42 110 17900 N, 22 360 31600 E; 1636 m; 1710 m; 7.vii.2013, MP & KO TRI 302 Bulgaria; Rhodopi Mts, Golyam Perelik Mt; 41 360 00800 N, 24 340 48500 E; 2140 m; 11.vii.2013, MP & KO Cyanus velenovskyi GR (Greek) genetic sublineage TRI 144 Greece; Pella, E slope of Mt Kaimaktsalan; 40 550 0300 N, 21 480 4100 E; 2160 m; 13.vi.2010, MO & KO TRI 241 Greece; Florina, Vernon Mts, W slope of Mt Vitsi; 40 380 53000 N, 21 220 47100 E; 1921 m; (A) 22.v.2012, MO & KO; (B) 13.vii.2013, MP & KO TRI 305 Greece; Florina, Varnous Mts, Mt Vigla; 40 480 35300 N, 21 150 07000 E; 1990 m; 13.vii.2013, MP & KO Cyanus velenovskyi RS (Serbian) genetic sublineage TRI 258 Serbia; Stara Planina Mts, Mt Midzor; (A) near the ski resort Babin zub; 43 220 45600 N, 22 370 3900 E; 1592 m; 15.vi.2012, KO; (B) S slope; 43 220 4400 N, 22 400 4800 E; 1662 m; 15.vi.2012, CL Cyanus nissanus RS genetic sublineage TRI 256 Serbia; Svrljiska Planina Mts, Mt Popova glava, summit area; 43 210 23400 N, 22 050 30200 E; 1008 m; (A) 14.vi.2012, CL & KO; (B) 19.v.2013, MO, JO & KO Cyanus orbelicus Vich (Vichren) genetic sublineage TRI 151 Bulgaria; Pirin Mts, Mt Vihren, above Vihren cottage; 41 450 2700 N, 23 240 1700 E; 2440 m; (A) 17.vi.2010; (B) 10.viii.2010; (C) Bakushevata Mura; 41 450 5400 N, 23 240 4400 E; 2217 m; 10.viii.2011, KO Cyanus tuberosus GR/MK (Greek–Macedonian) genetic sublineage TRI 240 Greece; Larisa, Deskati, E of Loutro; 39 560 55400 N, 21 550 40800 E; 732 m; 21.v.2012, MO & KO TRI 170 Greece; Grevena, Vourinos Mts, S from Exarchos, 40 090 0200 N, 21 370 3000 E, 684 m; 28.v.2011, KO TRI 149 Greece; Kozani, Pieria Mts, SE of Velvendo; 40 130 5200 N, 22 060 3900 E; 1278 m; 16.vi.2010, KO TRI 162 Greece; Kilkis, Paiko Mts; (A) S of Livadhia; 40 590 4000 N, 22 170 4600 E; 1160 m; (B) NW of Livadhia; 41 010 N, 22 160 3000 E; 1180 m; 26.v.2011, KO TRI 238 Greece; Thesaloniki, Mt Hortiatis; 40 350 57800 N, 23 060 01200 E; 826 m; 20.v.2012, MO & KO TRI 237 Greece; Serres; (A) Lailas; 41 140 22700 N, 23 340 22000 E; 1393 m; (B) between Kato Vrondous and Eleonas; 41 130 57200 N, 23 390 45300 E; 915 m; 19.v.2012, MO & KO TRI 235 Greece; Drama, Voras Mts; (A) Falakro; 41 170 36600 N, 24 020 29700 E; 1467 m; (B) Volokas; 41 180 10900 N, 23 590 34900 E; 1012 m; 19.v.2012, MO & KO TRI 273 Macedonia/Greece; Kozuf Mts; (A) Mt Keci Kaj; 41 110 48000 N, 22 140 33800 E; 1715 m; 22.vi.2012, KO; (B) above ski centre; 41 110 01100 N, 22 120 59700 E; 1912 m; 22.vi.2012, CL; (C) Mt Mala Rupa; 41 090 52500 N, 22 140 37600 E; 1863 m; 22.vi.2012, CL Cyanus tuberosus MK (central Macedonian) genetic lineage
No. of plants analysed AFLP
Flow cytometry N/P‡ (ploidy; ST§)
Morphometric analyses
8
4/3 (2n 20; SL)
20
5
3/3 (2n 20; SL)
15
8
3/3 (2n 20; SL)
12
8
4/3 (2n 20; SL)
11
8
6/0 (2n 20; SL)
20
8 8
2/4 (2n 20; SL) 2/6 (2n 20; SL)
15 31
8 8
5/0 (2n 20; SL) 0/4 (2n 20; SL)
20 41
8
8/0 (2n 20; SL)
32
8
5/0 (2n 20; SL)
20
8
0/4 (2n 20; SL)
20
8
5/4 (2n 20; SL)
12
8
5/0 (2n 20; SL)
20
8
3/4 (2n 20; SL)
29
8
1/2 (2n 20; SL) 0/2 (2n 30; SL)
20
7
3/7 (2n 20; SL)
24
8
2/0 (2n 20; SL)
15
5
0/4 (2n 20; SL)
13
8
0/4 (2n 20; SL)
18
7
11/0 (2n 20; SL)
35
8
2/0 (2n 20; SL)
14
8
5/0 (2n 20; SL)
27
8
1/0 (2n 20; SL)
35
8
3/0 (2n 20; SL)
52
(continued)
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
1075
TABLE 1. Continued Population (voucher) code*
Locality and collection details†
Macedonia; Dren Mts, pass above Stavica; 41 150 03500 N, 21 350 41100 E; 1132 m; (A) 21.vi.2012, CL & KO; (B) 23.v.2013, MO, JO & KO TRI 271 Macedonia; Dren Mts; (A) Veprchani; 41 160 11900 N, 21 450 24400 E; 832 m; 22.vi.2012, KO; (B) 20.v.2013, MO, JO & KO; (B) above Dunje; 41 270 27200 N, 21 460 07400 E; 1479 m; 22.vi.2012, CL TRI 262 Macedonia; Ramno Mts, above Brodec; (A) 42 080 35800 N, 21 260 57500 E; 1100 m; 15.vi.2012, CL & KO; (C) 42 090 50900 N, 21 260 53000 E; 1386 m; 25.v.2013, MO, JO & KO TRI 264 Macedonia; Vodno Mts, Skopje, near Millenium cross; (A) 42 580 11700 N, 21 230 54600 E; 922 m; 16.vi.2012, CL; (B) 41 570 583400 N, 21 230 433400 E; 1022 m; 24.v.2013, MO, JO & KO TRI 266 Macedonia; Babuna Mts, Prisad pass; 41 260 54000 N, 21 370 50200 E 41 260 2400 N, 21 360 5500 E; 809–1100 m; (A) 19.vi.2012, CL & KO; (B) 24.v.2013, MO, JO & KO TRI 265 Macedonia; Babuna Mts, above Nezhilovo; 41 400 11800 N, 21 290 2900 E; 2064 m; 17.vi.2012, CL Cyanus tuberosus HR (Croatian) genetic lineage TRI 247 Croatia; Lika, Mt Postak, along a road to Rasticevo; 44 130 28400 N, 16 050 48800 E; 752 m; 5.vi.2012, MO, JO & KO TRI 252 Croatia; Split-Dalmatia County, above Sinj; (A) near Zelovo; 43 450 04400 N, 16 330 12200 E; 800 m; 7.vi.2012, MO, JO & KO; (B) along a road to Sutina; 43 420 04000 N, 16 360 52500 E; 402 m; 7.vi.2012, MO, JO & KO TRI 253 Croatia; Mosor Mts, near Rasca; 43 330 48800 N, 16 420 33800 E 43 330 22300 N, 16 420 22500 E; 312–361 m; 8.vi.2012; MO, JO & KO Cyanus pindicola (Griseb.) Sojak TRI 145 Greece; Pella, NE slope of Mt Piperitsa; 40 510 3800 N, 21 440 3200 E; 1812 m; 14.vi.2010, MO & KO TRI 148 Greece; Pieria, Oros Olympos Mts, below Zolotas cottage; 40 040 5600 N, 22 220 5000 E; 1709 m; 15.vi.2010, MO & KO TRI 166 Vermio Mts, Seli 40 320 0300 N, 22 010 3000 E; 1662 m; 27.v.2011, KO TRI 168 Vourinos Mts 40 110 56200 N, 21 390 27700 E; 1330 m; 28.v.2011, KO TRI 175 Greece; Larisa, Kato Olympos Mts, 2 km NE from Kallipefki; 39 580 2500 N, 22 280 4600 E; 1141 m; 29.v.2011, CL & KO TRI 267 Greece; Galicica Mts; (A) Galicica Mt; 40 570 43200 N, 20 490 41300 E; 1843 m; 20.vi.2012, KO; (B) Magaro Mt; 40 560 22100 N, 20 490 18700 E; 2 163 m; 20.vi.2012, CL Cyanus epirotus (Halacsy) Holub TRI 182 Greece; Ioannina, Tzoumerka Mts, Kryopigi; 39 270 0900 N, 21 070 3900 E; 1443 m; 2.vi.2011, CL & KO TRI 183 Greece; Ioannina, Tzoumerka Mts, close to Mparos pass; 39 360 36600 N, 21 080 32900 E; 1676 m; 2.vi.2011, CL & KO TRI 185 Greece; Ioannina, Mitsikeli Mts, above Ioannina; 39 430 01600 N, 20 530 0700 E; 1669 m; 3.vi.2011, CL & KO TRI 187 Greece; Ioannina, Timfi Mts, above the road to Vradeto; 39 550 0300 N, 20 480 2100 E; 1738 m; 3.vi.2011, CL & KO Cyanus thirkei (Sch.Bip.) Holub TRI 228 Bulgaria; Stara Planina Mts, Varbishki prohod, above Sadovo, 42 530 17600 N, 26 390 05800 E; 686 m; (A) 15.v.2012, MO & KO; (B) 20.v.2013, MO, JO & KO TRI 229 Bulgaria; Balgarova near Burgas; 42 370 58300 N, 27 160 48200 E; 99 m; (A) 16.v.2012, KO; (B) 21.v.2013, MO, JO & KO TRI 230 Bulgaria; Otmanli near Burgas; 42 250 49100 N, 27 320 59100 E; 52 m; 16.v.2012, MO & KO TRI 294 Bulgaria; Karnobatska Planina Mts, between Valchin and Lazarevo; 42 440 26600 N, 26 530 54700 E; 271 m; 21.v.2013, MO, JO & KO TRI 295 Bulgaria; Karnobatska Planina Mts, N of Aytos; 42 450 04600 N, 27 150 16700 E; 347 m; 21.v.2013, MO, JO & KO Cyanus pichleri TRI 179 Greece; Peloponnese, Mt Panachiko; 38 150 5200 N, 21 510 0600 E; 1146 m; 31.v.2011, CL & KO TRI 186 Greece; Ioannina, Timfi Mts, Vradeto; 39 530 53200 N, 20 460 23800 E; 1339 m; 3.vi.2011, CL & KO TRI 188 Greece; Kastoria, between Nea Kotyli and Kato Nestorio; 40 220 16600 N, 21 020 41400 E; 1268 m; 4.vi.2011, CL & KO TRI 231 Bulgaria; Strandzha Mts, Golyamo Bukovo, above road E79; 42 140 28300 N, 27 060 36600 E; 340 m; 17.v.2012, MO & KO TRI 239 Greece; Mt Chortiatis, above village of Chortiatis; 40 350 57800 N, 23 060 01200 E; 826 m; 20.v.2012, MO & KO Cyanus fuscomarginatus (K.Koch) Greuter TRI 270
No. of plants analysed AFLP
Flow cytometry N/P‡ (ploidy; ST§)
Morphometric analyses
8
4/2 (2n 20; SL)
20
8
5/4 (2n 20; SL)
10
8
8/4 (2n 20; SL)
19
8
5/1 (2n 20; SL)
20
8
5/3 (2n 20; SL)
20
8
1/4 (2n 20; SL)
16
8
2/4 (2n 20; SL)
17
8
6/8 (2n 20; SL)
11
8
3/4 (2n 20; SL)
8
6
0/3 (2n 20*; SL)
–
6
0/1 (2n 40*; SL)
–
5 5 6
3/0 (2n 40; SL) 3/0 (2n 40; SL) 0/5 (2n 40*; SL)
– – –
6
1/0 (2n 20; SL)
–
3
0/3 (2n 20; SL)
–
6
0/3 (2n 20; SL)
–
6
0/5 (2n 20; SL)
–
6
0/3 (2n 20; SL)
–
8
–
8
3/0 (2n 20; BP) 0/7 (2n 30; BP) 3/4 (2n 20; BP)
–
8 8
0/4 (2n 20; BP) 6/0 (2n 20; BP)
– –
8
6/0 (2n 20; BP)
–
6
0/2 (2n 40; SL)
–
5
0/3 (2n 40; SL)
–
5
0/2 (2n 40; SL)
–
6
0/1 (2n 40; SL)
–
6
0/4 (2n 40; SL)
–
(continued)
1076
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group TABLE 1. Continued
Population (voucher) code*
Locality and collection details†
No. of plants analysed AFLP
Ukraine; Crimea, Aj-Petrinskaja Jajla, Mt Aj-Petri; 44 270 2900 N, 34 030 3700 E; 1178 m; 7.vii.2012, KO TRI 282 Ukraine; Crimea, Catyr-Dag Jalja, near Marmornaja cave; 44 470 2500 N, 34 160 4800 E; 973 m; 9.vii02012, KO Cyanus lingulatus (Lag.) Holub JAV Spain; Sierra de Javalambre; 40 130 31 N, 00 200 56 W; 1853 m; 3.vi.2010; SS, KM & JZ PUR Spain; Sierra del Moncayo, Purujosa; 41 420 47 N, 1 440 54 W; 1280 m; 8.v.2011; SS, KM & JZ ˆ in-Leuch; 33 124210 N, TRI 289 Morocco; Moyen Atlas, between Sources de l0 Oum-er-Rbia and A 05 206370 W; 1734 m; 1052013; AG, JK, MS & KO TRI 291 Morocco; Rif Mts, Tleta Ketama; 34 524060 N, 04 336100 W; 1573 m; 11.v.2013; AG, JK, MS & KO Cyanus graminifolius (Lam.) Olsavska TRI 65 France; Hautes Alpes, Montagne de Ce´u¨se; 44 310 0200 N, 05 520 0800 E; 1572 m; 15.vii.2007, KO, MP & DD TRI 101 France; Alpes de Hte Provence, Montagne du Cheiron, Pic de Fourneuby; 43 490 29400 N, 06 530 32500 E; 1422 m; 6.vii.2008; KO, MP & VK TRI 196 Italy; Pollino Mts, Colle del Dragone; 39 540 46400 N, 16 070 30800 E; 1531 m; 14.vi.2011, KO (SAV: TRI 196) TRI 198 Italy; Calabria, La Sila Mts, 4 km SW from Moccone; 39 200 05700 N, 16 230 56500 E; 1613 m; 14.vi.2011, KO (SAV: TRI 198) TRI 278
Flow cytometry N/P‡ (ploidy; ST§)
Morphometric analyses
8
0/4 (2n 40; SL)
–
8
0/3 (2n 40; SL)
–
8 7
0/5 (2n 20; BP) 0/5 (2n 20*; BP)
– –
8
6/0 (2n 20; BP)
–
8
1/0 (2n 20; BP)
–
8
0/5 (2n 40*; BP)
–
6
0/7 (2n 40; BP)
–
6
0/3 (2n 40; BP)
–
6
0/5 (2n 40; BP)
–
Populations of C. napulifer group are assigned to genetic lineages and sublineages according to the results of AFLP data analyses. *Population codes of populations with presumably hybrid origin are in bold. Collectors: AG, A. Guttova; IH, I. Hodalova; VK, V. Kolarcik; CL, C. Lo¨ser; KM, K. Marhold; JK, J. Kucera; JO, J. Olsavska; MO, M. Olsavsk y; MS, M. Slov a k; KO, K. Ol s avsk a ; MP, M. Pern y ; SS, S. Spaniel; JZ, J. Zozomov a -Lihov a . ‡ Number of newly analysed plants/number of plants analysed in a previous study (Olsavska et al., 2013) § Standard used in flow cytometry analysed: BP, Bellis perennis; SL, Solanum lycopersicum. †
C. lingulatus and C. graminifolius. These were regarded as a single species to ensure the monophyly of the group. Due to the large computational demands of the algorithm, we analysed a subsample comprising two randomly selected individuals for each population (AFLP dataset 4). We ran two analyses with different h priors in order to allow for different current and ancestral population sizes, as applied in Kolar et al. (2016): (1) mean h prior of 0043 (c distribution, a ¼ 15, b ¼ 35) and (2) mean h prior ¼ 01 (c distribution, a ¼ 12, b ¼ 110 prior for large population sizes); the remaining parameters were left at their default values. The analyses were checked for convergence using Tracer v1.6, making sure that Bayesian runs reached an effective sample size exceeding 200 after burn-in. The posterior distributions of species trees were visualized using DensiTree v2.2.0 (Bouckaert and Heled, 2014). Genetic diversity within and among the taxa/populations/genetic (sub)lineages under study (AFLP datasets 1, 2 and 4) was estimated using the following statistical parameters: (1) total number of AFLP multilocus genotypes; (2) average number of AFLP fragments (6 s.d.); (3) average proportion of pairwise differences between individuals (Nei’s gene diversity, DNei; Nei and Li, 1979) using the R script AFLPdat (Ehrich, 2006); and (4) percentage of polymorphic markers (P%) using Pop-Gene 132 (Yeh et al., 1997). Genetic divergence among analysed taxa/populations/genetic (sub)lineages was explored by calculating the number of (5) private fixed fragments (diagnostic); (6) private fragments (exclusive); (7) rare fragments (present at a frequency of < 10% of investigated individuals); and (9)
unique fragments, using FAMD 1108 beta software (Schlu¨ter and Harris, 2006) and by calculating of the frequency downweighted marker value (Scho¨nswetter et al., 2005), as implemented in AFLPdat (DW1; Ehrich, 2006). Within the C. napulifer group (AFLP dataset 3), we also explored within- and between-group variance using analyses of molecular variance (AMOVA). AMOVAs were performed using Arlequin 311 (Excoffier et al., 2005) based on Euclidean pairwise distances and a significance test with 10 000 permutations.
cpDNA data analyses
To address the evolutionary relationships among perennial representatives of the genus Cyanus, we explored cpDNA sequences deposited in the GeneBank database (NCBI, http:// www.ncbi.nlm.nih.gov/nucleotide), particularly those generated by C. J. Lo¨ser, G. Akaydin and F. H. Hellwig (University of Jena, Germany). We used 114 sequences of trnC-ycf6 and 99 sequences of 30 rps16-50 trnK(UUU) cpDNA regions (Shaw et al., 2005, 2007) for 15 Cyanus taxa. This dataset included not only an overwhelming majority of species analysed in our study with x ¼ 10, but also species with x ¼ 11 (Supplementary Data Table S1). All cpDNA sequences were manually edited and aligned using BIOEDIT (v7.0.4.1; Hall, 1999). Ambiguously aligned regions were removed with Gblocks (Castresana, 2000; Talavera and Castresana, 2007). The method of statistical parsimony was
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group employed to generate a haplotype network (TCS v.1.18; Clement et al., 2000) depicting the general pattern of relationships among the haplotypes. TCS was run with a 95% connection limit and with indels coded as missing data.
Flow cytometric analyses
The dataset resulting from relative nuclear DNA content analyses comprised both newly generated (136 individuals) and previously published data (194 individuals; Olsavska et al., 2013); between one and 14 plants per population were analysed by flow cytometry using the fluorochrome 40 ,6-diamidino-2phenylindole (DAPI, specific for A–T base pairs) (Table 1). DAPI flow cytometry was chosen because a previous study, which was focused on the C. triumfetti and C. montanus groups (Olsavska et al., 2012), proved a strong correlation between absolute DNA content (estimated using the intercalating fluorochrome propidium iodide) and relative DNA content (estimated using DAPI) (r ¼ 099, P5% in genome size were performed to confirm the reliability of the estimated values. The relationship between chromosome numbers and DNA content values was verified using previously published chromosome counts (Table 1; Olsavska et al., 2013). Relative DNA content was calculated as the ratio of the position of the G1 peak of the standard [Solanum lycopersicum ‘Stupicke´ polnı rane´’, 2C DNA ¼ 196pg (Dolezel et al., 1992) or Bellis perennis, 2C DNA ¼ 338pg (Scho¨nswetter et al., 2007)] to that of the sample. When B. perennis was used as the internal standard, the resulting relative DNA content was recalculated relative to the S. lycopersicum standard based on the ratio between B. perennis and S. lycopersicum (0578) calculated from eight estimations performed at different times. To compare samples differing in ploidy level, relative DNA content per monoploid genome was used. The Tukey–Kramer test (Tukey’s test for unequal sample size; Tukey, 1977) was used to test for differences in relative DNA content between populations of the C. napulifer group (triploids from population 256, investigated only in the previous study, were not included) and in monoploid relative DNA content between taxa or genetic (sub)lineages. All analyses were carried out using STATISTICA 91 (StatSoft Inc., 2013). For populations of C. nissanus and C. thirkei (48 plants in total), ploidy levels were estimated by flow cytometry for each individual subjected to AFLP because in a previous study (Olsavska et al., 2013), two ploidy levels (2x and 3x) were identified in some populations of these taxa. Sample preparation followed the procedure described above, but silica-dried material was used, and one to three plants were analysed simultaneously; such flow cytometric estimations were not used to calculate relative DNA content.
1077
Morphometric analyses
Morphological variability was studied using multivariate morphometric analyses based on 710 individuals from 34 populations of the C. napulifer group (Table 1). A total of 24 quantitative and nine qualitative characters were measured or scored on herbarium material, and three ratios were subsequently computed (Supplementary Data Table S2). Depending on the number of states recorded (Table S2), multistate qualitative characters were binary-coded for further analyses as two (ROOT, SI, LBWS, LBAS, LLI, LUI, BUL, BBL), three (FECO) or five (SECR) dummy binary variables. Each multistate qualitative character with s states was decomposed into s– 1 dummy variables (to avoid complete linear dependency of the sth variable). First, all quantitative characters were tested for normality of distribution (Shapiro–Wilk test). Subsequently, Spearman’s (non-parametric) correlation coefficient was computed for all characters because several characters departed from a normal distribution. No pair of characters showed a strong correlation (arbitrary value >095) that would potentially distort further computations. To gain insight into the overall phenetic relationships and to estimate the level of congruence among morphological, genome size and genetic (AFLP) variation, non-metric multidimensional scaling (NMDS, based on Euclidean distance and standardization by variable range; Kruskal, 1964) and canonical discriminant analysis (CDA; Krzanowski, 1990) were performed. Population samples used in the NMDS analyses were characterized by mean values of all 40 characters. The CDA analyses were based on individual plants and 26 quantitative and ratio characters. Qualitative characters were excluded because they were uniform in some groups of individuals. The analyses were based both on the whole set of data (NMDS 1) and on the partial datasets representing genetic (sub)lineages as revealed by Bayesian clustering (NMDS 2–6, CDA 1–5). The analyses were carried out using STATISTICA 91 (StatSoft Inc., 2013), Canoco 50 (ter Braak and Smilauer 2012) and SAS 93 (SAS Institute Inc., 2011).
Ecological and phenological trait analyses
To reveal the life strategies of populations of the C. napulifer group, the following ecological and phenological data were collected for each population: altitude, type of bedrock, type of biotope and flowering time (Table S2). These data, together with selected morphological characters describing the root system (Table S2), were subjected to NMDS analyses (NMDS 7) in the same way as morphological characters alone. Nei’s gene diversity, frequency down-weighted marker values and mean DNA content of particular populations were entered as supplementary characters. RESULTS AFLP analyses
The overall reproducibility of AFLP dataset 1 (entire dataset) was high because the error rate was minute, being only 019%. The analysis yielded 442 unambiguously scorable fragments
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
1078
Bayesian clustering based on AFLP dataset 2 (the C. napulifer group not including C. thirkei; Fig. 3B) yielded stable results (with a similarity coefficient between pairs of runs equal to 1) only for K ¼ 5. The results support the recognition of five genetic lineages, the composition of which largely contradicts recent taxonomic concepts: (1) C. napulifer and some accessions of C. orbelicus from Bulgaria (C. napulifer þ C. orbelicus BU; coloured pink in Fig. 3B); (2) C. orbelicus from Serbia and from Bulgaria, and C. velenovskyi from Greece (coloured blue in Fig. 3B); (3) C. tuberosus from Greece and the border region of Greece and the Former Yugoslav Republic of Macedonia, C. velenovskyi, C. nissanus from Serbia, and C. orbelicus from Mt Vichren (Bulgaria) (coloured yellow in Fig. 3B); (4) C. tuberosus from the central Former Yugoslav Republic of Macedonia (C. tuberosus MK; coloured orange in Fig. 3B); and (5) C. tuberosus from Croatia (C. tuberosus HR; coloured red in Fig. 3B). The last lineage was genetically most distinct from all the others. Obvious traces of genetic admixture were found in four populations (TRI 151, TRI 153, TRI 265 and TRI 273). Further NJ (not shown), NN (Fig. 3A) and PCoA analyses (Supplementary Data Fig. S2) based on AFLP dataset 2 showed deeper structuring within the second and third genetic lineages that correlated with geographical origin. Within the second
from 460 profiles and four selective primer combinations. Fragment length ranged from 51 to 488bp, and 436 (986%) of fragments were polymorphic. The overwhelming majority of profiles (453, 991%) were unique, and identical ones were fairly rare (Supplementary Data Table S3). The average (6 s.d.) number of fragments per population varied from 60648 to 91644 in diploids and from 75674 to 110616 in tetraploids (Table S3). All three genetic analyses, namely NJ (not shown), NN (Fig. 2) and PCoA (Supplementary Data Fig. S1), all based on AFLP dataset 1, revealed essentially similar patterns, with the ‘core’ genetic lineage encompassing solely diploids from the Balkan Peninsula, assignable to the C. napulifer group. This core lineage indicates that the C. napulifer group represents a genetically well-defined unit. Samples outside the core genetic lineage formed well-delimited and speciesspecific genetic lineages with high bootstrap support (Fig. 2, Fig. S1): C. epirotus (2x), C thirkei (2x) and C. pichleri (4x) from the Balkan Peninsula, and C. fuscomarginatus (4x) from Crimea. Exceptions were populations of C. graminifolius (4x) from south-eastern France and Calabria, which clustered together with populations of C. lingulatus (2x) from Spain and Morocco.
C. pindicola C. epirotus
57
100
C. lingulatus
99 89
99 65
99 71
66 66 99
77 99
C. graminifolius 100
98 83
100
99 100
C. fuscomarginatus
100
C. pichlerii
Cyanus napulifer group
C. thirkei
10·0
FIG. 2. Neighbour-net diagram based on AFLP profiles of 460 Cyanus individuals. The numbers indicate bootstrap probabilities >50 % obtained by neighbour-joining tree analysis of the same dataset. The symbols and colours are same as in Fig. 1; symbols with thick outlines are used for individuals from populations with presumably hybrid origin.
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group genetic lineage (coloured blue in Fig. 3B), samples of C. velenovskyi from Greece formed their own well-supported cluster (C. velenovskyi GR sublineage) separated from samples of C. orbelicus from Serbia and Bulgaria (C. orbelicus BU/RS sublineage). The most heterogeneous third lineage (coloured yellow in Fig. 3B) was split into several mostly well-supported sublineages: (1) C. tuberosus from Greece (C. tuberosus GR/MK sublineage), being genetically closest to the C. tuberosus MK lineage; (2) C. velenovskyi from Serbia (C. velenovskyi RS sublineage, TRI 258); (3) C. nissanus from Serbia (C. nissanus RS A
1079
sublineage, TRI 256); and (4) C. orbelicus from Mt Vichren (C. orbelicus Vich; TRI 151). Whereas C. velenovskyi RS and C. nissanus RS were resolved in close proximity of C. tuberosus GR/MK and C. tuberosus HR in all analyses; the position of C. orbelicus Vich was uncertain and varied across analyses. In summary, C. orbelicus Vich showed genetic affinities not only to C. orbelicus BU/RS and C. napulifer þ orbelicus BU (Fig. 3A), but also to populations of C. tuberosus GR/MK (Fig. S2). The results of Bayesian analysis based on the dataset including the C. napulifer group without the four genetically
C.orbelicus BU/RS
C. orbelicus Vich TRI 151
C. velenovskyi GR
TRI 153 C.napulifer + C. orbelicus BU
99
TRI 265 90 65
82
C. tuberosus MK
63
10·0
100 100
99
C. tuberosus HR
TRI 273 C. nissanus RS C. velenovskyi RS
C. tuberosus HR
TRI 265
C. tuberosus MK
TRI 273
C. tuberosus GR/MK
C. orbelicus Vich
C. nissanus RS
C. velenovskyi RS
C. velenovskyi GR
C. napulifer + orbelicus BU
B
C. orbelicus BU/RS
TRI 153
TRI 151
C. tuberosus GR/MK
FIG. 3. (A) Neighbour-net diagram based on AFLP profiles of 264 individuals of the Cyanus napulifer group. Numbers indicate bootstrap probabilities >50 % obtained by neighbour-joining tree analysis of the same dataset (Table S1). The symbols and colours are the same as in Fig. 1; symbols with thick outlines are used for individuals from populations with presumably hybrid origin. (B) Results of Bayesian clustering with K¼5 (BAPS software) of AFLP profiles for 264 individuals of the C. napulifer group. Each individual is represented by a vertical bar coloured according to its cluster assignment. Genetic lineages (coloured solid lines) and sublineages (coloured dashed lines) are indicated below the vertical bars.
1080
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
TABLE 2. Distribution of AFLP fragments across taxa of the genus Cyanus (A) and genetic (sub)lineages of the C. napulifer group (B) Taxon/genetic (sub)lineage
NIND
NPHEN
NFRAGM
(A) Taxa of the genus Cyanus (AFLP dataset 1) C. napulifer group 264 258 67 6 561 C. pindicola 34 34 69 6 74 C. epirotus 21 21 74 6 38 C. thirkei 40 40 76 6 458 C. pichleri 28 28 84 6 385 C. fuscomarginatus 16 15 99 6 328 C. lingulatus 45 45 81 6 712 C. graminifolius 12 12 104 6 71 (B) Genetic lineages of the C. napulifer group (AFLP dataset 3) C. napulifer þ orbelicus BU 45 44 71 6 451 C. orbelicus BU/RS 32 30 66 6 429 C. velenovskyi GR 24 23 65 6 401 C. velenovskyi RS 8 8 68 6 55 C. nissanus RS 8 8 64 6 463 C. orbelicus Vich – – – C. tuberosus GR/MK 52 52 64 6 49 C. tuberosus MK 40 40 66 6 63 C. tuberosus HR 24 24 72 6 436
P ( %)
NDG
NPR
NRARE
NUNI
7488 4263 2535 3341 3502 2396 3641 2097
0 0 0 1 0 1 1 2
41 4 7 6 3 6 9 1
287 86 54 70 63 39 58 25
9 0 0 1 0 0 2 0
– – – – – – – –
011 011 007 009 011 01 008 009
3935 3225 2751 1095 1716 – 4882 5355 2751
1 0 0 0 1 – 0 1 5
19 6 3 2 4 – 113 33 12
44 25 29 11 22 – 62 71 41
5 2 0 1 1 – 5 16 4
264 83 91 20 38 – 150 243 243
009 01 009 005 008 – 011 011 007
DW
Nei
NIND, number of individuals analysed; NPHEN, number of AFLP multilocus phenotypes; NFRAGM, average number of AFLP fragments per individual 6 standard deviation; P (%), percentage of polymorphic markers; NDG, number of private fixed (diagnostic) fragments; NPR, number of private (exclusive) fragments; NRARE, number of rare fragments (present at a frequency of
On the origins of Balkan endemics: the complex evolutionary history of the Cyanus napulifer group (Asteraceae) n ova1 and Jaromır Kucera1 Katarına Olsavska1,*, Marek Slovak1, Karol Marhold1,2, Eliska Stub 1
Institute of Botany, Slovak Academy of Sciences, D ubravsk a cesta 9, SK-84523 Bratislava, Slovakia and 2Department of Botany, Faculty of Science, Charles University, Ben atsk a 2, CZ-128 01 Praha 2, Czech Republic *For correspondence. E-mail [email protected] Received: 20 January 2016 Returned for revision: 2 March 2016 Accepted: 3 June 2016 Published electronically: 21 July 2016
Background and Aims The Balkan Peninsula is one of the most important centres of plant diversity in Europe. Here we aim to fill the gap in the current knowledge of the evolutionary processes and factors modelling this astonishing biological richness by applying multiple approaches to the Cyanus napulifer group. Methods To reconstruct the mode of diversification within the C. napulifer group and to uncover its relationships with potential relatives with x ¼ 10 from Europe and Northern Africa, we examined variation in genetic markers (amplified fragment length polymorphisms [AFLPs]; 460 individuals), relative DNA content (40 ,6-diamidino-2phenylindole [DAPI] flow cytometry, 330 individuals) and morphology (multivariate morphometrics, 40 morphological characters, 710 individuals). To elucidate its evolutionary history, we analysed chloroplast DNA (cpDNA) sequences of the genus Cyanus deposited in the GenBank database. Key Results The AFLPs revealed a suite of closely related entities with variable levels of differentiation. The C. napulifer group formed a genetically well-defined unit. Samples outside the group formed strongly diversified and mostly species-specific genetic lineages with no further geographical patterns, often characterized also by a different DNA content. AFLP analysis of the C. napulifer group revealed extensive radiation and split it into nine allopatric (sub)lineages with varying degrees of congruence among genetic, DNA-content and morphological patterns. Genetic admixture was usually detected in contact zones between genetic lineages. Plastid data indicated extensive maintenance of ancestral variation across Cyanus perennials. Conclusion The C. napulifer group is an example of a rapidly and recently diversified plant group whose genetic lineages have evolved in spatio-temporal isolation on the topographically complex Balkan Peninsula. Adaptive radiation, accompanied in some cases by long-term isolation and hybridization, has contributed to the formation of this species complex and its mosaic pattern. Key words: Cyanus napulifer group, Centaurea; Compositae, Balkan Peninsula, AFLP, cpDNA, flow cytometry, multivariate morphometrics, endemics, allopatric speciation, hybridization, incomplete lineage sorting.
INTRODUCTION The Balkan Peninsula is one of the most important European centres of plant diversity. Compared with any other region of Europe, it harbours the richest flora, including a high percentage of endemics (e.g. Eastwood, 2004; Krystufek and Reed, 2004; Stevanovic et al., 2007; Hewitt, 2011). Its astonishing diversity was likely generated by a combination of factors and processes. One of the most important is undoubtedly the refugial character of the Balkan Peninsula. The paleo-environmental conditions in this region allowed the survival of relic taxa, including some of tertiary origin (Tzedakis, 2004; Peev and Delcheva, 2007). Furthermore, the territory of today’s Balkan Peninsula served as a crossroad for Asian lineages and taxa during east–west colonization of Europe since the Early Oligocene (32 million years ago; Hewitt, 2011; Manafzadeh et al., 2014). The heterogeneous climate and complex physical geography of the Balkan Peninsula provided a wide range of conditions for different species and vegetation types to evolve in. The patchy, topographically varied landscape of the Balkans, with its mosaic of
mountain ranges, islands and peninsulas that acted as migration barriers for most organisms, significantly contributed to the formation of Balkan biodiversity (Stevanovic, 1996; Reed et al., 2004; Meshinev, 2007; Peev and Delcheva, 2007; Velikov and Stoyanova, 2007; Radford and Ode´, 2009). The Balkan flora is nevertheless still largely unexplored and calls for further investigation (but see Perny et al., 2005; Frajman and Oxelman, 2007; Stefanovic et al., 2008; Frajman and Schneeweiss, 2009; Scho¨nswetter and Schneeweiss, 2009; Bardy et al., 2010; Kolarcik et al., 2010; Kucera et al., 2010; Surina et al., 2011; Kuzmanovic et al., 2013; Lakusic et al., 2013; Ronikier and Zalewska-Gałosz 2014). The Balkan Peninsula accommodates almost 27% (18 out of 67 species) of perennial taxa of Cyanus sect. Protocyanus (Dobrocz.) Olsavska (Olsavska et al., 2013) and is considered its evolutionary centre (Hellwig, 2004). A large proportion of Cyanus perennials from the Balkans belong to the taxonomically extremely complex group of C. napulifer (Olsavska et al., 2013), which is the focus of the present study. The C. napulifer group, as defined by Olsavska et al. (2013), comprises five
C The Author 2016. Published by Oxford University Press on behalf of the Annals of Botany Company. V
All rights reserved. For Permissions, please email: [email protected]
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
1072
Balkan endemics – C. napulifer (Rochel) Sojak (s.str.), C. orbelicus (Velen.) Sojak, C. velenovskyi (Adamovic) Wagenitz & Greuter, C. nissanus (Petrovic) Sojak and C. tuberosus (Vis.) Sojak; one species occurring in the Balkan Peninsula and Turkey [C. thirkei (Sch. Bip.) Holub]; and one recently described Turkish species [C. eflanensis Kaya & Bancheva (Greuter, 2006–2009; Kaya and Bancheva, 2009)]. Although their mutual phylogenetic relationships remain unresolved, previous genetic studies indicate their close evolutionary relationships. However, C. thirkei seems to occupy a separate position, and no information is available on C. eflanensis (Borsic et al., 2011; Lo¨ser, 2012). Morphologically, the C. napulifer group is clearly characterized by a suite of common features, mainly by rhizomatous and/or tuberous roots and narrow stem leaves. All the taxa are diploids and share the base chromosome number of x ¼ 10 (no data are available for C. eflanensis), save for a few triploids (2n3x30) found in C. nissanus and C. thirkei (Olsavska et al., 2013). Two other Balkan species, C. diospolitanus Bancheva & Stoyanov and C. pseudoaxillaris (Stef. & T. Georgiev) Holub, previously included in Cyanus sect. Napuliferi (Stef. & T.Georgiev) Bancheva & Raimondo (Bancheva and Raimondo, 2003; Bancheva and Stoyanov, 2009), were proved to be genetically, karyologically and morphologically unrelated to the C. napulifer group (Borsic et al., 2011; Olsavska et al., 2013) and so were not included in the present study. The C. napulifer group represents a promising model for the study of homoploid speciation promoted mainly by geographical isolation and adaptation to new ecological niches. However, previous efforts to infer the evolutionary relationships of Cyanus perennials and the mechanisms of their diversification largely failed because of the insufficient or conflicting phylogenetic signals inferred from genetic sequence markers (Borsic et al., 2011; Lo¨ser, 2012). To uncover the evolutionary relationships among members of this group and the factors responsible for their diversification, we applied a combination of approaches that has proved to be effective in similar studies (e.g. Scho¨nswetter et al., 2009; Bardy et al., 2010; Kucera et al., 2010; Olsavska et al., 2011). Specifically, we employed fast-evolving multilocus amplified fragment length polymorphism (AFLP) markers, chloroplast DNA (cpDNA) data, flowcytometric data and multivariate morphometrics. The aims of our current study were: (1) to reveal which populations/taxa constitute the C. napulifer group and to uncover their relationships with other Cyanus taxa with x ¼ 10 occurring in Europe and Northern Africa; (2) to infer the overall genetic, karyological and morphological patterns of Balkan members of the C. napulifer group and to examine their congruence; and (3) to identify the evolutionary processes and factors that triggered speciation and divergence within the C. napulifer group and to infer biogeographical implications. MATERIALS AND METHODS Sampling design
The sampling strategy was focused on members of the Cyanus napulifer group occurring predominantly in the Balkan Peninsula. A total of 38 population samples of C. napulifer (s.str.), C. orbelicus, C. velenovskyi, C. nissanus, C. thirkei and
C. tuberosus were collected from natural habitats in 2009–13 (May/August) throughout the Balkan Peninsula [Bulgaria (BU), Croatia (HR), Greece (GR), Former Yugoslav Republic of Macedonia (MK), Serbia (RS)] (Fig. 1, Table 1). To elucidate the evolutionary relationships of the C. napulifer group, material from 30 populations of other morphologically similar European and North African species of the section Protocyanus, with the base chromosome number of x ¼ 10 , was collected for AFLP and DNA content analyses (Table 1, Fig. 1): C. pindicola (Griseb.) Sojak, C. epirotus (Halacsy) Holub, C. lingulatus (Lag.) Holub, C. graminifolius (Lam.) Olsavska, C. fuscomarginatus (K.Koch) Greuter, and C. pichleri (Olsavska et al., 2013). Plant material
Intact fresh leaves were collected from three to eight plants per locality and dried in silica gel for AFLP analyses. Roots/rhizomes from up to 14 plants were transferred to an experimental garden in Banska Bystrica, Slovakia (48 450 08900 N, 19 090 29000 E, 390m a.s.l.) for flow cytometry estimation. Depending on the abundance of flowering individuals, 10–52 plants of the C. napulifer group (not including C. thirkei) were collected at original localities for morphometric analyses. The type of bedrock and biotope were recorded during the field research. Each population of the species under investigation had at least several hundreds of individuals, and the taxa under study are neither endangered nor protected by law. Voucher specimens were deposited in the herbarium of the Institute of Botany of the Slovak Academy of Sciences (SAV). DNA extraction and AFLP generation
Total genomic DNA was extracted using the DNeasy Plant Mini Kit (Qiagen, Du¨sseldorf, Germany). The AFLP procedure (Vos et al., 1995) followed a modified protocol described by Mered’a et al. (2008). To obtain highquality AFLP profiles, we tested 38 selective primer combinations using four individuals from different, geographically distant taxa. Final analyses were performed using four primer combinations labelled with different fluorescent dyes and producing clear, unambiguous and polymorphic profiles: EcoRIACA-(6-FAM)/MseI-CTG, EcoRI-AAG-(VIC)/MseI-CAG, EcoRI-AAC-(NED)/MseI-CTG and EcoRI-ACA-(PET)/MseICAC (Applied Biosystems). Amplification products were pooled and fragment analysis was performed on an ABI 3100 Avant capillary sequencer of the BITCET Consortium, Department of Molecular Biology, Comenius University, Bratislava. The internal size standard GeneScan –500 LIZ (Applied Biosystems) was employed for size calibration of the AFLP profiles obtained. AFLP trace files were read and analysed using DAX software (Van Mierlo Software Consultancy, Netherlands). Only unambiguously scorable markers in the size range of 50–500 base pairs (bp) were scored manually as present (1) or absent (0). The reproducibility of AFLP profiles was determined using 35 replicates (76% of the final dataset). This dataset comprised randomly chosen re-extracted samples of all taxa under study, which were analysed and scored independently. The error rate (Bonin
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group Cyans napulifer group C. napulifer + orbelicus BU
C. orbelicus BU/RS C. velenovskyi GR C. velenovskyi RS C. nissanus RS C. orbelicus Vich C. tuberosus GR/MK C. tuberosus MK C. tuberosus HR
1073
other Cyanus taxa: C. pindicola C. epirotus C. thirkei C. pichleri C. fuscomarginatus C. lingulatus C. graminifolius
AL - Albania, BA - Bosnia and Herzegovina, BG - Bulgaria, HR - Croatia, GR - Greece, ME - Montenegro, MK - Former Yugoslav Republic of Macedonia, RS - Serbia
FIG. 1. Distribution map of Cyanus samples used in molecular, DNA-content and morphometric analyses. Samples of C. napulifer group are assigned to genetic lineages (coloured solid lines) and sublineages (coloured dashed lines) according to the results of AFLP data analyses. Symbols with thick lines are used for populations of presumably hybrid origin.
et al., 2004) was calculated as the ratio of mismatches (1 versus 0) to matches (1 versus 1) between replicated samples.
AFLP data analyses
To reveal the genetic pattern of the taxa under study, we assembled four AFLP datasets. The first dataset, AFLP dataset 1, comprised all analysed taxa (460 individuals); the second one, AFLP-dataset 2 (264 individuals), comprised only accessions of the C. napulifer group excluding C. thirkei, which early on turned out to be genetically clearly distinct from the rest of the group (see the Results section). AFLP dataset 3 (233 individuals) comprised accessions of the C. napulifer group excluding four genetically admixed populations (TRI 151, TRI 153, TRI 265 and TRI 273; see the Results section). For the purposes of multispecies coalescent analyses, we created AFLP dataset 4 (120 individuals), composed of two randomly selected individuals for each population analysed except four genetically admixed ones (TRI 151, TRI 153, TRI 265 and TRI 273; see also below in this section). The overall genetic pattern of relationships among populations and individuals was investigated using the following approaches based on AFLP datasets 1–3: (1) neighbour-joining tree (NJ; Saitou and Nei, 1987) based on Sørensen’s ( ¼ Dice) similarities transformed into a distance matrix with bootstrap
support for clades and subclades computed using 5000 replicates; (2) split network (neighbour-net, NN) analysis based on Sørensen’s similarities transformed into a distance matrix (Huson and Bryant, 2006); and (3) principal coordinates analysis (PCoA; Krzanowski, 1990) based on Jaccard’s similarity matrix. The NJ tree analysis and PCoA were performed using Famd 1108 beta software (Schlu¨ter and Harris, 2006), while NN was conducted using SplitsTree 4 (Huson and Bryant, 2006). To search for genetically homogeneous groups within the C. napulifer group encompassing only diploid taxa (AFLP dataset 2), Bayesian model-based clustering as implemented in Structure 223 (Falush et al., 2007) was employed. The settings for the analysis were the same as those used in Slovak et al. (2012a, b). The R script Structure-sum-2009 (part of AFLPdat; Ehrich, 2006) was used to calculate the statistics necessary for determining the correct K value (Evanno et al., 2005). Graphical outputs of analyses were generated using Clumpp v11.1 (Jakobsson and Rosenberg, 2007) and Distruct (Rosenberg, 2004). To infer hierarchical relationships between genetic lineages of the C. napulifer group and related species with x ¼ 10, we computed a species trees in a multispecies coalescent framework using SNAPP v2.2.0 (Bryant et al., 2012). A priori ‘species’ delimitation was done with respect to the results of NN, NJ, PCoA and structure analyses, with the exceptions of
1074
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
TABLE 1. Locality details, including geographical coordinates, altitude, the date of collection and the collector(s) of plants of Cyanus species included in the molecular (AFLP), DNA content flow cytometric and morphometric analyses Population (voucher) code*
Locality and collection details†
Cyanus napulifer group Cyanus napulifer þ orbelicus BU (central Bulgarian) genetic lineage TRI 214 Bulgaria; Stara Planina Mts, Korduna peak; 42 450 1900 N, 24 040 3600 E; 1605 m; (A) 2.viii.2011; (B) 8.vii.2013, MP & KO TRI 215 Bulgaria; Stara Planina Mts, 2 km SE from Vezhen cottage; 42 450 47500 N, 24 240 10800 E; 1869 m; 2.viii.2011, KO TRI 216 Bulgaria; Stara Planina Mts, Troyanski prokhod, Beclemento; 42 460 3700 N, 24 360 5200 E; 1601 m; 3.viii.2011, KO TRI 217 Bulgaria; Stara Planina Mts, Troyanski prokhod, Mt Viloto; 42 460 4700 N, 24 390 5600 E; 1603 m; 3.viii.2011, KO TRI 300 Bulgaria; Stara Planina Mts, Buzludzha; 42 440 07000 N, 25 230 51800 E; 1400 m; 9.vii.2013, MP & KO TRI 222 Bulgaria; Rila Mts, Belmeken dam; 42 090 5100 N, 23 460 5700 E; 1963 m; 9.viii.2011, KO TRI 153 Bulgaria; Rila Mts; (A) 25 km NE from Granchar Hut; 42 080 0800 N, 23 360 5400 E; 2294 m; 19.vi.2010, KO; (B) 12 km N from Granchar Hut; 42 070 5100 N, 23 350 4500 E; 2277 m; 20.vi.2010, KO Cyanus orbelicus BU/RS (Bulgarian–Serbian) genetic sublineage TRI 296 Serbia; Kopaonik Mts Karaman; 43 170 34000 N, 20 490 47700 E; 1910; 6.vii.2013, MP & KO TRI 260 Serbia; Mt Besna Kobila; (A) below the peak; 42 320 57900 N, 22 130 56200 E; 1662 m; (B) N of the peak; 42 330 29900 N, 22 140 09400 E; 1559 m; 15.vi.2012, CL & KO TRI 298 Bulgaria; Ossogovska planina; (A) Kulin Kamak; 42 110 41400 N, 22 350 54300 E; (B) Gramadite; 42 110 17900 N, 22 360 31600 E; 1636 m; 1710 m; 7.vii.2013, MP & KO TRI 302 Bulgaria; Rhodopi Mts, Golyam Perelik Mt; 41 360 00800 N, 24 340 48500 E; 2140 m; 11.vii.2013, MP & KO Cyanus velenovskyi GR (Greek) genetic sublineage TRI 144 Greece; Pella, E slope of Mt Kaimaktsalan; 40 550 0300 N, 21 480 4100 E; 2160 m; 13.vi.2010, MO & KO TRI 241 Greece; Florina, Vernon Mts, W slope of Mt Vitsi; 40 380 53000 N, 21 220 47100 E; 1921 m; (A) 22.v.2012, MO & KO; (B) 13.vii.2013, MP & KO TRI 305 Greece; Florina, Varnous Mts, Mt Vigla; 40 480 35300 N, 21 150 07000 E; 1990 m; 13.vii.2013, MP & KO Cyanus velenovskyi RS (Serbian) genetic sublineage TRI 258 Serbia; Stara Planina Mts, Mt Midzor; (A) near the ski resort Babin zub; 43 220 45600 N, 22 370 3900 E; 1592 m; 15.vi.2012, KO; (B) S slope; 43 220 4400 N, 22 400 4800 E; 1662 m; 15.vi.2012, CL Cyanus nissanus RS genetic sublineage TRI 256 Serbia; Svrljiska Planina Mts, Mt Popova glava, summit area; 43 210 23400 N, 22 050 30200 E; 1008 m; (A) 14.vi.2012, CL & KO; (B) 19.v.2013, MO, JO & KO Cyanus orbelicus Vich (Vichren) genetic sublineage TRI 151 Bulgaria; Pirin Mts, Mt Vihren, above Vihren cottage; 41 450 2700 N, 23 240 1700 E; 2440 m; (A) 17.vi.2010; (B) 10.viii.2010; (C) Bakushevata Mura; 41 450 5400 N, 23 240 4400 E; 2217 m; 10.viii.2011, KO Cyanus tuberosus GR/MK (Greek–Macedonian) genetic sublineage TRI 240 Greece; Larisa, Deskati, E of Loutro; 39 560 55400 N, 21 550 40800 E; 732 m; 21.v.2012, MO & KO TRI 170 Greece; Grevena, Vourinos Mts, S from Exarchos, 40 090 0200 N, 21 370 3000 E, 684 m; 28.v.2011, KO TRI 149 Greece; Kozani, Pieria Mts, SE of Velvendo; 40 130 5200 N, 22 060 3900 E; 1278 m; 16.vi.2010, KO TRI 162 Greece; Kilkis, Paiko Mts; (A) S of Livadhia; 40 590 4000 N, 22 170 4600 E; 1160 m; (B) NW of Livadhia; 41 010 N, 22 160 3000 E; 1180 m; 26.v.2011, KO TRI 238 Greece; Thesaloniki, Mt Hortiatis; 40 350 57800 N, 23 060 01200 E; 826 m; 20.v.2012, MO & KO TRI 237 Greece; Serres; (A) Lailas; 41 140 22700 N, 23 340 22000 E; 1393 m; (B) between Kato Vrondous and Eleonas; 41 130 57200 N, 23 390 45300 E; 915 m; 19.v.2012, MO & KO TRI 235 Greece; Drama, Voras Mts; (A) Falakro; 41 170 36600 N, 24 020 29700 E; 1467 m; (B) Volokas; 41 180 10900 N, 23 590 34900 E; 1012 m; 19.v.2012, MO & KO TRI 273 Macedonia/Greece; Kozuf Mts; (A) Mt Keci Kaj; 41 110 48000 N, 22 140 33800 E; 1715 m; 22.vi.2012, KO; (B) above ski centre; 41 110 01100 N, 22 120 59700 E; 1912 m; 22.vi.2012, CL; (C) Mt Mala Rupa; 41 090 52500 N, 22 140 37600 E; 1863 m; 22.vi.2012, CL Cyanus tuberosus MK (central Macedonian) genetic lineage
No. of plants analysed AFLP
Flow cytometry N/P‡ (ploidy; ST§)
Morphometric analyses
8
4/3 (2n 20; SL)
20
5
3/3 (2n 20; SL)
15
8
3/3 (2n 20; SL)
12
8
4/3 (2n 20; SL)
11
8
6/0 (2n 20; SL)
20
8 8
2/4 (2n 20; SL) 2/6 (2n 20; SL)
15 31
8 8
5/0 (2n 20; SL) 0/4 (2n 20; SL)
20 41
8
8/0 (2n 20; SL)
32
8
5/0 (2n 20; SL)
20
8
0/4 (2n 20; SL)
20
8
5/4 (2n 20; SL)
12
8
5/0 (2n 20; SL)
20
8
3/4 (2n 20; SL)
29
8
1/2 (2n 20; SL) 0/2 (2n 30; SL)
20
7
3/7 (2n 20; SL)
24
8
2/0 (2n 20; SL)
15
5
0/4 (2n 20; SL)
13
8
0/4 (2n 20; SL)
18
7
11/0 (2n 20; SL)
35
8
2/0 (2n 20; SL)
14
8
5/0 (2n 20; SL)
27
8
1/0 (2n 20; SL)
35
8
3/0 (2n 20; SL)
52
(continued)
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
1075
TABLE 1. Continued Population (voucher) code*
Locality and collection details†
Macedonia; Dren Mts, pass above Stavica; 41 150 03500 N, 21 350 41100 E; 1132 m; (A) 21.vi.2012, CL & KO; (B) 23.v.2013, MO, JO & KO TRI 271 Macedonia; Dren Mts; (A) Veprchani; 41 160 11900 N, 21 450 24400 E; 832 m; 22.vi.2012, KO; (B) 20.v.2013, MO, JO & KO; (B) above Dunje; 41 270 27200 N, 21 460 07400 E; 1479 m; 22.vi.2012, CL TRI 262 Macedonia; Ramno Mts, above Brodec; (A) 42 080 35800 N, 21 260 57500 E; 1100 m; 15.vi.2012, CL & KO; (C) 42 090 50900 N, 21 260 53000 E; 1386 m; 25.v.2013, MO, JO & KO TRI 264 Macedonia; Vodno Mts, Skopje, near Millenium cross; (A) 42 580 11700 N, 21 230 54600 E; 922 m; 16.vi.2012, CL; (B) 41 570 583400 N, 21 230 433400 E; 1022 m; 24.v.2013, MO, JO & KO TRI 266 Macedonia; Babuna Mts, Prisad pass; 41 260 54000 N, 21 370 50200 E 41 260 2400 N, 21 360 5500 E; 809–1100 m; (A) 19.vi.2012, CL & KO; (B) 24.v.2013, MO, JO & KO TRI 265 Macedonia; Babuna Mts, above Nezhilovo; 41 400 11800 N, 21 290 2900 E; 2064 m; 17.vi.2012, CL Cyanus tuberosus HR (Croatian) genetic lineage TRI 247 Croatia; Lika, Mt Postak, along a road to Rasticevo; 44 130 28400 N, 16 050 48800 E; 752 m; 5.vi.2012, MO, JO & KO TRI 252 Croatia; Split-Dalmatia County, above Sinj; (A) near Zelovo; 43 450 04400 N, 16 330 12200 E; 800 m; 7.vi.2012, MO, JO & KO; (B) along a road to Sutina; 43 420 04000 N, 16 360 52500 E; 402 m; 7.vi.2012, MO, JO & KO TRI 253 Croatia; Mosor Mts, near Rasca; 43 330 48800 N, 16 420 33800 E 43 330 22300 N, 16 420 22500 E; 312–361 m; 8.vi.2012; MO, JO & KO Cyanus pindicola (Griseb.) Sojak TRI 145 Greece; Pella, NE slope of Mt Piperitsa; 40 510 3800 N, 21 440 3200 E; 1812 m; 14.vi.2010, MO & KO TRI 148 Greece; Pieria, Oros Olympos Mts, below Zolotas cottage; 40 040 5600 N, 22 220 5000 E; 1709 m; 15.vi.2010, MO & KO TRI 166 Vermio Mts, Seli 40 320 0300 N, 22 010 3000 E; 1662 m; 27.v.2011, KO TRI 168 Vourinos Mts 40 110 56200 N, 21 390 27700 E; 1330 m; 28.v.2011, KO TRI 175 Greece; Larisa, Kato Olympos Mts, 2 km NE from Kallipefki; 39 580 2500 N, 22 280 4600 E; 1141 m; 29.v.2011, CL & KO TRI 267 Greece; Galicica Mts; (A) Galicica Mt; 40 570 43200 N, 20 490 41300 E; 1843 m; 20.vi.2012, KO; (B) Magaro Mt; 40 560 22100 N, 20 490 18700 E; 2 163 m; 20.vi.2012, CL Cyanus epirotus (Halacsy) Holub TRI 182 Greece; Ioannina, Tzoumerka Mts, Kryopigi; 39 270 0900 N, 21 070 3900 E; 1443 m; 2.vi.2011, CL & KO TRI 183 Greece; Ioannina, Tzoumerka Mts, close to Mparos pass; 39 360 36600 N, 21 080 32900 E; 1676 m; 2.vi.2011, CL & KO TRI 185 Greece; Ioannina, Mitsikeli Mts, above Ioannina; 39 430 01600 N, 20 530 0700 E; 1669 m; 3.vi.2011, CL & KO TRI 187 Greece; Ioannina, Timfi Mts, above the road to Vradeto; 39 550 0300 N, 20 480 2100 E; 1738 m; 3.vi.2011, CL & KO Cyanus thirkei (Sch.Bip.) Holub TRI 228 Bulgaria; Stara Planina Mts, Varbishki prohod, above Sadovo, 42 530 17600 N, 26 390 05800 E; 686 m; (A) 15.v.2012, MO & KO; (B) 20.v.2013, MO, JO & KO TRI 229 Bulgaria; Balgarova near Burgas; 42 370 58300 N, 27 160 48200 E; 99 m; (A) 16.v.2012, KO; (B) 21.v.2013, MO, JO & KO TRI 230 Bulgaria; Otmanli near Burgas; 42 250 49100 N, 27 320 59100 E; 52 m; 16.v.2012, MO & KO TRI 294 Bulgaria; Karnobatska Planina Mts, between Valchin and Lazarevo; 42 440 26600 N, 26 530 54700 E; 271 m; 21.v.2013, MO, JO & KO TRI 295 Bulgaria; Karnobatska Planina Mts, N of Aytos; 42 450 04600 N, 27 150 16700 E; 347 m; 21.v.2013, MO, JO & KO Cyanus pichleri TRI 179 Greece; Peloponnese, Mt Panachiko; 38 150 5200 N, 21 510 0600 E; 1146 m; 31.v.2011, CL & KO TRI 186 Greece; Ioannina, Timfi Mts, Vradeto; 39 530 53200 N, 20 460 23800 E; 1339 m; 3.vi.2011, CL & KO TRI 188 Greece; Kastoria, between Nea Kotyli and Kato Nestorio; 40 220 16600 N, 21 020 41400 E; 1268 m; 4.vi.2011, CL & KO TRI 231 Bulgaria; Strandzha Mts, Golyamo Bukovo, above road E79; 42 140 28300 N, 27 060 36600 E; 340 m; 17.v.2012, MO & KO TRI 239 Greece; Mt Chortiatis, above village of Chortiatis; 40 350 57800 N, 23 060 01200 E; 826 m; 20.v.2012, MO & KO Cyanus fuscomarginatus (K.Koch) Greuter TRI 270
No. of plants analysed AFLP
Flow cytometry N/P‡ (ploidy; ST§)
Morphometric analyses
8
4/2 (2n 20; SL)
20
8
5/4 (2n 20; SL)
10
8
8/4 (2n 20; SL)
19
8
5/1 (2n 20; SL)
20
8
5/3 (2n 20; SL)
20
8
1/4 (2n 20; SL)
16
8
2/4 (2n 20; SL)
17
8
6/8 (2n 20; SL)
11
8
3/4 (2n 20; SL)
8
6
0/3 (2n 20*; SL)
–
6
0/1 (2n 40*; SL)
–
5 5 6
3/0 (2n 40; SL) 3/0 (2n 40; SL) 0/5 (2n 40*; SL)
– – –
6
1/0 (2n 20; SL)
–
3
0/3 (2n 20; SL)
–
6
0/3 (2n 20; SL)
–
6
0/5 (2n 20; SL)
–
6
0/3 (2n 20; SL)
–
8
–
8
3/0 (2n 20; BP) 0/7 (2n 30; BP) 3/4 (2n 20; BP)
–
8 8
0/4 (2n 20; BP) 6/0 (2n 20; BP)
– –
8
6/0 (2n 20; BP)
–
6
0/2 (2n 40; SL)
–
5
0/3 (2n 40; SL)
–
5
0/2 (2n 40; SL)
–
6
0/1 (2n 40; SL)
–
6
0/4 (2n 40; SL)
–
(continued)
1076
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group TABLE 1. Continued
Population (voucher) code*
Locality and collection details†
No. of plants analysed AFLP
Ukraine; Crimea, Aj-Petrinskaja Jajla, Mt Aj-Petri; 44 270 2900 N, 34 030 3700 E; 1178 m; 7.vii.2012, KO TRI 282 Ukraine; Crimea, Catyr-Dag Jalja, near Marmornaja cave; 44 470 2500 N, 34 160 4800 E; 973 m; 9.vii02012, KO Cyanus lingulatus (Lag.) Holub JAV Spain; Sierra de Javalambre; 40 130 31 N, 00 200 56 W; 1853 m; 3.vi.2010; SS, KM & JZ PUR Spain; Sierra del Moncayo, Purujosa; 41 420 47 N, 1 440 54 W; 1280 m; 8.v.2011; SS, KM & JZ ˆ in-Leuch; 33 124210 N, TRI 289 Morocco; Moyen Atlas, between Sources de l0 Oum-er-Rbia and A 05 206370 W; 1734 m; 1052013; AG, JK, MS & KO TRI 291 Morocco; Rif Mts, Tleta Ketama; 34 524060 N, 04 336100 W; 1573 m; 11.v.2013; AG, JK, MS & KO Cyanus graminifolius (Lam.) Olsavska TRI 65 France; Hautes Alpes, Montagne de Ce´u¨se; 44 310 0200 N, 05 520 0800 E; 1572 m; 15.vii.2007, KO, MP & DD TRI 101 France; Alpes de Hte Provence, Montagne du Cheiron, Pic de Fourneuby; 43 490 29400 N, 06 530 32500 E; 1422 m; 6.vii.2008; KO, MP & VK TRI 196 Italy; Pollino Mts, Colle del Dragone; 39 540 46400 N, 16 070 30800 E; 1531 m; 14.vi.2011, KO (SAV: TRI 196) TRI 198 Italy; Calabria, La Sila Mts, 4 km SW from Moccone; 39 200 05700 N, 16 230 56500 E; 1613 m; 14.vi.2011, KO (SAV: TRI 198) TRI 278
Flow cytometry N/P‡ (ploidy; ST§)
Morphometric analyses
8
0/4 (2n 40; SL)
–
8
0/3 (2n 40; SL)
–
8 7
0/5 (2n 20; BP) 0/5 (2n 20*; BP)
– –
8
6/0 (2n 20; BP)
–
8
1/0 (2n 20; BP)
–
8
0/5 (2n 40*; BP)
–
6
0/7 (2n 40; BP)
–
6
0/3 (2n 40; BP)
–
6
0/5 (2n 40; BP)
–
Populations of C. napulifer group are assigned to genetic lineages and sublineages according to the results of AFLP data analyses. *Population codes of populations with presumably hybrid origin are in bold. Collectors: AG, A. Guttova; IH, I. Hodalova; VK, V. Kolarcik; CL, C. Lo¨ser; KM, K. Marhold; JK, J. Kucera; JO, J. Olsavska; MO, M. Olsavsk y; MS, M. Slov a k; KO, K. Ol s avsk a ; MP, M. Pern y ; SS, S. Spaniel; JZ, J. Zozomov a -Lihov a . ‡ Number of newly analysed plants/number of plants analysed in a previous study (Olsavska et al., 2013) § Standard used in flow cytometry analysed: BP, Bellis perennis; SL, Solanum lycopersicum. †
C. lingulatus and C. graminifolius. These were regarded as a single species to ensure the monophyly of the group. Due to the large computational demands of the algorithm, we analysed a subsample comprising two randomly selected individuals for each population (AFLP dataset 4). We ran two analyses with different h priors in order to allow for different current and ancestral population sizes, as applied in Kolar et al. (2016): (1) mean h prior of 0043 (c distribution, a ¼ 15, b ¼ 35) and (2) mean h prior ¼ 01 (c distribution, a ¼ 12, b ¼ 110 prior for large population sizes); the remaining parameters were left at their default values. The analyses were checked for convergence using Tracer v1.6, making sure that Bayesian runs reached an effective sample size exceeding 200 after burn-in. The posterior distributions of species trees were visualized using DensiTree v2.2.0 (Bouckaert and Heled, 2014). Genetic diversity within and among the taxa/populations/genetic (sub)lineages under study (AFLP datasets 1, 2 and 4) was estimated using the following statistical parameters: (1) total number of AFLP multilocus genotypes; (2) average number of AFLP fragments (6 s.d.); (3) average proportion of pairwise differences between individuals (Nei’s gene diversity, DNei; Nei and Li, 1979) using the R script AFLPdat (Ehrich, 2006); and (4) percentage of polymorphic markers (P%) using Pop-Gene 132 (Yeh et al., 1997). Genetic divergence among analysed taxa/populations/genetic (sub)lineages was explored by calculating the number of (5) private fixed fragments (diagnostic); (6) private fragments (exclusive); (7) rare fragments (present at a frequency of < 10% of investigated individuals); and (9)
unique fragments, using FAMD 1108 beta software (Schlu¨ter and Harris, 2006) and by calculating of the frequency downweighted marker value (Scho¨nswetter et al., 2005), as implemented in AFLPdat (DW1; Ehrich, 2006). Within the C. napulifer group (AFLP dataset 3), we also explored within- and between-group variance using analyses of molecular variance (AMOVA). AMOVAs were performed using Arlequin 311 (Excoffier et al., 2005) based on Euclidean pairwise distances and a significance test with 10 000 permutations.
cpDNA data analyses
To address the evolutionary relationships among perennial representatives of the genus Cyanus, we explored cpDNA sequences deposited in the GeneBank database (NCBI, http:// www.ncbi.nlm.nih.gov/nucleotide), particularly those generated by C. J. Lo¨ser, G. Akaydin and F. H. Hellwig (University of Jena, Germany). We used 114 sequences of trnC-ycf6 and 99 sequences of 30 rps16-50 trnK(UUU) cpDNA regions (Shaw et al., 2005, 2007) for 15 Cyanus taxa. This dataset included not only an overwhelming majority of species analysed in our study with x ¼ 10, but also species with x ¼ 11 (Supplementary Data Table S1). All cpDNA sequences were manually edited and aligned using BIOEDIT (v7.0.4.1; Hall, 1999). Ambiguously aligned regions were removed with Gblocks (Castresana, 2000; Talavera and Castresana, 2007). The method of statistical parsimony was
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group employed to generate a haplotype network (TCS v.1.18; Clement et al., 2000) depicting the general pattern of relationships among the haplotypes. TCS was run with a 95% connection limit and with indels coded as missing data.
Flow cytometric analyses
The dataset resulting from relative nuclear DNA content analyses comprised both newly generated (136 individuals) and previously published data (194 individuals; Olsavska et al., 2013); between one and 14 plants per population were analysed by flow cytometry using the fluorochrome 40 ,6-diamidino-2phenylindole (DAPI, specific for A–T base pairs) (Table 1). DAPI flow cytometry was chosen because a previous study, which was focused on the C. triumfetti and C. montanus groups (Olsavska et al., 2012), proved a strong correlation between absolute DNA content (estimated using the intercalating fluorochrome propidium iodide) and relative DNA content (estimated using DAPI) (r ¼ 099, P5% in genome size were performed to confirm the reliability of the estimated values. The relationship between chromosome numbers and DNA content values was verified using previously published chromosome counts (Table 1; Olsavska et al., 2013). Relative DNA content was calculated as the ratio of the position of the G1 peak of the standard [Solanum lycopersicum ‘Stupicke´ polnı rane´’, 2C DNA ¼ 196pg (Dolezel et al., 1992) or Bellis perennis, 2C DNA ¼ 338pg (Scho¨nswetter et al., 2007)] to that of the sample. When B. perennis was used as the internal standard, the resulting relative DNA content was recalculated relative to the S. lycopersicum standard based on the ratio between B. perennis and S. lycopersicum (0578) calculated from eight estimations performed at different times. To compare samples differing in ploidy level, relative DNA content per monoploid genome was used. The Tukey–Kramer test (Tukey’s test for unequal sample size; Tukey, 1977) was used to test for differences in relative DNA content between populations of the C. napulifer group (triploids from population 256, investigated only in the previous study, were not included) and in monoploid relative DNA content between taxa or genetic (sub)lineages. All analyses were carried out using STATISTICA 91 (StatSoft Inc., 2013). For populations of C. nissanus and C. thirkei (48 plants in total), ploidy levels were estimated by flow cytometry for each individual subjected to AFLP because in a previous study (Olsavska et al., 2013), two ploidy levels (2x and 3x) were identified in some populations of these taxa. Sample preparation followed the procedure described above, but silica-dried material was used, and one to three plants were analysed simultaneously; such flow cytometric estimations were not used to calculate relative DNA content.
1077
Morphometric analyses
Morphological variability was studied using multivariate morphometric analyses based on 710 individuals from 34 populations of the C. napulifer group (Table 1). A total of 24 quantitative and nine qualitative characters were measured or scored on herbarium material, and three ratios were subsequently computed (Supplementary Data Table S2). Depending on the number of states recorded (Table S2), multistate qualitative characters were binary-coded for further analyses as two (ROOT, SI, LBWS, LBAS, LLI, LUI, BUL, BBL), three (FECO) or five (SECR) dummy binary variables. Each multistate qualitative character with s states was decomposed into s– 1 dummy variables (to avoid complete linear dependency of the sth variable). First, all quantitative characters were tested for normality of distribution (Shapiro–Wilk test). Subsequently, Spearman’s (non-parametric) correlation coefficient was computed for all characters because several characters departed from a normal distribution. No pair of characters showed a strong correlation (arbitrary value >095) that would potentially distort further computations. To gain insight into the overall phenetic relationships and to estimate the level of congruence among morphological, genome size and genetic (AFLP) variation, non-metric multidimensional scaling (NMDS, based on Euclidean distance and standardization by variable range; Kruskal, 1964) and canonical discriminant analysis (CDA; Krzanowski, 1990) were performed. Population samples used in the NMDS analyses were characterized by mean values of all 40 characters. The CDA analyses were based on individual plants and 26 quantitative and ratio characters. Qualitative characters were excluded because they were uniform in some groups of individuals. The analyses were based both on the whole set of data (NMDS 1) and on the partial datasets representing genetic (sub)lineages as revealed by Bayesian clustering (NMDS 2–6, CDA 1–5). The analyses were carried out using STATISTICA 91 (StatSoft Inc., 2013), Canoco 50 (ter Braak and Smilauer 2012) and SAS 93 (SAS Institute Inc., 2011).
Ecological and phenological trait analyses
To reveal the life strategies of populations of the C. napulifer group, the following ecological and phenological data were collected for each population: altitude, type of bedrock, type of biotope and flowering time (Table S2). These data, together with selected morphological characters describing the root system (Table S2), were subjected to NMDS analyses (NMDS 7) in the same way as morphological characters alone. Nei’s gene diversity, frequency down-weighted marker values and mean DNA content of particular populations were entered as supplementary characters. RESULTS AFLP analyses
The overall reproducibility of AFLP dataset 1 (entire dataset) was high because the error rate was minute, being only 019%. The analysis yielded 442 unambiguously scorable fragments
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
1078
Bayesian clustering based on AFLP dataset 2 (the C. napulifer group not including C. thirkei; Fig. 3B) yielded stable results (with a similarity coefficient between pairs of runs equal to 1) only for K ¼ 5. The results support the recognition of five genetic lineages, the composition of which largely contradicts recent taxonomic concepts: (1) C. napulifer and some accessions of C. orbelicus from Bulgaria (C. napulifer þ C. orbelicus BU; coloured pink in Fig. 3B); (2) C. orbelicus from Serbia and from Bulgaria, and C. velenovskyi from Greece (coloured blue in Fig. 3B); (3) C. tuberosus from Greece and the border region of Greece and the Former Yugoslav Republic of Macedonia, C. velenovskyi, C. nissanus from Serbia, and C. orbelicus from Mt Vichren (Bulgaria) (coloured yellow in Fig. 3B); (4) C. tuberosus from the central Former Yugoslav Republic of Macedonia (C. tuberosus MK; coloured orange in Fig. 3B); and (5) C. tuberosus from Croatia (C. tuberosus HR; coloured red in Fig. 3B). The last lineage was genetically most distinct from all the others. Obvious traces of genetic admixture were found in four populations (TRI 151, TRI 153, TRI 265 and TRI 273). Further NJ (not shown), NN (Fig. 3A) and PCoA analyses (Supplementary Data Fig. S2) based on AFLP dataset 2 showed deeper structuring within the second and third genetic lineages that correlated with geographical origin. Within the second
from 460 profiles and four selective primer combinations. Fragment length ranged from 51 to 488bp, and 436 (986%) of fragments were polymorphic. The overwhelming majority of profiles (453, 991%) were unique, and identical ones were fairly rare (Supplementary Data Table S3). The average (6 s.d.) number of fragments per population varied from 60648 to 91644 in diploids and from 75674 to 110616 in tetraploids (Table S3). All three genetic analyses, namely NJ (not shown), NN (Fig. 2) and PCoA (Supplementary Data Fig. S1), all based on AFLP dataset 1, revealed essentially similar patterns, with the ‘core’ genetic lineage encompassing solely diploids from the Balkan Peninsula, assignable to the C. napulifer group. This core lineage indicates that the C. napulifer group represents a genetically well-defined unit. Samples outside the core genetic lineage formed well-delimited and speciesspecific genetic lineages with high bootstrap support (Fig. 2, Fig. S1): C. epirotus (2x), C thirkei (2x) and C. pichleri (4x) from the Balkan Peninsula, and C. fuscomarginatus (4x) from Crimea. Exceptions were populations of C. graminifolius (4x) from south-eastern France and Calabria, which clustered together with populations of C. lingulatus (2x) from Spain and Morocco.
C. pindicola C. epirotus
57
100
C. lingulatus
99 89
99 65
99 71
66 66 99
77 99
C. graminifolius 100
98 83
100
99 100
C. fuscomarginatus
100
C. pichlerii
Cyanus napulifer group
C. thirkei
10·0
FIG. 2. Neighbour-net diagram based on AFLP profiles of 460 Cyanus individuals. The numbers indicate bootstrap probabilities >50 % obtained by neighbour-joining tree analysis of the same dataset. The symbols and colours are same as in Fig. 1; symbols with thick outlines are used for individuals from populations with presumably hybrid origin.
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group genetic lineage (coloured blue in Fig. 3B), samples of C. velenovskyi from Greece formed their own well-supported cluster (C. velenovskyi GR sublineage) separated from samples of C. orbelicus from Serbia and Bulgaria (C. orbelicus BU/RS sublineage). The most heterogeneous third lineage (coloured yellow in Fig. 3B) was split into several mostly well-supported sublineages: (1) C. tuberosus from Greece (C. tuberosus GR/MK sublineage), being genetically closest to the C. tuberosus MK lineage; (2) C. velenovskyi from Serbia (C. velenovskyi RS sublineage, TRI 258); (3) C. nissanus from Serbia (C. nissanus RS A
1079
sublineage, TRI 256); and (4) C. orbelicus from Mt Vichren (C. orbelicus Vich; TRI 151). Whereas C. velenovskyi RS and C. nissanus RS were resolved in close proximity of C. tuberosus GR/MK and C. tuberosus HR in all analyses; the position of C. orbelicus Vich was uncertain and varied across analyses. In summary, C. orbelicus Vich showed genetic affinities not only to C. orbelicus BU/RS and C. napulifer þ orbelicus BU (Fig. 3A), but also to populations of C. tuberosus GR/MK (Fig. S2). The results of Bayesian analysis based on the dataset including the C. napulifer group without the four genetically
C.orbelicus BU/RS
C. orbelicus Vich TRI 151
C. velenovskyi GR
TRI 153 C.napulifer + C. orbelicus BU
99
TRI 265 90 65
82
C. tuberosus MK
63
10·0
100 100
99
C. tuberosus HR
TRI 273 C. nissanus RS C. velenovskyi RS
C. tuberosus HR
TRI 265
C. tuberosus MK
TRI 273
C. tuberosus GR/MK
C. orbelicus Vich
C. nissanus RS
C. velenovskyi RS
C. velenovskyi GR
C. napulifer + orbelicus BU
B
C. orbelicus BU/RS
TRI 153
TRI 151
C. tuberosus GR/MK
FIG. 3. (A) Neighbour-net diagram based on AFLP profiles of 264 individuals of the Cyanus napulifer group. Numbers indicate bootstrap probabilities >50 % obtained by neighbour-joining tree analysis of the same dataset (Table S1). The symbols and colours are the same as in Fig. 1; symbols with thick outlines are used for individuals from populations with presumably hybrid origin. (B) Results of Bayesian clustering with K¼5 (BAPS software) of AFLP profiles for 264 individuals of the C. napulifer group. Each individual is represented by a vertical bar coloured according to its cluster assignment. Genetic lineages (coloured solid lines) and sublineages (coloured dashed lines) are indicated below the vertical bars.
1080
Ol savsk a et al. — The complex evolutionary history of the Cyanus napulifer group
TABLE 2. Distribution of AFLP fragments across taxa of the genus Cyanus (A) and genetic (sub)lineages of the C. napulifer group (B) Taxon/genetic (sub)lineage
NIND
NPHEN
NFRAGM
(A) Taxa of the genus Cyanus (AFLP dataset 1) C. napulifer group 264 258 67 6 561 C. pindicola 34 34 69 6 74 C. epirotus 21 21 74 6 38 C. thirkei 40 40 76 6 458 C. pichleri 28 28 84 6 385 C. fuscomarginatus 16 15 99 6 328 C. lingulatus 45 45 81 6 712 C. graminifolius 12 12 104 6 71 (B) Genetic lineages of the C. napulifer group (AFLP dataset 3) C. napulifer þ orbelicus BU 45 44 71 6 451 C. orbelicus BU/RS 32 30 66 6 429 C. velenovskyi GR 24 23 65 6 401 C. velenovskyi RS 8 8 68 6 55 C. nissanus RS 8 8 64 6 463 C. orbelicus Vich – – – C. tuberosus GR/MK 52 52 64 6 49 C. tuberosus MK 40 40 66 6 63 C. tuberosus HR 24 24 72 6 436
P ( %)
NDG
NPR
NRARE
NUNI
7488 4263 2535 3341 3502 2396 3641 2097
0 0 0 1 0 1 1 2
41 4 7 6 3 6 9 1
287 86 54 70 63 39 58 25
9 0 0 1 0 0 2 0
– – – – – – – –
011 011 007 009 011 01 008 009
3935 3225 2751 1095 1716 – 4882 5355 2751
1 0 0 0 1 – 0 1 5
19 6 3 2 4 – 113 33 12
44 25 29 11 22 – 62 71 41
5 2 0 1 1 – 5 16 4
264 83 91 20 38 – 150 243 243
009 01 009 005 008 – 011 011 007
DW
Nei
NIND, number of individuals analysed; NPHEN, number of AFLP multilocus phenotypes; NFRAGM, average number of AFLP fragments per individual 6 standard deviation; P (%), percentage of polymorphic markers; NDG, number of private fixed (diagnostic) fragments; NPR, number of private (exclusive) fragments; NRARE, number of rare fragments (present at a frequency of
