Introduction

The Balkan Peninsula belongs to the major European biodiversity centers and, at the same time, to the most mountainous regions of the continent (Turrill 1929; Kryštufek and Reed 2004). More than two-thirds of this area are covered by mountains (Reed et al. 2004). While its southern part climatically and floristically pertains to the oro-mediterranean system, the central and northern parts have a more pronounced affinity to the temperate Central European mountains (Horvat et al. 1974). Hence, the flora of the latter area, while being part of the Balkan endemism hotspot (Stevanović et al. 2007), is also strongly related to more northerly mountain ranges, especially the Alps and the Carpathians, with quite similar high-mountain communities (Horvat 1953; Stevanović et al. 2009). This mountainous region consists of four main groups: (1) Dinaric Alps or Dinaric Mountain System, (2) Scardo-Pindhic Mountain System, (3) Rhodope-Rila Mountain System, and (4) Balkan (Stara Planina) Mountain System (Fig. 1). The low mountain range in northern Serbia, where there is no alpine belt, represents a continuation of the Balkan chain toward the Carpathians and is geologically referred to as the Carpatho-Balkanides (Stevanović et al. 2009). These main and internally further divided areas form a topographically, geologically, and environmentally diverse system with alpine habitats characterized by varying extent of spatial isolation and altitudinal gradients.

Fig. 1
figure 1

A Distribution of Campanula orbelica based on available literature, herbarium and field sources (UTM grid zones 34 T and 35 T; dots correspond to basic MGRS squares of 10 × 10 km); dots with question mark indicate localities which require verification. B Sampled populations of Campanula orbelica (circles) against a background of regional topography (elevation classes indicated in the legend). Different colors of circles and black pictogram inserts depict detected genetic lineages: red circle with square insert—Pirin lineage, yellow circle with dot insert—Šar-Korab lineage and orange circle with triangle insert—Rila lineage. All symbols are consistently used throughout the figures (see Fig. 2). Solid and dashed lines show the main hierarchical phylogeographic structure based on a summary of results. Numbering of populations corresponds to Table 1

Full size image

The rich documentation of European phylogeography led to define major southerly Quaternary refugia on the Iberian, Apennine, and Balkan Peninsulas (Taberlet et al. 1998; Hewitt 1999). Regional phylogeographic studies focusing on these areas pointed not only to their role as a source of postglacial northward re-colonization but also to their heterogeneity and internally complex history of biota, preconditioned by diverse landscapes and a relatively low influence of glaciations (Nieto Feliner 2014; Pabijan et al. 2015). Recent studies have increasingly addressed the high evolutionary rates and phylogeographic structure of mountain plants in the Balkan Peninsula (for a review, see Španiel and Rešetnik 2022), but insight into phylogeographic history and divergence of regional intraspecific lineages remains incomplete and limited with respect to region and bedrock coverage. Intraspecific differentiation patterns of the high-mountain flora have been in particular advanced on the example of calciphilous plant species in the western part of the area, the Dinaric Mountains stretching along the Adriatic coast (Surina et al. 2011, 2014; Kutnjak et al. 2014; Đurović et al. 2021). Taxonomic and evolutionary studies of several plant groups further highlighted the importance of allopatric speciation within the Balkan Peninsula (Kučera et al. 2010; Kuzmanović et al. 2013; Lakušić et al. 2013; Janković et al. 2016; Španiel et al. 2017; Rešetnik et al. 2022). However, the intraspecific structure of alpine plant populations across the mountain systems within the central Balkan Peninsula, especially in its eastern part and its relationship to the western ranges, remains a major gap in phylogeographic studies (Španiel and Rešetnik 2022; but see, e.g., Kuzmanović et al. 2013; Petrova et al. 2015; Rešetnik et al. 2022). We thus lack a more detailed understanding of important aspects in how the mountain plant biodiversity has evolved in this region.

Here, we explore the high-mountain phylogeography in the central Balkan Peninsula using Campanula orbelica Pančić (Campanulaceae) as an example, a common species ecologically characteristic of siliceous alpine grasslands of the order Caricetalia curvulae (Puşcaş and Choler 2012), less frequently occurring also at lower sites in subalpine belt. This taxon provides an appropriate case study for several aspects. First, it occurs in all three well-delimited, main mountain systems of the central Balkan Peninsula: Scardo-Pindhic range in the west, Rhodopean (Rhodopes-Rila) and Balkan (Stara Planina) ranges in the east. Second, recent genetic and morphological evidence confirmed that it represents an evolutionarily distinct endemic species limited to the Balkan Peninsula. It is thus strictly allopatric with respect to its presumable close relative Campanula alpina Jacq., which occurs more to the north in the Alps and the Carpathians (Ronikier and Zalewska-Gałosz 2014). Hence, C. orbelica is well suited to study the evolutionary relationships among isolated population groups in the Balkan mountains as its late Quaternary genetic structure within the area is not influenced by gene flow from external lineages. Third, an earlier, phylogeny-oriented study suggested a significant intraspecific genetic structure in this species including two distinct plastid DNA clades (Ronikier and Zalewska-Gałosz 2014).

