Evolutionary response of cold‐adapted chasmophytic plants to Quaternary climatic oscillations in the Mountains of Central Asia (a world hotspot of biodiversity)

Past climatic oscillations are the main driving force of evolutionary changes in alpine species. Species' response to paleoclimatic oscillations is crucial in forecasting their future response in face of climate warming. The aim of this research is to explore the effect of climatic fluctuations on the evolutionary history, demography, and distribution of high‐mountain bellflowers (Campanula lehmanniana complex), the flagship and taxonomically problematic members of chasmophytic vegetation within an underexplored biodiversity hotspot, the Mountains of Central Asia.


| INTRODUC TI ON
Heterogeneous topography, distinct environmental gradients, and complex geological history are the main factors that shape the distribution and genetic patterns of alpine species and trigger unusually high biodiversity of high-mountains regions (AEgisdóttir et al., 2009;Cain et al., 2000;Frei et al., 2012b;Körner, 2003;Ortego & Knowles, 2022;Rahbek, Borregaard, Antonelli, et al., 2019).A combination of dramatic past climate changes with subsequent habitat shifts in high-mountain regions promoted speciation processes reflected currently in a high richness of narrow-ranged endemic species observed in alpine biomes throughout the globe (Flantua et al., 2019;Mosbrugger et al., 2018;Muellner-Riehl, 2019).
Understanding how species have responded to past climate oscillations can help forecast their response to current and future changes, which is particularly urgent in the most species-richest and at the same time highly vulnerable to climate warming mountain areas (Opała-Owczarek, 2019;Rahbek, Borregaard, Colwell, et al., 2019;Raxworthy et al., 2008;Razgour et al., 2013;Vintsek et al., 2022;Yu et al., 2019).Nevertheless, despite recent technological advances, the fine-scale evolutionary history of species from underexplored regions remains largely unknown, impairing our understanding of biodiversity patterns.
Environmental changes driven by Pleistocene climatic oscillations altered the distribution of suitable habitats for many taxa (Carnaval et al., 2009;He et al., 2013;Ortego et al., 2021;Tonzo & Ortego, 2021).The range dynamics along latitude and elevation gradients (down and upslope shifts) during glacial and interglacial periods caused recurrent population fragmentation and expansion, leading alternately to the processes of so-called mixing-isolationmixing (Abbott, 2019).The separation of species populations in glacial or interglacial refugia (Bennett & Provan, 2008;Hewitt, 2000; Muellner-Riehl, 2019) led many of them to divergence, and in consequence to speciation processes.Subsequent climate changes (alternate cooling and warming) contributed to range extension, enabling the geographical contact of previously isolated gene pools.Thus, the genetic picture of present-day populations is often the result of complex (contemporary and ancient) processes that have impacted the evolutionary history and demography of species (He et al., 2013;Hewitt, 2000;Muellner-Riehl, 2019;Tonzo & Ortego, 2021;Vintsek et al., 2022;Wróbel et al., 2023;Xue et al., 2021).
The complexity of alpine landscapes can favour adaptive radiation, enabling better exploitation of niches in 'island-like' environments, as well as neutral diversification processes connected with geographical isolation (Perrigo et al., 2020).With one of the largest elevation amplitudes in the world, diverse geology and relief, considerable glacier cover, extreme precipitation and temperature fluctuations reflecting both Mediterranean and temperate patterns, the Mountains of Central Asia harbour a great number of plant species and plant communities (Abdulina, 1999;Lazkov & Sultanova, 2014;Nobis et al., 2013Nobis et al., , 2020;;Nowak et al., 2022;Nowak, Nobis, et al., 2020;Nowak, Świerszcz, et al., 2020) and have been defined as one of the world's biodiversity hotspots (Mittermeier et al., 2011).
All that makes the region an exciting research area to explore the eco-evolutionary patterns of alpine species.The presence of vast areas with isolated, large valleys free of permanent ice during the coldest stages of the Pleistocene made this region an important glacial refugium for warm-temperate Pleistocene or Tertiary taxa, as well as an important interglacial refugium for several cold-adapted organisms currently persisting in fragmented populations in highly elevated habitat patches (Muellner-Riehl, 2019;Yu et al., 2019).A unique combination of environmental conditions has promoted an extraordinarily high rate of endemism in this area-for instance, nearly 1500 (30%) out of ca.4300 vascular plants occurring naturally in the central Pamir-Alai Mts are endemics (Nowak et al., 2011;Raduła et al., 2021).A recent analysis of the degree of endangerment to the flora also shows considerable threats to its richness, where 1627 (ca.38%) of all native species are threatened and 23 species (0.54%) already extinct (Nowak, Świerszcz, et al., 2020).For all the above reasons, the Mountains of Central Asia are also considered one of the most sensitive regions to climate change and biodiversity loss in the world (Ajjur & Al-Ghamdi, 2021;Hartmann et al., 2013).
Species of the genus Campanula s.l.(bellflower) are among the flagship contributors to the floristic diversity of mountain regions throughout the Northern Hemisphere.The highest species diversity of bellflowers is observed in the Old World (Jones et al., 2017;Lammers, 2007;Mansion et al., 2012) and a particularly high concentration of species, often with narrow geographical ranges resulting from adaptation to specific habitat conditions, is found in mountain regions (Cano-Maqueda et al., 2008;Castroviejo et al., 2010;Oganesian, 2008).One of the most interesting representatives of bellflowers, associated with rocky habitats (fissures, clefts, and ledges) of the Central Asian Mountains are species of the C. lehmanniana complex, flag representatives of chasmophytic vegetation.The wide range of morphological variability within and between populations has prompted the recognition of two to six endemic taxa with limited distribution restricted to particular mountain ranges (Fedorov, 1957;Rasulova, 1988;Victorov, 2002) and makes the complex an interesting model for phylogeographic research.
Campanula lehmanniana Bunge was described from the Zreravshan Range (Karatau Mts) by Bunge (1852) and is distributed in the western Pamir-Alai and western Tian Shan Mts (Figure S1).
Later, Franchet (1883), based on differences in stem length, bellflowers, Campanula lehmanniana, Campanulaceae, chasmophytes, demographic history, ecological niche modelling, high-throughput sequencing, molecular clock, phylogeography distinguished C. lehmanniana var.capusii Franch., which was subsequently raised to the rank of species by Fedorov (1957) or the rank of subspecies by Victorov (2002).Fedorov (1957) also distinguished from the complex a species with extremely long leafpetals, single and thin stems, describing it as C. eugeniae Fed., an endemic taxon of the western Tian Shan Mts (Figure S1).Although Kamelin (in Rasulova, 1988)  Fed.mainly based on stem length and the corolla (glabrous vs. shortly pubescent).All the mentioned above taxa belonging to the C. lehmanniana complex are considered microendemics, occurring at elevations from 900 to 3500 m, and are restricted to calcareous rocky-islands (steep rock faces, ledges, fissures ;Rasulova, 1988).
Due to insufficient knowledge about the origin and evolutionary relationship within the C. lehmanniana species complex, we integrated here a comprehensive set of phylogenomic approaches with reconstruction of past and present species distribution models to understand the effects of Q uaternary climatic oscillations on the genetic diversification of chasmophitic species in the Central Asian mountains.We attempt to answer the following questions: (i) do genomic data corroborate hybridisation events or introgression across these taxa, (ii) do the patterns of genomic variations reflect signals of historical population connectivity, e.g., gene flow during glacial periods followed by population isolations in the interglacials?(iii) whether representatives of the studied complex refer to one polymorphic taxon or represent two or more independently evolved lineages, (iv) what are the divergence time and demographic history of particular lineages of the analysed species complex?

