Genetic Variability in Balkan Paleoendemic Resurrection Plants Ramonda serbica and R. nathaliae Across Their Range and in the Zone of Sympatry

The genus Ramonda includes three Paleoendemic and Tertiary relict species that survived in refugial habitats of the Balkan Peninsula (R. nathaliae and R. serbica) and the Iberian Peninsula (R. myconi). They are all “resurrection plants,” a rare phenomenon among flowering plants in Europe. Ramonda myconi and R. nathaliae are diploids (2n = 2x = 48), while R. serbica is a hexaploid (2n = 6x = 144). The two Balkan species occur in sympatry in only two localities in eastern Serbia, where tetraploid potential hybrids (2n = 4x = 96) were found. This observation raised questions about the existence of gene flow between the two species and, more generally, about the evolutionary processes shaping their genetic diversity. To address this question, genetic markers (AFLP) and an estimate of genome size variation were used in a much larger sample and at a larger geographic scale than previously. The combination of AFLP markers and genome size results suggested ongoing processes of interspecific and interploidy hybridization in the two sites of sympatry. The data also showed that interspecific gene flow was strictly confined to sympatry. Elsewhere, both Ramonda species were characterized by low genetic diversity within populations and high population differentiation. This is consistent with the fact that the two species are highly fragmented into small and isolated populations, likely a consequence of their postglacial history. Within sympatry, enormous variability in cytotypes was observed, exceeding most reported cases of mixed ploidy in complex plant species (from 2x to >8x). The AFLP profiles of non-canonical ploidy levels indicated a diversity of origin pathways and that backcrosses probably occur between tetraploid interspecific hybrids and parental species. The question arises whether this diversity of cytotypes corresponds to a transient situation. If not, the question arises as to the genetic and ecological mechanisms that allow this diversity to be maintained over time.

Introduction

The Mediterranean basin is one of the richest regions in terms of biodiversity, with 23,000 to 25,000 plant species (Thompson, 2020). The Balkan Peninsula, as part of the Mediterranean Basin, one of the global hotspots of biodiversity (Myers et al., 2000), is one of the centers of species richness and endemism in Europe with more than 8,000 plant species, of which about 30% are endemic (Turill, 1929; Stevanović et al., 2007). Several geological, climatic, and biogeographical factors have been identified that have contributed to the richness and endemism of the flora of this region (Hewitt, 1996, 1999; Stevanović et al., 1999, 2007; Krstić et al., 2012; Kougioumoutzis et al., 2021). Polyploidization and interspecific hybridization have played an important role especially in the evolution of many plant taxa in the Balkan Peninsula, e.g., Alyssum (Španiel et al., 2017), Cardamine (Šlenker et al., 2021), Centaurium (Banjanac et al., 2017), Cerastium (Niketić et al., 2013, 2021), Cotoneaster (Bogunić et al., 2021), Edraianthus (Lakušić et al., 2009), Knautia (Frajman et al., 2015), Luzula (Bačič et al., 2016), Sesleria (Kuzmanović et al., 2013; Lazarević et al., 2015), Sorbus (Hajrudinović et al., 2015), Veronica (Padilla-García et al., 2018).

Polyploidization and hybridization are considered important forces in plant evolution (Wendel, 2000; Mallet, 2007; Soltis et al., 2009; Jiao et al., 2011; Van de Peer et al., 2021). Polyploids may differ from their diploid relatives in their morphology, phenology, physiology, biotic interactions, distribution, preferences for ecological niches, and invasiveness (e.g., Levin, 2002; Raabová et al., 2008; Maherali et al., 2009; te Beest et al., 2012; Weiss-Schneeweiss et al., 2013; Stutz et al., 2016; Decanter et al., 2020; Ulum et al., 2021). Many studies suggest better fitness of polyploids (Comai, 2005; Maherali et al., 2009; te Beest et al., 2012; Rey et al., 2017; Van de Peer et al., 2021), but opposite results can also be found in the literature (Nunvářová Kabátová et al., 2019). Polyploid individuals may be favored under stressful conditions, including major climatic or geological changes, and may provide an advantage in surviving and adapting to such conditions (Otto and Whitton, 2000; Comai, 2005; Weiss-Schneeweiss et al., 2013; Van de Peer et al., 2017, 2021; Decanter et al., 2020). Climate change may also lead to shifts in species distributions that bring previously separate species into contact and allow hybridization, including related species that have evolved to different ploidy levels (Vallejo-Marín and Hiscock, 2016).

There is increasing interest in the genus Ramonda as a model to study the role of interploidy gene flow in enhancing genetic and phenotypic diversity and the ecological factors that may contribute to the maintenance of cytotype diversity. The genus Ramonda Rich. (Gesneriaceae) includes three species—R. myconi (L.) Rchb., R. nathaliae Pančić & Petrović and R. serbica Pančić. All three species are considered relicts of the Tertiary (Košanin, 1921; Turill, 1929; Petrova et al., 2015), remnants of the tropical-subtropical flora from the time when the climate in southern Europe was warmer and more humid. Apart from being able to survive almost complete desiccation of their aboveground tissues in a state of anabiosis when habitat conditions are dry, and rehydrate when humidity becomes favorable again (Živković et al., 2005; Rakić et al., 2014), the three species also share similar habitats. They occur exclusively on north-facing rocky slopes of gorges, canyons, and mountain ravines, where conditions are wetter and more sheltered, allowing them longer active periods (Košanin, 1921; Rakić et al., 2014). However, R. myconi and R. serbica are found only on limestone, while R. nathaliae also inhabits ultramafic soils (Košanin, 1921; Stevanović and Stevanović, 1985; Lazarević et al., 2013; Rakić et al., 2013). Ramonda myconi and R. nathaliae are diploids with 2n = 2x = 48 chromosomes and mean 2C values of 2.59 and 2.32 pg, respectively, while R. serbica is hexaploid with 2n = 6x = 144 and mean 2C = 7.91 pg (Siljak-Yakovlev et al., 2008). Relatively high basic chromosome number (x = 24) indicates that the three species may be paleopolyploids (Rakić et al., 2014). Also, based on the similarity of monoploid genome size in R. myconi and R. serbica, the same authors hypothesized past existence of a common diploid ancestor with 2n = 24 chromosomes from which later originated R. nathaliae on one side, and R. myconi and R. serbica on the other.

Despite ecological similarities, the ranges of three Ramonda species differ, with R. myconi inhabiting the NE part of the Iberian Peninsula, while R. nathaliae and R. serbica are endemics of the Balkan Peninsula. Ramonda nathaliae is distributed mainly in the Aegean watershed area (North Macedonia, N Greece and SW Serbia—one locality in SC Kosovo at the foothill of Mt. Šara) with small enclaves in E Serbia (the Black Sea watershed). In contrast, R. serbica is mainly present in the areas of the Adriatic watershed (Albania, S and SE Montenegro, W North Macedonia, SW Serbia—including W-SW Kosovo, NW Greece) and to a smaller extent in the Black Sea (SW, E and NE Serbia, NW Bulgaria) and the Aegean watersheds (one population in NW part of North Macedonia; Figure 1A, detailed map in Rakić et al., 2014). The range of the two species in the Balkans overlaps only in E Serbia in a narrow area of sympatry (Figure 1A) discovered only a few decades ago (Stevanović et al., 1986). In only two places in the nearby gorges of the Jelašnica and Nišava rivers (6–7 km apart), R. nathaliae and R. serbica can be found in the same macrohabitat, in the close proximity to each other opening possibility for hybridization. At the time of discovery of these two sites of sympatry, the authors did not notice any intermediate phenotypes that might indicate hybridization between the two species. However, thanks to a preliminary cytogenetic study at the two sites of sympatry, several individuals with intermediate chromosome number (2n = 4x = 96) and genome size (2C = 5.14 pg) were found at both locations (Siljak-Yakovlev et al., 2008). The presence of tetraploid plants has suggested the possibility of interploidy and thus interspecific hybridization. While two Ramonda species from the Balkan Peninsula can be clearly distinguished on the basis of leaf and flower characteristics (Pančić, 1874, 1884; Petrović, 1882), field and laboratory observations of tetraploid individuals have shown that it is very difficult to distinguish them from the parent species on the basis of morphological traits. Preliminary studies of flower morphology suggest greater similarity of tetraploid hybrids to R. serbica (Figure 1B) and very careful use of certain flower characters (the color of the anthers and the angle formed by the lines connecting the petal base and the points of maximum petal width) when attempting to distinguish tetraploid individuals from the canonical representatives of diploid R. nathaliae and hexaploid R. serbica (Lazarević et al., 2014).

FIGURE 1

Figure 1. (A) Geographic range of R. nathaliae and R. serbica and distribution map of investigated populations. Green color and green dots represent range and investigated populations of R. nathaliae. Red color and red dots represent range and investigated populations of R. serbica. Detailed distribution of two species in the zone of sympatry is presented in the upper left inset. Geographical position of the investigated area in Europe is presented in the upper right inset. SO—sympatry at Oblik, SRK—sympatry at Radovanski Kamen. See Table 1 for detailed population legend. (B) R. nathaliae from population Divljana (a); R. serbica from population Radovanski Kamen (b); two species side by side in the sympatry Radovanski Kamen (c); tetraploid hybrid from Radovanski Kamen (d) and hybrid of higher ploidy level [probably backcross; (e)] showing some flowers with six petals (f).

Thus, the genetic status of the tetraploid individuals is unclear. Genetic analyses are needed to advance the study of possible hybridization between these two species. However, cytotype diversity in these two contact zones has not yet been thoroughly evaluated because initial studies were based on only a few individuals (Siljak-Yakovlev et al., 2008). Furthermore, if tetraploids are true hybrids, the question arises whether backcrosses between tetraploids and the two parental cytotypes are possible. Finally, the potential ability of these two species to hybridize raises the more general question of how genetically distinct they are from each other and what ecological and historical conditions have led to their differentiation. However, their genetic diversity has never been studied throughout their range up to date.

The aims of the present study were: (1) to investigate the variation in genome size within the two locations of sympatry between R. nathaliae and R. serbica, but also in several monospecific populations, to infer cytotype diversity; (2) to generate genetic profiles (using AFLP markers) of individuals within these two zones to test whether there is gene flow between the two species and to gain further insight into the genetic mechanisms that might explain the observed cytotype diversity; (3) to assess the genetic diversity and population structure of both species and their genetic differentiation using AFLP markers on a sample of monospecific populations covering most of their range in the Balkans. Knowledge of these genetic patterns is essential to provide a reference situation for the analysis of gene flow between the two species and especially in the two local areas of sympatry.

Materials and Methods

Plant Material

Analyses included two species of the family Gesneriaceae from the Balkan Peninsula, Ramonda nathaliae and R. serbica from monospecific populations, and a significant sample of individuals from the only two sites of sympatry known to date, both located in E Serbia (Table 1; Figure 1A). For chromosome counting, live plants were collected in the field and grown at the Institute of Botany and Botanical Garden “Jevremovac” in Belgrade and at the Laboratoire d’Ecologie, Systématique et Evolution à l’Université Paris-Saclay. For genome size estimation (GS), fresh leaves were collected in the field, brought to the laboratory in zip bags and stored in the refrigerator at 4°C until analysis. Leaves used for DNA extraction were collected directly in the field and immediately placed in silica gel-filled bags. Voucher specimens are deposited in the herbarium of the Institute of Botany and Botanical Garden “Jevremovac,” Faculty of Biology, University of Belgrade (BEOU).