In an alpine plant range formed as disjunct environmental islands, it is mainly the low-elevation areas among the mountain systems that create obvious barriers to dispersal and gene flow, triggering intraspecific differentiation. Availability of suitable bedrock also places spatial barriers to species migrations and establishment. The tectonic history of geographical ranges is further overlaid by the dynamic impact of changing environments on species in a spatial context (Serra-Varela et al. 2015; Cao et al. 2018). During the Pleistocene, species’ ranges were chiefly impacted by the recurrent climatic fluctuations, which directly influenced the species’ niche availability through time. Past isolation in discrete refugia may also have pushed for adaptations and shifts in ecological niches, which could have further contributed to genetic isolation of lineages. It can generally be expected that the Balkan Peninsula, an area of high landscape complexity and moderate direct impact of glaciations in the past, offers a complex biogeographical history of species and patterns of glacial refugia (see also Španiel and Rešetnik 2022).

We aim at unraveling the phylogeographic pattern of a common high-mountain species of siliceous habitats in the central part of the Balkan Peninsula, and discuss possible drivers underlying the genetic diversity and divergence patterns. The genetic analysis of Campanula orbelica was based on a range-wide sampling, AFLP (Amplified Fragment Length Polymorphism; Vos et al. 1995) genome-wide fingerprinting, and sequencing of several non-coding plastid DNA regions. We further inferred and modeled the climatic niches specific for the entire species and for its main genetic lineages to test how the dynamics of their availability over time might have influenced the contemporary genetic patterns and whether the observed genetic differentiation might be associated with environmental (climatic niche) shifts. Based on these data, we addressed the following questions: (i) does the genetic structure of C. orbelica indicate intraspecific diversification implying a complex Pleistocene history with isolation in local refugia; (ii) if divergent groups are observed in extant populations, what are the likely drivers of diversification? More generally, we aim to answer whether intraspecific diversification patterns correlate with isolation in separate mountain systems (hence, with geographical distances), with landscape barriers (low-elevation formations, incompatible bedrock), or whether they are possibly driven by other factors, such as shifts in the environmental niche?

Materials and methods

Plant material

The study is based on the population sampling and genetic data earlier used for the general analysis of phylogenetic relationships within C. alpina s.l. (Ronikier and Zalewska-Gałosz 2014), and hence, the coherence of populations at the species level is confirmed. One to several populations were collected in each mountain range where the species is reported (14 populations in total) to cover the distribution range of Campanula orbelica (Fig. 1). This included both the western part of the range (chiefly the area of the Šar-Planina–Korab massifs), and the mountain ranges of the eastern Balkan Peninsula; extent of sampling in particular ranges depended on commonness of the species. All samples were used for the AFLP analysis, while a subset of nine populations representing all separate mountain ranges was used for the plastid DNA sequence survey (with three samples per population; Table 1).

Table 1 Sampled populations of Campanula orbelica: numbering, geographic origin (country: BG—Bulgaria; MK—North Macedonia), major mountain units, coordinates and altitudes, population sampling for AFLP analysis and plastid DNA sequencing (in parentheses), AFLP-based genetic parameters of populations (Poly—proportion of polymorphic bands; Div—gene diversity)
Full size table

DNA isolation, AFLP fingerprinting, and plastid DNA sequencing

A detailed description of molecular methods was included in Ronikier et al. (2008) and Ronikier and Zalewska-Gałosz (2014). In brief, total DNA was isolated from 10 to 15 mg of dried plant material using DNeasy 96 Plant Mini Kit (Qiagen), according to the manufacturer’s protocol. AFLP analysis was performed according to Vos et al. (1995) with modifications and selective primers selection as described by Ronikier et al. (2008). Reproducibility of AFLP markers was assessed according to Bonin et al. (2004) using 21 within-plate and nine between-plates replicates (in a joint AFLP analysis of C.orbelica and C. alpina).

Three non-coding plastid DNA regions were used to assess DNA sequence variation in C. orbelica: psbD-trnT, 3′ rps16-5′ trnK and trnS-trnfM, using universal primers of Demesure et al. (1995) and Shaw et al. (2007).

AFLP and plastid DNA sequencing data analysis

