Maja Lazarević

Sonja Siljak-Yakovlev
email: sonia.yakovlev@universite-paris-saclay.fr

Agathe Sanino

Marjan Niketić

Françoise Lamy

Damien D Hinsinger

Gordana Tomović

Branka Stevanović

Vladimir Stevanović

Thierry Robert
email: thierry.robert@universite-paris-saclay.fr

Real Gonzalo Nieto Feliner

Marcial Escudero

Eva Maria Temsch

Maja Lazarević
email: majalazarevićmajat@bio.bg.ac.rs

Marjan Niketić

Gordana Tomović

Branka Stevanović

Vladimir Stevanović

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

Keywords: AFLP, gene flow, genetic diversity, genome size, Gesneriaceae, interploidy hybridization, mixed-ploidy zones, sympatry

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Thompson | Plant Evolution in the Mediterranean: Insights for Conservation[END_REF]. The Balkan Peninsula, as part of the Mediterranean Basin, one of the global hotspots of biodiversity [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF], 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 [START_REF] Turill | The Plant-Life of the Balkan Peninsula. A Phytogeographical Study[END_REF][START_REF] Stevanović | Mapping the endemic flora of the Balkans -a progress report[END_REF]. Several geological, climatic, and biogeographical factors have been identified that have contributed to the richness and endemism of the flora of this region [START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF][START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF][START_REF] Stevanović | Characteristics of the flora of Serbia and its phytogeographical division" in Red Data Book of Flora of Serbia 1 -Extinct and Critically Endangered Taxa[END_REF][START_REF] Stevanović | Mapping the endemic flora of the Balkans -a progress report[END_REF][START_REF] Krstić | The Neogene lakes on the Balkan land[END_REF]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 [START_REF] Banjanac | Chemodiversity of two closely related tetraploid Centaurium species and their hexaploid hybrid: metabolomic search for high-resolution taxonomic classifiers[END_REF], Cerastium [START_REF] Niketić | Towards resolving the systematics of Cerastium subsection Cerastium (Caryophyllaceae): a cytogenetic approach[END_REF][START_REF] Niketić | Diversification within ploidy-variable Balkan endemic Cerastium decalvans (Caryophyllaceae) reconstructed based on genetic, morphological and ecological evidence[END_REF], Cotoneaster [START_REF] Bogunić | Article 873471 reproductive mode variation of Cotoneaster integerrimus (Rosaceae) from the Balkans. Plan[END_REF], Edraianthus [START_REF] Lakušić | Edraianthus × lakusicii (Campanulaceae) a new intersectional natural hybrid: morphological and molecular evidence[END_REF], Knautia [START_REF] Frajman | Cytotype diversity and genome size variation in Knautia (Caprifoliaceae, Dipscoideae)[END_REF], Luzula [START_REF] Bačič | Diversification of Luzula sect. Luzula (Juncaceae) on the Balkan Peninsula: a cytogenetic approach[END_REF], Sesleria [START_REF] Kuzmanović | Genetic, cytological and morphological differentiation within the Balkan -Carpathian Sesleria rigida sensu Fl. Eur. (Poaceae): a taxonomical intricate tetraploid -octoploid complex[END_REF][START_REF] Lazarević | Patterns of cytotype distribution and genome size variation in the genus Sesleria Scop. (Poaceae)[END_REF], Sorbus [START_REF] Hajrudinović | Towards a better understanding of polyploidy Sorbus (Rosaceae) from Bosnia and Herzegovina (Balkan Peninsula), including description of a novel, tetraploid apomictic species[END_REF], Veronica [START_REF] Padilla-García | The challenge of species delimitation in the diploid-polyploid complex Veronica subsection Pentasepalae[END_REF].