| DNA extraction
Genomic DNA was extracted using the Genomic Mini AXE Plant kit (A&A Biotechnology).At the beginning, dried plant leaf tissues were ground using an MM400 (Retsch) mixer mill and 3-5 mm glass beads to a fine powder.The NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific) were used to assess the purity and concentration of DNA, whereas 1% agarose gel electrophoresis was used to check the quality.To meet the DArTseq genotyping protocol requirements, the concentration of each sample was adjusted to 70 ng/μL, while for PCR for direct ITS and cpDNA sequencing, up to 10 ng/μL.
In both cases, agarose gel electrophoresis was used to detect the presence or absence of the target sequence and to check the length of the fragment.Prior to sequencing, PCR products were purified using the Exo-BAP Kit (EURx).For paired-end Sanger sequencing of ITS and cpDNA, the products of PCR were sent to an external company (Genomed).Resulting sequences were manually verified and aligned using BioEdit ver.7.0.5.3 (Hall, 1999).

| Genomic library preparation, DArT sequencing and DArT data filtering
DArTseq is a genome complexity reduction method widely used in ecological, evolutionary, and phylogenetic studies (Al-Beyroutiová et al., 2016;Baiakhmetov et al., 2020;Edet et al., 2018;Ketema et al., 2020;Melville et al., 2017;Smýkal et al., 2018).The method is based on the Next Generation Sequencing (NGS) platform, optimised for each organism and application in order to select the most appropriate complexity reduction method (Cruz et al., 2013;Kilian et al., 2012;Sansaloni et al., 2011).This section was performed according to the procedures previously published (Baiakhmetov et al., 2020).Briefly, genome complexity reduction using restriction enzymes and high-throughput polymorphism detection (Kilian et al., 2012) was performed by Diversity Arrays Technology Pty Ltd.
Based on testing several enzyme combinations, the PstI-MseI were selected.For a detailed description of the DArTseq genomic library preparation, see the Supplementary text in Appendix S1.
For the downstream analyses, we applied co-dominant single nucleotide polymorphisms (SNPs) markers processed in R v. 4.0.3(R Core Team, 2020) with the additional dartR v.2.3.3 package (Gruber et al., 2018).Unless otherwise stated, we applied the following filtration steps: (1) monomorphic loci were removed, (2) loci identified (=called) in greater than 95% of all individuals were kept, (3) loci with a scoring reproducibility of 100% were kept, (4) sequence tags contained more than one SNP were filtered to keep randomly selected one of them, (5) SNPs were filtered based on the criteria of a minor allele frequency (MAF) threshold of 1%.

| Phylogenomic inference, divergence time estimation, and species delimitation
We estimated species trees using data from the direct sequencing of the ITS region, six cpDNA regions, and DArTseq derived SNPs.
Four approaches were applied to analyse the above-mentioned data.

| Bayesian inference
First, we conducted a phylogenetic analysis based on ITS and cpDNA regions by using Bayesian Inference (BI) analysis performed in MrBayes 3.2.6 (Ronquist et al., 2011).Markov chain Monte Carlo (MCMC) simulation was set for 4,000,000 generations, sampling one of every 1000 generations.The first 12.5% of iterations were discarded as a burn-in fraction.The results of the Bayesian analysis were placed into visual form using FigTree v. 1.4.0 software (Rambaut, 2018).BI analysis was performed for ITS concatenated cpDNA regions.Alignment was based on 67 individuals from the C. lehmanniana complex and 3 individuals of C. incanescens as an outgroup.

| Neighbour joining
The phylogenetic tree based on DArTseq derived SNPs was constructed using the neighbour joining (NJ) method based on the Euclidean distance between individuals by using R package dartR (Gruber et al., 2018).The analysis was performed based on 9158 SNPs derived from 247 individuals from the C. lehmanniana complex.
For publication purposes, the phylogenetic tree was visualised and edited using FigTree v. 1.4.4 (Rambaut, 2018).

| Molecular clock based on cpDNA and ITS regions
We performed a molecular clock analysis in BEAST2 (Bouckaert et al., 2019).We used the model averaging method by running the analysis with bModelTest (Bouckaert & Drummond, 2017) as the site model.We applied a one-time calibration by setting a log-normal distribution with a mean of 47.35 Mya and a standard deviation (SD) of 0.08 for the potential divergence of Campanula lehmanniana and Campanula incanescens.We used the results of previous research, which assumed that the divergence of the evolutionary lineages containing the above-mentioned taxa took place  highest posterior density (HPD)] (Jones et al., 2017), and on this basis, we estimated the SD used in our research.Then, divergence time was estimated using the strict clock model and Calibrated Yule prior.In total, we ran the analysis three times independently, with 50 million MCMC generations for each run.LogCombiner v. 2.6.3 (BEAST package) was used for log and tree files combining.In both cases, 10% of generations were discarded as burn-in from each run.
In the next step, Tracer v.1.7.1 (Rambaut et al., 2018) was used to check the effective sample size (ESS) values for log files.As all ESSs exceeded 200, we summarised the final maximum clade credibility tree in TreeAnnotator v. 2.6.3 (part of the BEAST package).The visualisation and editing of the final tree with Bayesian posterior probabilities and 95% credibility intervals (CI) were made using FigTree v. 1.4.0 (Rambaut, 2018).The analysis was performed for the same individuals as for the BI.