AFLP data files were aligned (using ROX-500 internal standard; Applied Biosystems) and manually scored in Genographer 2.1 (http://sourceforge.net/projects/genographer). Dominant AFLP markers were scored in the size range 50–500 bp and coded as present (1) or absent (0) in a binary data matrix. Only markers unambiguously scoreable (well separated) and repeatable among duplicates were considered in the scoring. The error rate was calculated as the share of mismatches in the scoring of AFLP profiles of replicated individuals. Basic statistics were calculated for each population using the R-script AFLPdat (R Development Core Team 2004; Ehrich 2006): the total number of AFLP bands, percentage of polymorphic markers, and Nei’s gene diversity in populations (Nei 1987). The number of private markers restricted to defined groups was noted for geographical groups (separate mountain massifs) and for phylogeographic groups revealed in the Balkan Peninsula. General relationships of individuals in C. orbelica were assessed based on AFLP markers by (i) a Neighbor-Net to analyze the clustering pattern of individuals enabling multiple connections, calculated using SplitsTree 4.12.6 (Huson and Bryant 2006) and based on a matrix of uncorrected P-distances, with a bootstrap analysis based on 1000 replicates computed as implemented in the program, and (ii) a Principal Coordinate Analysis (PCoA) based on Euclidean distances and calculated in MVSP 3.13p (Kovach 1999). The phylogeographic structure (presence of genetic groups) was further tested using a Bayesian analysis of population structure in STRUCTURE 2.2, applying the recessive allele model for dominant markers (Pritchard et al. 2000; Falush et al. 2007). The admixture model with uncorrelated allele frequencies was used. Numbers of K ranging from 1 to 10 were tested with ten replicates per K. 1 × 106 Markov Chain Monte Carlo repetitions were applied with a burn-in period of 200.000. The computations were carried out at the Bioportal of the University of Oslo (http://www.bioportal.uio.no/). Outputs of all STRUCTURE runs were analyzed using the R-script Structure-sum (Ehrich 2006). Average similarity coefficients among runs were calculated for each K to verify the consistency of replicated runs. The following values were further observed to assess the most appropriate number of clusters: (i) the lnP(D) values, estimates of posterior probabilities provided in STRUCTURE outputs, examined as a function of increasing K; (ii) similarity among runs for each K; (iii) ΔK values, estimating the change of the likelihood function with respect to K and treated as an indicator of the most reliable clustering structure (Evanno et al. 2005). Hierarchical analyses of molecular variance (AMOVA) partitioning were performed using Arlequin 3.5.1.3 (Excoffier et al. 2005) with definition of groups based on the phylogeographic inference; this program was also used to calculate pairwise FST values among populations.

To support reconstruction of putative ancestral areas and range history, a continuous phylogeographic analysis using relaxed random walks (RRW; Lemey et al. 2010) was performed on the AFLP data set using the program BEAST v.1.8.4 (Drummond et al. 2012). A simple substitution model with estimated base frequencies was used for phylogeny inference, a strict clock was applied, and a Bayesian skyline coalescent prior (Drummond et al. 2005) with a piecewise-linear skyline model was set to model population growth. The diffusion process was modeled by a relaxed random walk (log-normal RRW) process, which is an extension of a phylogenetic Brownian motion process that rescales the precision matrix by a branch-specific scalar that is drawn independently from an identical distribution (Nylinder et al. 2014), in our case, a log-normal distribution centered on 1. We specified a prior exponential distribution on the standard deviation of the log-normal distribution with a mean of 5. We added random jitter with a window size 1.0 to the tips, as more individuals were sampled from the same location. Two independent analyses of the diffusion inference were run for 10,000,000 generations, logging parameters every 10,000 generations. The performance of the analysis, mixing and effective sample size values for all parameters were checked in Tracer 1.6.0 (Rambaut et al. 2014). The trees from both runs were combined after removing 10% burn-in using LogCombiner, and a maximum clade credibility (MCC) tree was produced and annotated with Tree Annotator (both programs are part of the BEAST package) and visualized with FigTree 1.4.2 (Rambaut et al. 2014). The diffused MCC tree with annotated diffusion estimates was visualized in SPREAD v.1.0.7 (Bielejec et al. 2011) and projected together with polygons representing ancestral areas using ArcGIS 10.8.1 (ESRI, Redlands, CA, USA).

The plastid DNA sequences were manually assembled from forward and reverse direction sequencing results and aligned using BioEdit 7.1.11 (Hall 1999). Sequences of the three plastid DNA regions were concatenated for analysis, assuming linked inheritance. Relationships between plastid haplotypes were analyzed using TCS 1.21 (Clement et al. 2000), which constructs networks by implementing the statistical parsimony algorithm described by Templeton et al. (1992). Assuming that structural mutations represent single evolutionary events, indels longer than one base pair were coded to represent single base-pair changes. A poly-A/T stretch was excluded from the analysis as mononucleotide repeat characters are prone to homoplasy at large geographic and temporal scales (Ingvarsson et al. 2003). The analyses were run with the default parsimony connection limit of 95% and treating gaps as a fifth character state.

GenBank accession numbers of all DNA sequences are available in Table S1.

Assessing the climatic niches of the species and its genetic lineages

The analyses were based on 337 occurrences of Campanula orbelica (Table S2). The list of occurrences was compiled based on recent field studies, analysis of herbarium material deposited at BEO, BEOU, HMMNH, KRAM, and SOM (acronyms following Thiers 2023), and literature data. Field occurrence data were recorded using a GPS device (Garmin eTrex Legend HCx and Garmin eTrex Vista C). All other distribution data were georeferenced using the OziExplorer 3.95 4 s program. To present the geographical distribution of Campanula orbelica (Fig. 1), chorological data are presented on a raster map with squares of approximately 10 km × 10 km based on the Universal Transverse Mercator (UTM) projection (Lampinen 2001). Latitude and longitude are given with respect to World Geodetic System 84 (WGS84).

Following the genetic analyses, each occurrence was directly classified (based on genetic data) or extrapolated to one of the three allopatric lineages revealed by AFLP markers and plastid DNA haplotypes (see Fig. 2). Contemporary climatic data (monthly temperature and precipitation, and derived bioclimatic variables) were taken from WorldClim (Hijmans et al. 2005) at a 30arc-sec (ca. 1 × 1 km). Past climate data reaching back to the LGM (21 kya) in 1 ky steps downscaled originally by Maiorano et al. (2013) to a 0.5° resolution were further downscaled using geographically weighted regressions to match the spatial resolution of current climate data. Bioclimatic variables were calculated on the downscaled historical climate data following the WorldClim procedure. Climatic niche overlap among all clades was compared following Broennimann et al. (2012) using the same bioclimatic variables as were used for modeling the geographic lineage distributions (see below). Following Warren et al. (2008), we tested for niche equivalency using Schoener’s D, and by setting up the test for niche equivalency using randomizations (here 1000 random draws).

Fig. 2
figure 2

Phylogeographic structure within Campanula orbelica. Numbering of populations indicated on the figures refers to Table 1. A Principal Coordinate Analysis (PCoA) diagram of AFLP data based on Euclidean distances showing individual genetic relationships; different symbols depict populations as indicated in the legend. B Neighbor-Net based on AFLP data showing individual genetic relationships; bootstrap values above 75 for main clusters are indicated. C statistical parsimony network of plastid DNA haplotypes (95% parsimony connection limit). Haplotypes are marked with capital letters A–I, circle size reflects the sample size for a given haplotype. Black dots represent intermediate haplotypes missing in the sample set. Lines represent links among haplotypes; solid and pointed lines indicate nucleotide substitutions and structural mutations (indels), respectively

Full size image

As potentially suitable habitats in the past may have had a different geographic pattern, which may explain today’s observed genetic structure, the potential niche was also hindcasted into the past. As only limited occurrence data were available for populations representing the Pirin lineage, which impedes constructing a distribution model for this lineage, we modeled the species as one entity to find past geographic patterns in the potential distribution of C. orbelica to shed light on the genetic structure observed today. A set of bioclimatic variables (as available from www.worldclim.org) was selected and additionally calculated on the downscaled past monthly temperature and precipitation maps, such that the variables are not strongly correlated (|r|< 0.7) and are ecologically interpretable. The selected bioclim variables are: mean annual temperature (bio1), mean temperature of wettest (bio8) and driest (bio9) quarter, annual precipitation (bio12), and precipitation of the driest month (bio14). We stratified the entire study region into ten climatic strata by binning these five bioclimate layers regularly and we randomly sampled ca. 150 pseudo-absence points from each climatic stratum.

A generalized additive model (GAM; binomial family with logit-link) was used to model climatic suitability for C. orbelica in the study region. Pseudo-absences were weighted, such that their summed weight equaled the summed weight of all known occurrences. We allowed for rather low complexity in response curves by limiting the degrees of freedom in the GAM smoothers to a maximum of 3. The fitted model performed well in a split-sample (70% training; 30% testing) cross-validation procedure; means across 100 split-sample replicates: AUC = 0.998, TSS = 0.981, kappa = 0.948. Partial response curves (Suppl. Fig. 1) indicate how C. orbelica responds along the climatic gradients used in the model.

Results

Genetic data and phylogeographic patterns

Two samples were excluded from the analysis due to lack of amplification. In total, 201 AFLP fragments were scored, of which 83.1% were polymorphic. The estimated error rate of the analysis was ca. 2% and 5% for within- and between-plates replicates, respectively (3.5% over all replicates). The proportion of polymorphic fragments in populations varied from about 6% in Stara Planina (population 1), Rhodopes (2) and Vitoša (9) to 14% and 16.49% in the Rila Mountains (populations 3 and 5, respectively). The gene diversity varied from 0.03 (populations 1, 2, and 9) to 0.06–0.07 mainly in the Rila Mountains (populations 3, 4, 5, 7, 12; Table 1). Average values for the proportion of polymorphic markers and gene diversity in the entire data set were 10.30% ± 2.92% and 0.043 ± 0.012, respectively.

A clear phylogeographic structure was detected in C. orbelica based on AFLP data. Both Neighbor-Net and PCoA analyses primarily identified the existence of three distinct lineages (Figs. 1, 2); one of them (hereafter referred to as Pirin lineage) comprised all populations from the Pirin mountains (10–12), the second one (Šar-Korab lineage) all populations from the westernmost, disjunct part of the range in the Šar-Planina and Korab Mountains (13, 14), and the largest genetic group (Rila lineage) contained the remaining populations from the eastern Balkan mountain ranges (Rila, Vitoša, Stara Planina and Rhodopes; populations 1–9). The consistency of these groups was supported by the Bayesian inference in STRUCTURE. K = 3 was designated as the optimal partition based on the ΔK parameter (Fig. 3). It was also the only K value with a similarity among runs equaling 1. The groups were strictly congruent with those found in the clustering and ordination analyses, and the individual assignments showed virtually no admixture (Fig. 3). In a hierarchical division of groups, the strongest genetic break appears to isolate the three populations located across the Pirin Mts. from the remaining samples. These populations are the most distant along the first PCoA axis, and in the STRUCTURE analysis for K = 2, the Pirin populations were segregated as a separate genetic pool against populations from all remaining mountain ranges (data not shown). The significance of these phylogeographic groups was confirmed by a hierarchical AMOVA, which attributed 25.99% of variation to the among-groups level, with an FST value of 0.46 (Table 2). Within the mountain ranges (excluding Rhodopes and Stara Planina represented only by single populations), the lowest average pairwise FST value characterized the Rila Mountains (0.118 ± 0.082) despite the highest number of populations being available here, while this value was 0.247 for Šar-Planina/Korab and 0.277 ± 0.137 for Pirin. Average pairwise FST values calculated for populations among groups were higher for Pirin (0.472 ± 0.082 with Rila and 0.550 ± 0.040 with Šar-Planina) than for Rila with Šar-Planina (0.442 ± 0.072). The highest number of regional private markers was found for Pirin (23) followed by Šar-Planina/Korab (14). In the third phylogeographic lineage, populations from the Rila mountains had ten private markers, while other mountain ranges only three (Rhodopes, Stara Planina) or one (Vitoša) each; however, the third group treated together was characterized by 22 private markers.

Fig. 3
figure 3

Overview of the Bayesian analysis of population structure results (STRUCTURE) for populations of Campanula orbelica. A assignment of individuals to groups at K = 3. Colors and graphical symbols refer to Fig. 1 and numbering of populations to Table 1. B average similarity coefficients among pairs of runs for each K (with standard deviations). C ln P(Data) values in function of K. D delta K values (according to Evanno et al. 2005) for each K

Full size image
Table 2 Analysis of molecular variance (AMOVA) of Campanula orbelica in the Balkan Peninsula without and with using the genetic groups detected
Full size table

The MCC tree based on AFLP data used for the relaxed random walk analysis displayed groups congruent with grouping of individuals shown by other analyses. The RRW-based continuous phylogeographic analysis indicated that the ancestral area of C. orbelica was in the eastern part of its distribution range, likely in the Pirin or Rila Mts. (Fig. 4). From there, it gradually expanded into the neighboring areas and it also dispersed at an early stage to the western part of its current range (Šar-Planina/Korab). Following this step, extension of regional ranges occurred both in the east and west with dispersal to the Stara Planina mountains as the most recent notable event inferred.

Fig. 4
figure 4

Selected outputs of the continuous phylogeographic analysis in BEAST based on the MCC tree of C. orbelica from AFLP data, visualized using SPREAD. Estimated ancestral node areas (with 80% highest posterior density) are marked as shaded polygons). The white star in the upper left map indicates the inferred starting point of diversification. Numbers indicate relative timescale of diversification (from 1 = origin to 0 = present)