Polyploidization and hybridization are considered important forces in plant evolution [START_REF] Wendel | Genome evolution in polyploids[END_REF][START_REF] Mallet | Hybrid speciation[END_REF][START_REF] Soltis | Polyploidy and angiosperm diversification[END_REF][START_REF] Jiao | Ancestral polyploidy in seed plants and angiosperms[END_REF][START_REF] Van De Peer | Polyploidy: an evolutionary and ecological force in stressful times[END_REF]. Polyploids may differ from their diploid relatives in their morphology, phenology, physiology, biotic interactions, distribution, preferences for ecological niches, and invasiveness (e.g., [START_REF] Levin | The Role of Chromosomal Change in Plant Evolution[END_REF][START_REF] Raabová | Niche differentiation between diploid and hexaploid Aster amellus[END_REF][START_REF] Maherali | Genome duplication and the evolution of physiological responses to water stress[END_REF][START_REF] Xte Te Beest | In support of Winge's theory of hybridization followed by chromosome doubling[END_REF][START_REF] Weiss-Schneeweiss | Evolutionary consequences, constraints and potential of polyploidy in plants[END_REF][START_REF] Stutz | Ploidy level in the genus Leucanthemum correlates with resistance to a specialst herbivore[END_REF][START_REF] Decanter | Ecological niche differences between two polyploid cytotypes of Saxifraga rosacea[END_REF][START_REF] Ulum | Polyploidy improves photosynthesis regulation within the Ranunculus auricomus complex (Ranunculaceae)[END_REF]. Many studies suggest better fitness of polyploids [START_REF] Comai | The advantages and disadvantages of being polyploidy[END_REF][START_REF] Maherali | Genome duplication and the evolution of physiological responses to water stress[END_REF][START_REF] Xte Te Beest | In support of Winge's theory of hybridization followed by chromosome doubling[END_REF][START_REF] Rey | The interplay between aridity and competition determines colonization ability, exclusion and ecological segregation in the heteroploid Brachypodium distachyon species complex[END_REF][START_REF] Van De Peer | Polyploidy: an evolutionary and ecological force in stressful times[END_REF], but opposite results can also be found in the literature [START_REF] Nunvářová Kabátová | Does geography, evolutionary history or ecology drive ploidy and genome size variation in the Minuartia verna group (Caryophyllaceae) across Europe?[END_REF]. 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 [START_REF] Otto | Polyploid incidence and evolution[END_REF][START_REF] Comai | The advantages and disadvantages of being polyploidy[END_REF][START_REF] Weiss-Schneeweiss | Evolutionary consequences, constraints and potential of polyploidy in plants[END_REF][START_REF] Van De Peer | The evolutionary significance of polyploidy[END_REF][START_REF] Košanin | La distribution geographique des deux especes de Ramondia du Balkan (in French)[END_REF][START_REF] Decanter | Ecological niche differences between two polyploid cytotypes of Saxifraga rosacea[END_REF]. 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 [START_REF] Vallejo-Marín | Hybridization and hybrid speciation under global change[END_REF].

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 [START_REF] Košanin | La distribution geographique des deux especes de Ramondia du Balkan (in French)[END_REF][START_REF] Turill | The Plant-Life of the Balkan Peninsula. A Phytogeographical Study[END_REF][START_REF] Petrova | 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[END_REF], 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 [START_REF] Živković | Low-molecular weight substances in the poikilohydric plant Ramonda serbica during dehydration and rehydration[END_REF][START_REF] Rakić | Resurrection plants of the genus Ramonda: prospective survival strategies -unlock further capacity of adaptation, or embark on the path of evolution? Front[END_REF], 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 [START_REF] Košanin | La distribution geographique des deux especes de Ramondia du Balkan (in French)[END_REF][START_REF] Rakić | Resurrection plants of the genus Ramonda: prospective survival strategies -unlock further capacity of adaptation, or embark on the path of evolution? Front[END_REF]. However, R. myconi and R. serbica are found only on limestone, while R. nathaliae also inhabits ultramafic soils [START_REF] Košanin | La distribution geographique des deux especes de Ramondia du Balkan (in French)[END_REF][START_REF] Stevanović | Asplenio cuneifolii-Ramondaetum nathaliae -new chasmophytic community on serpentine rocks in Macedonia[END_REF][START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF][START_REF] Rakić | The resurrection flowering plant Ramonda nathaliae on serpentine soil -coping with extreme mineral element stress[END_REF]. 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 [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF]. Relatively high basic chromosome number (x = 24) indicates that the three species may be paleopolyploids [START_REF] Rakić | Resurrection plants of the genus Ramonda: prospective survival strategies -unlock further capacity of adaptation, or embark on the path of evolution? Front[END_REF]. 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 [START_REF] Rakić | Resurrection plants of the genus Ramonda: prospective survival strategies -unlock further capacity of adaptation, or embark on the path of evolution? Front[END_REF]. 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 [START_REF] Stevanović | Sympatric area of the sibling and endemo-relict species Ramonda serbica Panč. and R. nathaliae Panč. Et Petrov. (Gesneriaceae) in Southeast Serbia (Yugoslavia)[END_REF]. 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 [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF]. 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 [START_REF] Pančić | Flora Kneževine Srbije[END_REF][START_REF] Pančić | Dodatak flori Kneževine Srbije[END_REF][START_REF] Petrović | Kraljevsko-srpska državna štamparija[END_REF], 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 [START_REF] Lazarević | Morphological differences between the flowers of Ramonda serbica, R. nathaliae and their hybrid[END_REF]. 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 [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF]. 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).

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 [START_REF] Feulgen | Mikroskopisch-chemischer Nachweis einer Nukleinsäure vom Typus der Thymonukleinsäure und die darauf beruhende elektive Färbung von Zellkernen in mikroskopischen Präparaten[END_REF] 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 [START_REF] Marie | A cytometric exercise in plant DNA histograms, with 2C values for 70 species[END_REF] 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 [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF].

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 [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF].

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 [START_REF] Vos | AFLP: a new technique for DNA fingerprinting[END_REF] 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 [START_REF] Vekemans | AFLP-SURV Version 1.0. Distributed by the author[END_REF] 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 [START_REF] Zhivotovsky | Estimating population structure in diploids with multilocus dominant DNA markers[END_REF] with a non-uniform prior distribution of allele frequencies.

