Phylogenomics unravels Quaternary vicariance and allopatric speciation patterns in temperate‐montane plant species: A case study on the Ranunculus auricomus species complex

The time frame and geographical patterns of diversification processes in European temperate‐montane herbs are still not well understood. We used the sexual species of the Ranunculus auricomus complex as a model system to understand how vicariance versus dispersal processes in the context of Pleistocene climatic fluctuations have triggered speciation in temperate‐montane plant species. We used target enrichment sequence data from about 600 nuclear genes and coalescent‐based species tree inference methods to resolve phylogenetic relationships among the sexual taxa of the complex. We estimated absolute divergence times and, using ancestral range reconstruction, we tested if speciation was enhanced by vicariance or by dispersal processes. Phylogenetic relationships among taxa were fully resolved with some incongruence in the position of the tetraploid R. marsicus. Speciation events took place in a very short time at the end of the Mid‐Pleistocene Transition (830–580 thousand years ago [ka]). A second wave of intraspecific geographical differentiation occurred at the end of the Riss glaciation or during the Eemian interglacial between 200 and 100 ka. Ancestral range reconstruction suggests a widespread European ancestor of the R. auricomus complex. Vicariance has triggered allopatric speciation in temperate‐montane plant species during the climatic deterioration that occurred during the last phase of the Mid‐Pleistocene Transition. Vegetation restructuring from forest into tundra could have confined these forest species into isolated glacial macro‐ and microrefugia. During subsequent warming periods, range expansions of these species could have been hampered by apomictic derivatives and by other congeneric competitors in the same habitat.