Full size image

A complete coverage of the studied plastid regions was obtained for all 27 samples. The alignments of the psbD-trnT, 3′ rps16-5′ trnK, and trnS-trnfM regions were 1333 bp, 759 bp, and 1077 bp long, respectively, and the concatenated alignment was 3169 bp long. The whole plastid sequence data set included four nucleotide substitutions and nine structural mutations (insertions/deletions from one to eight base-pair long). Nine haplotypes (A–I) were identified in the phylogeographic analysis. Two haplotypes with the highest frequency were A (0.30) present in several eastern Balkan ranges and E (0.22), which was unique and fixed for both western Balkan populations (13, 14). The TCS network of haplotypes displayed a sequence of closely related haplotypes, most of which were delimited based on single indel differences (Fig. 2). Populations from Šar-Planina and Korab (13, 14) and Stara Planina (1) were fixed for unique, local haplotypes (E and F, respectively). The remaining populations, except one from the Rila Mountains (3), had two or three haplotypes. Regional distribution of cpDNA variation was mostly consistent with the AFLP groups. In particular, populations from the Rila, Rhodopes, and Vitoša mountains contained only closely related haplotypes (A–C) located at one edge of the network, while plants from the Pirin mountains contained haplotypes grouped in the opposite part of the network (G–I).