| SNAPP
For SNPs, we used SNAPP v.1.5.1 (Stange et al., 2018), as it utilises the multi-species coalescent approach and is well-suited for analyses of genome-wide data deducing the species tree and divergence times directly from SNPs.To focus on the separation of individual evolutionary lineages within the studied complex, we used 3345 DArTseq-derived SNPs.We applied a one-time calibration using our own estimation from the molecular clock based on cpDNA and ITS regions, setting a log-normal distribution with a mean of 5.8 Mya and SD of 0.15 for the split between C. lehmanniana and C. eugeniae.The analysis was performed twice independently, with 2 million MCMC generations for each run using BEAST2 (Bouckaert et al., 2019).The final maximum clade credibility tree was validating and summarising as described above.Lastly, we visualised a pattern across all the posterior trees using DensiTree v. 2.01 (Bouckaert, 2010).Due to its long computation time, SNAPP was performed for only 20 individuals from the C. lehmanniana complex (2-4 individuals per clade separated in NJ).

| Inference of population genetic structure and admixture
To infer the population genetic structure and admixture, we followed the above-mentioned sequence of filtering steps with exception of the called loci (100% instead of 95%).Additionally, genomic regions exhibiting high linkage disequilibrium (r 2 > 0.5) were determined and removed based on 'the sliding window' approach using plink 1.9 (Chang et al., 2015).We consider a window of 50 SNPs, remove one of a pair of SNPs within window, and shift the window by 5 SNPs forward repeating the procedure.
We used two methods for quantifying genetic structure and admixture.First, we used fastSTRUCTURE software, which quickly and resourcefully implements the Bayesian clustering algorithm STRUCTURE (Raj et al., 2014).We tested a number of clusters (Kvalues) ranging from 1 to 10 using the default convergence criterion of 10e −6 and default prior (simple).The estimation of most likely K value was used with the best choice function implemented in fast-STRUCTURE.Second, genetic clustering was performed using the maximum likelihood method ADMIXTURE, testing for 1 up to 10 K (Alexander et al., 2009).Choosing of best K value was assessed by determining the replicate with the lowest cross-validation error.In both cases, the resulting ancestry coefficient matrices (Q-matrices) were imported into STRUCTURE PLOT (Ramasamy et al., 2014) to visualise bar plots.

| D-statistics and the f4-ratio
Since phylogenetic bifurcating trees are often insufficient to capture the evolutionary history of taxa showing admixture (Durand et al., 2011), we used ABBA/BABA tests based on the D statistic and f4-ratio implemented in D-suite software (Malinsky et al., 2021) to examine introgression within the C. lehmanniana complex.We performed the ABBA/BABA test for 247 individuals assigned to 8 genetic clades identified in NJ.As an outgroup, we used 3 individuals of C. incanescens.We used the 5149 SNPs that were left after all stages of filtration for this particular species set.We calculated D-statistics and f4-ratio for all combinations of trios among the studied genotypes.f4-ratio estimates the admixture proportion in percentages when multiplied by 100.For the graphical presentation of ABBA/ BABA test results, we used NJ tree with the percentage admix between clades displayed manually on a tree.For an additional description, see the Supplementary text in Appendix S1.

| Treemix graph
We analysed the potential presence of introgression and determined the direction of gene flow using Treemix v. 1.13 (Pickrell & Pritchard, 2012).Treemix uses genome-wide allele frequency data and Gaussian approximation to genetic drift to construct models which, unlike a bifurcate tree, not only allow the determination of population splits but also of gene flow.We ran Treemix analyses on the same 8 genetic clusters used for ABBA/BABA tests, using an estimated ML tree and assuming the independence of all SNPs with a window size of one SNP (k = 1).As Treemix is vulnerable to the lack of data, we filtered for sites without missing data by specifying a 1.0 threshold during the 'call rate' filtration step, resulting in a final analysis based on 1300 SNPs.We tested a range of migration events (m from 1 to 10).We estimated the optimal number of migration edges by using an ad hoc statistic based on the second-order rate of change in log-likelihood implemented in OptM v. 0.1.6(Fitak, 2021).

| Models of interspecific gene flow
We used fastsimcoal2 (Excoffier et al., 2021) to evaluate a hypothesis about the presence or absence of gene flow between populations from three studied mountain regions (Alai Mts, Tian Shan Mts, and Zeravshan-Hissar Mts) (compare the Supplementary text in Appendix S1).We used 45 individuals (15 from each of the 3 abovementioned groups) and 11,067 SNPs for the analysis.We used easySFS.pyscript to generate the site frequency spectrum (SFS).To construct a complete data matrix, we 'projected down' our sample size up to 13 individuals per group, maximising the number of segregating sites.Each model was run 100 replicated times considering 200,000 coalescent situations to calculate composite likelihood, and 40 optimisation (ECM) cycles to estimate parameters by using maximum composite likelihood (Excoffier et al., 2021).To choose the best-supported model, we used Akaike Information Criterion (AIC).To determine the best parameter estimates, we selected the 10 most likely replicate runs (with the smallest difference between the observed and estimated likelihood) and calculated the mean for all parameters.We calculated their 95% confidence intervals (90 percentiles from 100 parametric bootstraps SFS) by using R package 'jmuOutlier' (Garren, 2022).For additional information, see the Supplementary text in Appendix S1.Estimated models assumed a mutation rate of 7 × 10 −9 per generation per locus assumed for the model plant Arabidopsis thaliana (Belfield et al., 2020;Ossowski et al., 2010).Time estimation of taxa split was based on SNAPP results, assuming a 2-year generation time-the shortest possible time after which perennial bellflowers produce seeds (Shulkina et al., 2003).