| INTRODUC TI ON
Pleistocene climatic fluctuations have caused periodic range shifts (e.g., north-south) and/or range expansions and contractions in northern hemisphere plant species. The alternation of cold and warm periods has characterized the entire Quaternary (the last 1.8 million years). In Europe, cold periods have produced southward shifts in the distribution ranges of temperate species (in the so-called "Mediterranean refugia"; Brewer, Cheddadi, de Beaulieu, Reille, & Contributors, 2002;Taberlet, Fumagalli, Wust-Saucy, & Cosson, 1998) and/or the survival of these in few, isolated central European refugial areas (e.g., Magri et al., 2006;Naydenov, Senneville, Beaulieu, Tremblay, & Bousquet, 2007). On the other hand, warmer interglacials have favoured range expansion towards the north and the formation of contact zones among diverged populations of the same species (Hewitt, 1999;Magri, 2008). As a result, the survival of formerly coherent population groups in isolated refugia promoted allopatric speciation. Range expansions and the consequent formation of secondary contact zones resulted in hybridization, and eventually in the formation of new taxa through homoploid hybrid speciation or allopolyploidization (Abbott et al., 2013;Stebbins, 1984).
Much less information is available on montane and/or subalpine plant species of the northern hemisphere. Research of forest plants has focused on the dominating tree species, such as beech, oaks and pine (e.g., Brewer et al., 2002;Magri, 2008;Naydenov et al., 2007).
In cold periods of the Pleistocene, large parts of Central Europe were dominated by steppe vegetation, and deciduous forests were mostly restricted to southern and main refugia in the Balkan and Iberian peninsulas (Stuchlik & Wójcik, 2001). However, more distant local microrefugia outside the main southern glacial refugial areas should also be considered (Rull, 2009). For the boreomontane species Polygonatum verticillatum, putative glacial refugia might have existed in the foothills of heavily glaciated mountain ranges, such as the Alps and the Carpathians (Kramp, Huck, Niketić, Tomović, & Schmitt, 2008). In warmer periods, the ecoclimatic barriers for forests disappeared, leaving only the highest chains of the Alps and Carpathians as a geographical barrier. Hence, reforestation of the lower elevational zones around these mountain areas and in the lowlands has provided ample opportunities for herbaceous understorey plants to expand their range via geodispersal. In contrast to "jump" dispersal over a persisting barrier, geodispersal describes a more or less continuous range expansion after disappearance of a barrier, and would result in a single, widespread species (Ronquist & Sanmartín, 2011;Sanmartín, 2012). Geodispersal was found to be important for boreotropical forest plants (Couvreur et al., 2011), but little is known for temperate forest plants. Zhao, Wang, Ma, Liang, and He (2013), for example, found vicariance and allopatric speciation in a temperate deciduous forest plant in China.
The cosmopolitan genus Ranunculus L. with ~600 species originated and started diversifying in the early Miocene . Emadzade, Gehrke, Linder, and Hörandl (2011) have elucidated how both vicariance and long-distance dispersal have to be invoked to explain the present distribution patters of lowland and montane species in the genus. Intercontinental and transoceanic dispersal of forest species occurred mostly in the Ranunculus polyanthemos clade . Other clades comprising forest species, such as R. sect. Auricomus, have not yet been fully resolved. The Eurasian Ranunculus auricomus complex is nested as a monophyletic group within this large clade (R. sect. Auricomus) that comprises otherwise Asian, Arctic and North American species (Emadzade et al., 2015). Whereas apomictic polyploids of the R. auricomus complex are widespread and abundant in various temperate and boreal biomes, by contrast sexuals are montane and subalpine species of Mediterranean and temperate European mountains, inhabiting deciduous forest ecosystems or natural and anthropogenic meadows. Although these forests are widespread over temperate Europe, these taxa have restricted or disjunct distribution ranges, and usually occur at the foothills of the main European mountain ranges ( Figure 1). Both vicariance or jump dispersal could have shaped this pattern. It is also unclear why speciation occurred despite the lowland areas being connected during interglacial periods and providing habitats for geodispersal.
The whole R. auricomus complex shows a pattern of geographical parthenogenesis: sexual progenitors have much smaller and more southern distribution areas than the numerous polyploid apomictic derivatives (Bierzychudek, 1985;Hörandl, 2006). Apomictic lineages can expand their range rapidly due to uniparental reproduction, both via long-distance dispersal and geodispersal (Baker, 1967;Cosendai, Wagner, Ladinig, Rosche, & Hörandl, 2013;Hörandl, 2006;Kirchheimer et al., 2018). However, why sexual progenitors of such complexes do not manage to expand their distribution areas and fail to fill their potential range remain poorly understood (see e.g., simulations in Kirchheimer et al., 2018). Hence, it is important to disentangle the biogeographical histories and speciation processes of sexual species from their apomictic derivatives to understand the phenomenon of geographical parthenogenesis.
In the present study, we address three major questions concerning the spatiotemporal diversification of temperate-montane species in Europe, using as study object the sexual species of the R. auricomus complex. (a) Did colder periods during glaciation cycles enhance speciation processes in temperate-montane species? (b) Is vicariance, dispersal or geodispersal the most probable scenario to explain the present distribution of sexual taxa of R. auricomus? (c) Did the species survive glacial periods in southern macrorefugia or also in microrefugia, and (if the latter scenario is the most plausible) why are they still restricted to these areas? Therefore, we: (a) applied target enrichment of nuclear genes to resolve phylogenetic relationships among taxa of the complex; (b) use coalescent-based methods to estimate the species divergence times in the context of past climatic oscillations; (c) apply ancestral range reconstruction methods to test vicariance versus dispersal speciation scenarios and identify potential macro-and microrefugia; and (d) discuss the relevance of our findings for understanding the geographical parthenogenesis pattern.

| Plant material
We included samples from all the described sexual species of the Ranunculus auricomus complex. Plant material used for DNA extraction consisted of silica-gel dried leaves or herbarium specimens (Table 1). We follow the taxonomy of Karbstein et al. (2020) (see synonymy in Table 1). Two additional Ranunculus species not belonging to the R. auricomus complex were included in the analyses as an outgroup: R. pygmaeus from the clade of Ranunculus sect. Auricomus, and R. sceleratus from the next sister clade .