Current and historical climatic niches of lineages

Our genetic data showed three distinct, allopatric lineages (see above; Figs. 1, 2). The niche overlap between the Pirin and Šar-Korab lineages, as well as Rila and Šar-Korab lineages is very low (D ≤ 0.001); consequently, their climatic niches differ significantly (p = 0.001; Table 3). The Pirin and Rila lineages overlap considerably in terms of their niche (D = 0.520; orange vs. red dots in Fig. 5). In summary, these findings parallel the geographic distances. The Pirin lineage occupies a climatic niche that is close to that of the adjacent Rila lineage, but differs from the climatic niche of the geographically more distant Šar-Korab lineage. The latter occupies a climatic niche that does not overlap with either of the eastern lineages.

Table 3 Niche overlap among three genetic lineages of Campanula orbelica
Full size table
Fig. 5
figure 5

Genetic lineages of Campanula orbelica in climatic space. Grey dots represent the background climate of the entire study region. Red dots show the Pirin lineage, orange dots—the Rila lineage, and yellow dots—the Šar Planina/Korab lineage

Full size image

When projected in space, the model predicts highly suitable climates in congruence with areas of known C. orbelica occurrences (Fig. 6; current). The hindcasting reveals that the area with suitable climate extends spatially to lower elevations when going backward in time toward the LGM (Fig. 6; 20 k). Hence, in colder climate time periods, C. orbelica had ecological niche extent allowing a more widespread distribution compared to its current distribution that is restricted to high-mountain regions, yet a contiguous distribution range may not have existed across the three regions at the LGM.