| Inference of demographic history
We reconstructed the demographic history of populations representing the C. lehmanniana complex.We used Stairway Plot v.
2.1.1,which is capable of handling folded SNP frequency spectra (without knowing the ancestral alleles).The program does not need the reference genome; therefore, it is particularly suitable for nonmodel organisms (Liu & Fu, 2020).Since the analysis requires a total number of observed nucleic sites, including polymorphic and monomorphic, we remove only loci scored all missing (NA) loci with repeatability <1.We filter data for each of the studied groups separately.We used the databases of 119,281 up to 191,089 SNPs depending on studied group.We choose populations with at least 5 individuals.We grouped populations into groups corresponding to clusters obtained from ADMIXTURE and computed the SFS as described in the previous section.We ran Stairway Plot fitting a flexible multi-demographic model, considering the same generation time and mutation rate as for Models of interspecific gene flow.
We performed 200 bootstrap replicates to estimate 95% confidence intervals.

| Ecological niche modelling
We performed ecological niche modelling (ENM) to predict the geographical distribution of climatically suitable niches for the studied C. lehmanniana complex taxa in the present day and during the last glacial maximum (LGM, ~21,000 years BP) and also to determine whether it supports our hypothesis of secondary contact.
For this purpose, we used the maximum entropy algorithm implemented in Maxent software v. 3.3.3(Phillips et al., 2006;Phillips & Dudík, 2008).Distribution data were obtained from our own field studies and from the revision of the herbarium collections.Based on the results of our phylogenetic and species-delimited analysis, we decided to run separate models for C. lehmanniana and C. eugeniae.
For the preliminary models, we used all 19 bioclimatic variables derived from the WorldClim v. 1.4 database (available online: http:// www.world clim.org/)(Hijmans et al., 2005), from which we selected for the final model 6 uncorrelated: bio2, bio3, bio10, bio11, bio17, bio19.For the LGM, we tested 2 general circulation models (GCM), CCSM4 and MPI-ESM-P.A detailed description of modelling methodology is described in the Supplementary text in Appendix S1.

| Inference of population genetic structure and admixture
FastSTRUCTURE analysis of SNPs revealed the most likely number of clusters at a K value of 6 (Figure 2a).For ADMIXTURE, 'best' K was inferred as K = 10 (Figure 2b).Both analyses revealed a strong population genetic structure, with most genotypes consistent with the clades identified in NJ (Figure 1).Both analyses defined C2/3 as a genotype of admixture origin sharing over 40% of its markers with C3, and the remaining shared in half with C6 and C2.In the ADMIXTURE analysis, <10% is shared with C1B, which may suggest historical events of secondary contact between C. lehmanniana and C. eugeniae (Figure 2b).The results of fastSTRUCTURE and ADMIX-TURE analyses revealed limited evidence of secondary contact be-  S1.Pictures show the variations in the morphology of the studied bellflowers.
NJ results, confirming the existence of two geographically wellseparated clusters, with no evidence of admixture (Figure 2).

| D-statistics and f4-ratio
A statistically significant excess of ABBA patterns (D > 0) supports the complex history of post-speciation admixture within the studied group (Table S3).Statistically significant results were obtained between C1 and C1B clades of C. eugeniae and C4, C5, and C6 clades of C. lehmanniana.In all these cases, the proportions of gene flow are rather low, <10% (6.5%-9.8%),but statistically significant.The analysis also revealed statistically significant signatures of gene flow with varying proportions between the C. lehmanniana clusters.The highest proportion of gene flow has been detected between C5 and C2/3 as well as C4 and C5 with an intensity of 74.7% and 33.4%, respectively (Figure 3a and Table S3).

| Treemix graph
Treemix analyses support the model with four admixture events between and within the examined genetic clades.The results only partially coincide with the results of the ABBA-BABA test and show weak directional gene flow from the common ancestor of the C1 and C1B subclades of C. eugeniae to C. lehmanniana (C4), as well as much stronger evidence of gene flow between different clusters (and/or their common ancestor) within C. lehmanniana (Figure 3b).S1.Clades C1-C2 marked on the yellow line refer to populations of Campanula eugeniae, whereas clades C2/3-C6 marked on the blue line refer to populations of C. lehmanniana (clades numbered C1-C6 according to Figure 1). of C. eugeniae from Thian Shan Mts (C1B including populations 30 and 31) as well as two neighbouring groups of populations from Seven Lakes Valley (C6A and C6B including populations 9,10,11 and 8,7,4).

| Inference of past demographic history
We observed more or less stable levels of Ne in most northern populations of C. lehmanniana located between the Syr-Darya and Zeravshan River valleys as well as from the Pastruddarya and neighbouring Iskanderdarya River valleys (C5A including populations 2, 3, 6 as well as C5B including populations 16, 17, 18) (Figure 5).

| Environmental niche modelling
Maxent models show high performance in predicting current suitable conditions for both C. lehmanniana and C. eugeniae, with average AUC of 0.981 and 0.907, respectively, for models with 30 arc-seconds grid, as well as 0.971 and 0.878 for 2.5 arc-minutes under CCSM4 respectively and 0.970 and 0.894 for 2. with the highest contributions to the model are bio19, bio17, bio10, bio2, bio3, and bio11 (Table 1; Figure 6).The climate-suitable areas model shows distribution patterns roughly congruent with observed distributions for both species studied (Figures 7, S1, and S3).However, potentially suitable areas are much wider than those actually occupied, especially in the case of C. eugeniae, which is a narrow endemic currently known only from several localities in the Talass, Fergana and Chatkal ridges within the Tian Shan Mts.In general, current distribution models based on maps with a spatial resolution of 2.5 arc-minutes show wider areas of potential suitability than models based on maps with a resolution of 30 arc-seconds, which show greater accuracy and specificity, and therefore a narrower overall area of potential suitability.The model under CCSM4 and MPI-ESM-P atmospheric general circulation models showed a increase of potentially suitable area during the LGM period (especially in the case of highly suitable areas) for both studied species (Figures 7c-f and S4; Table S4), when compared with the present distribution models (Figures 7a,b and S3).
According to both CCSM4 and MPI-ESM-P atmospheric general circulation models, during the LGM, ranges of the studied species overlapped in a relatively large area within the Alai and western Tian Shan Mts (Figure S5), which may confirm the hypothesis about secondary historical contact between the studied species, and is consistent with the results of introgression tests described above.The analysis showed that during LGM, the potential range of both species shifted to lower-elevated mountain locations, however, this distribution change was much less pronounced in case of C. eugeniae compared with C. lehmanniana.

| DISCUSS ION
In the genomics era, genome-wide dataset availability has become crucial to understanding the evolutionary history of any organism.Our genomic analyses support the taxonomic division of the Campanula lehmanniana complex into two main species, C. eugeniae and C. lehmanniana (Li, 1987).Furthermore, the analyses suggest that across the Pleistocene, the examined populations distributed within the Pamir-Alai and Tian Shan Mts experienced flickering connectivity reflected by strong genetic fragmentation during interglacial periods, followed by glacial admixtures (Cano-Maqueda et al., 2008;Ortego & Knowles, 2022;Perrigo et al., 2020;Tonzo & Ortego, 2021).

| Species delimitation and divergence time estimation
The India-Eurasia collision taking place between ~55-52 Mya (initial contact) and ~34-20 Mya (final collision) led to the uplift of Central and South Asian mountain massifs, which later experienced further geo-climatic changes favouring unusual radiation of plants (Aitchison et al., 2007;Luo & Li, 2022;Sehsah et al., 2022;Van Hinsbergen et al., 2012).The beginning of this collision can probably be estimated as the starting point of the evolution of Asian Campanula s.l., estimated at ~57 Mya, and coinciding with the onset of global climatic transitions when temperature fluctuations and mean temperatures decreased until the late Oligocene (Jones et al., 2017;Zachos et al., 2008).The Miocene (~23-5 Mya) appeared to be the most 'active' of all geological periods for the evolution of many flowering plant lineages (Cervantes et al., 2016;Löhne et al., 2007), including Campanula s.l.(Jones et al., 2017).Our results are in line with previous findings, suggesting that the split of C. lehmanniana s.l.into two main genetic clades that refer to C. eugeniae and C. lehmanniana occurred in the late Miocene (~5.8 or 5.64 Mya) or early Pliocene (Figure 1).At this time, the Sino-Himalayan region was undergoing climate cooling (Zhisheng et al., 2001).Subsequent long-term declines or increases in temperature during Quaternary glacial-interglacial cycles were responsible for range fluctuations and the so-called 'flickering connectivity' promoting the extensive diversification of alpine plants (Flantua et al., 2019;Mosbrugger et al., 2018;Raduła et al., 2021;Vintsek et al., 2022;Wróbel et al., 2023).These processes also appear to play an important role in splitting both above-mentioned Campanula species.According to our molecular clock based on concatenated ITS and cpDNA as well as SNAPP analyses, C. lehmanniana split ~2.7-1.68Mya and C. eugeniae ~ 1.5-1 Mya (Figure 1) into several further genetic clades currently distributed in isolated geographical regions within the Hissar-Zeravshan and Tian Shan Mts, respectively (Figure 2).
Specimens belonging to subclade C2/3 are characterised by fairly tall stems, but they mainly differ in the character of the calyx and corolla (glabrous vs. pubescent).However, some other specimens of C. hissarica determined using the height of stems were located in subclades with typical representatives of C. lehmanniana (clade C5, Figure 1).A similar issue is observed for specimens corresponding to the description of C. capusii (Fedorov, 1957;Rasulova, 1988).
The taxon was originally described as C. lehmanniana var.capusii Franch.based on specimens collected in the Zeravshan Mts and characterised by short stems (Franchet, 1883).However, the samples that fall into the description of the taxon are segregated in our analyses into four different subclades (C2, C3, C4, and C6).Thus, the independence of this taxon (Fedorov, 1957;Rasulova, 1988) is rather doubtful.

| 'Mixing-Isolation-Mixing' and flickering population connectivity during past climatic oscillations
Both fastSTRUCTURE and ADMIXTURE analyses revealed a strong divergence of particular genetic clades within the analysed C. lehmanniana complex, which is in line with the studies on other endemic taxa associated with island-type rocky refugia within deep valleys of the Central Asian mountains (Vintsek et al., 2022).
Both above-mentioned analyses defined C2/3 subclade as having admixture origin, which may suggest historical events of secondary contact between C2, C3, and even C1B genetic subclades, Lakes Valley, Figure 2b).This suggests that pollen dispersal, similarly as in other plants related to such habitat, is relatively low (Frei et al., 2012a;Kuss et al., 2007;Vintsek et al., 2022).