| Gene selection and probe design
To find single-copy genes and select target regions for phylogenetic analyses, we used transcriptomes from two diploid species (R. notabilis and R. cassubicifolius, as synonym "R. carpaticola" in Pellino et al., 2013), a tetraploid sexual accession of R. cassubicifolius (Pellino et al., 2013), and R. brotherusii (Chen, Zhao, Wang, & Moody, 2015), an Asian species F I G U R E 1 Distribution map of the accessions included in the present study (a) and detailed map of specimens collected in Slovenia (b). Numbers correspond to the map IDs in Table 1. Different symbols are for different species, following the treatment given in Karbstein et al. (2020) of R. sect. Auricomus. markerminer version 1.0 (Chamala et al., 2015) was used to identify putative single-copy orthologous loci for probe design with Arabidopsis thaliana as reference. We selected exons found in at least two of the four transcriptomes, with lengths ranging from 120 to 960 bp and a minimum variability of two single nucleotide polymorphisms (SNPs)/120 bp, using a custom python script (available at: https://github.com/Claud iaPae tzold /Marke rMine rFilt er.git). We identified 2,628 exonic regions belonging to 736 target genes (Table S1).
The selected target regions (genes) ranged from 121 to 5,291 bp and included at least one exonic fragment each.
Arbor Biosciences produced a MYbaits Target Enrichment kit with 20,000, 120-bp-long, in-solution, biotinylated baits based on target sequence information. The final bait panel consisted of 17,988 probes, 14,632 of which are unique, and tiling at 2 × density.

| Library preparation
Sequencing libraries were prepared using either the "NEBNext for details on library preparation).

| Read processing und alignment
We checked the quality of raw reads with fastqc (available at: http:// www.bioin forma tics.bbsrc.ac.uk/proje cts/fastqc). Further processing of the raw reads was done using the pipeline hybphylomaker (Fér & Schmickl, 2018). Sequence adapters were removed, and reads were quality-trimmed using trimmomatic version 0.32 (Bolger, Lohse, & Usade, 2014), with the default settings used in hybphylomaker. Duplicated reads were removed with fastuniq version 1.1 (Xu et al., 2012). Quality-trimmed individual raw reads were mapped to a reference sequence and then merged into contigs that are aligned for each gene separately. As pseudoreference for read mapping, we used a sequence consisting of the concatenated target exonic sequences, separated by stretches of 800 Ns. Mapping to the pseudoreference genome was done with bwa (Li & Durbin, 2009), and consensus sequences of the mapped reads were produced with consensusfixer (available at: https://github.com/cbg-ethz/Conse nsusF ixer), because this is the only approach available in hybphylomaker able to call ambiguity DNA codes in case of multiple bases per site in the mapped reads. For consensusfixer we used the following settings: minimum relative abundance of the alternative base ("plurality" in the setting file of hybphylomaker) of 0.2 and a minimum read coverage for ambiguity calling ("mincov") of 5.
Consensus sequences were matched to sequences of the target exons to produce PSLX files using blat (Kent, 2002 Note: Accepted names following Karbstein et al. (2020) are in bold type, and younger synonyms are in inverted commas. Countries are abbreviated according to ISO code 3166-2. Altitude is given in metre above sea level (m a.s.l.). Range codes are those used in biogeobears for the ancestral range reconstruction and correspond to those in Figure 1. Map IDs refer to the samples IDs given in Figure  Samples from type material or collected from "locus classicus" and nearby.

TA B L E 1 (Continued)
therefore combined to produce exon-wise matrices with "assem-bled_exons_to_fastas.py" (Weitemier et al., 2014). Matrices were aligned with mafft version 7.029 (Katoh & Standley, 2013), using the default program settings, and then gene-wise concatenated with amas (Borowiec, 2016). We filtered potential paralogues by running the script "HybPhyloMaker4a2_selectNonHet.sh" and setting the maximum number of heterozygous sites per locus ("maxhet" in the hybphylomaker settings file) to 10. We continued the pipeline with both data sets, filtered and unfiltered. hybphylomaker performs two consecutive steps of missing-data filtering. First, sequences with more than a certain percentage of Ns in an alignment ("missingpercent" in the settings file) are deleted. We set this option to 40.

| Phasing
Because the importance of retrieving allele information for a correct phylogeny estimation (especially in a recently diverged group) has been emphasized in recent studies (Andermann et al., 2019;Eriksson et al., 2018), we tested the impact of allele phasing on our phylogenetic analyses. To avoid the loss of allelic information during the process of allele mapping and consensus sequence production, we extracted the.bam files produced after mapping the reads to the pseudoreference sequence from the hybphylomaker pipeline.
These were phased with samtools version 0.1.19 ), using the combination of the commands "samtools sort" and "samtools phase." The phased.bam and.bai files were then placed in a new hybphylomaker working directory for further processing within the pipelines workflow. A new "/10rawreads/" folder was also placed in the hybphylomaker working directory, containing only the file with the modified sample names. The pipeline was therefore resumed for computing of the consensus sequences (this time allele-wise consensus sequences) by calling the script "HybPhyloMaker2_readmapping.sh" but specifying "mapping = no" in the setting file. The alignments obtained are hereafter called allele alignments. We ran the pipeline further as explained above, with the only differences that (a) matrices of exons belonging to the same gene were not concatenated; and (b) exon trees (instead of gene trees) were inferred and used for further analyses.

| Position filtering
We evaluated the effect of filtering alignments for positions that could add phylogenetic noise to the phylogenetic reconstructions, the so-called "phantom" spike positions (ambiguous positions, usually situated close to indel-rich regions of the alignment, with abnormally high substitution rates). We followed the procedure illustrated by Fragoso-Martínez et al. (2017), calculating first the Phylogenetic Informativeness (Townsend, 2007) and net Phylogenetic Informativeness profiles using hyphy (Pond, Frost, & Muse, 2005) in the web portal PhyDesign (López-Giráldez & Townsend, 2011). As tree inputs, we used trees (different trees for different data sets; consensus or allele alignments, paralogure filtering or not, etc.) obtained from hybphylomaker by concatenating gene alignments and fasttree (Price, Dehal, & Arkin, 2010). The phylograms were then transformed to ultrametric trees using Penalized Likelihood methods (Sanderson, 2002) as implemented in the r package "ape" (Paradis & Schliep, 2018), with correlated rates and smoothing parameter (λ) set to 1. The relative timescale was set to 1 at the root.
Second, we used the estimated substitution rate per locus from  Table S2. The cleaned alignments were then reintroduced in the hybphylomaker pipeline, just before the missing-data filtering step.

| Gene tree and species tree/network analyses
For all different data sets, species trees were estimated (a) using maximum likelihood (ML) as implemented in raxml version 8.2.4 (Stamatakis, 2014) on the concatenated data sets (concatML) and (b) applying the coalescent-based method astral-iii (Chao, Rabiee, Sayyari, & Mirarab, 2018). For the concatML analyses, concatenated alignments were produced in hybphylomaker using amas and saved in phylip format. Analyses were run on the Scientific Computer Cluster of the GWDG (https://www.gwdg.de/) using raxml, with alignments partitioned by genes, the GTRGAMMA model and 100 bootstrap replicates. Concatenated analyses were only performed for the consensus alignments, because there was no way of knowing how the phased alleles of a sample combined across loci. We also quantified branch support and conflict on the concatML trees with the Quartet Sampling method (Pease, Brown, Walker, Hinchliff, & Smith, 2018).
For the coalescent-based analyses, first gene trees (or exon trees in the case of the allele alignments) were reconstructed with raxml. Analyses were run with 100 standard bootstrap replicates, with the GTRGAMMA model and partitioning by exons. Gene tree characteristics (average bootstrap support, average branch length, etc.) and correlations among alignments and gene tree characteristics were calculated and plotted in r version 3.5.2 (R Core Team, 2018), using modified scripts from Borowiec (2016) as implemented in hybphylomaker. Gene trees were rooted, collapsing branches with bootstrap values <50, and combined into single newick files using Newick Utilities (Junier & Zdobnov, 2010). Species trees were inferred applying the coalescent-based algorithm implemented in astral-iii version 5.6.3 (Chao et al., 2018) with 100 multilocus bootstrap replicates. To assess the amount of gene tree conflict on branches, we measured the quartet support on the astral trees (Sayyari & Mirarab, 2016) for the main, the first alternative and the second alternative topologies ("-t 8" option in astral).
In addition, and to test for the extent of reticulate evolution among the sexual representatives of the R. auricomus complex, we inferred neighbor-net networks on the concatenated data sets using splitstree version 4.14.6 (Huson & Bryant, 2006

| Divergence time estimation
To obtain absolute divergence times, we used two secondary calibration points for the analyses. The crown age of the R. auricomus complex was set according to Hörandl (2004 The clock rate was derived from standard substitution rates for plant nuclear genes as in Pellino et al. (2013). Assuming a standard substitution rate of 5e-9 (Wolfe, Gouy, Yang, Sharp, & Li, 1989) and a generation time of 3 years for members of the R. auricomus complex (Pellino et al., 2013), the beast clock rate translates to 1.66e-6 mutations per thousand years. Formerly described species were treated as independent lineages. Tetraploid accessions of R. cassubicifolius were treated separately from the diploid representatives of the species.
Because the main goal of the analyses was to estimate the time frame of speciation events, and because these analyses were based only on a reduced number of loci, we used a soft-constrained approach on the species tree topology in order to reach faster convergence and to make the topology of the *beast tree similar to the topology of the species trees obtained with the other methods. In particular, the monophyly of the clade with the Illyrian taxa was enforced.
Substitution models for each locus were adjusted to those found with jmodeltest, and the Birth-Death model was used as prior on the species tree. Two independent analyses were run for 2 × 10 9 generations, sampling every 100,000th generation. Effective sample size (ESS) values and convergence between independent analyses were checked in tracer version 1.6 (Rambaut & Drummond, 2007).
Results of the two analyses were merged using logcombiner version 2.5.2 (Bouckaert et al., 2019) applying a burn-in of 10%. Finally, the remaining 18,000 trees were used to construct a maximum clade credibility tree with a posterior probability limit set to 0.5 and "Mean Heights" for node heights using treeannotator version 2.5.2 (Bouckaert et al., 2019). Another set of analyses was performed as described above, but including only the ingroup taxa, as many authors have highlighted that using outgroups in coalescent-based analyses violates some of the model assumptions (Agudo, Tomasello, Àlvarez, & Oberprieler, 2017;Drummond & Bouckaert, 2015).
An additional set of analyses was performed for the above-mentioned data sets (including and excluding the outgroup) using only clock rates or calibrations (and then applying an exponential prior with mean equal to 0.01 for clock rates) to obtain absolute time estimates. This was done to evaluate the effect of using both calibration and clock rate priors in the same analysis, and to test for congruence of the resulting time estimates.

| Ancestral area reconstruction
To evaluate hypotheses about past distribution patterns and dispersal and/or vicariance events in the evolutionary history of the sexual representatives of the R. auricomus complex, we performed an ancestral area reconstruction with the r package biogeobears version 1.1.2 (Matzke, 2018). This approach allows several models of biogeographical evolution to be compared via likelihood-based model selection (Matzke, 2013). Because sexual R. auricomus species grow within or near the main European mountain ranges, we assigned them to the seven following areas: Northern Carpathian

| Gene capture and sequencing
The average number of obtained reads per sample was 1,373,158 (range 867,000-to 2,178,952). Read quality was relatively high, and on average only 3.27% of reads per sample were discarded after quality check (

| Alignments and the impact of the different filtering schemes
The number of alignments used for the phylogenetic analyses varied according to the filtering scheme applied (Table 2). For the "consensus data set," about 15% of loci were excluded after missing data filtering and more than 25% after paralogue filtering (more than 50% of the total alignment length). Concerning the "allele alignments," working with exons made it possible to retrieve more loci (more than 80% and 50% of the total alignment length before and after paralogue filtering, respectively). More details are given in Table 2.
Allele phasing generated alignments with a considerably higher proportion of PI sites compared to the consensus alignments (5.82% and 1.4% on average, respectively). Paralogue filtering resulted in a loss of a fraction of PI sites in both the consensus and the allele data sets. When selecting the 50 loci for the beast analyses, 14.02% (78,458 bp) of the total matrix length was retrieved, comprising >23% of the total PI sites (1,971). Lists of the alignments used in the different filtering schemes, their characteristics and substitution models of the loci selected for the beast analyses are given in Table S5.

| Phylogenetic analyses and species trees
Species trees obtained from the concatML and coalescent-based analyses of the consensus data set showed topological incongruences across most of the filtering scenarios. When using allele alignments, coalescent-based species trees obtained from different filtering treatments always showed the same topology, differing only slightly in support values and branch lengths. When taking the average local posterior probability (LPP) of the astral bootstrap trees as a criterion to compare different treatments, the paralogue filtering tool offered in hybphylomaker significantly worsened the results (in a post-hoc Tukey's test only differences between treatments with and without paralogue filtering were significant; Table S8), whereas position filtering improved support values slightly in some cases ( Figure S1). No significant differences were noticed between "consensus" and "allele" data sets. Accordingly, and as a matter of simplicity, only results of the "consensus" and "allele" analyses, without

| Age estimation
Results from the two *beast analyses (including and excluding the outgroup) were congruent to a large extent, varying only slightly on the age TA B L E 2 Alignment characteristics and performance of different filtering scenarios for the consensus and the allele data set  estimates. As for the concatML and the astral species tree, relationships within R. notabilis s.l. and within R. cassubicifolius s.l. are blurred and phylogenetic relationships receive low support (Figure 4; Figure S9).
The crown age of the R. auricomus complex is estimated to be ~733.5 thousand years ago (ka) ( Table 3), somewhat older (~750 ka) in the analyses including outgroups ( Figure S9). Most of the speciation events occurred between 830 and 580 ka. The only exception is the Apenninian R. marsicus, which diverged from R. notabilis s.l. more recently (266-107 ka). However, the latter result does not seem to be supported by the concatML and by the species tree analyses, in which R. marsicus is rather found to be related to R. envalirensis (but see support values in Figures 2 and 3), and its divergence seems to have occurred much earlier. Results of analyses including only calibrations or clock rates were congruent to those of the main analyses (Table S9).
Differentiation within R. notabilis s.l., R. envalirensis s.l. and R. cassubicifolius s.l. occurred in the last 200 thousand years, with the crown ages of these three species estimated as 176, 147 and 130 ka, respectively (166, 128 and 119 ka in the analysis including outgroups; Figure S9).

| Ancestral area
The LRT shows that results of models including the x parameter were significantly better (p < .05) than those without this parameter (Table S10). The j parameter did not improve the model fit for DEC and DIVA-LIKE, whereas it considerably ameliorated the BAYAREA-LIKE model fit. Overall, the best-fitting model was the DIVA-like + x (Table 4).

| Phylogenetic relationships among sexual species of the R. auricomus complex
Disentangling phylogenetic relationships among recently and rapidly diverged lineages remains a challenge (Knowles & Chan, 2008). In the present study, we resolved the phylogeny of the sexual representatives of the less than 1-Myr-old Ranunculus auricomus complex.
To this end, we analysed Target Enrichment data using concatenated ML and coalescent-based species-tree inference methods. Target Enrichment is considered to be a very suitable approach for resolving deep phylogenies, as the baits are usually designed for highly conserved regions (McCormack et al., 2013;McKain, Johnson, Uribe-Convers, Eaton, & Yang, 2018). Nevertheless, the potential of this method for resolving shallow phylogenies and application in population genetic studies has been highlighted recently (Harvey, Smith, Glenn, Faircloth, & Brumfield, 2016). Until now, this sequencing approach has been applied a few times to resolve phylogenetic relationships in recently radiating plant and animal groups (e.g., the 7.2-Myr-old genus Heuchera [Folk, Mandel, & Freudenstein, 2015]; the ~3-Myr-old carnivorous plant Sarracenia [Stephens et al., 2015]; the rapid island radiation of the Philippine Shrews [<4. In this sense, our study represents the most outstanding example of employing Target Enrichment data in recently and rapidly radiating species complexes. The main problems when working with rapidly radiating groups are: (a) rampant incomplete lineage sorting (ILS) due to the fact that ancestral polymorphisms are carried through a series of nearly simultaneous speciation events; and (b) the lack of information due to the slow mutation rate of the (most often) coding regions used in studies employing a Target Enrichment approach (Giarla & Esselstyn, 2015).
Using information from several (hundreds of) independent loci from F I G U R E 4 Maximum clade credibility tree estimated in *beast using the 50 most informative loci from the consensus data set. Blue bars indicate 95% highest posterior density (HPD) intervals of the age estimate (see Table 3). Numbers above branches indicate Bayesian posterior probabilities; only values above 0.7 are shown. On the right, bars indicate the current nomenclature according to Karbstein et al. (2020). The curve below the tree represents the δ 18 O modified from Lisiecki and Raymo (2005). Valleys correlate with cold periods (glaciations) and peaks with warm periods ( multiple individuals per species, and applying coalescent-based species tree methods should improve phylogenetic inference and overcome the above-mentioned problems (Kumar, Filipski, Battistuzzi, Pond, & Tamura, 2012). However, in Grummer, Morando, Avila, Sites, and Leaché (2018) even 580 loci were not enough to provide significant support for interspecific relationships in the L. fitzingerii species group. In R. auricomus, 500-600 loci were able to resolve phylogenetic relationships among the main lineages of the complex.
Incongruence was found mostly within the Illyrian group and within the clade of R. cassubicifolius s.l. There, relationships among formerly described species were inconsistent across analyses and clades did not receive significant support, confirming the broader species delimitation provided by Karbstein et al. (2020). Diversification processes within these two species (i.e., R. notabilis s.l. and R. cassubicifolius s.l.) took place recently and in a very short time frame (i.e., the last 176 and 136 kyr, respectively). We believe that the lack of phylogenetic divergence might be the cause of the poorly resolved intraspecific relationships in these two species. As already observed in comparable complexes (Grummer et al., 2018;Stephens et al., 2015  Note: Estimates for the parameters dispersal parameter (d), extinction (e), x (dispersal probability as a function of distance) and j (jump speciation events at cladogenesis). The favoured model is in bold type.
(ILS or reticulation producing incongruence in the data set) are not necessarily mutually exclusive and a definitive and unequivocal reconstruction of the processes giving rise to this polyploid taxon, which requires a different experimental design and more specific analyses (e.g., coalescent-based species network inference), is far beyond the scope of the present study.
Concerning differences among tree inference methods (concatenated vs. coalescence-based) and performance of data sets from different filtering settings, we note that concatenation usually resulted in higher bootstrap values (compared to equivalent coalescent-based analyses) even though topologies were more discordant among trees. Similar patterns have already been observed in empirical and simulation studies (Herrando-Moraira & The Cardueae Radiations Group, 2018; Kubatko & Degnan, 2007;Weisrock et al., 2012). Remarkably, trees obtained by applying coalescent-based analyses to the phased data set were topologically identical, and minor differences concerned only branch lengths and support values ( Figures S5-S7). This is further evidence that using allelic information together with coalescent-based species tree approaches can considerably increase the trustworthiness of phylogenetic analyses (Andermann et al., 2019). Concerning the performance of different phasing and filtering settings, the paralogue filtering implemented in hybphylomaker significantly worsened branch supports. Applying this option to diploids possibly resulted in the elimination of many nonparalogue (but variable) alignments.
Moreover, paralogy might not be an issue in our data set, because the baits kit was specifically designed for members of the R. auricomus complex. Therefore, it is rather unlikely that the selected single-copy regions went through duplication in some of these very closely related taxa. No significant amelioration in support values of the bootstrap astral trees was produced when applying phasing (but see above considerations on tree topologies) or position filtering ( Figure S1).

| Time frame of diversification processes
The crown age of the R. auricomus complex was estimated to be ~730 ka (Table 3). Hörandl (2004) estimated the divergence time between three sexual species using isozyme allelic frequencies and genetic distances after Nei (1975). Ages of these divergence events were estimated to 914,000 and 535,000 years, respectively, even though they describe the same split in the phylogeny of the complex (the crown age of the whole group). Therefore, using a broad prior for the crown age of the complex and a relatively informative prior for the clock rate (according to the mutation rates estimated for the complex by Pellino et al., 2013), we were able to estimate the age of this node more precisely. Moreover, we obtained the time frame of all speciation events among sexual representatives of R. auricomus.
Interestingly, almost all these events have taken place in a very short time between 830 and 580 ka (Figure 4). The only exception is the tetraploid species R. marsicus, which seems to have diverged much more recently (but see considerations about uncertainty above).
The divergence estimates for all these speciation events coincide with the Mid-Pleistocene Transition (MPT; Clark et al., 2006;Lisiecki & Raymo, 2005), one of the most important climatic changes occurring in the Pleistocene. During this interval (1.2 Ma to 650 ka), the relatively low-amplitude 41-kyr climate cycles of the earlier Pleistocene were progressively replaced by high-amplitude 100kyr cycles (Lisiecki & Raymo, 2005). Glaciations became longer and more severe, and the average global ice volume increased consistently, with larger increases at high and middle latitudes of Eurasia and North America (Elderfield et al., 2012;Raymo, Lisiecki, & Nisancioglu, 2006). This process became progressively stronger, and evidence from all northern continents indicates that Marine Isotope Stage (MIS) 22 (~870-880 ka) was the first of the major glaciation events characterizing the later Pleistocene (Head & Gibbard, 2005).
At the end of the MPT (MIS 16; ~650 ka), the most substantial glaciation yet experienced in the northern hemisphere took place (Head & Gibbard, 2005). This tremendous climatic transition was probably the consequence of the increase in atmospheric moisture Diversification processes within species took place more recently, most probably within the last two glaciation cycles. Diverging populations of R. envalirensis s.l., R. notabilis s.l. and R. cassubicifolius s.l. differentiated allopatrically (in the first case) or parapatrically (in the latter two cases) within the last 160 ka ( Figure 4). As for the alpine taxa, and probably also for montane plant species, the last two glaciation cycles were a triggering force shaping intraspecific genetic patterns, but were not sufficiently long to complete speciation processes. Overall, our age estimates are concordant with those found in previous studies Hörandl, 2004;Pellino et al., 2013).
In Hörandl (2004) and , diversification processes within R. cassubicifolius were found to be slightly older, 317 and 640 ka, respectively. In our analyses, sexual autotetraploids in R. cassubicifolius (Hörandl & Greilhuber, 2002) were found to be ~95,000 years old. This result is in accordance with previous studies estimating the maximum age of hexaploid apomictic derivatives as 80,000 years (Pellino et al., 2013), and confirms the general assumption that the origin of polyploids is associated with the last glaciation cycle.

| Biogeographical history
In our biogeographical analyses, AIC preferred models supporting the existence of a widespread ancestor for the R. auricomus complex, stressing the importance of vicariance in triggering speciation events among sexual representatives of the complex (Table 4). The ancestor of the whole complex was reconstructed to be distributed in the whole of Europe, although with some uncertainty ( Figure 5).
A first vicariance event separated populations with a northern dis-  (Qiu, Guan, Fu, & Comes, 2009). Surprisingly, the sexual species of the R. auricomus complex did not show long-distance dispersal, which has been observed in many other plant genera (Knapp et al., 2005), and also in other forest species from other clades of the genus Ranunculus . The achenes as diaspores exhibit no obvious specialized structures potentially suitable for long-distance dispersal via anemo-or zoochory. The absence of long-distance dispersal in the sexual R. auricomus species may also reflect limited abilities to adapt to new environments or novel niches. Outcrossing and self-sterility in the diploids may be other reasons for low colonization abilities (Hörandl, 2008). Nevertheless, moderate range expansions in warmer periods between 600 and 200 ka should be considered, especially in the R. cassubicifolius, R. envalirensis and R. notabilis clades, as their respective crown group areas are larger than their stem group areas (Figure 5a). Such geodispersal during interglacials might have potentially resulted in secondary contact hybridization and an allopolyploid origin of the apomictic taxa. However, range expansions remained at a regional scale and did not allow sexual taxa to cover again the whole ancestral area of the R. auricomus complex. In fact, dispersal and vicariance scenarios are not mutually exclusive (Sanmartín, 2009).
A second wave of intraspecific vicariance occurred ~200-100 ka in the Alpine-Carpathian system (R. cassubicifolius s.l.), between the Pyrenees and the Massif Central (R. envalirensis s.l.) and, to a lesser extent, within the Illyrian clade (R. notabilis s.l.). Most taxa occur today disjunctly in very restricted areas with few populations, although appropriate habitats would be broadly available. The disjunct patterns suggest that survival of the last glacial periods occurred also in microrefugia where local conditions protected plant species from unfavourable regional conditions (Rull, 2009). A distribution pattern suggesting survival in proximal microrefugia is specifically apparent in R. cassubicifolius in the northern foothills and foreland of the Alps. Here, the species occurs today only in patchy areas that remained ice-free between the large glacier tongues of the Last Gacial Maximum (LGM; see map in Paun, Stuessy, & Hörandl, 2006).
Strikingly, the species are still distributed just in their very restricted areas, despite post-glacial reforestation of temperate Europe. The failure of the sexual species in post-glacial geodispersal may be explained by biotic barriers, namely congeneric competitors.
The origin of polyploid apomictic derivatives is presumably younger than 100 thousand years (Pellino et al., 2013), and they rapidly colonized various habitats following the LGM (Paun et al., 2006). The rapid range expansion of younger, polyploid apomictic derivatives of the R. auricomus complex could have blocked range expansion of sexual progenitors (Hörandl, 2006). Polyploid apomicts could be more competitive, but they could also benefit from having no minority cytotype disadvantage (Levin, 1975): sexual outcrossers need a conspecific pollinator for reproduction, while the apomicts can reproduce as single individuals (Kirchheimer et al., 2018). Fertilization of sexual individuals by sympatric apomictic pollen donors can even result in the formation of aneuploid, introgressed offspring. Such processes were observed in diploid R. notabilis populations mixed with tetraploid apomicts (Hörandl, Greilhuber, & Dobes, 2000) and may also occur between sexual and apomictic cytotypes of R. marsicus (Dunkel, 2011

ACK N OWLED G EM ENTS
We thank the German Research Foundation for project funding (DFG, Ho4395/10-1) to E.H., within the priority programme "Taxon-Omics: New Approaches for Discovering and Naming Biodiversity" (SPP 1991), and the Botanical Museum, University of Oslo (O) for the material of R. pygmaeus. We also thank Mr F. G. Dunkel for providing herbarium vouchers and living material for the recently discovered diploid populations from the Illyrian region and France.

AUTH O R CO NTR I B UTI O N S
S.T. and E.H. conceived the ideas; S.T., K.K., L.H. and E.H. collected materials and data; S.T., K.K. and C.P. analysed the data; and S.T.
wrote the paper with contributions from all co-authors.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw reads used in this study are deposited in the Sequence Read Archive (SRA) under the BioProject PRJNA628081.