Fig. 6
figure 6

Potential distribution of Campanula orbelica according to climatic niche modeling, projected to today’s (current) and past (1 kya, 5 ka, 10 kya, and 20 kya) climate. See the color legend for projected suitability

Full size image

Discussion

Divergent lineages indicate several discrete Pleistocene refugia

The distribution range of Campanula orbelica extends over three well-delimited areas geographically defined by the mountain systems of the Balkan Peninsula, namely the Scardo-Pindhic range in the west, the Rhodopean (Rhodopes-Rila) and Balkan (Stara Planina) ranges in the east. In the Rhodopes-Rila system the species occurs in several distinct yet not very distant mountain ranges, such as Pirin, Rola, Vitoša, and Rhodopes. The data its occurence in the Metohijske Prokletije mts. in the southeastern Dinarides and in the western parts of the Stara Planina mts. in Serbia require confirmation in the field (Fig. 1A).

Three well-defined genetic lineages were revealed (Figs. 2, 3), with a significant impact on the species’ genetic variation (Table 2). They were, however, only partially correlated with the above main mountain systems. The highest divergence characterized the Pirin lineage populations (10–12), as supported by the primary division against the remaining populations in the structure analysis at K = 2, highest pairwise Fst values with other groups of populations, and by the highest number of AFLP private markers restricted to a single mountain range. The plastid lineage containing the Pirin populations was resolved as sister to the rest of the Balkan range in the phylogeny of C. alpina s.l. (Ronikier and Zalewska-Gałosz 2014). While a major phylogeographic break could a priori be expected between populations separated by the main ca. 200 km-wide east–west geographical disjunction, the Pirin divergence did not correlate with geographical range discontinuities. This range forms a direct southward prolongation of the Rila massif, separated only by the lowering of the Predel pass (1140 m), and hence, populations from these two ranges are isolated by the smallest spatial geographical distance (Fig. 1). The general congruence of AFLP and plastid data sets indicates that the isolation of the Pirin populations from the rest of the range, including the nearby Rila Mts., may have persisted for a long period. As no other detailed intraspecific case studies in the central Balkan Peninsula extend over the entire (both east and west) regional range of a high-mountain species to cover all important parts of the island-like distribution, more data from other species are needed to test whether the pattern found for the Pirin Mts. in C. orbelica represents an idiosyncratic species history, or whether it unravels a more significant biogeographical pattern. A recent study of the siliceous alpine snowbed species R. crenatus Waldst. & Kit. (Kuzmanović et al. 2021) suggested colonization of the Rila Mts. from Šar-Planina by a lineage that had diverged earlier. This study focused on the western Balkan Peninsula and no larger sampling across the eastern range was included, which would have allowed a comparison of Rila/Pirin populations. However, another intraspecific study of a silicicolous high-mountain plant (Carex curvula All., a dominant, wind-pollinated sedge), which included samples from Rila and Pirin, indicated that they formed a single, well-supported group in a European range-wide analysis, with a rather low among-population variation (Puşcaş et al. 2008). No significant differentiation was found among single local populations of Cerastium decalvans subsp. orbelicum from the Rila, Pirin, and Rhodopes mountains either (Niketić et al. 2022). These available scarce comparative data do not support at the moment a strong biogeographical pattern repeated across species.