| Implication for conservation
A demographic history of many high-mountain, cold-adapted plants and animals from different regions worldwide indicates negative impact of paleoclimate warming on their range and effective population size (Guerrina et (Gao et al., 2021;Li et al., 2020;Pepin et al., 2015).Due to the climate warming in the Mountains of Central Asia, recognised as biodiversity hotspots (Mittermeier et al., 2011), is stronger than the global mean warming trend (Hartmann et al., 2013), we can forecast that the demographic population decline of cold-adapted species will proceed.
Given the insular distribution of rock habitats, the elevation shifts and/or dispersion of chasmophytes to more suitable areas may be difficult even under the optimistic scenario (Auld et al., 2022;Vintsek et al., 2022).
The patterns of genetic variation within the C. lehmanniana complex allowed us to identify within the studied area deep genetic structure of the populations and to show the occurrence of several (2-8) genetic lines that can be referred to as Evolutionary Significant Units (Canestrelli et al., 2014;Moritz, 1994).Particular attention should be paid to the Zeravshan-Hissar mountain range, which harbours a number of unique intraspecific lineages of the studied alpine chasmophytes.Therefore, the region can be defined as a hotspot of genetic diversity which occurrence is often connected with past climatic oscillations and presence of both glacial and interglacial refugia (Canestrelli et al., 2014;He et al., 2015).
The south-western part of the Central Asian biodiversity hotspot, including the Zeravshan-Hissar range, harbour a number of coldadapted endemic species (Nowak et al., 2011;Raduła et al., 2021;Tojibaev et al., 2022).The region is also most sensitive to future climate warming (Vintsek et al., 2022), and if the trend continues, many of these taxa will likely run out of suitable habitat in the near future.Therefore, identifying and protecting of micro-hotspots within alpine biodiversity hotspots is contemporary challenge for conservation management (Canestrelli et al., 2014;Ludovicy et al., 2022;Noroozi et al., 2022).

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
The SNP dataset derived from the DArTseq pipeline in the genlight format as well as alignments from ITS and cpDNA in the fasta format are available via Figshare repository, https://doi.org/10.6084/m9.figshare.21710519.The ITS sequences used for described C. hissarica Kamelin ex Rasulova and two additional subspecies within the C. lehmanniana complex, i.e., C. lehmanniana subsp.integerrima Kamelin ex Rasulova and C. lehmanniana subsp.pseudohissarica Kamelin ex Rasulova, distinguished from C. lehmanniana and C. capusii (Franch.) for the preliminary identification of bellflowers, revised and collected during field studies.Fully developed specimens representing all taxa distinguished within the Campanula lehmanniana complex were studied under a light microscope.Plant material for population genetic analyses was collected from the Tian Shan (Talassky Alatau, Fergana Mts, Chatkalsky Alatau) and Pamir-Alai Mts (Alai Mts and Zeravshan-Hissar Mts including the Hissar Mts, Zeravshan Mts, and Turkestanian Mts).Leaf tissues from 262 specimens were used for molecular analyses.Of this number, 70 samples were used for the phylogenetic analyses, including 67 of C. lehmanniana s.l. as well as 3 of C. incanescens (outgroup).The remaining 247 samples from 31 localities, including 47 C. eugeniae samples (5 populations) and 200 C. lehmanniana s.l.samples (26 populations), were used for population molecular analyses.All voucher specimens used in molecular analyses are listed in Table Phylogenetic trees were largely consistent in topology regardless of analysed DNA markers or the method used for reconstruction.Trees based on cpDNA and DArT SNPs revealed two strongly supported clades with samples from the Tian Shan Mts, morphologically corresponding to the description of C. eugeniae, and also individuals from the Alai Mts and Zeravshan-Hissar Mts, morphologically identified as C. lehmanniana (Figures 1 and S2b; Appendix S1).The tree derived solely from the ITS region (Figure S2a) revealed four polytomic clades corresponding to the geographical distribution (particular regions or river valleys) of analysed samples (Tian Shan and Alai-Zeravshan-Hissar Mts).The Tianshanian clade of C. eugeniae includes two subclades reflecting a geographic pattern of population distribution.Subclade C1A comprises C. eugeniae populations from the Fergana Mts, while subclade C1B comprises populations from the Talassky Alatau, both isolated by the fairly wide Naryn River valley.The largest clade from the Alai Mts and Zeravshan-Hissar Mts refers to C. lechmanniana s.l. and consists of six subclades (C2, C2/3, C3, C4, C5, C6) with morphologically highly variable individuals, partially reflecting the geographical pattern of distribution.Populations within subclade C2 are distributed in the Alai Mts, subclade C2/3 consists of individuals from the Hissar Mts (mainly the Varzob River valley), subclade C3 contains populations from the Turkestan Mts located in the Yangiaryk River valley, whereas subclades C4-C6 comprises populations from the Zeravshan Mts (C4 from the Yagnob River valley; C5 from the Pastruddarya and neighbouring Iskanderdarya River valleys, and C6 from Kshtut, Urech and neighbouring Seven Lakes Valley; Figure 1).The molecular clock analysis based on concatenated cpDNA and ITS data, as well as the SNAPP analysis based on SNP data, suggest that C. lehmanniana and C. eugeniae split from a common ancestor during the late Miocene ca.5.8 or 5.64 Mya (based on cpDNA and ITS or DArTseq-based SNP, respectively), or even later, in early Pliocene if we consider 95% HPD intervals (4.16-7.41Mya).The split of the C. eugeniae clade within Tian Shan occurred in the Pleistocene, ~1.52 Mya (according to SNAPP) or ~1 Mya (according to the molecular clock based on ITS and cpDNA), while the C. lehmanniana clade split within the Zeravshan-Hissar and Alai Mts took place in the early Pleistocene ~2.69 or 1.68 Mya according to molecular clock based on ITS and cpDNA and SNAPP, respectively (Figure 1b,c).All subclades within the Zeravshan-Hissar Mts clade (C2/3-C6) diversified in the Pleistocene (Figure1b,c).Based on SNAPP, we estimate that this event occurred during interglacial periods of the Middle Pleistocene (from 1.1 to 0.07 Mya), whereas molecular clock based on ITS and cpDNA results suggest that most of the clades split in the early Pleistocene, ~2.28-0.71Mya (Figure1b,c).In SNAPP, the three topologies contained in the 95% HPD tree set differed only between subclade C5 and C6 relationships inferred for C. lehmanniana (Figure1c).
tween C4 and C5, as well as C5 with C6.These admixtures support unresolved population relationships within C. lehmanniana revealed by SNAPP among C5 and C6 clades located in the same geographic region.ADMIXTURE suggests the existence of three geographically clustered genotypes within C6, and two within C5, with signs of admixture between the geographically closest groups of populations, whereas within the C. eugeniae clade, the fastSTRUCTURE showed only one cluster, and ADMIXTURE was consistent with F I G U R E 1 Phylogenetic relationship estimated from (a) Neighbour Joining phylogenetic tree based on Euclidean distance (DArT, 9158 SNPs, 247 samples), (b) molecular clock rooted with Campanula incanescens (concatenated ITS and cpDNA, 6331 bp, 70 samples) and (c) SNAPP (3345 loci, 20 samples).Campanula lehmanniana in blue, C. eugeniae in yellow.Codes for specimens are explained in Table

Fastsimcoal2
analyses performed for 3 groups of populations from 3 distinct mountain ranges (Alai Mts, Zeravshan-Hissar Mts and Tian Shan Mts) support scenario assuming intra-and interspecific gene flow after the split due to secondary contact between groups of populations from Alai and Zeravshan-Hissar, Alai and Tian Shan, as well as Zeravshan-Hissar and Tian Shan Mts.(Figure 4).Taking into account the 2-year generation time of these species, it is estimated that the admixture between C. lehmanniana from the Alai Mts and Zeravshan-Hissar Mts took place relatively recently, over 3000 years BP (Figure 4).Same as the admixture between C. lehmanniana from the Alai Mts and C. eugeniae from the Tian Shan mountains which occurred ca.1000 years BP.The oldest admixture event is estimated between the most geographically distant populations of C. eugeniae from the Tian Shan Mts and C. lehmanniana from the Zeravshan-Hissar Mts, which took place over 125,000 years BP when climate began to cool (Figure 4).
Stairway plot revealed that 4 groups of populations (C1A including populations27, 28, 29; C2 including populations 23, 24, 26; C4 including populations 12, 13, 14; and C6C including populations 19, 15,   20)  shared highly similar pattern of historical fluctuations in effective population size (Ne) over time that roughly mirrored known events of climatic upheaval.Those groups of populations underwent long-term steady expansion during the Pleistocene and noticeable decline in Ne after the LGM with a large population shrinking just after the LGM or at the beginning of the Holocene.A different pattern, characterised by an increase of Ne, is shown by the most northern group of populations F I G U R E 2 Genetic structure estimated from (a) fastSTRUCTURE and (b) ADMIXTURE (9054 SNPs, 247 samples).Codes for specimens are explained in Table 5 arcminutes under MPI-ESM-P, respectively.The climatic variables F I G U R E 3 Analysis of introgression based on (a) ABBA/BABA tests shown on a Neighbour Joining tree (5149 SNPs, 250 samples, including Campanula incanescens as an outgroup).Arrows are coloured according to the admixture proportion in percentage based on the f4-ratio (multiplied by 100; see Section 2); (b) Maximum likelihood tree inferred with treemix (1300 SNPs and 247 samples) showing the most likely migration event (m = 4).The direction of gene flow is represented with an arrow coloured according to the percentage of alleles (weight) originating from the source.Codes of clades according to Figure 1.F I G U R E 4 Best-supported model of inter and intraspecific gene flow inferred by using fastsimcoal2.Parameters include ancestral (θ A-ZH , θ A-TS , θ TS-ZH ) and contemporary (θ A , θ ZH , θ TS ) effective population sizes, the timing of admixture (Tm 1 , Tm 2 , Tm 3 ) and migration rates (m 1 , m 2 , m 3 ) between different pairs.The table presents point estimates as well as lower and upper 95% confidence intervals (in brackets) for each parameter for the most likely species divergence model inferred with fastsimcoal2; A, Alai Mts; m, migration rate; Tm, migration time (measured in generations); TS, Tian Shan Mts; ZH, Zeravshan-Hissar Mts; θ, population effective sizes.F I G U R E 5 Stairway plots showing changes in effective population size (Ne) through relative time; demographic history inference with folded SNP frequency spectra.Estimates are based on the assumption of 7 × 10 −9 per site per generation per locus and a generation time of 2 years.Clade numerations (C1-C6) according to Figure 1 with populations of Campanula lehmanniana (in blue) and C. eugeniae (in yellow).

F I G U R E 6
Diagrams of bioclimatic variables oscillations for current and Last Glacial Maximum (medium of CCSM4 and MPI-ESM-P) climatic conditions; yellow shows the populations of Campanula eugeniae, whereas blue shading-C.lehmanniana.
rather than long-distance pollen dispersal or hypothetical corridors of gene flow.The low-dispersal capacity of the rocky-related bellflowers and the considerable topographic complexity (i.e., high ridges and lack of suitable habitats consisting of solid calcareous rocky substrates with small ledges, crevices, cracks, and fissures on steep walls) within the Tian Shan and Pamir-Alai Mts might explain the low or absent genetic admixture both among geographically separated populations as well as between populations located even within the same river valley (e.g., within Seven

F I G U R E 7
Maps of current (a, b) and Last Glacial Maximum (LGM) (c, d; e, f) climate-suitable areas for Campanula lehmanniana (a, c, e) and C. eugeniae (b, d, f); predicted distributions for the LGM based on the CCSM4 (c, d) and MPI-ESM-P (e, f) general atmospheric circulation models.Models with a spatial resolution of 2.5 arc-minutes (approximately 5 × 5 km).The high genetic differentiation and observed regional differences in phenotypic traits prove that the Quaternary climatic oscillation and repeated cycles of mixing-isolation-mixing(Abbott, 2019;Vintsek et al., 2022) greatly influenced the populations of the examined bellflowers.Long-lasting isolation within high-mountain rocky-islands of different regions during the interglacials induces only limited seed and insect-mediated pollen dispersal between populations, and very limited or even completely absent gene exchange between particular populations from different valleys.However, our analyses revealed that some population connectivity with rare hybridisation events occurred during the last glaciations.The colder climate during these periods forced high-elevated glaciers to advance hundreds of meters below their present locations(Abramowski et al., 2006;Owen et al., 2002).Consequently, populations of C. eugeniae and C. lehmanniana went downslope.During the LGM, the potential distribution range of Tianshanian species was much wider than is observed today (Figures 7, S3 and S4), and partially overlapped with the range of Hissaro-Zeravshanian species, leading to gene flow between them.Based on our analyses (Admixture, ABBA/BABA), population C2/3 of C. lehmanniana refers to the morphological description of both C. hissarica and C. lehmanniana subsp.pseudohissarica, arising from hybridisation between different C. lehmanniana genetic clades as well as a small admixture of C1 population of C. eugeniae (fast-STRUCTURE, Figure 2), which had a much wider distribution range during the LGM (Figure 7).Furthermore, we detected that the C1 clade of C. eugeniae also had a secondary, but weak (<10%) gene flow with C4, C5, and C6 genetic clades during the LGM.Similar gene flow was detected between currently isolated populations of Stipa zeravshanica occurring in rocky habitats of the Zeravshan Mts (Vintsek et al., 2022).Based on Fastsimcoal2 analysis, it is estimated that the most ancient gene flow occurred between C. eugeniae and C. lehmanniana almost 130,000 years BP, during the secondary contact of the ranges of the two taxa (Figure S5) when the temperature started to cool just before the last glaciation.It is also worth noting that contrary to ABBA/BABA results, Fast-simcoal2 analysis shows bidirectional gene flow from the Alaian population of C. lehmanniana to both C. eugeniae from the Tian Shan Mts and C. lehmanniana from the Zeravshan-Hissar Mts (C2 subclade).It is estimated that such gene flow took place ca.1000 and 3000 years BP, respectively, possibly during cooling associated with temperature fluctuations in the Holocene (Figure4).These results are in line with findings of other authors(Hewitt, 2000;Muellner-Riehl, 2019; Tonzo & Ortego, 2021), highlighting that the response of alpine species from different mountain regions of the temperate zone are similar and reflect both past climatic oscillations and contemporary population fragmentation in island-like habitats.A similar phylogeographic pattern could be applicable to other co-distributed and ecologically similar alpine species from Central Asia, supporting the general paradigm that cold-adapted species thrive during cooler intervals and are limited during climate warming(Vintsek et al., 2022).

4. 3 |
Demographic history and paleo and present niche modelsEstimating long-term demographic history is essential to clarify the genetic background of a species.Demographic history dynamics inferred from genetic data provide useful information to understand how climatic oscillations affected population fitness and viability(Selwood et al., 2015).There are two main demographic responses to climate changes of cold-adapted species: (a) increase in population sizes during glaciations and a decrease in interglacials mainly due to temperature rise, the so-called 'interglacial refugia model'(Bennett & Provan, 2008;He et al., 2015;Ortego et al., 2021;Ortego & Knowles, 2022; Tonzo & Ortego, 2021;Vintsek et al., 2022;Yang et al., 2022), and (b) decrease in population sizes during glaciations and an increase in population sizes during interglacials (expansion from glacial refugia).The second one is mainly due to the rise in temperature and subsequent increase of suitable habitat areas, typical for lower elevation species(Fernandez et al., 2020;Guerrina et al., 2022;Wielstra et al., 2013),andmore exceptionally, for plants from high-mountain plateaus (Muellner-Riehl, 2019; Wróbel et al., 2023).Our results for the demographic population reconstruction are in line with the first mentioned above response and show that many of analysed populations of C. eugeniae and C. lehmanniana underwent parallel and substantial changes in effective population size with similar trajectories, and right after LGM, they experienced an extraordinary decline of effective population size.These results are linked with those obtained in analyses of past and current distribution models, where the distribution ranges of both examined species (Figures 7, S3 and S4) clearly decreased, affected by changes in climatic conditions.The temperature increase after glaciation (Figure 6) resulted in the disappearance of potential distribution ranges in the western Pamir for C. eugeniae and in Darvaz and the southern Hissar Mts for C. lehmanniana.According to our analyses, precipitation in the coldest and driest quarters (bio17, bio19) has the strongest impact on the past and present distribution of the examined species (Figure 6), which is in line with findings of other studies supporting the strong effect of precipitation on vascular plant distribution in the Central Asian mountains (Nowak et al., 2022; Vintsek et al., 2022).The bioclimatic factors affecting the distribution of Campanula lehmanniana and C. eugeniae potentially underlie also the occurrence of many other closely related and endemic plant species as well as their plant associations exhibiting similar patterns of geographical division between the Tian Shan and western Pamir-Alai Mts.For instance, such geographical vicariants also include: Asperula laevis, Carex litvinovii, Scutellaria andrachnoides, Silene litophila, Spiraea ferganensis, Stipa gracilis distributed in the Tian Shan, and respectively, Asperula chukavinae, A. oppositifolia, Carex koshewnikowii, Scutellaria megalodonta, Silene samarcandensis, Spiraea baldshuanica, Stipa zeravshanica distributed in the western Pamir-Alai Mts.However, in the Alai Mts (northern Pamir-Alai region), the precipitation in the coldest and driest quarter is more similar to that of the Tian Shan rather than to the Zeravshan-Hissar region.Thus, these specific environmental conditions most probably forced populations of C. lehmanniana (C2) occurring in that region to genetic adaptation.

5
| CON CLUS IONS Our comprehensive study, including phylogenomic, demographic, and distributional analyses for C. lehmanniana complex, shows that past climatic oscillations and contemporary population fragmentation in island-like habitats may rule the phylogeography of alpine species.The Quaternary climatic oscillations resulting in cyclically recurring episodes were the main driving force for currently observed genetic patterns.The presence of strong genetic structures within currently existing populations reflects the strong geographical isolation within both C. eugeniae and C. lehmanniana populations, with very limited gene flow between them.The parallel demographic response of both analysed species to paleoclimatic changes may outline the general response of other cold-adapted chasmophytes to ongoing climate warming.The contemporary challenge for nature conservation in the alpine biodiversity hotspots should be to establish new protected areas, including micro-hotspots of genetic variation as well as to create a map of regions with potentially replacement habitats for transferring selected populations of plant species, most sensitive to climate change.Particular emphasis should also be placed on preserving the genetic resources of such species in gene banks(Díez et al., 2018;Peres, 2016) as well as establishing ex situ breeding collections in local botanical gardens.ACK N O WLE D G E M ENTSTo the curators of AA, B, K, KRA, FRU, TAD, OPUN, LE, P and MW herbaria, for making available the collection of the genus Campanula and its allies for our studies.To Dr. Polina D. Gudkova for the help during the fieldwork, to Georgy Lazkov for giving us samples of Campanula for phylogenetic studies.To Dr. Anna Rucińska for the discussion of molecular analysis used in the study.We are grateful also to the Editor Dr. Christopher Burridge and the anonymous reviewers for providing valuable suggestions and improvements to the previous version of the manuscript.The data utilised in the studies have been lawfully acquired in accordance with The Nagoya Protocol.FU N D I N G I N FO R M ATI O NThe research was supported by National Science Centre, Poland, Grant Number: 2018/29/B/NZ9/00313.
(Hartmann et al., 2013)wles, 2022; Tonzo & Ortego, 2021).Our results show that C. lehmanniana and C. eugeniae also demonstrate a parallel negative response to past climatic oscillations, and since a similar population decline pattern was observed also in Stipa zeravshanica and S. gracillis (Vintsek et al., 2022), we may extrapolate it for other Tianshanian and western Pamir-Alaian chasmophytes.Currently, high-mountainareas experience more rapid changes in temperature than lower elevations(Hartmann et al., 2013).The phenomenon, in which the warming rate of air temperature is intensified with elevation, is called elevation-dependent warming (EDW); however, the rate of EDW is diverse in different mountain areas