TABLE 1

Table 1. List of Ramonda populations included in the study. Dashed line represents the boundary between populations: used only for genome size estimation, or for two different molecular studies. In the last three columns, “1” indicates which analyses the individual populations were used for. Number of individuals included in molecular analyses is given in parentheses in the column AFLP.

Chromosome Preparation and Estimation of Nuclear DNA Content by Flow Cytometry

Mitotic chromosome plates were prepared for two populations (one of R. nathaliae and one of R. serbica) as an addition to previously published chromosome numbers. This was done using pretreatment with 0.002 M 8-hydroxyquinoline for 4.5 h at 16°C, followed by cold fixation in 3/1 (v/v) ethanol/acetic acid for 24–48 h, hydrolysis for 14 min at 60°C, staining with Schiff reagent (Feulgen and Rossenbeck, 1924) and applying classical squash technique. Two to five individuals per population were used for chromosome counting.

Genome size was estimated for the first time in ten monospecific populations (four of R. nathaliae and six of R. serbica) as an addition to previously published results, as well as in two localities where these species grow in sympatry. Total nuclear DNA amount was determined by flow cytometry according to Marie and Brown (1993) on fresh leaves from at least five individuals per monospecific population (or three individuals in only two populations) and from 213 individuals from two hybrid zones. Pisum sativum cv. Long Express (2C = 8.37 pg) and Lycopersicon esculentum cv. Roma (2C = 1.99 pg) were used as internal standards. Nuclei were stained using propidium iodide, and measurements were done by Elite ESP flow cytometer (Beckman-Coulter, Roissy, France) equipped with argon laser. Both methods are described in detail by Siljak-Yakovlev et al. (2008).

DNA Extraction and AFLP Analyses

AFLP analyses were performed by two independent molecular biology assays. A first experiment was designed to evaluate the genetic structure of populations of both species based on an initial sample of monospecific populations (16 populations in total—five of R. nathaliae and 11 of R. serbica, Table 1). This first analysis is referred to as the large-scale geographic study. The second AFLP experiment was above all designed to assess the hybridization potential of the two species when they occur in sympatry. Ten populations were studied: four monospecific populations of each species were included in this second analysis, as well as individuals collected at two sites of sympatry, Radovanski Kamen and Oblik (Table 1). Each individual from sites of sympatry was first identified as R. nathaliae, R. serbica, or “potential hybrid” based on its morphological characteristics (leaves and flowers). The material collected from individuals in sympatry consisted of leaves from 26 individuals potentially corresponding to R. serbica, 29 to R. nathaliae, and 40 individuals with ambiguous phenotypes at Oblik, and 11 individuals potentially corresponding to R. serbica, nine to R. nathaliae, and 35 individuals with ambiguous phenotypes at Radovanski Kamen. At each location of sympatry, sampling was conducted along a transect across the entire population. We systematically avoided sampling adjacent plants to eliminate the possibility of sampling clonal individuals since vegetative reproduction also occurs (Lazarević et al., 2013).

Fragments of leaves dried in silica gel weighting 21–26 mg were used for DNA extraction, which was performed using the NucleoSpin^®96 Plant Kit (Macherey–Nagel) according to the manufacturer’s instructions. AFLP analyses were performed according to the method described by Vos et al. (1995) with some modifications (see Supplementary Material 1 for details).

Three primer combinations were selected for the full analysis (fluorescent dyes in parentheses): in both studies, EcoRI-AAC/MseI-CTC (FAM) and EcoRI-CAG /MseI-CTC (VIC) were used, whereas the third primer pair was different: EcoRI-ACA/MseI-CTC (NED) was used in the first AFLP experiment (monospecific populations only), and EcoRI-CAG/MseI-CTT (VIC) in the second experiment (hybridization study).

Amplified fragments were detected using the ABI Prism 3,100 sequencer on the Gentyane platform (INRAE, Clermont-Ferrand, France) and analyzed by GENESCAN 1.1 software (Applied Biosystems). DNA fragment detection results were indicated as present (1) or absent (0) by GeneMapper software version 4.0. Only unambiguous peaks in the size range of 50 to 500 bp were scored, with the minimum peak height set at 30. Each peak from each profile was checked by eye and its status changed accordingly. Finally, singletons (fragments that occur only in a single individual in the entire dataset) were removed to increase the reliability of the AFLP data.

To ensure the best possible reliability of AFLP markers, all samples from each AFLP experiment were managed in a single manipulation (same digestion, ligation, and PCR amplification mixtures and conditions for all samples). Two types of repeatability tests were performed. Four individuals, two from R. nathaliae and two from R. serbica, were used for the preliminary repeatability test. They were duplicated from the same amplification and then, genotyped with the same mixture of formamide and molecular weight marker. In this way, the reproducibility of the detection of the AFLP profiles was estimated (97.15%). In the second test, twenty individuals were amplified twice with two independent mixtures and the percentage of amplification reproducibility was 99.07%.

Because the samples from the two AFLP experiments were managed separately, the two datasets were not combined for the population genetic analyses to ensure better reliability of the analyses.

Statistical Analyses and Population Genetics

Statistical analyses of genome size were performed in R4.1.2. For monospecific populations, each species was analyzed separately because genome size values differ greatly between the two Ramonda species. For small sample sizes, nonparametric tests were used (Fligner-Killeen test for homogeneity of rank variation and Kruskal–Wallis and Nemenyi tests for median comparisons). For pairwise comparisons, p values were corrected according to the Holm correction for multiple testing. Student t-test and Welch test (when homoscedasticity was rejected) were used for mean comparisons for larger sample sizes (n > 15).

AFLP-SURV 1.0 (Vekemans, 2002) was used to infer genetic diversity within populations. Allele frequencies were estimated from the presence or absence of fragment phenotypes using the Bayesian method described in Zhivotovsky (1999) with a non-uniform prior distribution of allele frequencies.

Genetic dissimilarity between individuals and populations was estimated using Shared Allele Distance (Chakraborty and Jin, 1993) with POPULATIONS 1.2.30 (Langella, 1999). Principal coordinate analysis (PCoA) was performed on the individual distance matrix. Neighbor-joining dendrograms were built using the population and individual distance matrices.