The second phylogeographic break separates the geographically most distant Scardo-Pindhic range (populations of the Šar Planina and Korab), characterized by their unique and fixed haplotype E. The continuous phylogeographic analysis placed the onset of diversification for C. orbelica in the east followed by early dispersal to the Scardo-Pindhic range, which formed the current major disjunction (Fig. 4). Clear differentiation of lineages indicates that gene flow between the western and eastern ranges was restricted for a long time.

Finally, the geographically largest group comprises populations from several mountain ranges in the eastern part of the range. Although distributed in two distinct mountain systems, Rhodopean (Rila, Vitoša, Rhodopes) and Balkan (Stara Planina), this group is genetically coherent and, according to the continuous phylogeographic analysis, it dispersed relatively late in the species’ range expansion history with Stara Planina being the most recently colonized range. Highest population genetic diversity values and highest numbers of private markers support the Rila mountains to be the primary evolutionary center of this group, from which it spread to other areas with nowadays lower genetic diversity of populations (Table 1). The latter areas also have much lower numbers of private markers, although this indicator may be biased by single populations only available in Stara Planina and Rhodopes. In the Neighbor-Net analysis, minor branches of peripheral populations (Rhodopes, Stara Planina) are nested within the Rila mountains, while—again challenging the geographical correlation—the least remote Vitoša has a more distinct position (Fig. 2). Also, a separate plastid haplotype in Stara Planina (F) contrasts to some extent with the AFLP results, although short distances among haplotypes allow to mostly refer to a large pattern only. In any case, our data indicate a large-scale, recent dispersal and/or efficient gene flow within the central Bulgarian mountains, which could be correlated with the projected Last Glacial Maximum extension of suitable climatic niche space (Fig. 5). In a study of Haberlea rhodopensis Friv., an eastern Balkan palaeoendemic, some trend of latitudinal genetic diversification was found with separation of the populations in the Bulgarian Rhodopes and Stara Planina mountains (Petrova et al. 2015). Somewhat higher diversification of a rock-dwelling H. rhodopensis compared to C. orbelica growing in grasslands may be linked with differences in long-term habitat isolation; however, in the context of the Tertiary age of H. rhodopensis, its moderate regional diversification is considered to reflect a recent, late Pleistocene gene flow restriction and it does not support a long-term isolation.

Overall, the phylogeography of the Balkan populations reveals the presence of well-delimited groups which probably evolved through geographic isolation. Puşcaş et al. (2008) suggested that the relatively pronounced genetic isolation of Carex curvula populations in the Balkan peninsula compared to other European mountain regions was due to a smaller impact of the last glaciation and in consequence a longer lasting isolation with restricted possibilities for gene flow among mountain ranges. This would contrast with more northerly mountains such as the Carpathians where the significant lowering of the alpine belt may have facilitated gene flow in populations of high-mountain taxa in the past (Puşcaş et al. 2008; Ronikier 2011). In Campanula orbelica, however, the Balkan Peninsula populations have similar characteristics as those of the putative vicariant species C. alpina in the Carpathians, with several groups, similar diversity levels, and influence of cryptic barriers (Ronikier et al. 2008). Phylogeographic lineages likely constitute genetic traces of multiple refugia, and this pattern further confirms the general concept of ‘refugia within refugia’ (Gómez and Lunt 2007) for survival of cold-adapted species within the macrorefugium of the Balkan Peninsula (Surina et al. 2011; Pabijan et al. 2015; Đurović et al. 2021).

Seeking correlates for genetic divergence patterns

The divergence of phylogeographic lineages is only in part explicitly correlated with spatial distance and obvious geographical/ecological barriers (e.g., large valleys and low-elevation plains). Genetic delimitation of the disjunct Šar-Korab lineages from the Scardo-Pindhic system was expected due to its spatial isolation. Modeling the climatic space of the delimited intraspecific groups clearly parallels their geographic distribution (Fig. 5) and the main east–west disjunction is also correlated with a climatic niche separation. Thus, geographical distance can be considered to have acted here, after erly dispersal, as primary driver of long-term isolation, significant ecological niche diversification, and genetic divergence. Hindcasting of suitable niche space for the entire species supported the maintenance of long-term isolation between the eastern and western parts of the range (Fig. 6).