Genetic dissimilarity between individuals and populations was estimated using Shared Allele Distance [START_REF] Chakraborty | Determination of relatedness between individuals by DNA fingerprinting[END_REF] with POPULATIONS 1.2.30 [START_REF] Langella | POPULATIONS version 1.2.30: population genetic software[END_REF]. 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 [START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and windows[END_REF]; 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 modelbased 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 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. 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 [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] implemented in Structure Harvester v.0.6.94 [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF]. Then, 50 new replicates of the MCMC method were performed to determine the best K value. CLUMPP v.1.1.2 [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF] 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 [START_REF] Ramasamy | STRUCTURE PLOT: a program for drawing elegant STRUCTURE bar plots in user friendly interface[END_REF].

In addition, and to evaluate the robustness of the classification used, a Discriminant Analysis of Principal Components (DAPC, [START_REF] Jombart | Discriminant analysis of principal components: a new method for the analysis of genetically structured populations[END_REF] 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.

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. [START_REF] Doležel | Nuclear DNA content and genome size of trout and human[END_REF]. 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).

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 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 firstgeneration 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 [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF][START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF][START_REF] Lazarević | Morphological differences between the flowers of Ramonda serbica, R. nathaliae and their hybrid[END_REF]. 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 largescale 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.

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.

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 [START_REF] Glišić | Development of the female x-generation and embryo in Ramonda[END_REF] found n = 36 to be the most frequent haploid chromosome number in R. nathaliae, later both [START_REF] Ratter | Some chromosome numbers in the Gesneriaceae[END_REF] and [START_REF] Contandriopoulos | Contribution à l' étude caryologique des Gésneriacées d'Europe et de leur germination. 91[END_REF] showed that R. nathaliae is a diploid with 2n = 2x = 48 chromosomes, and this is confirmed by broad and detailed cytogenetic studies by [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF] and [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF]. 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 [START_REF] Glišić | Development of the female x-generation and embryo in Ramonda[END_REF], 2n = 96 in plants from Mt. Timfi-Vikos Gorge in Greece [START_REF] Contandriopoulos | Contribution à l' étude caryologique des Gésneriacées d'Europe et de leur germination. 91[END_REF], 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 [START_REF] Daskalova | Initial determination of polymorphism and in vitro conservation of some Ramonda serbica and Ramonda nathaliae populations from Albania, Macedonia and Bulgaria[END_REF]. 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 [START_REF] Daskalova | Initial determination of polymorphism and in vitro conservation of some Ramonda serbica and Ramonda nathaliae populations from Albania, Macedonia and Bulgaria[END_REF], and only 6x individuals were found. In the last 15 years, [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF] and [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF], 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 [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF]. 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, [START_REF] Daskalova | Initial determination of polymorphism and in vitro conservation of some Ramonda serbica and Ramonda nathaliae populations from Albania, Macedonia and Bulgaria[END_REF] could not clearly distinguish the two species. More recently, [START_REF] Petrova | Low genetic diversity in small leading adge populations of the European paleoendemic Ramonda serbica (Gesneriacae) in Bulgaria[END_REF] 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 [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF][START_REF] Lazarević | Morphological differences between the flowers of Ramonda serbica, R. nathaliae and their hybrid[END_REF]. 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 [START_REF] Hughes | Quaternary glacial history of the Mediterranean mountains[END_REF], raising the possibility of existence of multiple and scattered refugia for plant and animal species survival in the Balkan region [START_REF] Hewitt | Some genetic consequences of ice ages, and their role in divergence and speciation[END_REF]. 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 [START_REF] Menkovic | Glacial morphology of Serbia Yugoslavia, with comments on the Pleistocene glaciation of Montenegro, Macedonia and Albania[END_REF][START_REF] Hughes | Quaternary glacial history of the Mediterranean mountains[END_REF], to 270 m in some coastal parts of the Dinarides (Marjanac and [START_REF] Marie | A cytometric exercise in plant DNA histograms, with 2C values for 70 species[END_REF]. 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 [START_REF] Košanin | La distribution geographique des deux especes de Ramondia du Balkan (in French)[END_REF][START_REF] Stevanović | Sympatric area of the sibling and endemo-relict species Ramonda serbica Panč. and R. nathaliae Panč. Et Petrov. (Gesneriaceae) in Southeast Serbia (Yugoslavia)[END_REF][START_REF] Lazarević | Cytogenetics, palynology and phylogeography of genus Ramonda (Gesneriaceae) in the Balkan Peninsula[END_REF]. Also, the occurrence of a series of lakes during the Miocene and Pliocene, which spread in different parts of the Balkan Peninsula [START_REF] Krstić | The Neogene lakes on the Balkan land[END_REF], 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 [START_REF] Krstić | The Neogene lakes on the Balkan land[END_REF], 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 [START_REF] Dubreuil | Genetic structure and diversity in Ramonda myconi (Gesneriaceae): effects of historical climate change on a preglacial relict species[END_REF]. 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 [START_REF] Petrova | Morphological variation, genetic diversity and genome size of critically endangered Haberlea (Gesneriaceae) populations in Bulgaria do not support the recognition of two different species[END_REF]. Finally, [START_REF] Petrova | Low genetic diversity in small leading adge populations of the European paleoendemic Ramonda serbica (Gesneriacae) in Bulgaria[END_REF] 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 [START_REF] Picó | Regional-scale demography of Ramonda myconi: remnant population dynamics in a preglacial relict species[END_REF][START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF][START_REF] Petrova | Morphological variation, genetic diversity and genome size of critically endangered Haberlea (Gesneriaceae) populations in Bulgaria do not support the recognition of two different species[END_REF]. The overall genetic diversity of the two studied species is low, similar to what was found for R. myconi in the Pyrenees [START_REF] Dubreuil | Genetic structure and diversity in Ramonda myconi (Gesneriaceae): effects of historical climate change on a preglacial relict species[END_REF] and H. rhodopensis in Bulgaria and Greece [START_REF] Petrova | Morphological variation, genetic diversity and genome size of critically endangered Haberlea (Gesneriaceae) populations in Bulgaria do not support the recognition of two different species[END_REF][START_REF] Petrova | 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[END_REF]. 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 [START_REF] Stevanović | Sympatric area of the sibling and endemo-relict species Ramonda serbica Panč. and R. nathaliae Panč. Et Petrov. (Gesneriaceae) in Southeast Serbia (Yugoslavia)[END_REF]. 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 [START_REF] Petrova | 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[END_REF] 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 [START_REF] Brownfield | Unreduced gamete formation in plants: mechanisms and prospects[END_REF][START_REF] Mason | Unreduced gametes: meiotic mishap or evolutionary mechanism?[END_REF][START_REF] Kreiner | Evolutionary dynamics of unreduced gametes[END_REF]. 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 [START_REF] Kolář | Mixed-ploidy species: progress and opportunities in polyploidy research[END_REF] 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 [START_REF] Castro | Cytotype distribution at a diploid -hexaploid contact zone in Aster amellus (Asteraceae)[END_REF], as well as niche differentiation and shifts in flowering times between cytotypes [START_REF] Raabová | Niche differentiation between diploid and hexaploid Aster amellus[END_REF]. 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 [START_REF] Sonnleitner | Distribution and habitat segregation on different spatial scales among diploid, tetraploid and hexaploid cytotypes of Senecio carniolicus (Asteraceae) in the eastern Alps[END_REF]. 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 [START_REF] Sonnleitner | Parental ploidy strongly affects offspring fitness in heteroploid crosses among three cytotypes of autopolyploid Jacobaea carniolica (Asteraceae)[END_REF].

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 mixedploidy 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 [START_REF] Glombik | Competition of parental genomes in plant hybrids[END_REF]. 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 [START_REF] Siljak-Yakovlev | Genome size variation and polyploidy in the resurrection plant genus Ramonda: cytogeography of living fossils[END_REF]. 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 [START_REF] Szadkowski | Polyploid formation pathways have an impact on genetic rearrangements in resynthesized Brassica napus[END_REF][START_REF] Liu | The Brassica oleracea genome reveals the asymmetrical evolution of polyploid genomes[END_REF][START_REF] Mason | Homoeologous exchanges, segmental allopolyploidy, and polyploid genome evolution[END_REF]. 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 [START_REF] Kolář | Mixed-ploidy species: progress and opportunities in polyploidy research[END_REF]. 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 [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF]. Our results showed therefore that pre-and/or post-zygotic barriers to hybridization that usually exist between different cytotypes and/or plant species [START_REF] Vallejo-Marín | Hybridization and hybrid speciation under global change[END_REF] are not operative in two sites of sympatry [START_REF] Mallet | Hybrid speciation[END_REF][START_REF] Abdelaziz | Asymmetric reproductive barriers and gene flow promote the rise of a stable hybrid zone in the Mediterranean high mountains[END_REF].

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 [START_REF] Taylor | The genetic architecture of reproductive isolation in Louisiana irises: hybrid fitness in nature[END_REF], limiting therefore genetic introgression into restricted regions of genomes (adaptive introgression, e.g., [START_REF] Marburger | Interspecific introgression mediates adaptation to whole genome duplication[END_REF]. 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 [START_REF] De Storme | Plant speciation through chromosome instability and ploidy change: cellular mechanisms, molecular factors and evolutionary relevance[END_REF]. In addition, mainly because of chromosome pairing during meiosis, triploids are often sterile or at least with greatly reduced fertility [START_REF] Comai | The advantages and disadvantages of being polyploidy[END_REF] 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 [START_REF] Husband | The role of triploid hybrids in the evolutionary dynamics of mixed-ploidy populations[END_REF]. 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;[START_REF] Peskoller | Do pentaploid hybrids mediate gene flow between tetraploi Senecio disjunctus and hexaploid S. carniolicus s. str. (S. carniolicus aggregate, Asteraceae)[END_REF].

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 [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF] 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 [START_REF] Lihová | Hybridization between two polyploid Cardamine (Brassicaceae) species in North-Western Spain: discordance between morphological and genetic variation patterns[END_REF][START_REF] Čertner | Does hybridization with a widespread congener threaten the long-term persistence of the eastern alpine rare local endemic Knautia carinthiaca?[END_REF][START_REF] Tendal | 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[END_REF][START_REF] Beirinckx | Hidden hybridization and habitat differentiation in a Mediterranean macrophyte, the euryhaline genus Ruppia[END_REF], as seems to be the case with both sites of sympatry for Ramonda [START_REF] Lazarević | Morphological differences between the flowers of Ramonda serbica, R. nathaliae and their hybrid[END_REF].

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 [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF]. On the other hand, seeds of tetraploid individuals are drastically smaller than those of parental species and have very low or no germination capacity [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF], which is often the case in interspecific hybrids [START_REF] Mallet | Hybridization, ecological arces and the nature of species: empirical evidence for the ease of speciation[END_REF][START_REF] Schmickl | Adaptive introgression: how polyploidy reshapes gene flow landscapes[END_REF]. Interestingly, with the exception of R. nathaliae in Oblik, seeds of both parental cytotypes have shown lower germination rates compared to monospecific populations [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF], 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 [START_REF] Jacquemyn | The impact of hybridization on long-term persistence of polyploid Dactylorhiza species[END_REF]. 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 [START_REF] Abdelaziz | Asymmetric reproductive barriers and gene flow promote the rise of a stable hybrid zone in the Mediterranean high mountains[END_REF] 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 [START_REF] Lazarević | Pollen and seed morphology of resurrection plants from the genus Ramonda (Gesneriaceae): relationship with ploidy level and relevance to their ecology and identification[END_REF]. 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 [START_REF] Schmickl | Adaptive introgression: how polyploidy reshapes gene flow landscapes[END_REF].

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.

FIGURE 1 |

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

FIGURE 2 |

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

FrontiersFIGURE 4 |FIGURE 5 |

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

Frontiers

  in Plant Science | www.frontiersin.org 13 April 2022 | Volume 13 | Article 873471

Frontiers

  in Plant Science | www.frontiersin.org 14 April 2022 | Volume 13 | Article 873471

Frontiers

  in Plant Science | www.frontiersin.org 15 April 2022 | Volume 13 | Article 873471

Frontiers

  in Plant Science | www.frontiersin.org 16 April 2022 | Volume 13 | Article 873471

TABLE 1 |

 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.

	Pop. code	Species	Locality	Altitude	Latitude	Geological	Voucher	Collectors	2n	2C	AFLP
				(m)	longitude	substrate	number
	N7	R. nathaliae Serbia (Kosovo): Mt.	1,100	N 42°14′6.46″	Limestone	21843	Lazarevic, M.,
			Šara-Gotovuša		E 21°4′30.14″		BEOU	Lazarevic, P.
	N8	R. nathaliae North Macedonia: Kozjak	1,610	N 41°23′56.72″	Limestone	29886	Niketic, M., Tomovic, G.
					E 21°41′13.17″		BEOU
	S11	R. serbica	Albania: Kruje	574	N 41°30′30.06″	Limestone	32469	Lakušic, D.,
					E 19°47′45.13″		BEOU	Kuzmanovic, N.,
								Lazarevic, P., Alegro, A.
	A	R. nathaliae Serbia: Mt. Suva planina	1750	N 43°11′6.8″	Limestone	20019	Lakušić D.			1 (26)
					E 22°10′0.7″		BEOU
	B	R. nathaliae North Macedonia: gorge of	311.5	N 41°57′9.6″	Limestone	31669	Lazarević M., Sanino A.,			1 (26)
			r. Treska, near Matka		E 21°17′51.5″		BEOU	Siljak-Yakovlev S.,
			village					Matevski V., Stevanović B.,
								Stevanović V.
	C	R. nathaliae North Macedonia: gorge of	266.1	N 41°49′18.9″	Serpentine	31670	Idem	1		1 (26)
			r. Pcǐnja, near the road to		E 21°41′4.4″		BEOU
			Veles
	D	R. nathaliae North Macedonia: gorge of	160	N 41°24′32.8″	Limestone	31671	Tomović G., Stevanović B.,			1 (26)
			r. Vardar, near Demir Kapija		E 22°15′42.1″		BEOU	Stevanović V.
	E	R. nathaliae Greece: Mt. Vermion,	ca. 1,000 N 40°35′5.4″	Limestone	31668	Niketić M.	1		1 (26)
			gorge in the vicinity of		E 22°1′28.5″		BEOU
			village Sella, west of
			Noussa
	F	R. serbica	Serbia: Radovanski	564	N 43°17′2.5″	Limestone	31672	Lazarević M., Lazarević P.,	1		1 (30)
			Kamen, above the village		E 22°4′16.4″		BEOU	Tomović G., Niketić M.,
			C ̌ukljenik					Jović D., Petrović B.
	G	R. serbica	Serbia: Mt. Rtanj	985	N 43°45′36.55″	Limestone	31673	Lazarević M., Lazarević P.,			1 (8)
					E 21°55′40.01″		BEOU	Nešić D.
	H	R. serbica	North Macedonia: gorge of	807	N 41°39′57.1″	Limestone	Personal	Stevanović V., Tomović G.			1 (26)
			r. Radika, Baricǩa klisura		E 20°37′3.6″		herbarium
	I	R. serbica	Montenegro: gorge of r.	420	N 42°8′32.5″	Limestone	19673	Stevanović B.,	1		1 (26)
			Kroni e Murici		E 19°13′9.9″		BEOU	Stevanović V.
	J	R. serbica	Montenegro: gorge of r.	250	N 42°9′19.7″	Limestone	9880	Stevanović B.,	1		1 (26)
			Kroni e Besit		E 19°12′48.4″		BEOU	Stevanović V.
	K	R. serbica	North Macedonia: gorge of	620	N 41°22′11.1″	Limestone	2081	Stevanović V., Tomović G.			1 (26)
			r. Crni Drim, near the		E 20°35′57.6″		BEOU
			village Lukovo
	L	R. serbica	Greece: Mt. Timfi, gorge of	ca. 650	N 40° 00′00″	Limestone	31674	Niketić M., Tomović G.	1		1 (26)
			r. Vikos, near the village		E 20° 54′39.6″		BEOU
			Vikos
	M	R. serbica	Serbia: gorge of r.	920	N 42°55′39.6″	Limestone	31675	Lazarević M., Lazarević P.,			1 (26)
			Mojstirska suhovara		E 20°25′58.4″		BEOU	Pavicévić D.
	N	R. serbica	Serbia: gorge of r.	310	N 44°1′38.1″	Limestone	370/90	Lazarević M., Tomović G.,			1 (26)
			Lazareva reka		E 21°57′16.2"			Lakušić D.
	O	R. serbica	Serbia: C ̌iflik	370	N 43°13′3″	Limestone	31676	Lazarević M., Lazarević P.,	1		1 (26)
					E 22°25′6.2″		BEOU	Tomović G., Niketić M.,
								Jović D., Petrović B.
	P	R. serbica	Serbia: gorge of r. Godulja	750	N 42°55′27.4″	Limestone	31677	Lazarević M., Lazarević P.,			1 (5)
					E 20°20′0.4″		BEOU	Pavicévić D.
	N3	R. nathaliae North Macedonia: gorge of	165.8	N 41°41′15.2″		33186	Lazarevic, M., Siljak-			1 (25)
			r. Vardar near Veles		E 21°47′56.2″		BEOU	Yakovlev, S., Sanino, A.,
								Matevski, V., Stevanovic, B.,
								Stevanovic, V.
	N5	R. nathaliae Serbia: slopes of Mt. Suva	346.6	N 43°10′19.3″	Limestone	20638	Lazarevic, M.,	1		1 (17)
			planina near village Divljana		E 22°18′10"		BEOU	Zlatkovic, B.,
			in the vicinity of Bela					Stevanovic, V.
			Palanka

(Continued) Frontiers in Plant Science | www.frontiersin.org 6 April 2022 | Volume 13 | Article 873471

TABLE 1 |

 1 Continued

	Pop. code	Species	Locality	Altitude	Latitude	Geological	Voucher	Collectors	2n	2C	AFLP
				(m)	longitude	substrate	number
	N6	R. nathaliae North Macedonia: gorge of	295.7	N 41°26′14.1″	Limestone	2063	Lazarevic, M., Šiljak-		1	1 (25)
			r. Raec near Kavadarci		E 21°49′6.27″		BEOU	Yakovlev, S., Sanino, A.,
								Matevski, V., Stevanovic,
								B., Stevanovic, V.
	Nrk	R. nathaliae Serbia: pure subpopulation	576.7	N 43°17′1.2″	Limestone	20641	Stevanovic, V.,	1	1	1 (16)
			from Radovanski Kamen		E 22°4′23.5″		BEOU	Lazarevic, M.
	S1	R. serbica	North Macedonia: Galicǐca,	950	N 40°57′02.63″	Limestone	Personal	Stevanovic, B., Tomovic, G.,			1 (24)
			Zli do		E 20°48′02.73″		herbarium	Stevanovic, V.,
	S7	R. serbica	Serbia: gorge of river	295.4	N 43°16′56″	Limestone	20642	Lazarevic, M.,	1	1	1 (24)
			Jelašnica in the vicinity of		E 22°3′42″		BEOU	Stevanovic, V.
			Niš
	S9	R. serbica	Serbia: gorge of r. Nišava,	308.5	N 43°18′38.1″	Limestone	20639	Lazarevic, M., Zlatkovic,	1	1	1 (23)
			20 km east from Niš		E 22°10′53.3″		BEOU	B., Stevanovic, V.
	S10	R. serbica	Serbia: gorge of r. Crna	1,100	N 42°56′46.3″	Limestone	20623	Lazarevic, M.,	1	1	1 (18)
			Reka in the vicinity of		E 20°28′09.7″		BEOU	Lazarevic, P.,
			monastery St. Peter c.					Pavicévic, D.
			8 km southwards from
			village Ribaricí
	Radovanski	Sympatry	Serbia: Radovanski Kamen	574.5	N 43°17′5.1″	Limestone	20640	Lazarevic, M., Zlatkovic, B.,	1	1	1 (55)
	Kamen				E 22°4′16.6″		BEOU	Stevanovic, V.
	Oblik	Sympatry	Serbia: gorge of r. Nišava	418.2	N 43°19′47.4″	Limestone	20636	Idem	1	1	1 (95)
			in the vicinity of village		E 22°8′44.9″		BEOU
			Ostrovica-Oblik

Frontiers in Plant Science | www.frontiersin.org 7 April 2022 | Volume 13 | Article 873471

TABLE 2 |

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

	Pop. Code	Species	Locality	N	2n	Ploidy level (x) 2C in pg and	CV (%)	1Cx in pg and
							(Mbp) a		(Mbp)
	A	R. nathaliae	RS: Mt. Suva planina	5		2	2.28 (2229)*	2.34	1.14 (1115)
	B	R. nathaliae	MK: Matka	5		2	2.26 (2208)	1.50	1.13 (1104)
	C	R. nathaliae	MK: Pcǐnja	5	48**	2	2.27 (2217)**	1.54	1.13 (1108)
	E	R. nathaliae	GR: Mt. Vermion	5	48	2	2.28 (2231)*	1.93	1.14 (1116)
	N3	R. nathaliae	MK: Vardar -Veles	5		2	2.35 (2295)**	0.50	1.17 (1147)
	N5	R. nathaliae	RS: Divljana	5	48*	2	2.33 (2277)*	1.95	1.16 (1138)
	N6	R. nathaliae	MK: Raec -Kavadarci	5		2	2.37 (2318)	0.88	1.18 (1159)
	N7	R. nathaliae	RS (XK): Gotovuša	5		2	2.32 (2266)	2.17	1.16 (1133)
	N8	R. nathaliae	MK: Kozjak	5		2	2.50 (2445)	0.92	1.25 (1222)
		R. nathaliae	Mean value	45		2	2.33 (2278)	3.36	1.16 (1139)
	G	R. serbica	RS: Mt. Rtanj	5		6	7.76 (7587)	1.04	1.29 (1264)
	H	R. serbica	MK: Radika	5		6	7.88 (7706)	0.77	1.31 (1284)
	I	R. serbica	ME: Kroni e Murici	5	~150*	6	7.94 (7763)*	2.38	1.32 (1294)
	J	R. serbica	ME: Kroni e Besit	5	144*	6	7.77 (7601)*	0.50	1.30 (1267)
	K	R. serbica	MK: Crni Drim	5		6	7.95 (7772)	0.32	1.32 (1295)
	L	R. serbica	GR: Mt. Timfi	5	144**	6	7.58 (7416)**	2.29	1.26 (1236)
	M	R. serbica	RS: Mojstir	5		6	7.91 (7737)	2.49	1.32 (1289)
	N	R. serbica	RS: Lazareva reka	5		6	7.58 (7413)*	1.37	1.26 (1236)
	O	R. serbica	RS: C ̌iflik	5	144	6	7.79 (7619)	2.92	1.30 (1270)
	P	R. serbica	RS: Godulja	5		6	7.61 (7446)**	6.28	1.27 (1241)
	S7	R. serbica	RS: Jelašnica	5	144**	6	7.85 (7678)**	0.86	1.31 (1280)
	S9	R. serbica	RS: Nišava	5	> 144*	6	8.11 (7927)*	1.86	1.35 (1321)
	S10	R. serbica	RS: Crna Reka	3	144*	6	8.32 (8136)*	1.82	1.39 (1356)
	S11	R. serbica	AL: Kruje	3	-	6	7.76 (7592)	1.04	1.29 (1265)
		R. serbica	Mean value	66		6	7.83 (7660)	3.15	1.31 (1277)

a Values in Mbp calculated based on pg = 978 Mbp

TABLE 3 |

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

	Locality	Taxon	Ploidy level (x)	2n	2C (pg)	CV (%)	Min-Max of 2C	N
							(pg)
	RS: Oblik	R. nathaliae	2	-	2.34	2.00	2.27-2.45	30
	RS: Radovanski	R. nathaliae	2	48*	2.31	2.81	2.22-2.51	45
	Kamen
	RS: Oblik	R. serbica	6	-	7.95	2.01	7.72-8.29	23
	RS: Radovanski	R. serbica	6	144*	8.00	4.38	7.48-8.45	18
	Kamen
	RS: Oblik	Class I	4	96*	5.14	1.41	5.03-5.27	25
		Class II	-	-	6.12	-	-	1
		Class III	-	-	10.32	12.83	8.94-12.06	4
	RS: Radovanski	Class I	4	96*	5.18	1.98	4.96-5.39	45
	Kamen	Class II	-	-	6.05	5.25	5.52-6.55	14
		Class III	-	-	9.98	5.14	9.54-10.88	6
		Class IV	-	-	3.70	8.97	3.35-4.01	3

FIGURE 3 | Box-plot showing the distribution of genome size of each taxon and each genome class (according to

Table 3). Frontiers in Plant Science | www.frontiersin.org 10 April 2022 | Volume 13 | Article 873471

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.

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.

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.

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

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.