The distribution of genetic variation (between species, between populations, and within populations) was estimated by analysis of molecular variance (AMOVA) performed in Arlequin 3.5.1.2 (Excoffier and Lischer, 2010; http://cmpg.unibe.ch/software/arlequin35/Arlequin35.html). To test for correlations between geographic distance and genetic variance, a Mantel test was performed using the same software. The significance of Φ-statistics was tested based on the distributions of Φ-statistics calculated under the null hypothesis from 1,000 virtual samples obtained by random permutations of the original data set.

The genetic structure of the monospecific and the two populations of sympatry was studied using genetic model-based and non-genetic model-based methods. For the first type of analysis, Bayesian inference of genetic clusters and assignment of individuals was performed using STRUCTURE v.2.3.4 (Pritchard et al., 2000). STRUCTURE analyses were performed assuming admixture, correlated allele frequencies, and no prior population information. In a first step, 10 replicates were performed for each value of the genetic cluster (K), where K varied from 1 to 12, with a burn-in period of 50,000 iterations followed by a run length of 120,000 iterations of the Markov chain. The most probable number of groups was determined using the method of Evanno et al. (2005) implemented in Structure Harvester v.0.6.94 (Earl and vonHoldt, 2012). Then, 50 new replicates of the MCMC method were performed to determine the best K value. CLUMPP v.1.1.2 (Jakobsson and Rosenberg, 2007) was used to check for similarities between runs and average the membership probabilities. The final bar graphs showing the individual admixture coefficients were created using Structure Plot v.2.0 (Ramasamy et al., 2014).

In addition, and to evaluate the robustness of the classification used, a Discriminant Analysis of Principal Components (DAPC, Jombart et al., 2010) implemented in the adegenet package for R software was used for assigning each studied individual from the first study to a specific group based on its genetic characteristics.

Finally, in order to integrate cytogenetic data and genetic inferences, we provided distributions of ancestry coefficients values (Q-values) estimated by the Bayesian clustering analysis in sympatry and monospecific populations, according to ploidy levels inferred from genome size.

Results

Chromosome Counts and Genome Size in Monospecific Populations

As expected, chromosome number showed that all individuals sampled in monospecific populations of R. nathaliae were diploid with 2n = 2x = 48 chromosomes. The chromosome number of all individuals sampled in monospecific populations of R. serbica also confirmed that this species is hexaploid (2n = 6x = 144).

The average 2C DNA value was 2.33 pg in R. nathaliae and 7.83 pg in R. serbica (data from monospecific populations only; Table 2). Although differences in genome sizes between populations of R. nathaliae were significant (Kruskal–Wallis Chi-squared = 28.019, df = 8, p value < 1e-3), pairwise Nemenyi comparison tests yielded only three significant comparisons, all involving the population from Mt. Kozjak (N8) with slightly larger GS compared to other populations. There were also significant differences in genome size values between populations of R. serbica (Chi-squared = 40.405, df = 13, p value < 1e-3). However, this was the result of differences between populations with the smallest (Mt. Timfi and Lazareva Reka) and populations with the largest (Nisava and Crna Reka) genome size values.

TABLE 2

Table 2. Chromosome numbers and nuclear DNA contents of investigated monospecific populations of R. nathaliae and R. serbica (N = number of samples for 2C DNA value).

Large-Scale Geographic Study of Genetic Diversity in Their Entire Range

Genetic Diversity Revealed With AFLP Markers

In the first molecular study, a total of 1,077 AFLP fragments were obtained with the three primer pairs. Once singletons excluded 163 fragments were present only in R. nathaliae, and 235 fragments were private to R. serbica. However, none of these fragments was fixed for alternative alleles in the two species (fragments present in all individuals of one species and absent in all individuals of the other species). Low frequency fragments dominated in both species. The number of fragments detected was moderately higher in R. serbica than in R. nathaliae (938 in the former and 866 in the latter), a trend expected due to the larger genome size of R. serbica.

Within-population genetic diversity was moderately similar in R. serbica (mean Hj = 0.105) and R. nathaliae (mean Hj = 0.092; Supplementary Table 1). Very little variation from the expected heterozygosity was observed between populations of each species. The lowest genetic diversity was found in R. nathaliae from Mt. Suva Planina in E Serbia (A), and the highest values were found in R. nathaliae from Matka in North Macedonia (B) and R. serbica from Mt.Timfi in Greece (L).

Genetic Differentiation and Population Genetic Structure

AMOVA and PCoA analyses showed that most of the genetic diversity within the entire sample are explained by genetic differentiation between the two species, as would be expected in the case of truly independent evolutionary lineages (Φct = 0.718, see Supplementary Table 2; Figure 2A). The genetic differentiation of populations within species is also very high (Φsc = 0.318). However, AMOVA analyses for individual species showed that population genetic differentiation was much higher in R. nathaliae (Φst = 0.432) than in R. serbica (Φst = 0.244).

FIGURE 2

Figure 2. Principal Coordinate analysis performed on the large geographical scale AFLP dataset using, the Shared Allele Distance computed on: (A) individuals of both R. serbica and R. nathaliae species (the first two axes represented 48.8 and 4.0% of the total variation, respectively); (B) on populations of R. serbica alone (the first two axes represented 52.2 and 30.0% of the total variation, respectively).

Total genetic diversity (Ht) was equally low in both species with values of Ht = 0.116. Among the populations of R. nathaliae, the lowest genetic distance was observed between the populations of Mt. Vermion (E) and Demir Kapija (D) and the highest between Mt. Suva Planina (A) and Matka (B) (Supplementary Figure 1). In R. serbica, a population from Mt. Timfi (L), at the southern edge of the species’ range, and Godulja (P) in SW Serbia, were genetically the most distant from the others (Supplementary Figure 1). PCoA analysis of populations of R. serbica revealed a clear geographic pattern of genetic differentiation, contrasting northern populations from Serbia (G, F, P, N O, and M) on one side and southwestern populations from Greece, North Macedonia, and Montenegro (I, J, H, K, and L) on the other (Figures 1A, 2B). The relationship between pairwise geographic distances between populations and genetic distances was not significant for R. nathaliae (Z = r = 0.51; p = 0.131), but was significant in R. serbica (Z = r = 0.56; p < 1e-3). However, this relationship is mainly due to the differentiation between northern and southern populations of this species. Overall, geographic distance between populations explained their genetic distances poorly, as can be clearly seen for northern populations G and F on the one hand, and M, N, O, and P on the other.

DAPC analysis also showed that the vast majority of individuals of both species could be unambiguously assigned to their population of origin (Supplementary Figure 2). There were only a few cases of assignment uncertainty or assignment to a different population in R. serbica between two populations in E Serbia (Radovanski Kamen (F) and Mt. Rtanj (G)) that were only 50 km apart, and also between two geographically very close populations (I and J), both from Mt. Rumija in Montenegro.

Interestingly, a few R. serbica individuals from several populations (F, J, K, L, and N) showed a slightly intermediate genetic profile between the two species (Figure 2A), although they were well assigned to their own species in the DAPC analysis.

Evidence for Hybridization in Two Sites of Sympatry

Genome Size Variability in the Two Sites of Sympatry

Individuals from both locations of sympatry showed tremendous but discontinuous variability in genome size (Table 3; Figure 3). A first group of individuals had a mean genome size of 2.32 pg (range between 2.22 and 2.51 pg). The two mean GS values of this group in Oblik and in Radovanski Kamen were not significantly different (Student-t = 1.79, df = 73, p value = 0.078). The GS of this group was also not significantly different from the mean genome size of R. nathaliae estimated from monospecific populations (Student-t = 0.50, df = 118, p value = 0.615). Thus, the plants in this group certainly corresponded to diploid individuals.

TABLE 3

Table 3. Chromosome numbers, nuclear DNA contents and ploidy levels of individuals from sympatry between R. nathaliae and R. serbica (N = number of samples for 2C DNA value).

FIGURE 3

Figure 3. Box-plot showing the distribution of genome size of each taxon and each genome class (according to Table 3).

A second group of individuals had a mean genome size of 7.98 pg (range between 7.48 and 8.45 pg). This group showed no significant difference between Oblik and Radovanski Kamen (Welch-t = −0.067, df = 22.53, p value = 0.96). Although the mean GS of this group at the two sites was significantly higher than the mean of the monospecific populations of R. serbica (Student-t = −2.851, df = 105, p value = 0.005), it remained within the range of variation of R. serbica. Pairwise comparisons showed that only the comparison with the Lazareva Reka population was at the threshold of significance (data not shown). Therefore, this group of individuals was likely hexaploid at both sites of sympatry.

However, other classes of genome size values suggest that there are several other ploidy levels at both sites of sympatry with genome composition contributed by both species. Indeed, a first class of individuals had a mean GS value of 5.17 pg (range between 4.96 and 5.39) (cf. Class I in Table 3). The individuals of this class in Oblik and Radovanski Kamen did not show any significant difference in their mean GS (Student-t = −1.65, df = 68, p value = 0.104). The mean GS of this entire class was very close, although significantly higher (Student-t = 7.79, df = 69, p value = 4.8e-11) than the expected value for a first-generation tetraploid interspecific hybrid (5.08 pg) calculated from the mean GS values of the two Ramonda species. This supported that these individuals are interspecific tetraploid hybrids. It is noticeable that a few individuals exhibited a combination of morphological characteristics of the two species (Figure 1B), suggesting that they may be first-generation interspecific hybrids. These plants were very abundant in our sample, representing 70 out of 214 individuals analyzed in the two sites of sympatry.

Finally, two other classes of GS values were also observed in the two sites of sympatry. Fifteen analyzed individuals had a higher genome size than the hypothetical tetraploid hybrids, ranging from 5.52 pg to 6.55 pg (cf. Class II in Table 3, Figure 3). In some of these individuals, the GS was close to what would be expected in a 5x interspecific hybrid. For example, a 5x individual with two R. nathaliae genome complements and three R. serbica genome complements would be expected to have a GS of about 6.27 pg. Thus, the plants in this Class II could correspond to individuals with mosaic genomes in which the proportion of the two species is different than in a first-generation tetraploid hybrid.

In both populations, a few individuals (10 in total) had a higher genome size than the hexaploid R. serbica (from 8.94 to 12.06 pg, cf. Class III in Table 3).

Finally, three individuals sampled at Radovanski Kamen were found to have an intermediate genome size between R. nathaliae and tetraploid interspecific hybrids (Class IV, Table 3). Their GS values were close to what would be expected in triploid individuals of R. nathaliae.

Taxonomic identification of individuals in monospecific populations of the two Ramonda species based on leaf and flower morphological characters is always trivial and unambiguous (Siljak-Yakovlev et al., 2008; Lazarević et al., 2013, 2014). However, this was not the case for all individuals sampled at the two sites of sympatry (Figure 1B). For this reason, the GS values were also used to confirm or assign a taxonomic identification to be used in the analyses of the AFLP data (see below). In the following, all individuals in sympatry with the same genome size as R. nathaliae were classified as this species, and the same principle was applied to R. serbica. Other individuals are hereafter referred to as potential “hybrids” and assigned to their GS classes as indicated in Table 3.

Genetic Variability by AFLP Markers in the Two Sites of Sympatry

The following analyses included eight monospecific populations (four for each species) that were used as references in addition to the individuals sampled at the two sites of sympatry (Table 1). After removing singleton fragments, 1,263 AFLP fragments remained for analyses. None of the fragments were present in all individuals of one species and absent in all individuals of the other species, similarly to what was found in the large-scale geographic study. However, some fragments were highly abundant in one species and absent in the other, consistent with the data from the first set of populations (see above).

A strong contrast in AFLP fragment frequencies between R. serbica and R. nathaliae would have been expected because of the genome size and ploidy level of R. serbica. In fact, the distribution of the difference in fragment frequencies between the two species was clearly centered at zero (median = −8e-4; mean = 0.049) (Supplementary Figure 3). Thus, the dominant alleles at the AFLP markers (presence of the fragment) were almost equally distributed between the two species. Importantly, this result supports the assumption that in this study the dominance effect at AFLP markers was unlikely to cause bias in favor of either species in individual assignment in subsequent clustering analyses.

The strong genetic differentiation between the two species, already reported in the large-scale geographic study, was confirmed when only monospecific populations were considered (Φct = 0.548; Figure 4A). In contrast, it was noticeable that no genetic differentiation (Φct = −0.009) was observed between diploid and hexaploid individuals from the two sites in sympatry, as also evidenced by the strong mixing of individuals from both ploidy levels in the dendrogram (Figure 4B). The tetraploid plants also exhibited similar genetic variability as the diploid or hexaploid plants in sympatry as shown by the very strong mixing of the three ploidy levels in Figure 4B. However, clear and even strong genetic differentiation was observed between individuals from the two sites of sympatry for each ploidy level (Φst = 0.22 for diploids and Φst = 0.10 for hexaploids). Consistent with this result, clusters of individuals from each location of sympatry, regardless of their ploidy level, are clearly visible in Figure 4B.

FIGURE 4

Figure 4. Dendrogram built on inter-individual Shared Allele Distance using: (A) individuals of R. serbica and R. nathaliae sampled in allopatric monospecific populations (each color indicates a specific population whose code names are referred to in Table 1); (B) individuals from the two sites of sympatry [Oblik (Obl) and Radovanski Kamen (RK)]. All those individuals contributed to the interspecific gene flow study.

Bayesian clustering yielded two genetic clusters as the solution that best fit the data, according to the Evanno test (Supplementary Figure 4) and the very clear plate obtained by the likelihood of the data at K = 2 (data not shown). This result was to be expected, since species differentiation is the factor that contributes most to the genetic structure of the population in the entire data set. Indeed, from the data of monospecific populations, it appears that the two genetic clusters correspond to one of each species (Supplementary Figure 5). DAPC analysis confirmed this result. Indeed, although K = 3 was the optimal solution (Supplementary Figure 6), one of those two clusters was composed only with one individual of the S7 monospecific population of R. serbica (Supplementary Table 3). As for the Bayesian clustering analysis, the two main clusters corresponded to R. nathaliae and R. serbica. We will therefore refer to “nathaliae” and “serbica” for these two clusters below. The distribution of individual ancestry coefficients showed two very different patterns between the monospecific populations and sites of sympatry (Figure 5, Supplementary Figure 5). Most individuals in monospecific populations had high ancestry values of the species where they belonged (mean ancestry of 0.93 for individuals of R. serbica to the genetic cluster “serbica” and 0.96 for individuals of R. nathaliae to the genetic cluster “nathaliae”; Figure 5A). Therefore, the AFLP data allowed unambiguous assignment of nearly all individuals sampled in monospecific populations to their own species, as already observed in the large-scale geographic study. Only three individuals of R. serbica, all from the gorge of Nišava river (S9), had an intermediate ancestry coefficient (Figure 5A), with the ancestry coefficient to the “serbica” genetic cluster ranging from 0.43 to 0.59. In comparison, the lowest value of ancestry to the “nathaliae” genetic cluster was up to 0.75 for individuals sampled in monospecific populations of R. nathaliae.

FIGURE 5

Figure 5. Distribution of individual ancestry coefficients (inferences from the Bayesian clustering analysis) to the “serbica” genetic cluster for: (A) R. nathaliae (red) and R. serbica (blue) allopatric populations; (B) for individuals sampled in the two sites of sympatry, represented according to their genome size. “Other” refers to individuals that could not be undoubtedly assigned to a given ploidy level from their genome size.

The two sites of sympatry showed a very different and more complex genetic pattern, similar to what was observed for genome size. The main difference with the monospecific populations was that there was much less agreement between the genome size values (and thus the inferred ploidy level) and the AFLP profiles. In fact, when 2x and 6x individuals were considered, the distribution of ancestry coefficients was bimodal for both cytotypes (Figure 5B). Surprisingly, 17 of the 35 diploid individuals had high ancestry coefficients to the hexaploid R. serbica (mean = 0.88; range: 0.68–0.99 for these 17 individuals). All but one of them were sampled in the Oblik site of sympatry. The remaining 18 diploid individuals had very high ancestry coefficients to R. nathaliae, as expected based on their assumed ploidy level (mean = 0.96; range: 0.83–1). In contrast, among hexaploid individuals, five of the total 30 individuals, all from Oblik, showed strong ancestry to R. nathaliae, contrary to what was expected based on their genome size.

The bimodal character of the distribution of the ancestry coefficient distribution was also observed for tetraploid individuals (Figure 5B), although the right side of the distribution was much flatter than that observed for monospecific R. serbica populations. This also testifies to the high genetic heterogeneity of tetraploid plants. The mean ancestry value of the 40 individuals belonging to the right side of the distribution was 0.80, very close to the value of 0.75 expected in the case of a first-generation hybrid resulting from a cross between hexaploid R. serbica and diploid R. nathaliae. However, ten individuals, all from Oblik (Supplementary Figure 5), formed the left side of the distribution (i.e., individuals unambiguously assigned to R. nathaliae based on AFLP markers). Their mean ancestry value to “serbica” genetic cluster was 0.05 (range: 0.01–0.2), very close to the value obtained for R. nathaliae individuals from monospecific populations.

Apart from these three groups of individuals, eight individuals from Radovanski Kamen had rather mixed genetic profiles with a GS lower than R. serbica but higher (between 5.61 and 6.55 pg) than potential 4x individuals (cf. Class II in Table 3), but with a mean ancestry coefficient of 0.78 (range between 0.56 and 0.99) much closer to R. serbica than to R. nathaliae. Of the nine individuals in the Class III of genome size (ranging from 8.94 to 10.5 pg, Table 3), three were strongly assigned to the R. serbica genetic cluster, while two had a very high ancestry coefficient to R. nathaliae. This last observation suggests the existence of autopolyploid individuals with high ploidy level descended from R. nathaliae. Finally, the three individuals with a genome size between 3.35 and 4.01 pg (i.e., a GS between 2x and 4x) were unambiguously assigned to R. nathaliae.

It is noticeable that individual assignment posterior probabilities inferred by the DAPC analysis and ancestry coefficients inferred through the Bayesian approach were highly correlated (0.95 for all individuals in sympatry and 0.96 for the whole sample).

Both GS and AFLP patterns showed that more complex genomic structures than first-generation interspecific tetraploid hybrids are present at both sites of sympatry.

Discussion

Chromosome Number, Genome Size, and Ploidy Level in Allopatric Populations

In the past, there were uncertainties concerning the chromosome number and ploidy levels of Balkan Ramonda species. Although Glišić (1924) found n = 36 to be the most frequent haploid chromosome number in R. nathaliae, later both Ratter (1963) and Contandriopoulos (1966) showed that R. nathaliae is a diploid with 2n = 2x = 48 chromosomes, and this is confirmed by broad and detailed cytogenetic studies by Siljak-Yakovlev et al. (2008) and Lazarević et al. (2013). Still, some uncertainties remained about chromosome number and ploidy level in R. serbica based on the following reports: 2n = 72 in R. serbica individuals collected in Sićevo Gorge in E Serbia (Glišić, 1924), 2n = 96 in plants from Mt. Timfi—Vikos Gorge in Greece (Contandriopoulos, 1966), polyploid forms as a mixture of several chromosome numbers, from less than 72 to more than 144 chromosomes, with 2n = 72 and 2n = 96 being most common in only two populations in Albania (Daskalova et al., 2012). In this last study, five individuals per population were analyzed and hexaploid individuals were not detected. However, one of the populations analyzed for genome size in the present study is Kruje, the same population as in Daskalova et al. (2012), and only 6x individuals were found. In the last 15 years, Siljak-Yakovlev et al. (2008) and Lazarević et al. (2013), as well as the present study, have analyzed the genome size and/or chromosome number of a total of 133 individuals sampled from 16 monospecific populations of R. serbica throughout the species range. The vast majority of them were hexaploid individuals with 2n = 6x = 144 chromosomes. Only in one population from Montenegro were two individuals found with a larger genome size, suggesting the possibility of the existence of 8x and even 10x individuals (Siljak-Yakovlev et al., 2008). All this indicates possible but rare alternative chromosome numbers in monospecific populations.

Genetic Diversity in Ramonda Species From the Balkan Peninsula

To date, very few studies have examined the genetic diversity of the two Balkan Ramonda species. In a study on the genetic variability of three populations of R. serbica from Albania and three populations of R. nathaliae from North Macedonia using RAPD markers, Daskalova et al. (2012) could not clearly distinguish the two species. More recently, Petrova et al. (2018) studied five populations of R. serbica in Bulgaria using ISSR markers. They showed that this species has low genetic variability and significant differentiation among populations. The present study is the first investigation of the genetic diversity of R. nathaliae and R. serbica as a whole, using the same molecular markers and including a large coverage of their respective ranges in the Balkan Peninsula. The number of AFLP markers is also high, and thus, we are confident that genetic diversity and structure patterns here reported are representative. The high level of genetic differentiation is in agreement with morphological data that have shown that the two species are broadly consistent across their range, with the exception of the two sites of sympatry (Lazarević et al., 2013, 2014). Therefore, it is clear that the two species actually correspond to two independent evolutionary lineages, suggesting a unique origin of hexaploid monospecific populations in their range. Moreover, both species showed a similar picture: a very high level of genetic differentiation between populations within species and very low signals of genetic admixture between populations of each species, except between a few pairs of geographically close populations.

What information about the evolutionary forces that shaped genetic diversity in these two species can be inferred from this pattern? There is no doubt that the dynamics and extent of the glacial periods in the mountains of the central part of the Balkan Peninsula contributed to the genetic and chorological differentiation of these two species. Glaciation in this region was much less pronounced than in other parts of Europe (Hughes et al., 2006), raising the possibility of existence of multiple and scattered refugia for plant and animal species survival in the Balkan region (Hewitt, 1996). Available data show that the snowline in different parts of the Mediterranean region, and especially in the Balkan Peninsula, extended to different heights, from down to 1,250–1,600 m in the highest mountains of the region (Menkovic et al., 2004; Hughes et al., 2006), to 270 m in some coastal parts of the Dinarides (Marjanac and Marjanac, 2004). There is a possibility that the Balkan Ramonda species were mountain plants that descended to refugia in the surrounding gorges and canyons during the Ice Age, where they are still predominantly found today (Košanin, 1921; Stevanović et al., 1986; Lazarević, 2012). Also, the occurrence of a series of lakes during the Miocene and Pliocene, which spread in different parts of the Balkan Peninsula (Krstić et al., 2012), may have had a significant impact on the distribution of plant species and contributed to the isolation of Ramonda populations from each other.

Thus, the strong genetic differentiation within species and between populations exhibited by the two Ramonda species is likely the result of long reproductive isolation of populations and strong drift effects during glacial periods, but possibly before as well (Krstić et al., 2012), as the geographic distribution of populations is highly fragmented. A similar pattern was found for R. myconi in N Spain and the central Pyrenees, corresponding to its entire range, with high genetic differentiation between regions (Dubreuil et al., 2008). This pattern led the authors to suggest that this species occurred in several different but genetically unrelated refugia in this area. Similarly, in another Tertiary relict species from the family Gesneriaceae Haberlea rhodopensis Friv. two distinct genetic clusters were identified in Bulgaria, and the existence of which was suggested to be a consequence of isolation between refugial populations (Petrova et al., 2014). Finally, Petrova et al. (2018) have shown that the genetic diversity of R. serbica populations in Bulgaria is very low. The authors hypothesize that this is the result of both post-glacial fragmentation and a particularly pronounced lack of genetic connectivity for the Bulgarian populations, which are on the periphery of the species range.

This pattern of variability is expected in long-living plants with geographically isolated populations and probably with limited seed dispersal, as is the case in Ramonda species (Picó and Riba, 2002; Lazarević et al., 2013; Petrova et al., 2014). The overall genetic diversity of the two studied species is low, similar to what was found for R. myconi in the Pyrenees (Dubreuil et al., 2008) and H. rhodopensis in Bulgaria and Greece (Petrova et al., 2014, 2015). In the post-glacial period, eventual dispersal outside refugial areas may also have led to repeated demographic bottlenecks, which should enhance or at least maintain genetic differentiation of populations. All the more, these two species still have a very patchy distribution and most populations are still spatially isolated. It is noticeable that the population of R. nathaliae with the lowest genetic diversity is found on Mt. Suva Planina, which is at the highest altitude of all populations of this species studied here. Living conditions are harsh there, and rosettes are much smaller than elsewhere (authors’ personal observation). It is therefore likely that the very low genetic diversity of this population is the result of a very small effective size, due to the combined effects of small size and potential selection.

The current ranges of these two species do not overlap, and only two narrow areas of sympatry have been observed so far (Stevanović et al., 1986). This observation and the marked genetic differentiation between these two species may indicate both that the two species were separated prior to the ice ages and that their refugia were distant from each other. However, despite the high genetic differentiation between the two species, no specific markers (alternative alleles fixed in either species) were found among the hundreds of AFLP fragments revealed in this study. However, this result should be confirmed on the basis of orthologous sequence data. The absence of specific markers may be evidence of a relatively recent divergence between the two species, which is consistent with the results of Petrova et al. (2015) on Gesneriaceae species. Based on chloroplastic and ITS sequence polymorphisms, these authors inferred that the timing of the divergence between R. nathaliae and R. serbica may have occurred in the middle or last Pleistocene. This divergence is one of the youngest divergence times among Gesneriaceae species estimated in this study. This could also mean that the genus Ramonda is a relict of the Tertiary, but that the individual species are actually younger.

Interploidy Gene Flow

Polyploidy is common in plant species and is mainly favored by meiosis defects leading to the production of unreduced gametes (Brownfield and Köhler, 2011; Mason and Pires, 2015; Kreiner et al., 2017). Intrageneric variability in ploidy level is widespread in the literature, and situations have been described where different cytotypes coexist in large contact zones as well as in geographically much narrower areas. Therefore, the possibility of interploidy gene flow within the same species or between related species has been a focus of attention for many years. However, the vast majority of cases reported in the literature in natural species involved diploid and tetraploid cytotypes, whereas cases of contact zones between diploid and hexaploid plants have been reported less frequently (see Kolář et al., 2017 for a comprehensive review of reported cases). In Aster amellus L., for example, diploid and hexaploid cytotypes are usually found in separate locations, while a few parapatric and one mixed-ploidy population have been observed. Based on fine-scale cytotype screening, the authors concluded that there is no evidence for the existence of tetraploid plants (i.e., first-generation hybrids) probably as a consequence of the exclusion of minority cytotypes (Castro et al., 2012), as well as niche differentiation and shifts in flowering times between cytotypes (Raabová et al., 2008). One of the most comparable situations to that observed in Balkan Ramonda is the case of Senecio carniolicus Willd. (syn. Jacobaea carniolica (Willd.) Schkrank), an autopolyploid complex with diploid and hexaploid plants as the most widespread ploidy levels. In contrast to Ramonda species, a high number of mixed ploidy populations were found in this species (Sonnleitner et al., 2010). Tetraploid individuals were observed in restricted areas, which the authors believe may correspond to glacial refugia in the eastern Alps. Other ploidy levels (3x, 5x, 7x−9x) are also observed, but at very low frequency. The authors noted that there should be barriers limiting the production and/or maintenance of interploid hybrids, to which the distinction between 2x and 6x habitats contributes. Lower fitness of progeny from crosses between diploids and polyploids also clearly contributes to limited gene flow between cytotypes (Sonnleitner et al., 2013).

It is necessary to emphasize that in the case of Ramonda, in contrast to the examples cited above, it is a question of two cytotypes (2x and 6x) belonging not to one but to two different species. This makes hybridization more difficult and the genome constitution of the hybrids extremely diverse and sometimes difficult to explain. The case of Ramonda documented in this article therefore stands out from most studies of mixed-ploidy populations. It is interesting because, first, it shows enormous variability in genome size at micro-level in the only sites of sympatry between two species. This variation in GS allows us to hypothesize that the observed ploidy levels are also diverse, ranging from at least 2x to 8x (or possibly up to 10x), including odd ploidy levels such as at least 3x and 5x. For example, the class III of genome size (range between 8.94 and 12.06 pg, Table 3) could include 7x, 8x, but also higher ploidy levels. Based on the mean 1Cx values of R. nathaliae and R. serbica, one would expect GS values between 8.15 and 9.24 pg for 7x individuals. For 8x individuals, the expected range of variation would be between 9.32 and 10.56 pg. Therefore, we could not exclude the possibility that 9x and/or 10x ploidy values were also present in our sample.

Octoploid individuals may have arisen from whole genome duplication (WGD) of tetraploid individuals, resulting in potential stabilization of the hybrid genome, a phenomenon already suggested by several authors (reviewed in Glombik et al., 2020). In contrast, the extended study of genome size in the two Balkan Ramonda species reported here showed that there was no cytotype variability within each species in monospecific populations, except for two individuals of R. serbica from Montenegro that were presumably 8x and 10x (Siljak-Yakovlev et al., 2008). Aneuploidy could also contribute to the observed variation in GS in both sites of sympatry. However, these hypotheses remain to be confirmed by cytogenetic studies. Finally, genome rearrangements and transposition activation have been documented in artificial interspecific hybridization experiments (Szadkowski et al., 2011; Liu et al., 2014; Mason and Wendel, 2020). We cannot rule out that such phenomena also contributed to the large variation in genome size observed in both sites of sympatry.

Nevertheless, the proportion of cytotypes distinct from the canonical ploidy levels 2x and 6x (a total of 33% of 4x plants) appears to be very high in these two sites of sympatry compared with most previously reported studies (Kolář et al., 2017). Other classes of cytotypes were also not rare [for example, 6% of plants suspected of being 5x (referred as class II in Table 3)].

The existence of hybridization between the two species, viz. between 2x and 6x cytotypes, in each of the two sites of sympatry was clearly demonstrated by: (i) the presence of individuals with an intermediate genome size compatible with the expected value for F1 tetraploid hybrids given the genome sizes of the two parental species; (ii) the presence of tetraploid individuals combining morphological features of the two species; and (iii) the observation that diploid and hexaploid individuals at the two sites of sympatry differed genetically much less (based on AFLP markers) than in monospecific populations. Several factors make this hybridization between the two species in sympatry possible. Both species flower at the same time. Although the flowers are slightly different, they are most likely pollinated by the same pollinators (Lazarević et al., 2013). Our results showed therefore that pre- and/or post-zygotic barriers to hybridization that usually exist between different cytotypes and/or plant species (Vallejo-Marín and Hiscock, 2016) are not operative in two sites of sympatry (Mallet, 2007; Abdelaziz et al., 2021).

The joint distribution of cytotypes and ancestry coefficients estimated from AFLP data in the two sympatry sites revealed a very singular pattern, namely a bimodality of these coefficients for diploids and hexaploids, compared to that observed in allopatric monospecific populations, but also for tetraploids. This striking observation may help to better understand the evolutionary dynamics that took place at these two sites by suggesting hypotheses on the pathways by which cytotype diversity was achieved and on the processes that might have been involved in the maintenance of tetraploid individuals and that would explain their genetic constitution.

The high number of tetraploid plants and their high AFLP diversity, similar to that of diploid and hexaploid plants in both sites of sympatry, suggests that hybridization between the two species is not an uncommon event in these sites. Alternatively, it is possible that only a few founders of tetraploid hybrids are at the origin of all currently observed tetraploid plants and that preferential breeding among plants of this cytotype has provided the opportunity for genetic recombination between the genomes of R. serbica and R. nathaliae. This would result in a wide variety of genetic combinations from the two parental genomes, including the existence of genetic profiles similar to the parental types (i.e., R. nathaliae or R. serbica). This latter possibility would imply that synapsis occurs between the chromosomes of the two species, which remains to be proven. However, in that case, one would expect to observe a continuous distribution of ancestry coefficients within the tetraploid group, from “nathaliae” to “serbica” and all possibilities in between. Data reported here showed a different pattern with a clear bimodality in this distribution. The coexistence of divergent genomes in interspecific hybrids is considered as a major mechanism of molecular and phenotypic novelty. This is explained by the huge structural and functional variation caused both by the change in ploidy level, if any, and by novel genetic interaction effects generated between the divergent genomes that merge (see Nieto Feliner et al., 2020, for a review). This variation is a raw material on which selection can operate. For example, negative epistasis can contribute to post-zygotic barriers to interspecific gene flow through a decrease in fitness of interspecific F1 hybrids and/or their progenies (Taylor et al., 2009), limiting therefore genetic introgression into restricted regions of genomes (adaptive introgression, e.g., Marburger et al., 2019). It is therefore possible that, among tetraploid Ramonda interspecific hybrids, selection in favor of “parental” combinations of co-adapted genes has taken place, because of weak survival or fertility of individuals having a mixed genome composition. The AFLP patterns of the different classes of each cytotype also strongly suggest that different pathways were involved in the production of individuals of the same ploidy level, the 4x first-generation hybrids, but also other more unusual and rare ploidy levels. Three putative triploid plants have been observed that show strong ancestry from R. nathaliae. This suggested that these triploids could have been produced by fertilization between an unreduced 2x gamete of R. nathaliae and either a haploid gamete of R. nathaliae or a 1x gamete with complete R. nathaliae genome complement from an interspecific tetraploid hybrid. The latter case would imply that triploids are restricted to sympatry. In the first case, the presence of some triploid individuals would also be expected in R. nathaliae populations, but this has not been reported so far. However, the cytotype variability of monospecific populations of R. nathaliae has always been based on only a few individuals, and it cannot be excluded that such triploid individuals exist but have not yet been identified. Indeed, despite the high sample size of the present study, only three potential triploids at sympatry were found, all at Radovanski Kamen. However, it is now known that the production of non-reduced gametes is highly variable between individuals within populations and that this ability is also influenced by environmental conditions (De Storme and Mason, 2014). In addition, mainly because of chromosome pairing during meiosis, triploids are often sterile or at least with greatly reduced fertility (Comai, 2005) and therefore rarely remain in natural populations. The above-described pathway to triploid plants opens another way to explain the existence of tetraploid plants with high ancestry coefficient to “nathaliae.” One could hypothesize at least two explanations. First, it is theoretically possible to obtain such profiles in progeny resulting from fertilization between a reduced gamete of a diploid R. nathaliae and a triploid gamete of a triploid R. nathaliae. The triploid bridge is considered to be very important for interploid gene flow and polyploidization (Husband, 2004). Secondly, unreduced (2x) gametes of R. nathaliae may also be involved in the production of tetraploid individuals that possess only chromosomes of R. nathaliae by autopolyploidy. The existence of such tetraploid individuals could also open the possibility of obtaining 6x progeny with only R. nathaliae genomes, the presence of which is suggested by the AFLP data. This hypothesis deserves attention because such tetraploid “nathaliae” would have inherited the ability to produce unreduced gametes (i.e., 4x gametes) from their R. nathaliae diploid parent.

Therefore, the AFLP profiles of the tetraploid individuals may indicate the existence of an autopolyploid-allopolyploid mixture at both sites of sympatry. They also suggest that there could be some barriers to reproduction between the two tetraploid hypothetical types, since otherwise tetraploid individuals with a much higher genome admixture would arise in the progeny from such crosses, which was not observed. However, the few individuals in our sample that exhibit such an intermediate AFLP profile do not allow us to rule out this last possibility. If this is indeed the case, the tetraploid level could also be a bridge between the two species, reinforcing interploidy gene flow in sympatry.

Pentaploid individuals were also found in both places. They may have originated in different ways. For example, they may simply have come from a cross between a tetraploid (of whatever type) and a hexaploid, or even from the union of an unreduced 2x gamete from R. nathaliae and a 3x gamete from a hexaploid. These pentaploid individuals could also play a role as mediators of gene flow, but without forming an independent evolutionary unit (Tel-Zur et al., 2020; Peskoller et al., 2021).

Obviously, other pathways could be invoked to explain all patterns of GS and AFLP profiles observed in this study. For example, it has been observed (Lazarević et al., 2013) that hybrid individuals produced pollen with highly variable sizes, a well-known proxy for irregular meioses in plants. It is therefore possible that 4x plants could produce haploid gametes with one R. serbica complement, a pathway for the production of 2x progenies with the genetic feature of R. serbica, an intriguing pattern we observed in our data.

Further studies are therefore needed to gain a comprehensive and detailed understanding of the biological mechanisms that explain the cytotype and molecular diversity observed at the two sites of sympatry. Further studies should focus on finding evidence for the production of unreduced gametes, especially in R. nathaliae, but also in triploids and tetraploids. In particular, a more detailed genetic study of the genetic profiles of the plants found at the two sites of sympatry, based on large-scale sequence data, should be useful to assess the extent of genome mosaicism with respect to each cytotype. This would help to assess the extent of genetic introgression between the two species and the potential for recombination between their genomes. However, it is already clear that hybridization between the two species is occurring and that, at the very least, backcrossing (in which tetraploid hybrids give rise to pentaploids) is possible. Unfortunately, backcrosses with one or both parental species often result in the appearance of individuals whose morphology is difficult to distinguish from the parental species, and identification of hybrids based on morphology becomes very difficult or impossible (Lihová et al., 2007; Čertner et al., 2015; Tendal et al., 2018; Beirinckx et al., 2020), as seems to be the case with both sites of sympatry for Ramonda (Lazarević et al., 2014).

Our data (both cytotypes and AFLP markers) therefore suggest that the first-generation interspecific hybrids are viable and possibly fertile. It has already been shown that pollen grains from tetraploid individuals are viable, although somewhat less so than pollen from parental cytotypes in sympatry (Lazarević et al., 2013). On the other hand, seeds of tetraploid individuals are drastically smaller than those of parental species and have very low or no germination capacity (Lazarević et al., 2013), which is often the case in interspecific hybrids (Mallet, 2008; Schmickl and Yant, 2021). Interestingly, with the exception of R. nathaliae in Oblik, seeds of both parental cytotypes have shown lower germination rates compared to monospecific populations (Lazarević et al., 2013), which could be the consequence of a reduction in fitness caused by various factors such as genetic introgression or inclusion of pollen from the hybrid individuals in pollination processes (Jacquemyn et al., 2016). The large number of hybrid individuals observed in sympatry suggests that they arise from regular gene flow between parental species in these sites, as shown by Abdelaziz et al. (2021) in two Erysimum species. If hybridization events between the two species were rare, competition between the cytotypes would have been expected to eliminate the hybrids due to their rarity (cytotype exclusion), unless there were counteracting factors (homogamy, enlargement of the hybrid niche, or heterosis). Vegetative reproduction may also contribute to their maintenance (Lazarević et al., 2013). It is also possible that the process of minority cytotype exclusion either does not work for some reason that should be investigated or has not yet eroded the diversity of cytotypes in the two sites of sympatry. This last argument could support the idea that the only two known secondary contact zones between the two species are recent and probably formed during repeated glaciations in the Pleistocene. Finally, the overall panorama of genetic and cytotype variability may indicate that a secondary contact zone between the two species emerged after a long period of isolation (Schmickl and Yant, 2021).

Conclusion

The present study shows the complex influence of various factors on the evolution of two representatives of the genus Ramonda in the Balkan Peninsula. The turbulent geological past of the Balkan Peninsula, together with climatic fluctuations throughout the Mediterranean, has influenced the distribution, ecology, speciation and genetic diversity of Ramonda species, leading to the formation of isolated populations in refugial habitats. During migrations favored by interglacial periods, a secondary contact zone between R. nathaliae and R. serbica formed in a narrow area in E Serbia. In this newly formed sympatry, an intense hybridization process takes place between the two species, leading to the occurrence of different ploidy levels (from 2x to > 8x) and individuals with mixed morphological characteristics, suggesting that not only crosses between the parental species occur, but also backcrosses and possibly duplication of the whole genome. These two sites of sympatry are places for a huge genetic mixing. Moderate pollen fertility, but low seed germination and low dispersal potential suggest that the effects of hybridization are limited to the sites of sympatry. However, the rich mixture of ploidy levels shows that these relict species, which have preserved a specific survival strategy, are not evolutionary “dead ends” and raises the question of possible ecological and evolutionary consequences of this phenomenon. These two locations constitute therefore an open-air laboratory for evolutionary and interspecific recombination studies.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

SS-Y, BS, VS, and TR designed the study. ML, SS-Y, BS, VS, TR, AS, MN, and GT collected plant material. ML and SS-Y performed cytometric measurements and chromosome counts. ML, AS, and FL performed molecular laboratory work and provided data for statistical analyses done by TR and DH. ML, TR, and SS-Y drafted the manuscript. All authors read, revised, and approved the final version of the manuscript.

Funding

This research was supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia (No. 451–03–68/2022–14/200178), the International Franco-Serbian “Pavle Savić” project, and the Federative Institute of Research (IFR 87 La Plante et son environnement) at Gif-Orsay. The l’Institut Français de Serbie provided grants to ML.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank all colleagues who helped us to collect plant material for this study: Vlado Matevski, Bojan Zlatković, Predrag Lazarević, Dmitar Lakušić, Nevena Kuzmanović, Antun Alegro, Danko Jović, and Dragan Nešić. Sequencing was performed using the Gentyane platform (INRAE, Clermont-Ferrand, France). We thank Mickaël Bourge and Nicolas Valentin, and Spencer C. Brown and Olivier Catrice for their assistance on the Cytometry Platform, Institute of Integrative Biology of the Cell (I2BC), Gif-sur-Yvette (http://www.imagif.cnrs.fr). Finally, we thank two reviewers and the editor Gonzalo Nieto Feliner for their valuable comments and suggestions that considerably improved our manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.873471/full#supplementary-material

References

Abdelaziz, M., Muñoz-Pajares, A. J., Berbel, M., Garcia-Muñoz, A., Gómez, J. M., and Perfectti, F. (2021). Asymmetric reproductive barriers and gene flow promote the rise of a stable hybrid zone in the Mediterranean high mountains. Front. Plant Sci. 12:687094. doi: 10.3389/fpls.2021.687094

PubMed Abstract | CrossRef Full Text | Google Scholar

Bačič, T., Frajman, B., and Dolenc Koce, J. (2016). Diversification of Luzula sect. Luzula (Juncaceae) on the Balkan Peninsula: a cytogenetic approach. Folia Geobot. 51, 51–63. doi: 10.1007/s12224-016-9235-2

CrossRef Full Text | Google Scholar

Banjanac, T., Dragićević, M., Šiler, B., Gašić, U., Bohanec, B., Nestorović Živković, J., et al. (2017). Chemodiversity of two closely related tetraploid Centaurium species and their hexaploid hybrid: metabolomic search for high-resolution taxonomic classifiers. Phytochemistry 140, 27–44. doi: 10.1016/j.phytochem.2017.04.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Beirinckx, L., Vanschoenwinkel, B., and Triest, L. (2020). Hidden hybridization and habitat differentiation in a Mediterranean macrophyte, the euryhaline genus Ruppia. Front. Plant Sci. 11:830. doi: 10.3389/fpls.2020.00830

PubMed Abstract | CrossRef Full Text | Google Scholar

Bogunić, F., Siljak-Yakovlev, S., Mahmutović-Dizdarević, I., Hajrudinović-Bogunić, A., Bourge, M., Brown, S. C., et al. (2021). Genome size, cytotype diversity and reproductive mode variation of Cotoneaster integerrimus (Rosaceae) from the Balkans. Plan. Theory 10:2798. doi: 10.3390/plants10122798

CrossRef Full Text | Google Scholar

Brownfield, L., and Köhler, C. (2011). Unreduced gamete formation in plants: mechanisms and prospects. J. Exp. Bot. 62, 1659–1668. doi: 10.1093/jxb/erq371

PubMed Abstract | CrossRef Full Text | Google Scholar

Castro, S., Loureiro, J., Procházka, T., and Münzbergová, Z. (2012). Cytotype distribution at a diploid – hexaploid contact zone in Aster amellus (Asteraceae). Ann. Bot. London 110, 1047–1055. doi: 10.1093/aob/mcs177

PubMed Abstract | CrossRef Full Text | Google Scholar

Čertner, M., Kolář, F., Schönswetter, P., and Frajman, B. (2015). Does hybridization with a widespread congener threaten the long-term persistence of the eastern alpine rare local endemic Knautia carinthiaca? Ecol. Evol. 5, 4263–4276. doi: 10.1002/ece3.1686

PubMed Abstract | CrossRef Full Text | Google Scholar

Chakraborty, R., and Jin, L. (1993). Determination of relatedness between individuals by DNA fingerprinting. Hum. Biol. 65, 875–895.

PubMed Abstract | Google Scholar

Comai, L. (2005). The advantages and disadvantages of being polyploidy. Nat. Rev. Genet. 6, 836–846. doi: 10.1038/nrg1711

PubMed Abstract | CrossRef Full Text | Google Scholar

Contandriopoulos, J. (1966). Contribution à l’étude caryologique des Gésneriacées d’Europe et de leur germination. 91 Congrès des Sociétés Savantes. Rennes 3, 271–280.

Google Scholar

Daskalova, E., Dontcheva, S., Zekaj, Z., Bacu, A., Sota, V., Abdullai, K., et al. (2012). Initial determination of polymorphism and in vitro conservation of some Ramonda serbica and Ramonda nathaliae populations from Albania, Macedonia and Bulgaria. Biotechnol. Biotec. Eq. 26, 16–25. doi: 10.5504/50YRTIMB.2011.0004

CrossRef Full Text | Google Scholar

De Storme, N., and Mason, A. (2014). Plant speciation through chromosome instability and ploidy change: cellular mechanisms, molecular factors and evolutionary relevance. Curr. Plant Biol. 1, 10–33. doi: 10.1016/j.cpb.2014.09.002

CrossRef Full Text | Google Scholar

Decanter, L., Colling, G., Elvinger, N., Heiðmarsson, S., and Matthies, D. (2020). Ecological niche differences between two polyploid cytotypes of Saxifraga rosacea. Am. J. Bot. 107, 423–435. doi: 10.1002/ajb2.1431

PubMed Abstract | CrossRef Full Text | Google Scholar

Doležel, J., Bartoš, J., Voglmayr, H., and Greilhuber, J. (2003). Nuclear DNA content and genome size of trout and human. Cytometry 51A, 127–128. doi: 10.1002/cyto.a.10013

CrossRef Full Text | Google Scholar

Dubreuil, M., Riba, M., and Mayol, M. (2008). Genetic structure and diversity in Ramonda myconi (Gesneriaceae): effects of historical climate change on a preglacial relict species. Am. J. Bot. 95, 577–587. doi: 10.3732/ajb.2007320

PubMed Abstract | CrossRef Full Text | Google Scholar

Earl, D. A., and vonHoldt, B. M. (2012). STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361. doi: 10.1007/s12686-011-9548-7

CrossRef Full Text | Google Scholar

Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., and Lischer, H. E. L. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Feulgen, R., and Rossenbeck, H. (1924). Mikroskopisch-chemischer Nachweis einer Nukleinsäure vom Typus der Thymonukleinsäure und die darauf beruhende elektive Färbung von Zellkernen in mikroskopischen Präparaten. Hoppe-Seyler’s Zeitschr. Physiol. Chem. 135, 203–248. doi: 10.1515/bchm2.1924.135.5-6.203

CrossRef Full Text | Google Scholar

Frajman, B., Rešetnik, I., Weiss-Schneeweiss, H., Ehrendorfer, F., and Schönswetter, P. (2015). Cytotype diversity and genome size variation in Knautia (Caprifoliaceae, Dipscoideae). BMC Evol. Biol. 15:140. doi: 10.1186/s12862-015-0425-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Glišić, L. J. (1924). Development of the female x-generation and embryo in Ramonda. dissertation thesis. Belgrade (Yugoslavia): University of Belgrade.

Google Scholar

Glombik, M., Bačovský, V., Hobza, R., and Kopecký, D. (2020). Competition of parental genomes in plant hybrids. Front. Plant Sci. 11:200. doi: 10.3389/fpls.2020.00200

PubMed Abstract | CrossRef Full Text | Google Scholar

Hajrudinović, A., Frajman, B., Schönswetter, P., Silajdžić, E., Siljak-Yakovlev, S., and Bogunić, F. (2015). Towards a better understanding of polyploidy Sorbus (Rosaceae) from Bosnia and Herzegovina (Balkan Peninsula), including description of a novel, tetraploid apomictic species. Bot. J. Linn. Soc. 178, 670–685. doi: 10.1111/boj.12289

CrossRef Full Text | Google Scholar

Hewitt, G. M. (1996). Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linn. Soc. 58, 247–276. doi: 10.1006/bijl.1996.0035

CrossRef Full Text | Google Scholar

Hewitt, G. M. (1999). Post-glacial re-colonization of European biota. Biol. J. Linn. Soc. 68, 87–112. doi: 10.1111/j.1095-8312.1999.tb01160.x

CrossRef Full Text | Google Scholar

Hughes, P. D., Woodward, J. C., and Gibbard, P. L. (2006). Quaternary glacial history of the Mediterranean mountains. Prog. Phys. Geography 30, 334–364. doi: 10.1191/0309133306pp481ra

CrossRef Full Text | Google Scholar

Husband, B. C. (2004). The role of triploid hybrids in the evolutionary dynamics of mixed-ploidy populations. Biol. J. Linn. Soc. 82, 537–546. doi: 10.1111/j.1095-8312.2004.00339.x

CrossRef Full Text | Google Scholar

Jacquemyn, H., van der Meer, S., and Brys, R. (2016). The impact of hybridization on long-term persistence of polyploid Dactylorhiza species. Am. J. Bot. 103, 1829–1837. doi: 10.3732/ajb.1600274

PubMed Abstract | CrossRef Full Text | Google Scholar

Jakobsson, M., and Rosenberg, N. A. (2007). CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23, 1801–1806. doi: 10.1093/bioinformatics/btm233

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiao, Y., Wickett, N., Ayyampalayam, S., Chanderbali, A. S., Landherr, L., Ralph, P. E., et al. (2011). Ancestral polyploidy in seed plants and angiosperms. Nature 473, 97–100. doi: 10.1038/nature09916

PubMed Abstract | CrossRef Full Text | Google Scholar

Jombart, T., Devillard, S., and Balloux, F. (2010). Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 11:94. doi: 10.1186/1471-2156-11-94

PubMed Abstract | CrossRef Full Text | Google Scholar

Kolář, F., Čertner, M., Suda, J., Schönswetter, P., and Husband, B. C. (2017). Mixed-ploidy species: progress and opportunities in polyploidy research. Trends Plant Sci. 22, 1041–1055. doi: 10.1016/j.tplants.2017.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Košanin, N. (1921). La distribution geographique des deux especes de Ramondia du Balkan (in French). Belgrade: Academie des Sciences et Arts

Google Scholar

Kougioumoutzis, K., Kokkoris, I. P., Panitsa, M., Kallimanis, A., Strid, A., and Dimopoulos, P. (2021). Plant endemism centres and biodiversity hotspots in Greece. Biology 10:72. doi: 10.3390/biology10020072

PubMed Abstract | CrossRef Full Text | Google Scholar

Kreiner, J. M., Kron, P., and Husband, B. C. (2017). Evolutionary dynamics of unreduced gametes. Trends Genet. 33, 583–593. doi: 10.1016/j.tig.2017.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Krstić, N., Savić, L., and Jovanović, G. (2012). The Neogene lakes on the Balkan land. Annales Géologiques de la Péninsule Balkanique 73, 37–60. doi: 10.2298/GABP1273037K

CrossRef Full Text | Google Scholar

Kuzmanović, N., Comanescu, P., Frajman, B., Lazarević, M., Paun, O., Schönswetter, P., et al. (2013). Genetic, cytological and morphological differentiation within the Balkan – Carpathian Sesleria rigida sensu Fl. Eur. (Poaceae): a taxonomical intricate tetraploid – octoploid complex. Taxon 62, 458–472. doi: 10.12705/623.13

CrossRef Full Text | Google Scholar

Lakušić, D., Rakić, T., Stefanović, S., Surina, B., and Stevanović, V. (2009). Edraianthus × lakusicii (Campanulaceae) a new intersectional natural hybrid: morphological and molecular evidence. Plant Syst. Evol. 280, 77–88. doi: 10.1007/s00606-009-0168-6

CrossRef Full Text | Google Scholar

Langella, O. (1999). POPULATIONS version 1.2.30: population genetic software. CNRS UPR9034.

Google Scholar

Lazarević, M. (2012). Cytogenetics, palynology and phylogeography of genus Ramonda (Gesneriaceae) in the Balkan Peninsula. dissertation thesis. [Belgrade (RS)]: University of Belgrade.

Google Scholar

Lazarević, M., Kuzmanović, N., Lakušić, D., Alegro, A., Schönswetter, P., and Frajman, B. (2015). Patterns of cytotype distribution and genome size variation in the genus Sesleria Scop. (Poaceae). Bot. J. Linn. Soc. 179, 126–143. doi: 10.1111/boj.12306

CrossRef Full Text | Google Scholar

Lazarević, M., Rakić, T., and Šinžar-Sekulić, J. (2014). Morphological differences between the flowers of Ramonda serbica, R. nathaliae and their hybrid. Bot. Serb. 38, 91–98.

Google Scholar

Lazarević, M., Siljak-Yakovlev, S., Lazarević, P., Stevanović, B., and Stevanović, V. (2013). Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification. Turk. J. Bot. 37, 872–885. doi: 10.3906/bot-1209-58

CrossRef Full Text | Google Scholar

Levin, D. A. (2002). The Role of Chromosomal Change in Plant Evolution. New York, NY: Oxford University Press.

Google Scholar

Lihová, J., Kučera, J., Perný, M., and Marhold, K. (2007). Hybridization between two polyploid Cardamine (Brassicaceae) species in North-Western Spain: discordance between morphological and genetic variation patterns. Ann. Bot. London 99, 1083–1096. doi: 10.1093/aob/mcm056

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Liu, Y., Yang, X., Tong, C., Edwards, D., Parkin, I. A. P., et al. (2014). The Brassica oleracea genome reveals the asymmetrical evolution of polyploid genomes. Nat. Commun. 5:3930. doi: 10.1038/ncomms4930

PubMed Abstract | CrossRef Full Text | Google Scholar

Maherali, H., Walden, A. E., and Husband, B. C. (2009). Genome duplication and the evolution of physiological responses to water stress. New Phytol. 184, 721–731. doi: 10.1111/j.1469-8137.2009.02997.x

CrossRef Full Text | Google Scholar

Mallet, J. (2007). Hybrid speciation. Nature 446, 279–283. doi: 10.1038/nature05706

CrossRef Full Text | Google Scholar

Mallet, J. (2008). Hybridization, ecological arces and the nature of species: empirical evidence for the ease of speciation. Phil. Trans. R. Soc. 363, 2971–2986. doi: 10.1098/rstb.2008.0081

PubMed Abstract | CrossRef Full Text | Google Scholar

Marburger, S., Monnahan, P., Seear, P. J., Martin, S. H., Koch, J., Paajanen, P., et al. (2019). Interspecific introgression mediates adaptation to whole genome duplication. Nat. Commun. 10:5218. doi: 10.1038/s41467-019-13159-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Marie, D., and Brown, S. (1993). A cytometric exercise in plant DNA histograms, with 2C values for 70 species. Biol. Cell. 78, 41–51. doi: 10.1016/0248-4900(93)90113-S

PubMed Abstract | CrossRef Full Text | Google Scholar

Marjanac, L., and Marjanac, T. (2004). “Glacial history of the Croatian Adriatic and coastal Dinarides” in Quaternary Glaciations — Extent and Chronology. Part I: Europe. eds. J. Ehlers and P. L. Gibbard (Amsterdam: Elsevier), 19–26.

Google Scholar

Mason, A. S., and Pires, J. C. (2015). Unreduced gametes: meiotic mishap or evolutionary mechanism? Trends Genet. 31, 5–10. doi: 10.1016/j.tig.2014.09.011

CrossRef Full Text | Google Scholar

Mason, A. S., and Wendel, J. F. (2020). Homoeologous exchanges, segmental allopolyploidy, and polyploid genome evolution. Front. Genet. 11:1014. doi: 10.3389/fgene.2020.01014

PubMed Abstract | CrossRef Full Text | Google Scholar

Menkovic, L., Markovic, M., Cupkovic, T., Pavlovic, R., Trivic, B., and Banjanac, N. (2004). “Glacial morphology of Serbia Yugoslavia, with comments on the Pleistocene glaciation of Montenegro, Macedonia and Albania” in Quaternary Glaciations – Extent and Chronology. Part I: Europe. eds. J. Ehlers and P. L. Gibbard (Amsterdam: Elsevier), 379–384.

Google Scholar

Myers, N., Mittermeier, R. A., Mittermeier, C. G., da Fonseca, G. A. B., and Kent, J. (2000). Biodiversity hotspots for conservation priorities. Nature 403, 853–858. doi: 10.1038/35002501

CrossRef Full Text | Google Scholar

Nieto Feliner, G., Casacuberta, J., and Wendel, J. F. (2020). Genomics of evolutionary novelty in hybrids and polyploids. Front. Genet. 11:792. doi: 10.3389/fgene.2020.00792

PubMed Abstract | CrossRef Full Text | Google Scholar

Niketić, M., Đurović, S. Z., Tomović, G., Schönswetter, P., and Frajman, B. (2021). Diversification within ploidy-variable Balkan endemic Cerastium decalvans (Caryophyllaceae) reconstructed based on genetic, morphological and ecological evidence. Bot. J. Linn. Soc. :boab037. doi: 10.1093/botlinnean/boab037

CrossRef Full Text | Google Scholar

Niketić, M., Siljak-Yakovlev, S., Frajman, B., Lazarević, M., Stevanović, B., Tomović, G., et al. (2013). Towards resolving the systematics of Cerastium subsection Cerastium (Caryophyllaceae): a cytogenetic approach. Bot. J. Linn. Soc. 172, 205–224. doi: 10.1111/boj.12050

CrossRef Full Text | Google Scholar

Nunvářová Kabátová, K., Kolář, F., Jarolímová, V., Krak, K., and Chrtek, J. (2019). Does geography, evolutionary history or ecology drive ploidy and genome size variation in the Minuartia verna group (Caryophyllaceae) across Europe? Plant Syst. Evol. 305, 1019–1040. doi: 10.1007/s00606-019-01621-2

CrossRef Full Text | Google Scholar

Otto, S. P., and Whitton, J. (2000). Polyploid incidence and evolution. Annu. Rev. Genet. 34, 401–437. doi: 10.1146/annurev.genet.34.1.401

CrossRef Full Text | Google Scholar

Padilla-García, N., Rojas-Andrés, B. M., López-González, N., Castro, M., Castro, S., Loureiro, J., et al. (2018). The challenge of species delimitation in the diploid-polyploid complex Veronica subsection Pentasepalae. Mol. Phylogenet. Evol. 119, 196–209. doi: 10.1016/j.ympev.2017.11.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Pančić, J. (1874). Flora Kneževine Srbije. Belgrade: Državna štamparija.

Google Scholar

Pančić, J. (1884). Dodatak flori Kneževine Srbije. Belgrade: Kraljevsko-srpska državna štamparija.

Google Scholar

Peskoller, A., Silbernagl, L., Hülber, K., Sonnleitner, M., and Schönswetter, P. (2021). Do pentaploid hybrids mediate gene flow between tetraploi Senecio disjunctus and hexaploid S. carniolicus s. str. (S. carniolicus aggregate, Asteraceae). Alpine Bot. 131, 151–160. doi: 10.1007/s00035-021-00254-x

CrossRef Full Text | Google Scholar

Petrova, G., Dzhambazova, T., Moyankova, D., Georgieva, D., Michova, A., Djilianov, D., et al. (2014). Morphological variation, genetic diversity and genome size of critically endangered Haberlea (Gesneriaceae) populations in Bulgaria do not support the recognition of two different species. Plant Evol. Syst. 300, 29–41. doi: 10.1007/s00606-013-0857-z

CrossRef Full Text | Google Scholar

Petrova, G., Moyankova, D., Nishii, K., Forrest, L., Tsiripidis, I., Drouzas, A. D., et al. (2015). The European paleoendemic Haberlea rhodopensis (Gesneriaceae) has an Oligocene origin and a Pleistocene diversification and occurs in a long-persisting refugial area in southeastern Europe. Int. J. Plant Sci. 176, 499–514. doi: 10.1086/681990

CrossRef Full Text | Google Scholar

Petrova, G., Petrov, S., and Möller, M. (2018). Low genetic diversity in small leading adge populations of the European paleoendemic Ramonda serbica (Gesneriacae) in Bulgaria. Nord. J. Bot. 36:01655. doi: 10.1111/njb.01655

CrossRef Full Text | Google Scholar

Petrović, S. (1882). Flora okoline Niša. Belgrade: Kraljevsko-srpska državna štamparija.

Google Scholar

Picó, F. X., and Riba, M. (2002). Regional-scale demography of Ramonda myconi: remnant population dynamics in a preglacial relict species. Plant Ecol. 161, 1–13. doi: 10.1023/A:1020310609348

CrossRef Full Text | Google Scholar

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. doi: 10.1093/genetics/155.2.945

CrossRef Full Text | Google Scholar

Raabová, J., Fischer, M., and Münzbergová, Z. (2008). Niche differentiation between diploid and hexaploid Aster amellus. Oecologia 158, 463–472. doi: 10.1007/s00442-008-1156-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Rakić, T., Ilijević, K., Lazarević, M., Gržetić, I., Stevanović, V., and Stevanović, B. (2013). The resurrection flowering plant Ramonda nathaliae on serpentine soil – coping with extreme mineral element stress. Flora 208, 618–625. doi: 10.1016/j.flora.2013.09.006

CrossRef Full Text | Google Scholar

Rakić, T., Lazarević, M., Jovanović, Ž. S., Radović, S., Siljak-Yakovlev, S., Stevanović, B., et al. (2014). Resurrection plants of the genus Ramonda: prospective survival strategies – unlock further capacity of adaptation, or embark on the path of evolution? Front. Plant Sci. 4:550. doi: 10.3389/fpls.2013.00550

CrossRef Full Text | Google Scholar

Ramasamy, R. K., Ramasamy, S., Bindroo, B. B., and Naik, V. G. (2014). STRUCTURE PLOT: a program for drawing elegant STRUCTURE bar plots in user friendly interface. Springerplus 3:431. doi: 10.1186/2193-1801-3-431

PubMed Abstract | CrossRef Full Text | Google Scholar

Ratter, J. A. (1963). Some chromosome numbers in the Gesneriaceae. Notes Roy. Bot. Gard. Edinburgh 24, 221–229.

Google Scholar

Rey, P. J., Manzaneda, A. J., and Alcántara, J. M. (2017). The interplay between aridity and competition determines colonization ability, exclusion and ecological segregation in the heteroploid Brachypodium distachyon species complex. New Phytol. 215, 85–96. doi: 10.1111/nph.14574

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmickl, R., and Yant, L. (2021). Adaptive introgression: how polyploidy reshapes gene flow landscapes. New Phytol. 230, 457–461. doi: 10.1111/nph.17204

PubMed Abstract | CrossRef Full Text | Google Scholar

Siljak-Yakovlev, S., Stevanovic, V., Tomasevic, M., Brown, S. C., and Stevanovic, B. (2008). Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils. Environ. Exp. Bot. 62, 101–112. doi: 10.1016/j.envexpbot.2007.07.017

CrossRef Full Text | Google Scholar

Šlenker, M., Kantor, A., Marhold, K., Schmickl, R., Mandáková, T., Lysak, M. A., et al. (2021). Allele sorting as a novel approach to resolving the origin of allotetraploids using Hyb-Seq data: a case study of the Balkan mountain endemic Cardamine barbaraeoides. Front. Plant Sci. 12:659275. doi: 10.3389/fpls.2021.659275

PubMed Abstract | CrossRef Full Text | Google Scholar

Soltis, D. E., Albert, V. A., Leebens-Mack, J., Bell, C. D., Paterson, A. H., Zheng, C., et al. (2009). Polyploidy and angiosperm diversification. Am. J. Bot. 96, 336–348. doi: 10.3732/ajb.0800079

CrossRef Full Text | Google Scholar

Sonnleitner, M., Flatscher, R., García, P. E., Rauchová, J., Suda, J., Schneeweiss, G. M., et al. (2010). Distribution and habitat segregation on different spatial scales among diploid, tetraploid and hexaploid cytotypes of Senecio carniolicus (Asteraceae) in the eastern Alps. Ann. Bot. London 106, 967–977. doi: 10.1093/aob/mcq192

PubMed Abstract | CrossRef Full Text | Google Scholar

Sonnleitner, M., Weis, B., Flatscher, R., García, P. E., Suda, J., Krejčiková, J., et al. (2013). Parental ploidy strongly affects offspring fitness in heteroploid crosses among three cytotypes of autopolyploid Jacobaea carniolica (Asteraceae). PLoS One 8:e78959. doi: 10.1371/journal.pone.0078959

PubMed Abstract | CrossRef Full Text | Google Scholar

Španiel, S., Marhold, K., and Zozomová-Lihová, J. (2017). The polyploid Alyssum montanum – A. repens complex in the Balkans: a hotspot of species and genetic diversity. Plant Syst. Evol. 303, 1443–1465. doi: 10.1007/s00606-017-1470-3

CrossRef Full Text | Google Scholar

Stevanović, V., Jovanović, S., Lakušić, D., and Niketić, M. (1999). “Characteristics of the flora of Serbia and its phytogeographical division” in Red Data Book of Flora of Serbia 1 – Extinct and Critically Endangered Taxa. ed. V. Stevanović (Belgrade: Ministry of Environment of the Republic of Serbia; Faculty of Biology, University of Belgrade; Institute for Protection of Nature of the Republic of Serbia), 393–399.

Google Scholar

Stevanović, V., Niketić, M., and Stevanović, B. (1986). Sympatric area of the sibling and endemo-relict species Ramonda serbica Panč. and R. nathaliae Panč. Et Petrov. (Gesneriaceae) in Southeast Serbia (Yugoslavia). Glasn. Inst. Bot. Bašte Univ. Beograd 20, 45–55.

Google Scholar

Stevanović, V., and Stevanović, B. (1985). Asplenio cuneifolii-Ramondaetum nathaliae - new chasmophytic community on serpentine rocks in Macedonia. Bull. Nat. Hist. Mus. Belgrade Ser. B 40, 75–87.

Google Scholar

Stevanović, V., Tan, K., and Petrova, A. (2007). Mapping the endemic flora of the Balkans – a progress report. Bocconea 21, 131–137.

Google Scholar

Stutz, S., Hinz, H. L., Konowalik, K., Müller-Schärer, H., Oberprieler, C., and Schaffner, U. (2016). Ploidy level in the genus Leucanthemum correlates with resistance to a specialst herbivore. Ecosphere 7:e01460. doi: 10.1002/ecs2.1460

CrossRef Full Text | Google Scholar

Szadkowski, E., Eber, F., Huteau, V., Lode, M., Coriton, O., Jenczewski, E., et al. (2011). Polyploid formation pathways have an impact on genetic rearrangements in resynthesized Brassica napus. New Phytol. 191, 884–894. doi: 10.1111/j.1469-8137.2011.03729.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, S. J., Arnold, M., and Martin, N. H. (2009). The genetic architecture of reproductive isolation in Louisiana irises: hybrid fitness in nature. Evolution 63, 2581–2594. doi: 10.1111/j.1558-5646.2009.00742.xte

CrossRef Full Text | Google Scholar

te Beest, M., Le Roux, J. J., Richardson, D. M., Brysting, A. K., Suda, J., Kubešová, M., et al. (2012). The more the better? The role of polyploidy in facilitating plant invasions. Ann. Bot. London 109, 19–45. doi: 10.1093/aob/mcr277

CrossRef Full Text | Google Scholar

Tel-Zur, N., Mouyal, J., Zurgil, U., and Mizrahi, Y. (2020). In support of Winge’s theory of hybridization followed by chromosome doubling. Front. Plant Sci. 11:954. doi: 10.3389/fpls.2020.00954

PubMed Abstract | CrossRef Full Text | Google Scholar

Tendal, K., Ørgaard, M., Larsen, B., and Pedersen, C. (2018). Recurrent hybridisation events between Primula vulgaris, P. veris and P. elatior (Primulaceae, Ericales) challenges the species boundaries: using molecular markers to re-evaluate morphological identifications. Nord. J. Bot. 36:e01778. doi: 10.1111/njb.01778

CrossRef Full Text | Google Scholar

Thompson, J. D. (2020). Plant Evolution in the Mediterranean: Insights for Conservation. New York, NY: Oxford University Press.

Google Scholar

Turill, W. B. (1929). The Plant-Life of the Balkan Peninsula. A Phytogeographical Study. Oxford: Clarendon Press.

Google Scholar

Ulum, F. B., Hadacek, F., and Hörandl, E. (2021). Polyploidy improves photosynthesis regulation within the Ranunculus auricomus complex (Ranunculaceae). Biology 10:811. doi: 10.3390/biology10080811

PubMed Abstract | CrossRef Full Text | Google Scholar

Vallejo-Marín, M., and Hiscock, S. J. (2016). Hybridization and hybrid speciation under global change. New Phytol. 211, 1170–1187. doi: 10.1111/nph.14004

PubMed Abstract | CrossRef Full Text | Google Scholar

Van de Peer, Y., Ashman, T.-L., Soltis, P. S., and Soltis, D. E. (2021). Polyploidy: an evolutionary and ecological force in stressful times. Plant Cell 33, 11–26. doi: 10.1093/plcell/koaa015

PubMed Abstract | CrossRef Full Text | Google Scholar

Van de Peer, Y., Mizrachi, E., and Marchal, K. (2017). The evolutionary significance of polyploidy. Nat. Rev. Genet. 18, 411–424. doi: 10.1038/nrg.2017.26

CrossRef Full Text | Google Scholar

Vekemans, X. (2002). AFLP-SURV Version 1.0. Distributed by the author. Laboratoire de Génétique et Ecologie Végétale, Université Libre de Bruxelles, Belgium.

Google Scholar

Vos, P., Hogers, R., Bleeker, M., Reijans, M., Lee, T. V. D., Hornes, M., et al. (1995). AFLP: a new technique for DNA fingerprinting. Nucleic Acids Res. 23, 4407–4414. doi: 10.1093/nar/23.21.4407

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss-Schneeweiss, H., Emadzade, K., Jang, T.-S., and Schneeweiss, G. M. (2013). Evolutionary consequences, constraints and potential of polyploidy in plants. Cytogenet. Genome Res. 140, 137–150. doi: 10.1159/000351727

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendel, J. F. (2000). Genome evolution in polyploids. Plant Mol. Biol. 42, 225–249. doi: 10.1023/A:1006392424384

CrossRef Full Text | Google Scholar

Zhivotovsky, L. A. (1999). Estimating population structure in diploids with multilocus dominant DNA markers. Mol. Ecol. 8, 907–913. doi: 10.1046/j.1365-294x.1999.00620.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Živković, T., Quartacci, M. F., Stevanović, B., Marinone, F., and Navari-Izzo, F. (2005). Low-molecular weight substances in the poikilohydric plant Ramonda serbica during dehydration and rehydration. Plant Sci. 168, 105–111. doi: 10.1016/j.plantsci.2004.07.018

CrossRef Full Text | Google Scholar