In contrast, no clear context was found to explain the phylogeographic structure observed within the eastern Balkan mountain ranges. Here, despite the presence of geographical disjunctions among massifs in the Rila lineage, efficient dispersal contrasts with the prominent genetic isolation of the Pirin populations, which is the most intriguing feature of our case study. Our sampling in the Rila and Pirin mountains spans various parts of these ranges including sites in the adjoining northern Pirin (Vihren) and southern Rila (Ribni Ezera). The crest lowering, which separates suitable high-mountain (ca. > 2000 m a.s.l.) habitats of the two ranges does not exceed ca. 10–15 km. While today’s upper forest limit, formed by Picea abies (L.) H.Karst. and Pinus peuce Griseb., reaches ca. 2000–2100 m, the palynological reconstruction suggests dominance of cold and dry open landscape vegetation in the mountains with a limited, patchy presence of tree species during the Younger Dryas period (Stefanova and Ammann 2003). The geographical border between the Pirin and Rila mountains thus did not constitute a strong and long-lasting ecological barrier for the high-mountain biota. This finding is confirmed by our hindcasted climatic niche for C. orbelica; currently populations occupy highly retracted potential habitats, but hindcasting the potential niche indicates a general increase in suitable climatic niche back in time to the Last Glacial Maximum-related climate cooling, a pattern congruent with findings for other alpine plant taxa (Alsos et al. 2009; Espíndola et al. 2012; but see Đurović et al. 2021). At the time where the area of potentially suitable climates is maximized (Last Glacial Maximum – 20 ka), there may have been a connection between areas now occupied by the Pirin-clade and remaining eastern populations, especially in the Rila Mts. (Fig. 6).

In an analysis of the current climatic niche of the Pirin and Rila groups, overlapping patterns of climatic niches of these eastern phylogeographic lineages do not concur with the strength of genetic differentiation (Fig. 5). Therefore, climatic niche shift at this scale can be rejected as a factor for this divergence. Edaphic conditions are another key ecological factor in shaping the biogeographical history of alpine plants (Alvarez et al. 2009) through their direct impact on habitat conditions. However, no edaphic barriers or major bedrock transitions occur that would establish a significant ecological barrier between the two ranges (Läßiger et al. 2008). Sharp phylogeographic breaks, which spatially segregated distinct genetic lineages, and sometimes likely predate several glacial episodes, are observed in the European mountains (e.g., Schönswetter et al. 2005; Ronikier et al. 2008, 2012) and they are not always correlated with significant geographical disjunctions or environmental transitions, thereby being challenging to interpret.

However, different bedrock conditions, although not apparent as a barrier among the lineages, may hypothetically have triggered adaptive divergence and lineage diversification within the geologically diverse Pirin area. Campanula orbelica is mainly confined there to typical siliceous habitats, but its populations also occur within the limestone area of the Vihren massif. Strict genetic partitioning was observed between populations of Campanula alpina from granitic and limestone massifs of the Eastern Alps, although it could be attributed either to bedrock-related differentiation or to historical isolation in distinct refugia (Ronikier et al. 2008). Bedrock-driven divergence was demonstrated for geographically adjacent siliceous and calcicolous populations of Ranunculus crenatus s.l. in the Balkan Peninsula mountains (Kuzmanović et al. 2021). Such a microspeciation mechanism could have also acted in C. orbelica, although subsequent local introgression between bedrock-differentiated populations may have limited the ecological speciation, while maintaining the genetic divergence of the local lineage. Notably, a slight differentiation between the limestone (no. 10, 11) and siliceous (12) populations was detected both in AFLP and plastid DNA data (Fig. 2). The Pirin populations also had the highest mean within-range pairwise Fst value which supports some degree of maintained diversification. A fine-scale genetic and morphological analysis of populations in the Pirin–Rila transition will be needed to further test the hypothesis of bedrock influence.

From the point of view of floristic composition, the Pirin range is characterized by a high number of endemic species. This is not a unique characteristic of this range as all high-elevation Bulgarian ranges reveal high levels of endemism (Peev and Delcheva 2007), but, notably, the regional stenoendemic elements are mostly linked to the limestone nucleus of the Vihren massif (Andreev et al. 2006), which supports the role of this area as a diversification center.

An alternative scenario to consider for explaining the divergence between the adjacent Rila and Pirin could be a recent dispersal to the Rila mountains from other massifs of the genetic Rila group (Rhodopes or Stara Planina) following earlier lineage diversification. However, it is not supported by the continuous phylogeographic analysis and also less plausible considering much higher genetic diversity of the Rila populations compared to those from Rhodopes and Stara Planina, and highest abundance of C. orbelica in the Rila mountains.

Clearly, a complex combination of processes was responsible for the genetic structure observed in our case study. Increased regional population sampling coupled with fine-scale analysis in the Pirin–Rila transition may help further pinpoint the most probable drivers of this pattern. Our results show the need for a detailed phylogeographic survey of the Balkan high-mountain species and may indicate patterns important for the conservation of evolutionarily important units in the high-mountain biota, which is one of the most important elements behind this major biodiversity center in Europe.