Hybrid zones in the European Alps impact the phylogeography of alpine vicariant willow species (Salix L.)

Introduction: In the European Alps, Pleistocene climate oscillations resulted in geographical range expansions and restrictions of species. Postglacial recolonizations often result in secondary contact hybridization of vicariant species, thereby creating hybrid zones with patterns of introgression. Here, we compare the genetic structure of two secondary contact hybrid zones between two vicariant willow species pairs occurring in the European Alpine System. Supplemented by morphological and ecological data, we try to understand the factors shaping the hybrid zones and their influence on geographical range filling patterns.

Methods: RAD sequencing and morphometric data were used to characterize biogeographical history, genetic diversity and the hybrid zone of each species pair. Vegetation relevés and species distribution models provided ecological context and support.

Key results: Results suggest that recolonization of the Alps happened from peripheral glacial refugia, resulting in broad secondary contact zones in the Eastern Alps in both species pairs. Both hybrid zones show introgression, but differ in symmetry and intensity of gene flow, in the type of introgressed loci, and in the geographical range. Habitat preferences and species distribution models do not indicate ecological barriers to recolonization.

Conclusions: Hybrid zones do not only affect the genetic structure of species by gene flow and introgression, but also appear to impact the biogeographical patterns of species.

1 Introduction

The European Alps are the highest mountain range within the European Alpine System and harbor a great species diversity (Ozenda, 1988; Kadereit, 2017). The Alps form an arc spanning over 1200 km in the east-west direction, and about 250 km in the north-south direction (Aeschimann et al., 2004). The western part comprises the highest peaks with mountain tops above 4000 m elevation above sea level (a.s.l.) and is geologically diverse. In the eastern part the mountain chains remain below 4000 m a.s.l. and are geologically more clearly structured. The Eastern Alps are divided into the Northern and Southern Alps with calcareous bedrock, and the Central Alps with predominant siliceous bedrock, but also local calcareous substrates (Supplementary Figure S1). The east-west differentiation and the geological substrate influenced the different biogeographical histories and ecological niches of alpine plants (Schönswetter et al., 2005; Alvarez et al., 2009). Also, the Pleistocene climate fluctuations impacted the evolution and geographic distribution of genetic lineages (Hewitt, 2004). In Central Europe, cold-adapted plant species have survived the last glaciations in different refugial areas. These include ice-free regions within the glaciated area (Nunatak refugia hypothesis (Stehlik, 2000; Schönswetter and Schneeweiss, 2019)), unglaciated areas at the periphery of the ice sheet (peripheral refugia hypothesis (Schönswetter et al., 2005; Holderegger and Thiel-Egenter, 2009; Carnicero et al., 2022)) or lowlands (lowland refugia hypothesis (Schönswetter et al., 2005; Holderegger and Thiel-Egenter, 2009; Guerrina et al., 2022)). Additionally, other mountains of the European Alpine System that were glaciated to a lower extent or remained unglaciated played an important role as refugial areas. Further glacial refugia for cold-adapted alpine plants are usually spatially restricted and occur in the Dinarids (Spaniel and Resetnik, 2022), the Apennines, the Carpathians and the mountains of the Iberian Peninsula (Gentili et al., 2015; Schmitt, 2017). Beside these large-scale factors, the climatic gradients with increasing elevation, and microhabitat conditions in the mosaic-like subalpine and alpine vegetation influence occurrences of species.

After the last glaciation period, species have started recolonizing newly available habitats and expanding their ranges from their refugia. The current range of cold-adapted plant species often represents an incomplete transient stage (Dullinger et al., 2012b). Every species has its own recolonization rate. Location of refugia, ecological niche, dispersal ability, geographical pattern of suitable habitats, and other evolutionary processes like hybridization or polyploidization modify the recolonization rate and influence the range filling (Hewitt, 1996; Tribsch and Schönswetter, 2003; Kirchheimer et al., 2016; Smycka et al., 2022). The impact of certain factors, such as hybridization, are complex and have been largely overlooked. Hybrid zones can arise between two closely related but previously isolated species that merge their ranges during the process of recolonization in a secondary contact zone. These zones provide opportunities to assess contemporary dynamics of selection and the evolution of reproductive isolation (Barton, 2001; Petit and Excoffier, 2009). Homoploid hybrids can create adaptive combinations in the parent species by introgression and potentially promote their range expansion (Abbott et al., 2013). Adaptive introgression is likely in loci that are associated to advantageous traits, whereas loci involved in reproductive isolation are less prone to introgression (Suarez-Gonzalez et al., 2018). Testing for adaptive introgression requires not only identification of introgressed loci, but also of signatures of selection, of phenotypic changes and fitness parameters (Suarez-Gonzalez et al., 2018). Hybrid zones also facilitate gene flow between different species and introduce new genetic variation that increases the genetic diversity and may be beneficial (Rieseberg and Carney, 1998). Genomic clines will also detect intrinsic genetic incompatibilities between hybridizing species due to reproductive isolation (Gompert and Buerkle, 2012), and can therefore help us to understand the role of selection in maintaining species barriers (Fitzpatrick, 2013). However, hybrid zones can also lead to competitive interactions between the hybrids and the parent species. Hybrids may outcompete the parent species and force them to retreat to their optimal niches, resulting in ecogeographical displacement of hybrids and parents (Abbott and Brennan, 2014; Kadereit, 2015).

Here we compare the spatial genetic structure of two hybrid zones in the European Alps by investigating two sister species pairs of willows (Salix L.). This genus comprises approximately 33 species in the European Alps, representing a broad range of geographic patterns, habitat preferences and ploidy levels (Hörandl, 1992; Aeschimann et al., 2004; Wagner et al., 2021). Hybridization is common in Salix and has even been reported between species of different taxonomic sections (Argus, 1974; Hörandl et al., 2012; Gramlich et al., 2016, 2018). The two well-defined sister species pairs S. alpina and S. breviserrata, as well as S. foetida and S. waldsteiniana, are endemic to Europe (Neumann and Polatschek, 1972; Büchler, 1986; Wagner et al., 2018, 2021). They show vicariant distributions following an east-west pattern in the European Alps and adjacent mountain systems, and share a contact area situated in the Eastern Alps (Hörandl, 1992). All four species investigated are diploid (2n = 38). They are dioecious, wind- and insect pollinated and produce numerous wind-dispersed seeds (Wagner et al., 2021).

Salix alpina occurs in the Eastern Alps, Carpathians, and the Balkan Peninsula. S. breviserrata is found from the Eastern to the Southwestern Alps and in the Cantabrian Mts. Records in the Apennines (Martini and Paiero, 1988) were examined in the field (E.H.) but could not be confirmed. Salix alpina grows on pure limestone and dolomite, as well as base-rich silicate while S. breviserrata colonizes basic silicate substrates but can also occur on limestone (Merxmüller, 1952; Hörandl, 1992). Both species are ascending to erect dwarf shrubs without runners or creeping rhizomes, i.e. without clonal growth; both growing on alpine screes, rocks and in snowbeds. Morphologically, the species differ mostly by leaf margin, which is entire and ciliate for S. alpina but densely denticulate or serrate for S. breviserrata; the flowers and fruits of the two species are not differentiated (Rechinger, 1964; Martini and Paiero, 1988; Hörandl, 1992; Hörandl et al., 2012). Intermediate forms with morphological characteristics of both species have been reported in the overlapping area of the two species ranges and were interpreted to represent natural homoploid hybrids (Hörandl, 1992).

Salix foetida is present from the Western Alps to Tyrol as well as in the Central Apennines. Doubtful occurrence records in the Pyrenees (Blanco, 1993) were examined in the field (L.P.) but turned out to be misidentifications. Salix waldsteiniana occurs in the Eastern Alps, Dinarids and Balkan Peninsula; some records from the Balkan Peninsula (Rila and Vitosha Mts.) turned out to belong to S. phylicifolia L. s.l. (Wagner et al., 2023). Salix foetida primarily inhabits moist environments such as swamps, snowbeds, and rivulets, thriving on acidic silicate substrates. In contrast, S. waldsteiniana predominantly grows on carbonate soils but can also occur on base-rich silicates. It is less sensitive to dry conditions and is commonly found on screes, in shrubberies, and along rivulets. Both species are medium-sized, erect shrubs without clonal growth, and occur in the subalpine zone. Morphologically, S. foetida is characterized by dense and prominent teeth on the leaf margin that are terminated with bright glands. Salix waldsteiniana has sparser and less pronounced teeth with darker glands. Flowers and fruits are not differentiated, but S. waldsteiniana tends to form larger catkins (Rechinger, 1964; Hörandl, 1992; Hörandl et al., 2012). In the contact area of the two species ranges, transitional forms with intermediate morphological characteristics have been confirmed to be hybrids (Marinček et al., 2023).

Our study combines restriction site-associated DNA sequencing (RAD-seq) data from samples out of the whole distribution ranges with ecological niche analyses, morphological analyses and species distribution models to reconstruct the phylogeography of the four species and characterize the contact zones. The hybrid zones are analyzed to identify their genetic structure, to quantify genomic patterns of introgression and to compare the species pairs with respect to degree and symmetry of introgression. With the combined information of genetic, morphological and ecological data, we try to understand the potential factors shaping the hybrid zones and their influence on geographical range filling patterns of species.

2 Materials and methods

2.1 Sampling and molecular data generation

Leaf material of 162 individuals of S. alpina and S. breviserrata was collected from 25 locations, and leaf material of 161 individuals of S. foetida and S. waldsteiniana was collected from 33 locations, representing the whole distribution ranges of all four species, including intermediate forms. The collection and determination of pure species and intermediate forms were based on Hörandl (1992), and the authors expertise. The dataset for S. foetida and S. waldsteiniana includes 102 individuals already included in a previous study (Marinček et al., 2023). Details about locations and phenotypes are given in Supplementary Table S1. Leaves were dried in silica gel and voucher specimens were deposited in the Herbarium of the University of Göttingen (GOET). The DNA was extracted using the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following a modified manufacturer’s protocol (Marinček et al., 2023). After quality check, the DNA was sent to Floragenex Inc. (Beaverton, USA) for library preparation following the protocol described in (Baird et al., 2008) and single-end RAD-sequencing. The restriction enzyme PstI was used for digestion and fragments between 300 bp and 500 bp were selected.

Each species pair was processed separately. First, a dataset including all individuals for each species pair and two accessions of Salix reticulata as outgroup was generated using Ipyrad v0.9.52 (Eaton and Overcast, 2020). Raw reads were demultiplexed and filtered for low-quality using default parameters and trimmed to 86 bp. Reads were clustered following the optimization as in Wagner et al. (2018), the minimum number of samples per locus was set to 50 (30%). This dataset was used to construct a phylogenetic tree that will allow us to sort out misidentifications or hybridization with other species among our samples, identify individuals with high genetic similarity and provide a phylogenetic framework for further analyses. Then, a second dataset for each species pair was generated using Stacks v2.65 (Catchen et al., 2013), which is more suitable for population genetic analysis. For this dataset, individuals were combined into groups of comparable size based on the previous phylogenetic analysis and on geographic proximity (Figures 1, 2; Supplementary Table S1). Six individuals were removed from the S. alpina and S. breviserrata dataset because they were found in a clade with less than 6 individuals, the minimum required sample size to recover within-population genetic diversity estimates (Nazareno et al., 2017). For S. foetida and S. waldsteiniana, fourteen individuals were removed (See red crosses in Figures 1, 2). A total of sixteen groups were defined for S. alpina and S. breviserrata, and sixteen groups for S. foetida and S. waldsteiniana (Figures 1, 2; Supplementary Table S1). The process_radtags program with default parameters was used to demultiplex raw reads, trim the reads to 86 bp and remove low-quality reads. The quality of the resulting reads was checked using FastQC v0.11.4 (Andrews, 2010). After parameter optimization using a subset of 30 individuals and the r80 method (Viricel et al., 2014; Paris et al., 2017; Rochette and Catchen, 2017), the denovo_map.pl program was used with the following parameters: -M 4, -n 4, -m 3 (-M maximum bp difference between two stacks within a sample; -m minimum coverage for a stack; -n maximum bp difference between stacks to be considered as orthologous across samples). At this point, individuals for which a large proportion of the reads represented PCR duplicates were removed. The populations program was further used to filter SNPs. Potential sequencing errors were removed by using a minimum minor allele count of 10 to process a SNP. Only SNPs with sequence data for at least 80% of the individuals were retained. Based on allelic balance at heterozygous sites, we checked for cross-contamination. We expected erroneous alleles coming from contamination to be much rarer than the correct allele. Individuals showing strongly unbalanced number of reads for alleles at biallelic sites were removed. Additionally, VCFtools was used to remove individuals with more than 15% missing genotypes, sites with more than 10% missing data and SNPs with excessive coverage (higher than two standard deviation above the mean) because they represented putative paralogous loci (Danecek et al., 2011). To remove linkage disequilibrium, only one SNP per locus was retained. An additional threshold for a minor allele frequency of 0.04 was applied for the dataset used with Structure, BGC, Treemix and adegenet (PCA). Further filters were applied before using Treemix and BGC to meet the analyses requirements. A summary of the pipelines, filtering steps and software used with the corresponding dataset is presented in Supplementary Figure S2.

Figure 1

Figure 1. RAxML phylogeny of Salix alpina, Salix breviserrata and the intermediate forms based on 65,497 loci (576,045 concatenated SNPs). Bootstrap support is color-coded for each branch. Crosses in front of the sample name indicate that the sample was not included in the further analyses because it did not group with enough other samples (red) or because it did not pass the filtering steps (blue). Samples sharing a high proportion of identical alleles pairs (>90%) are displayed in grey. Code on the right indicates the group ID used for genetic analyses. Countries where the samples were collected and phenotypes are indicated in the last two rightmost columns, respectively (AT, Austria; CH, Switzerland; ES, Spain; FR, France; IT, Italy; MK, North Macedonia; Sl, Slovenia, SK, Slovakia).

Figure 2

Figure 2. RAxML phylogeny of Salix foetida, Salix waldsteiniana and the intermediate forms based on 65,225 loci (575,194 concatenated SNPs). Bootstrap support is color-coded for each branch. Crosses in front of the sample name indicate that the sample was not included in the further analyses because it did not group with enough other samples (red), because it did not pass the filtering steps (blue), or due to potential cross-contamination (yellow). Samples sharing a high proportion of identical alleles pairs (>90%) are displayed in grey. Code on the right indicates the group ID used for genetic analyses. Countries where the samples were collected and phenotypes are indicated in the last two rightmost columns, respectively (AT, Austria; CH, Switzerland; ES, Spain; FR, France; IT, Italy; MK, North Macedonia; Sl, Slovenia, SK, Slovakia).

2.2 Phylogenetic and population genetic analyses

The relationship between all individuals of a species pair was investigated using a maximum-likelihood (ML) building approach with RAxML v8.2.4 (Stamatakis, 2014). For S. alpina and S. breviserrata, the analysis was made with 64,683 RAD loci (598,650 SNPs) in a concatenated alignment of 5,648,559 bp with 17.04% missing data. For S. foetida and S. waldsteiniana 65,225 RAD loci (575,194 SNPs) in a concatenated alignment of 5,704,214 bp with 19.45% missing data were used. The GTR + Γ model of nucleotide substitution was used and a bootstrapping analysis of 100 replicates was performed. The resulting trees were visualized with Figtree v1.4.4 (Rambaut and Drummond, 2012).

Based on 35,858 unlinked and quasi-neutral SNPs in 152 individuals for S. alpina and S. breviserrata, and 36,533 unlinked and quasi-neutral SNPs in 140 individuals for S. foetida and S. waldsteiniana, a principal component analysis (PCA) was performed using the R package adegenet (Jombart, 2008). Additionally, the genetic structure among individuals was also analyzed with the Bayesian clustering approach implemented in Structure v.2.3.4 (Pritchard et al., 2000). For the analysis, the admixture model was used with correlated allele frequencies and no other prior information (Falush et al., 2003). Eight repetitions for each number of genetic clusters (K) ranging from 1 to 10 were tested. Each run had 50,000 iterations as burn-in, followed by 150,000 MCMC iterations. The R package Pophelper 2.3.1 was used to estimate the optimal number of clusters K with the Evanno method, and plot the results (Evanno et al., 2005; Francis, 2017). To infer gene flow between the different groups of a species pair, the composite-likelihood approach implemented in Treemix v.1.13 was used (Pickrell and Pritchard, 2012). Maximum Likelihood trees allowing for a predefined number of migration events were computed. The number of migration events ranged from 0 to 30 and each was tested three times. The optimal number of migration events was determined using the Evanno method. Finally, node support was estimated for the ML tree by running 100 bootstrap replicates. The R package Bite v.2 was used to find the optimal number of migration events and for bootstrap analysis (Milanesi et al., 2017). Stacks (Catchen et al., 2013) was used to compute pairwise Wright’s F_st between groups, and other statistics (gene diversity, nucleotide diversity, observed heterozygosity and the number of private alleles). Finally, differential patterns of introgression in hybrid groups were investigated with the Bayesian genomic cline model implemented in the software BGC (Gompert and Buerkle, 2012). It quantifies locus-specific ancestry given a hybrid index and identifies outliers based on two parameters, α and β. The hybrid index represents the proportion of an individual’s genome that is inherited from a reference population. The genomic cline center (α) gives the direction of introgression. Outliers may be involved in directional selection towards homozygotes or selection against heterozygous genotypes (Gompert and Buerkle, 2012). The cline rate (β) gives SNP frequency changes. Negative values indicate a wider cline associated with a higher introgression rate than expected. These loci may be involved in adaptive introgression. Positive values are associated with a steeper cline and reproductive isolation (Gompert and Buerkle, 2012). For the analysis, reference groups were defined based on the Q-values of Structure for K=2. Individuals with Q < 0.1 or Q > 0.9 were assigned to the respective reference groups, and the remaining were considered hybrids. The hybrid index for each admixed individual was calculated using the software BGC (Gompert and Buerkle, 2012). The genomic cline analysis was run three time independently for 600,000 MCMC iterations after a 300,000 burn-in phase. A thinning interval of 100 was used. The runs were combined, and convergence was assessed visually. Outlier loci were detected when they met both criteria, that the 95% credible intervals for α and β excluded zero and that the median of posterior distribution was outside the interval between 1-0.975/2 and 0.975/2 quantiles of the probability distribution. The R package ClineHelpR was used to visualize the output (Martin et al., 2021).

2.3 Morphometric analysis

Morphometric measurements were made for 154 individuals from fifteen groups of S. alpina and S. breviserrata. For S. foetida and S. waldsteiniana, 70 individuals were measured and added to the 95 individuals already measured in Marinček et al. (2023), covering fifteen groups (Supplementary Table S2). Up to ten entire and healthy leaves per individual were scanned with a Canon CanoScan LiDE 220 (Canon, Tokyo, Japan). Leaves were randomly selected from the herbarium specimen. Each leaf was described by six measurements and four ratios (Supplementary Table S2) using the software ImageJ (Schneider et al., 2012). The characteristics were chosen based on previous morphological studies on Salix L (Kosiński et al., 2017; Marinček et al., 2023). Leaf metrics were averaged per individual. Correlations between characteristics were checked using Spearman’s correlation coefficients, and four uncorrelated characteristics (correlation ≤ 0.7) were further standardized and analyzed by Principal Component Analysis (PCA). The area of the leaves (A), the number of teeth (T), the width-to-length ratio (BLBW) and the ratio of the total length to the length of the widest part of the leaves (BLBWP) were retained. Individuals were assigned to a category based on Structure Q-values. Individuals with less than 10% admixture (Q < 0.1 and Q > 0.9) were assigned to the respective parent species, while the other were treated as hybrids. Tests and analyses were performed with the R package MorphoTools2 (Šlenker et al., 2022).

2.4 Ecological niche differentiation

To examine the ecological niches of the four species outside the contact zone, we gathered vegetation relevés from the EVA Database and from our collection sites in the contact zone (Supplementary Table S3). Unprecise relevés, duplicates and other relevés occurring inside the contact zone were removed to avoid including wrongly determined species. A total of 1214 relevés were kept (470 for S. alpina, 118 for S. breviserrata, 160 for S. foetida, and 466 for S. waldsteiniana). Four uncorrelated Landolt indicator values (Pearson r ≤ 0.7) were derived from accompanying species of the target species. Landolt indicator values classify species based on their realized ecological niche using expert knowledge (Landolt et al., 2010). We used temperature (T), soil moisture variability during the growing season (VSMGS), soil pH (SpH) and soil nutrient availability (SN). The indicators unweighted mean values for each relevé were analyzed by PCA from the R package stats. Vegetation relevés from our sampling (Supplementary Table S3) were also plotted to compare the ecology of the different groups. Additionally, Schoener’s D index of niche overlap was calculated for each species pair using the R package Ecospat (Di Cola et al., 2017). The index ranges from 0 (no overlap) to 1 (complete overlap). For this test, background values were obtained from relevés of other mountain Salix species.

2.5 Species distribution modeling

Species occurrences were obtained from the European Vegetation Archive EVA database (Chytrý et al., 2016), Infoflora (https://www.infoflora.ch/), GBIF database (https://www.gbif.org/) and our collections. Only occurrences with precise coordinates were kept, and occurrences falling within 100m from another occurrence were discarded using the R package rangeBuilder to limit the sampling effort differences (Davis Rabosky et al., 2016). A total of 3,370 occurrences were kept (488 for S. alpina, 489 for S. breviserrata, 2,013 for S. foetida, and 380 for S. waldsteiniana). Current and historic climatic conditions were derived respectively from the CHELSA database (Karger et al., 2017; Karger et al., 2018) and PaleoView (Fordham et al., 2017). Arithmetic means of the climate variables were used to calculate bioclimatic variables. Six uncorrelated (Pearson < 0.7) bioclimatic variables were used (BIO3: Isothermality; BIO5: Mean daily maximum air temperature of the warmest month; BIO6: Mean daily minimum air temperature of the coldest month; BIO7: Annual range of air temperature; BIO14: Precipitation amount of the driest month; BIO18: Mean monthly precipitation amount of the warmest quarter). To represent the topographic variation of climate in the mountain regions, the bioclimatic variables were downscaled to a 100 m x 100 m resolution as described in (Dullinger et al., 2012a) and using a kriging with elevation as covariable. SDMs (species distribution modelling) were calibrated using the R package Biomod2 (Thuiller et al., 2023). Presence data were supplemented with four randomly selected pseudo-absence datasets. Four techniques (Generalized Linear Model, Generalized Boosting Model, Random Forest, and Generalized Additive Model) were used to model the distribution of our species. Every model was repeated four times with every pseudo-absence dataset. The calibration of the models was done with 80% of the data and, the remaining 20% were used to evaluate the models. Models with a True Skill Statistic TSS score above 0.9 were put together to generate ensemble projections of potential species distribution based on the mean of single models. Projections were made in the recent and in the Last Glaciation Maximum (LGM; 21 Kya) climatic conditions. The results of the suitable habitats were visualized in QGIS (https://www.qgis.org/).

3 Results

3.1 Phylogenetic analyses

The ML phylogenies showed well-supported topologies (Figures 1, 2). The trees revealed distinct clades for all locations. For S. alpina and S. breviserrata, all locations were arranged in a nested pattern that followed a colonization from northeast to southwest. The Tatra group of S. alpina appeared in sister position to all the other locations (BS 100). All intermediate forms and individuals of S. breviserrata were derived from S. alpina. The clade of S. breviserrata was well supported (BS 100). For S. foetida and S. waldsteiniana, the individuals were divided into two well-supported clades (BS 100). Individuals from North Macedonia identified as S. waldsteiniana in the field were found within the clade of S. foetida. The clade comprising S. waldsteiniana individuals also comprised all phenotypically intermediate forms.

3.2 Population genetic analyses

For S. alpina and S. breviserrata, Structure identified K = 2 to be the optimal number of genetic clusters (Supplementary Figure S3). One cluster (blue) represented S. breviserrata, and the other cluster (red) S. alpina (Figure 3A). Groups of S. alpina located in the Tatra Mts. and North Macedonia showed pure ancestry, while groups in the Alps showed a moderate level of admixture with S. breviserrata. Intermediate forms were represented by admixture following a geographic cline. The genetic structure was confirmed with the PCA (Supplementary Figure S5). Treemix identified nine and 21 migration events as most likely. We chose m = 9 for clarity, under which 99.6% of the variation in the data was explained (Supplementary Figure S7). Overall, the topology suggested that S. breviserrata derived from intermediate forms, that derived themselves from S. alpina (Figure 3B). Most signals of gene flow were detected from Tatra Mts. (A_SK15) to alpine groups of S. alpina and intermediate forms. Three additional migration events were found between groups of the secondary contact zone. Group statistics identified most private alleles for S. alpina in the Tatra Mts. (A_SK15) and North Macedonia (A_MK16) (Table 1). The observed heterozygosity, gene diversity and nucleotide diversity were very similar for all groups, ranging from 0.185 to 0.227, from 0.177 to 0.225, and from 0.191 to 0.237, respectively. The inbreeding coefficient was slightly positive for all groups. Pairwise F_st values ranged between 0.046 and 0.147 (Supplementary Table S4). The individual hybrid index (Figure 4a) ranged from 0.17 to 0.82 and was strongly correlated with the longitude (Pearson’s r = 0.76, p < 0.001). Among the 526 loci used in BGC, 27 were negative α outliers and 51 negative β outliers (Figure 4b).

Figure 3

Figure 3. Summary of the phylogeographic analyses for 16 groups (23 locations) representing 152 individuals of Salix alpina and Salix breviserrata, based on 35,858 unlinked loci. (A) Map showing the results of Structure depicting genetic structure across the 23 locations for K=2. Only the average ancestry proportions of individuals for each location are displayed, providing an overview of the spatial distribution of genetic structure. Colors around the pie charts represent the group. Small dots in the background represent observations from GBIF (https://www.gbif.org, blue for Salix breviserrata, red for Salix alpina). The putative contact zone, based on phenotype only, is indicated by the dashed outline. (B) Gene flow between the 16 groups of Salix alpina and Salix breviserrata inferred by Treemix at the most likely number of migration m = 9. Arrows indicate the direction of the gene flow, and the color indicates the percentage of alleles that migrated from the source group (1% = yellow, 25% = orange, 50% = red). Bootstrap values for the nodes are color coded (dark green dot = 90-100%, green square = 75-90%, light green triangle = 50-75%). (C) Last Glaciation Maximum habitat suitability for Salix breviserrata and (D) for Salix alpina. Highly suitable areas are indicated in dark red. Areas covered by ice are indicated in grey. Arrows indicate potential recolonization routes. Circles indicate potential peripheral refugia. The European Commission, Eurostat (ESTAT), GISCO (2020). Countries, 2020 - Administrative Units - Dataset was used to draw the administrative boundaries ^© EuroGeographics.

Table 1

Table 1. Group statistics for Salix alpina and Salix breviserrata based on 115,347 loci.

Figure 4

Figure 4. (a) Histogram representing the distribution of the hybrid index of the admixed individuals for Salix alpina and Salix breviserrata. Hybrid index of 0 means that the individuals is pure Salix breviserrata, hybrid index of 1 means that the individual is pure Salix alpina. (b) Genomic cline generated by BGC depicting the probability of P1 (Salix alpina) alleles given the hybrid index for 526 SNPs. Outliers are color-coded. The black line gives the null expectation based on genome-wide admixture.

For S. foetida and S. waldsteiniana, K = 2 was identified to be optimal by Structure (Supplementary Figure S4). One cluster (green) was assigned to S. foetida, for which individuals outside the contact zone showed no or low level of admixture. The other cluster (orange) represented S. waldsteiniana and was pure only in Slovenia (Figure 5A). Individuals identified as S. waldsteiniana that were sampled in Italy, Austria, and North Macedonia, as well as the intermediate forms showed some level admixture. The main genetic clustering was confirmed by the PCA (Supplementary Figure S6). In the Treemix analysis, three migration events were detected to be optimal, and the resulting tree explained 99.5% of the total variation (Supplementary Figure S8). The overall topology suggested an early divergence between S. foetida and S. waldsteiniana (Figure 5B). Migration edges were found from an ancestral population of S. waldsteiniana to North Macedonia (W_MK16), from Slovenian S. waldsteiniana (W_SI13) to S. waldsteiniana in the secondary contact zone (W_IT10), and from the ancestry group of S. foetida to intermediate forms (H_IT8). Stacks found 650 private alleles for S. foetida in the secondary contact zone (F_IT7); other groups ranged from 78 to 0. The observed heterozygosity, gene diversity and nucleotide diversity were moderate, ranging from 0.162 to 0.222, from 0.113 to 0.215 and from 0.148 to 0.227, respectively (Table 2). Finally, the inbreeding coefficient was slightly negative for two groups and slightly positive for the rest. Pairwise F_st values ranged between 0.037 and 0.231 (Supplementary Table S5). Hybrid indices ranged from 0.45 to 0.98 (Figure 6A). Only two groups (F_IT7 and W_IT10) were identified as hybrids, the four other groups (H_IT8, H_IT9, W_IT11 and W_AT12) had hybrid indices higher than 0.97 and were identified as S. waldsteiniana. Among the 1167 loci used for the analysis, 1147 (98%) were identified as outliers, representing 173 positive α, 318 negative α, 443 positive β and 715 negative β outliers (Figure 6B). Noticeably, 171 positive β outliers were also negative α outliers, but only 10 were also positive α outliers.

Figure 5

Figure 5. Summary of the phylogeographic analyses for 16 groups (24 locations) representing 140 individuals of Salix foetida and Salix waldsteiniana, based on 36,533 unlinked loci. (A) Map showing the results of Structure depicting genetic structure across the 24 locations for K=2. Only the average ancestry proportions of individuals for each location are displayed, providing an overview of the spatial distribution of genetic structure. Colors around the pie charts represent the group. Small dots in the background represent observations from GBIF (https://www.gbif.org, green for Salix foetida, orange for Salix waldsteiniana). The putative contact zone, based on phenotype only, is indicated by the dashed outline. (B) Gene flow between the 16 groups of Salix foetida and Salix waldsteiniana inferred by Treemix at the most likely number of migration m = 3. Arrows indicate the direction of the gene flow, and the color indicates the percentage of alleles that migrated from the source group. (1% = yellow, 25% = orange, 50% = red). Bootstrap values for the nodes are color coded (dark green dot = 90-100%, green square = 75-90%, light green triangle = 50-75%). (C) Last Glaciation Maximum habitat suitability for Salix foetida and (D) for Salix waldsteiniana. Highly suitable areas are indicated in dark red. Areas covered by ice are indicated in grey. Arrows indicate potential recolonization routes. Circles indicate potential peripheral refugia. The European Commission, Eurostat (ESTAT), GISCO (2020). Countries, 2020 - Administrative Units - Dataset was used to draw the administrative boundaries ^© EuroGeographics.

Table 2

Table 2. Group statistics for Salix foetida and Salix waldsteiniana based on 128,200 loci.

Figure 6

Figure 6. (A) Histogram representing the distribution of the hybrid index of the admixed individuals for Salix foetida and Salix waldsteiniana. Hybrid index of 0 means that the individuals is pure Salix foetida, hybrid index of 1 means that the individual is pure Salix waldsteiniana. (B) Genomic cline generated by BGC depicting the probability of P1 (Salix waldsteiniana) alleles given the hybrid index for 1167 SNPs. Outliers are color-coded. The black line gives the null expectation based on genome-wide admixture.

3.3 Morphometric analysis

Salix alpina and S. breviserrata were mostly differentiated from each other (Supplementary Figure S9), with highest differentiation by number of teeth (T). The PCA showed that hybrids had intermediate forms and fell within the phenotypic range of both parent species.

For S. foetida and S. waldsteiniana, the PCA showed substantial overlap of the parent species (Supplementary Figure S10). The number of teeth (T) was again the character that contributed most to the species differentiation. Hybrids showed intermediate phenotypes.

3.4 Ecological niche differentiation

The PCA for S. alpina and S. breviserrata indicated overlapping ecology between the two species (Figure 7B). Salix breviserrata had the most variable environmental niche and was found in a broader range of soil moisture (VSMGS), whereas S. alpina appeared on a broader range of soil nutrients. The niche overlap was 63% between the two species.

Figure 7

Figure 7. (A) Current habitat suitability based on six uncorrelated high-resolution climatic layers for Salix breviserrata (up) and Salix alpina (down). High suitability is indicated in dark red. The current distributions are indicated by the thick dashed outline. Areas covered by ice during the Last Glacial Maximum are indicated in grey. The Alps are delimited by the thick continuous outline. (B) Principal component analysis of environmental variables for Salix breviserrata and Salix alpina. Ellipses represent 95% confidence intervals for Salix breviserrata (blue) and Salix alpina (red) outside the contact zone. Colored dots represent the groups used for genetic analyses. Eigenvectors represent the 4 uncorrelated Landolt indicator values: Temperature (T), variability of soil moisture during the growing season (VSMGS), soil pH (SpH) and soil nutrient availability (SN). The European Commission, Eurostat (ESTAT), GISCO (2020). Countries, 2020 - Administrative Units - Dataset was used to draw the administrative boundaries ^© EuroGeographics.

The PCA for S. foetida and S. waldsteiniana showed overlapping ecological niches (Figure 8B). Soil moisture during the growing season (VSMGS) contributed most to a differentiation between the two species. The niche overlap was 41%.

Figure 8

Figure 8. (A) Current habitat suitability based on six uncorrelated high-resolution climatic layers for Salix foetida (up) and Salix waldsteiniana (down). High suitability is indicated in dark red. The current distributions are indicated by the thick dashed outline. Areas covered by ice during the Last Glacial Maximum are indicated in grey. The Alps are delimited by the thick continuous outline. (B) Principal component analysis of environmental variables for Salix foetida and Salix waldsteiniana. Ellipses represent 95% confidence intervals for Salix foetida (green) and Salix waldsteiniana (orange) outside the contact zone. Colored dots represent the groups used for genetic analyses. Eigenvectors represent the 4 uncorrelated Landolt indicator values: Temperature (T), variability of soil moisture during the growing season (VSMGS), soil pH (SpH) and soil nutrient availability (SN). The European Commission, Eurostat (ESTAT), GISCO (2020). Countries, 2020 - Administrative Units - Dataset was used to draw the administrative boundaries ^© EuroGeographics.

3.5 Species distribution modelling

All single models were kept for the ensemble models. The ensemble model for each of the four species indicated high predictive accuracy with a TSS > 0.9. For S. breviserrata, projections in the LGM conditions identified potential glacial refugia at the foot of the Southwestern Alps, in Massif Central and in the Eastern Pyrenees (Figure 3C). For S. alpina, the projections found vast suboptimal areas in the periphery of the Eastern Alps and in the Carpathians (Figure 3D). Salix foetida had highly suitable conditions all along the Mediterranean coast, from the Pyrenees to the Alps, in the Southern Apennines (Italy) and in the Balkan Peninsula (Figure 5C). For S. waldsteiniana, suitable areas were mostly detected in the eastern part of the Alps, and in the Southern Carpathians (Figure 5D). According to current climate data, the projections of the highly suitable conditions were consistent with the current distribution for all four species (Figures 7A, 8A). However, SDM predicted additional suitable habitats beyond the secondary contact zones.

4 Discussion

4.1 Biogeographical history of parent species

Alpine plant species in the European Alpine System survived the last glaciations in refugia, which provided a source for post-glacial recolonization (Hewitt, 1996). Notable peripheral refugial areas of the Alps include the Southwestern, Southern and Eastern Alps (Schönswetter et al., 2005). Additionally, extra-alpine origins have already been found for several alpine taxa (Moore and Kadereit, 2013; Surina et al., 2014; Kuzmanović et al., 2021). The present survey analyzed the genetic composition of two vicariant species pairs from their entire distribution areas to assess their biogeographical history. Analyses of gene flow and species distribution modeling suggest recolonization of the Alps from different potential refugial areas.

For S. alpina, our results suggest that the Tatra Mts. and Balkans harbored refugial populations, characterized by high numbers of private alleles (Figures 1, 3A; Table 1). In accordance with other studies that covered the Carpathians and the Alps, we found some genetic differentiation between the two ranges (Mráz et al., 2007; Ronikier et al., 2008). However, biogeographical connections between the Eastern Alps and the Tatra Mts. are a common pattern in alpine plants (Schmitt, 2017). In S. alpina, gene flow from the Tatra Mts. towards the Alps highlights the connection between both regions (Figure 3B). The high number of Alpine-Carpathian endemics depicts the close relationship between both regions and supports our finding (Pawłowski, 1970). Colonization from the isolated North Macedonian group would have required long-distance dispersal (LDD) or the extinction of bridging populations in the northern Balkan Peninsula after recolonization, as suggested for Wulfenia carinthiaca (Surina et al., 2014). Indeed, to our knowledge, no other records exist for S. alpina in this region (Jalas and Suominen, 1988). A lack of phylogeographical differentiation within the Alps, as observed for S. alpina, was also found in Oxytropis campestris s.l. and Hypochaeris uniflora (Schönswetter et al., 2004; Patrik et al., 2007). Only one group of S. alpina in the Eastern Alps (A_AT11) has higher number of private alleles, lower gene diversity and lower nucleotide diversity compared to the other alpine groups (Table 1). This group is situated at the edge of the last glaciation ice sheet on an island-like mountain of carbonate rocks (Grebenzen), within an otherwise siliceous refugial area (Schönswetter et al., 2005) and could have survived the LGM in a marginal refugium. Its isolated position from other S. alpina populations allowed to keep its genetic distinctness. Salix breviserrata represents a monophyletic group deeply nested within the clade of S. alpina (Figures 1, 3B). We suggest that the phylogenetic tree failed to separate the two species as sister clades due to broad introgression between the species pairs (see below). Hybridization has been shown to distort tree topologies (e.g (McDade, 1992)). A recent dated phylogeny suggests a split of S. alpina and S. breviserrata in the Pliocene (Marinček et al., 2024) and rejects the interpretation that S. breviserrata evolved out of S. alpina during the colonization of the Alps. Based on genetic results only, the potential glacial refugium for S. breviserrata remains unresolved. However, the projection of suitable conditions in the LGM identified a potential southwestern alpine peripheral refugium. The Maritime Alps are a known potential refugium (Taberlet et al., 1998), and evidences were found, e.g. for Berardia subacaulis and multiple Gentiana species (Christe et al., 2014; Guerrina et al., 2022). The samples from the Cantabrian Mts. appear on terminal branches of the phylogeny, suggesting a western range expansion from the Alps rather than a refugium, similar as in Ranunculus parnassiifolius (Cires and Prieto, 2012).

The refugial areas of S. foetida are unclear because the basal branch of the clade (Figure 2) is represented by a strongly differentiated group (F_IT7) in the Southern Alps, showing a high number of private alleles but low gene and nucleotide diversity. Excess of heterozygotes in this group support an interpretation as rare hybrid × hybrid genotypes, as observed in other willow hybrid zones (Gramlich et al., 2018). The low genetic diversity of this group (F_IT7) can be best explained by a drastic recent anthropogenic reduction of population size (covering just 15m²; see Supplementary Table S3) which has removed other, introgressed genotypes (see Marinček et al. (2023)). Another potential factor influencing genetic diversity within populations could be biased sex ratios in dioecious plants (Rosche et al., 2018), which would require further studies for our populations. For the other groups, the results could not identify any potential glacial refugium based on genetic statistics but the projection of suitable conditions in the LGM suggest a southwestern alpine peripheral refugium. The lack of distinct genetic structure and the partially low genetic diversity in both S. breviserrata and S. foetida confirm our hypothesis of rapid postglacial recolonization of the Alps. After the last glaciation maximum, both species rapidly colonized the newly available habitats, leading to a large area sharing similar genetic structure. Willows are well-known pioneer species that rapidly colonize disturbed habitats and glacier forefields (Karrenberg et al., 2003; Gramlich et al., 2016). For S. waldsteiniana, results suggest a single refugial area. As for other species from the Eastern Alps, the refugium may have been located in the periphery or in the northern Dinaric Mts (Schönswetter et al., 2005; Durovic et al., 2017). Unexpectedly, the North Macedonian samples previously reported as S. waldsteiniana (Micevski, 1995) came out in genetic and morphometric analysis as hybrids with S. foetida. However, there is no other record of S. foetida in the Balkan Peninsula, but several ones for S. waldsteiniana (Jalas and Suominen, 1988). The genetic pattern might reflect a shared evolutionary history with S. foetida in the Central Apennines from Messinian or Pleistocene land bridge connections (Spaniel and Resetnik, 2022); similar distributions are also observed e.g. in Saxifraga prenja (Gerschwitz-Eidt et al., 2023). Our SDM suggests better habitat conditions for S. foetida than for S. waldsteiniana in the North Macedonia mountains, Figure 8A). Alternatively, LDD and hybridization events between S. foetida and other S. waldsteiniana populations (Jalas and Suominen, 1988) in the Balkan Peninsula cannot be excluded. Despite different potential refugia, postglacial recolonization of the sister species in the Alps happened mostly from eastern and southwestern marginal areas, with secondary contact hybridization in the Central Eastern and Southern Alps.

4.2 Secondary contact and hybrid zone

To understand the impact of hybridization on the biodiversity of the Alps, Parisod (2022) encouraged studies across suture zones where differentiated lineages come into sympatry. With the phylogeographic framework of the four species of interest (see above), we can follow the suggestion and investigate the two contact zones between the species pairs. Our study supplements other studies of hybrid zones in Salix (Hardig et al., 2000; Fogelqvist et al., 2015; Gramlich et al., 2016, 2018; Marinček et al., 2023) but explores for the first time hybrid zones over the whole distribution ranges. We reanalyzed the hybrid zone between S. foetida and S. waldsteiniana in a broader context than done previously (Marinček et al., 2023), and we newly characterized a homoploid hybrid zone of the species pair S. alpina and S. breviserrata. Our results suggest that secondary contact zones are species-specific, although the species pairs share identical, unspecialized pollination systems, efficient dispersal capabilities via wind-dispersed seeds, and similar biogeographical histories. In both cases we found introgressive hybridization which can be related to the low genetic divergence between the parental species pairs, as observed in other willows (Gramlich et al., 2018).

Both species pairs were represented by nonadmixed (parental) individuals dominating one marginal areas of the distribution, with hybrid individuals occurring in between. The morphometric data also confirm these findings. Other alpine plant species, as shown for Androsace helvetica and Androsace pubescens, likely originated in the northern Eastern Alps and Southwestern Alps, expanded their ranges toward the center of the Alps and hybridized after secondary contact (Schneeweiss et al., 2017). Typically, hybrid zones are described as “narrow regions” of hybridization (Barton and Hewitt, 1985; Harrison, 1993; Jiggins and Mallet, 2000). In this regard, we observed two distinct patterns. Specifically, S. alpina and S. breviserrata exhibited a broad hybrid zone, extending over more than 300 kilometers between the Zillertal Alps (Italy) and the Graz Mountains (Austria). Similar broad hybrid zones were found in Fagus sylvatica subsp. orientalis, which grows in the Caucasus and Asia Minor and in Tephroseris helenitis in the Eastern Alps (Pflugbeil et al., 2021; Hrivnák et al., 2024). The hybrid zone was larger than we expected based on morphology, because patterns of introgression into both parental species were found in individuals with phenotypes similar to pure parental individuals. Morphometric data confirmed that hybrids have a broad phenotypic variability and cover most of the parental variability (Supplementary Figure S9). Morphological variability is in willows strongly influenced by local habitat conditions (Neumann, 1981), and incongruences between phenotype and genetic admixture were already reported for S. foetida and S. waldsteiniana (Marinček et al., 2023). Overall, hybridization between S. alpina and S. breviserrata suggests the formation of hybrid swarms, as observed for other alpine species, e.g. in Saxifraga (Ebersbach et al., 2020). For S. foetida and S. waldsteiniana, the hybrid zone was narrower and situated between the Western and the Eastern Alps, on the line between Lake Garda and Innsbruck, known as the Brenner zone (Kerner, 1870; Thiel-Egenter et al., 2011). Our genetic data confirmed the presence of S. waldsteiniana and hybrids spread over 50 kilometers along the Brenner line, although it may extend to east Tyrol (Hörandl, 1992; not sampled here). Hybridization between S. foetida and S. waldsteiniana creates an asymmetric hybrid zone, where individuals belonging to one parental species co-occur with their hybrids.

The hybrid zones were further analyzed to characterize the hybrids and identify loci with exceptional patterns. For S. alpina and S. breviserrata, all individuals defined as hybrids in the analysis showed at least 16% admixture. The strong correlation between hybrid index and longitude confirmed the broad pattern of introgression rather than hybrid speciation or reinforcement, as found before in willows (Gramlich et al., 2018). All α outlier loci represented excess ancestry of S. alpina. This could be the result of a selective advantage towards homozygous genotypes from S. alpina in the secondary contact zone (Guo et al., 2023). Species distribution modelling in the present shows higher suitability in the secondary contact zone for S. alpina than S. breviserrata. Also, the ecology of the secondary contact zone is more similar to the ecological niche of S. alpina. All β outlier loci represented faster introgression than expected. No locus was identified as a positive β outlier that would potentially indicate reproductive isolation. These results confirmed the hypothesis of Hörandl that crossing barriers between the two species are low (Hörandl, 1992). The dominant pattern of introgression between the two closely related species and the high number of migration events found between groups suggest that gene flow might be sufficient to weaken species boundaries. In addition, there are no major differences in flowering times (Hörandl, 1992). As suggested by our ecological analysis, the habitat requirements for each species share similarities, with more than 60% of niche overlap. The availability of suitable habitats with both calcareous and base-rich siliceous substrates in the Eastern Alps may have enhanced the formation of a very broad and symmetric hybrid zone. Field observations and herbarium specimens within the hybrid zone (Hörandl, 1992) indicate normal catkin formation and seed set, i.e. no obvious postzygotic crossing barriers. If no postzygotic isolation can limit the interspecific hybridization, further hybridization and introgression might lead to the genetic homogenization of both species (Rieseberg and Carney, 1998).

In S. foetida and S. waldsteiniana, we observed an asymmetric hybrid zone. Asymmetry occurs when differences in gene flow are observed between the species or uneven barrier strength promotes the introgression in one direction more than in the other (Pickup et al., 2019). Both species have identical pollination and dispersal systems, and it is therefore very unlikely to observe asymmetric gene flow. Asymmetry could be due to differential introgression of adaptive loci versus loci involved in reproductive isolation (Suarez-Gonzalez et al., 2018). Asymmetric introgression could also arise from different abundance of parental species in the initial phases of the hybrid zone formation, as observed in a recently formed hybrid zone between S. helvetica and S. purpurea in the Western Alps (Gramlich et al., 2018). Our results rather suggest the presence of a stronger crossing barrier in S. waldsteiniana. Notably, 170 outlier loci indicating a slower introgression than expected also represented excess ancestry of S. waldsteiniana (against 10 representing excess of S. foetida). These loci are potentially indicative of reproductive isolation. In contrast, more loci indicating faster introgression and adaptive advantages also represented excess ancestry of S. foetida. Selection likely favored adaptation of S. foetida on habitats with more moist soils. Such habitats were in the Eastern Alps mostly available directly after the melting of the ice-shield on glacier forefields, where S. foetida is usually one of the first colonizers. In this time period, S. foetida was perhaps in the Eastern Alps more abundant than today and left introgressive signatures into S. waldsteiniana. After reforestation, S. waldsteiniana could have been the more efficient colonizer on drier sites like limestone screes, and in shrubberies of Pinus mugo and Alnus alnobetula, that are nowadays dominating the subalpine zone of the Northern and Southern ridges of the Eastern Alps but also are available more locally in the central parts (Figure 8A). More detailed studies on the recolonization process would be needed to support this hypothesis. The combination of data potentially explains the asymmetry of the hybrid zone and the presence of almost unadmixed S. waldsteiniana in the secondary contact zone (group W_IT11 detected by Structure, BGC and morphometrics as pure S. waldsteiniana).

4.3 Hypothetical impact of hybrid zones on the species’ recolonization of the Alps

The current projections for the four species revealed available suitable conditions beyond the hybrid zones. Only 30% of alpine plant species fill their ranges in the Alps almost completely (Dullinger et al., 2012b). Interestingly, all four alpine willow species included in this study (S. herbacea, S. reticulata, S. retusa and S. serpillifolia) filled their ranges between 88% and 100%. The range filling is positively correlated with an increase in dispersal capacity (Dullinger et al., 2012b). It is therefore not surprising that willows, known to be pioneer species that colonize open habitats showed almost complete range filling (Hörandl et al., 2012). Willows are specifically efficient as first colonizers of glacier forefields (Gramlich et al., 2016), and their adaptations to wind-dispersal of seeds also allows for long distance dispersal over mountain ridges. It is consequently unlikely that the topography of the European Alps impacts strongly the current recolonization of the four focal species and hinders complete range filling. Topographic barriers to recolonization were rather observed in genera that are not specifically adapted to wind-dispersal (e.g., in Ranunculus kuepferi, Kirchheimer et al. (2018)). Moreover, our ecological data suggest similar niches between sister species. Differentiation in soil preferences was often mentioned to characterize willow species (Martini and Paiero, 1988; Hörandl, 1992), and are in general a strong factor for explaining distributions of plants in the Eastern Alps (Alvarez et al., 2009). However, based on our ecological analysis, we believe that bedrock conditions alone cannot be regarded as a strong barrier to range filling. Accordingly, the hybrid zones do not fit a scenario of strong ecogeographical displacement of hybrids and parents, as suggested by Kadereit, but rather show clinal patterns of introgression (Kadereit, 2015). Other ecological factors like pollinators and phenology were not explored but are probably less important as the species of interest are both wind- and insect-pollinated, with unspecific flower morphology and main flowering time during June and July. The response to biotic and abiotic stress factors with secondary compounds, as a potential factor, was found to be very similar for the species pair S. foetida and S. waldsteiniana (Volf et al., 2023).

Eventually, the prevalence of admixed ancestry in the secondary contact zone could be interpreted as a higher fitness of hybrid populations in comparison with parent species (Arnold and Martin, 2010). Genetic patterns matching later-generation hybrids in both hybrid zones are expected when F1 hybrids are fertile and indicate no fitness reduction in hybrids (Gramlich et al., 2016, 2018). According to our field observations, hybrids of our species pairs do form fertile catkins and seeds, but Hörandl (1992) reported for intermediate forms sometimes aborted catkins in herbarium collections. However, these observations are for both species’ pairs insufficient to draw conclusions about hybrid fitness.

In conclusion, we suggest that the biogeographical history and ecological bedrock conditions shaped the distribution of the willow species pairs, as found in other alpine species (Alvarez et al., 2009). In the secondary contact zone, hybrid zones could have influenced the recolonization process by introgression. Further studies on reproductive interference would be needed to test this hypothesis, whereby changes in hybridization dynamics over time have to be considered (Kyogoku, 2020). Our present data suggest that the incomplete range filling could have been influenced by the high genetic similarity between sister species and the hybridization that ensues from it. We propose that parent species have not yet recolonized the Alps beyond the secondary contact zone because the genetically similar sister species has already occupied suitable habitats. Other ecological factors hindering range filling, such as continued reforestation, the up-moving of alpine species due to recent global warming (Dullinger et al., 2012a), and microclimatic conditions, have to be studied. However, our work provides a new perspective on the importance of hybrid zones as an evolutionary and biogeographical process and highlights the significance of considering and studying species across their entire distribution using multiple methods.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in https://www.ncbi.nlm.nih.gov, PRJNA873074.

Author contributions

LP: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. PM: Resources, Writing – review & editing. PK: Resources, Writing – review & editing. NW: Funding acquisition, Resources, Writing – review & editing. EH: Funding acquisition, Resources, Supervision, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This study was financially supported by the German research foundation Deutsche Forschungsgemeinschaft, DFG (Ho4395/7-1 to E.H. and Wa3684/2-1 to N.W.), and by the German Academic Exchange Grant Service (DAAD) for a grant (Personal ref. no.: 91795704, Funding ID: 57552336 to LP).

Acknowledgments

We thank Johannes Wessely (University of Vienna) for providing the high-resolution climatic layers and for his advice on species distribution modelling, B. Frajman (University of Innsbruck) for a sample from Croatia, and Margaux Tréguy and Samira Wachsmuth for help with morphometrics. We thank the referees for their valuable comments.

Conflict of interest

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

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

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

Supplementary material

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

Supplementary Figure 1 | Map of the study area with different mountain systems of Europe.

Supplementary Figure 2 | Bioinformatic workflow to extract and filter loci for downstream analyses.

Supplementary Figure 3 | Evanno analysis plots for Salix alpina and Salix breviserrata.

Supplementary Figure 4 | Evanno analysis plots for Salix foetida and Salix waldsteiniana.

Supplementary Figure 5 | Result of the principal component analysis for Salix breviserrata and Salix alpina.

Supplementary Figure 6 | Result of the principal component analysis for Salix foetida and Salix waldsteiniana.

Supplementary Figure 7 | Likelihood, variance explained, and Delta m values derived from TreeMix analysis for Salix alpina and Salix breviserrata.

Supplementary Figure 8 | Likelihood, variance explained, and Delta m values derived from TreeMix analysis for Salix foetida and Salix waldsteiniana.

Supplementary Figure 9 | Principal component analysis of morphological characters for Salix alpina and Salix breviserrata.

Supplementary Figure 10 | Principal component analysis of morphological characters for Salix foetida and Salix waldsteiniana.

Supplementary Table 1 | Location and sampling details of the 352 accessions included in the article.

Supplementary Table 2 | Leaves measurements used for morphometric analysis.

Supplementary Table 3 | Vegetation relevés from 19 sampling sites.

Supplementary Table 4 | Pairwise Fst values for the groups of Salix alpina and Salix breviserrata.

Supplementary Table 5 | Pairwise Fst values for the groups of Salix foetida and Salix waldsteiniana.

References

Abbott, R., Albach, D., Ansell, S., Arntzen, J. W., Baird, S. J., Bierne, N., et al. (2013). Hybridization and speciation. J. Evol. Biol. 26, 229–246. doi: 10.1111/j.1420-9101.2012.02599.x

PubMed Abstract | Crossref Full Text | Google Scholar

Abbott, R. J., Brennan, A. C. (2014). Altitudinal gradients, plant hybrid zones and evolutionary novelty. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 369 (1648), 20130346. doi: 10.1098/rstb.2013.0346

PubMed Abstract | Crossref Full Text | Google Scholar

Aeschimann, D., Lauber, K., Moser, D. M., Theurillat, J.-P. (2004). Flora alpina: atlas des 4500 plantes vasculaires des Alpes (Bern: Haupt Verlag).

Google Scholar

Alvarez, N., Thiel-Egenter, C., Tribsch, A., Holderegger, R., Manel, S., Schönswetter, P., et al. (2009). History or ecology? Substrate type as a major driver of spatial genetic structure in Alpine plants. Ecol. Lett. 12, 632–640. doi: 10.1111/j.1461-0248.2009.01312.x

PubMed Abstract | Crossref Full Text | Google Scholar

Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data (United Kingdom: Babraham Bioinformatics, Babraham Institute, Cambridge).

Google Scholar

Argus, G. W. (1974). An experimental study of hybridization and pollination in Salix (willow). Can. J. Bot. 52, 1613–1619. doi: 10.1139/b74-212

Crossref Full Text | Google Scholar

Arnold, M. L., Martin, N. H. (2010). Hybrid fitness across time and habitats. Trends Ecol. Evol. 25, 530–536. doi: 10.1016/j.tree.2010.06.005

PubMed Abstract | Crossref Full Text | Google Scholar

Baird, N. A., Etter, P. D., Atwood, T. S., Currey, M. C., Shiver, A. L., Lewis, Z. A., et al. (2008). Rapid SNP discovery and genetic mapping using sequenced RAD markers. PloS One 3, e3376. doi: 10.1371/journal.pone.0003376

PubMed Abstract | Crossref Full Text | Google Scholar

Barton, N. H. (2001). The role of hybridization in evolution. Mol. Ecol. 10, 551–568. doi: 10.1046/j.1365-294x.2001.01216.x

PubMed Abstract | Crossref Full Text | Google Scholar

Barton, N. H., Hewitt, G. M. (1985). Analysis of hybrid zones. Annu. Rev. Ecol. Syst. 16, 113–148. doi: 10.1146/annurev.es.16.110185.000553

Crossref Full Text | Google Scholar

Blanco, P. (1993). “Salix,” in Flora Iberica. Eds. Castroviejo, S., Cirujano, S., Montserrat, P., Munoz Garmendia, F., Paiva, J. (Real Jardin Botanico, Madrid), 477–517.

Google Scholar

Büchler, W. (1986). Neue Chromosomenzählungen in der Gattung Salix. II. Botanica Helv. 96, 135–143.

Google Scholar

Carnicero, P., Wessely, J., Moser, D., Font, X., Dullinger, S., Schönswetter, P. (2022). Postglacial range expansion of high-elevation plants is restricted by dispersal ability and habitat specialization. J. Biogeogr. 49, 1739–1752. doi: 10.1111/jbi.v49.10

PubMed Abstract | Crossref Full Text | Google Scholar

Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., Cresko, W. A. (2013). Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140. doi: 10.1111/mec.2013.22.issue-11

PubMed Abstract | Crossref Full Text | Google Scholar

Christe, C., Caetano, S., Aeschimann, D., Kropf, M., Diadema, K., Naciri, Y. (2014). The intraspecific genetic variability of siliceous and calcareous Gentiana species is shaped by contrasting demographic and re-colonization processes. Mol. Phylogenet. Evol. 70, 323–336. doi: 10.1016/j.ympev.2013.09.022

PubMed Abstract | Crossref Full Text | Google Scholar

Chytrý, M., Hennekens, S. M., Jiménez-Alfaro, B., Knollová, I., Dengler, J., Jansen, F., et al. (2016). European Vegetation Archive (EVA): an integrated database of European vegetation plots. Appl. Veget. Sci. 19, 173–180. doi: 10.1111/avsc.2016.19.issue-1

Crossref Full Text | Google Scholar

Cires, E., Prieto, J. A. F. (2012). The Iberian endemic species Ranunculus cabrerensis Rothm.: an intricate history in the Ranunculus parnassiifolius L. polyploid complex. Plant Syst. Evol. 298, 121–138. doi: 10.1007/s00606-011-0529-9

Crossref Full Text | Google Scholar

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., et al. (2011). The variant call format and VCFtools. Bioinformatics 27, 2156–2158. doi: 10.1093/bioinformatics/btr330

PubMed Abstract | Crossref Full Text | Google Scholar

Davis Rabosky, A. R., Cox, C. L., Rabosky, D. L., Title, P. O., Holmes, I. A., Feldman, A., et al. (2016). Coral snakes predict the evolution of mimicry across New World snakes. Nat. Commun. 7, 11484. doi: 10.1038/ncomms11484

PubMed Abstract | Crossref Full Text | Google Scholar

Di Cola, V., Broennimann, O., Petitpierre, B., Breiner, F. T., D’Amen, M., Randin, C., et al. (2017). ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography 40, 774–787. doi: 10.1111/ecog.2017.v40.i6

Crossref Full Text | Google Scholar

Dullinger, S., Gattringer, A., Thuiller, W., Moser, D., Zimmermann, N. E., Guisan, A., et al. (2012a). Extinction debt of high-mountain plants under twenty-first-century climate change. Nat. Climate Change 2, 619–622. doi: 10.1038/nclimate1514

Crossref Full Text | Google Scholar

Dullinger, S., Willner, W., Plutzar, C., Englisch, T., Schratt-Ehrendorfer, L., Moser, D., et al. (2012b). Post-glacial migration lag restricts range filling of plants in the European Alps. Global Ecol. Biogeogr. 21, 829–840. doi: 10.1111/j.1466-8238.2011.00732.x

Crossref Full Text | Google Scholar

Durovic, S., Schönswetter, P., Niketic, M., Tomovic, G., Frajman, B. (2017). Disentangling relationships among the members of the Silene saxífraga alliance (Caryophyllaceae): Phylogenetic structure is geographically rather than taxonomically segregated. TAXON 66, 343–364. doi: 10.12705/662.4

Crossref Full Text | Google Scholar

Eaton, D. A. R., Overcast, I. (2020). ipyrad: Interactive assembly and analysis of RADseq datasets. Bioinformatics 36, 2592–2594. doi: 10.1093/bioinformatics/btz966

PubMed Abstract | Crossref Full Text | Google Scholar

Ebersbach, J., Tkach, N., Röser, M., Favre, A. (2020). The role of hybridisation in the making of the species-rich arctic-alpine genus Saxifraga (Saxifragaceae). Diversity 12, 440. doi: 10.3390/d12110440

Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Falush, D., Stephens, M., Pritchard, J. K. (2003). Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164, 1567–1587. doi: 10.1093/genetics/164.4.1567

PubMed Abstract | Crossref Full Text | Google Scholar

Fitzpatrick, B. M. (2013). Alternative forms for genomic clines. Ecol. Evol. 3, 1951–1966. doi: 10.1002/ece3.2013.3.issue-7

PubMed Abstract | Crossref Full Text | Google Scholar

Fogelqvist, J., Verkhozina, A. V., Katyshev, A. I., Pucholt, P., Dixelius, C., Rönnberg-Wästljung, A. C., et al. (2015). Genetic and morphological evidence for introgression between three species of willows. BMC Evol. Biol. 15, 193. doi: 10.1186/s12862-015-0461-7

PubMed Abstract | Crossref Full Text | Google Scholar

Fordham, D. A., Saltré, F., Haythorne, S., Wigley, T. M. L., Otto-Bliesner, B. L., Chan, K. C., et al. (2017). PaleoView: a tool for generating continuous climate projections spanning the last 21 000 years at regional and global scales. Ecography 40, 1348–1358. doi: 10.1111/ecog.2017.v40.i11

Crossref Full Text | Google Scholar

Francis, R. M. (2017). pophelper: an R package and web app to analyse and visualize population structure. Mol. Ecol. Resour. 17, 27–32. doi: 10.1111/men.2017.17.issue-1

PubMed Abstract | Crossref Full Text | Google Scholar

Gentili, R., Bacchetta, G., Fenu, G., Cogoni, D., Abeli, T., Rossi, G., et al. (2015). From cold to warm-stage refugia for boreo-alpine plants in southern European and Mediterranean mountains: the last chance to survive or an opportunity for speciation? Biodiversity 16, 247–261. doi: 10.1080/14888386.2015.1116407

Crossref Full Text | Google Scholar

Gerschwitz-Eidt, M. A., Dillenberger, M. S., Kadereit, J. W. (2023). Phylogeny of Saxifraga section Saxifraga subsection Arachnoideae (Saxifragaceae) and the origin of low elevation shade-dwelling species. Ecol. Evol. 13, e9728. doi: 10.1002/ece3.9728

PubMed Abstract | Crossref Full Text | Google Scholar

Gompert, Z., Buerkle, C. A. (2012). bgc: Software for Bayesian estimation of genomic clines. Mol. Ecol. Resour. 12, 1168–1176. doi: 10.1111/men.2012.12.issue-6

PubMed Abstract | Crossref Full Text | Google Scholar

Gramlich, S., Sagmeister, P., Dullinger, S., Hadacek, F., Hörandl, E. (2016). Evolution in situ: hybrid origin and establishment of willows (Salix L.) on alpine glacier forefields. Heredity 116, 531–541. doi: 10.1038/hdy.2016.14

PubMed Abstract | Crossref Full Text | Google Scholar

Gramlich, S., Wagner, N. D., Hörandl, E. (2018). RAD-seq reveals genetic structure of the F-2-generation of natural willow hybrids (Salix L.) and a great potential for interspecific introgression. BMC Plant Biol. 18, 317. doi: 10.1186/s12870-018-1552-6

PubMed Abstract | Crossref Full Text | Google Scholar

Guerrina, M., Theodoridis, S., Minuto, L., Conti, E., Casazza, G. (2022). First evidence of post-glacial contraction of Alpine endemics: Insights from Berardia subacaulis in the European Alps. J. Biogeogr. 49, 79–93. doi: 10.1111/jbi.14282

Crossref Full Text | Google Scholar

Guo, J.-F., Zhao, W., Andersson, B., Mao, J.-F., Wang, X.-R. (2023). Genomic clines across the species boundary between a hybrid pine and its progenitor in the eastern Tibetan Plateau. Plant Commun. 4, 100574. doi: 10.1016/j.xplc.2023.100574

PubMed Abstract | Crossref Full Text | Google Scholar

Hardig, T. M., Brunsfeld, S. J., Fritz, R. S., Morgan, M., Orians, C. M. (2000). Morphological and molecular evidence for hybridization and introgression in a willow (Salix) hybrid zone. Mol. Ecol. 9, 9–24. doi: 10.1046/j.1365-294X.2000.00757.x

PubMed Abstract | Crossref Full Text | Google Scholar

Harrison, R. G. (1993). Hybrid zones and the evolutionary process (USA: Oxford University Press).

Google Scholar

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

Crossref Full Text | Google Scholar

Hewitt, G. M. (2004). Genetic consequences of climatic oscillations in the Quaternary. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 359, 183–195; discussion 195. doi: 10.1098/rstb.2003.1388

PubMed Abstract | Crossref Full Text | Google Scholar

Holderegger, R., Thiel-Egenter, C. (2009). A discussion of different types of glacial refugia used in mountain biogeography and phylogeography. J. Biogeogr. 36, 476–480. doi: 10.1111/j.1365-2699.2008.02027.x

Crossref Full Text | Google Scholar

Hörandl, E. (1992). Die Gattung Salix in Österreich. Abh. Zool.-Bot. Ges. Österreich. 27, 1–170.

Google Scholar

Hörandl, E., Florineth, F., Hadacek, F. (2012). Weiden in Österreich und angrenzenden Gebieten (willows in Austria and adjacent regions) (Vienna: University of Agriculture).

Google Scholar

Hrivnák, M., Krajmerová, D., Paule, L., Zhelev, P., Sevik, H., Ivanković, M., et al. (2024). Are there hybrid zones in Fagus sylvatica L. sensu lato? Eur. J. For. Res. 143, 451–464. doi: 10.1007/s10342-023-01634-0

Crossref Full Text | Google Scholar

Jalas, J., Suominen, J. (1988). Atlas Florae Europaeae: volume 3: distribution of vascular plants in Europe (Helsinki: The Committee for Mapping the Flora of Europe).

Google Scholar

Jiggins, C. D., Mallet, J. (2000). Bimodal hybrid zones and speciation. Trends Ecol. Evol. 15, 250–255. doi: 10.1016/S0169-5347(00)01873-5

PubMed Abstract | Crossref Full Text | Google Scholar

Jombart, T. (2008). adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics 24, 1403–1405. doi: 10.1093/bioinformatics/btn129

PubMed Abstract | Crossref Full Text | Google Scholar

Kadereit, J. W. (2015). The geography of hybrid speciation in plants. TAXON 64, 673–687. doi: 10.12705/644.1

Crossref Full Text | Google Scholar

Kadereit, J. W. (2017). The role of in situ species diversification for the evolution of high vascular plant species diversity in the European Alps—A review and interpretation of phylogenetic studies of the endemic flora of the Alps. Perspect. Plant Ecol. Evol. Syst. 26, 28–38. doi: 10.1016/j.ppees.2017.03.002

Crossref Full Text | Google Scholar

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., et al. (2017). Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4, 170122. doi: 10.1038/sdata.2017.122

PubMed Abstract | Crossref Full Text | Google Scholar

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., et al. (2018). Data from: Climatologies at high resolution for the earth’s land surface areas. EnviDat doi: 10.16904/envidat.228.v2.1

PubMed Abstract | Crossref Full Text | Google Scholar

Karrenberg, S., Kollmann, J., Edwards, P. J., Gurnell, A. M., Petts, G. E. (2003). Patterns in woody vegetation along the active zone of a near-natural Alpine river. Basic. Appl. Ecol. 4, 157–166. doi: 10.1078/1439-1791-00123

Crossref Full Text | Google Scholar

Kerner, A. J. (1870). Die natürlichen Floren im Gelände der deutschen Alpen: éditeur non identifié. Fromann, Jena.

Google Scholar

Kirchheimer, B., Schinkel, C. C. F., Dellinger, A. S., Klatt, S., Moser, D., Winkler, M., et al. (2016). A matter of scale: apparent niche differentiation of diploid and tetraploid plants may depend on extent and grain of analysis. J. Biogeogr. 43, 716–726. doi: 10.1111/jbi.2016.43.issue-4

PubMed Abstract | Crossref Full Text | Google Scholar

Kirchheimer, B., Wessely, J., Gattringer, A., Hülber, K., Moser, D., Schinkel, C. C. F., et al. (2018). Reconstructing geographical parthenogenesis: effects of niche differentiation and reproductive mode on Holocene range expansion of an alpine plant. Ecol. Lett. 21, 392–401. doi: 10.1111/ele.2018.21.issue-3

PubMed Abstract | Crossref Full Text | Google Scholar

Kosiński, P., Boratyński, A., Hilpold, A. (2017). Taxonomic differentiation of Salix retusa agg. (Salicaceae) based on leaf characteristics. Dendrobiology 78, 40–50. doi: 10.12657/denbio.078.005

Crossref Full Text | Google Scholar

Kuzmanović, N., Lakušić, D., Frajman, B., Stevanoski, I., Conti, F., Schönswetter, P. (2021). Long neglected diversity in the Accursed Mountains (western Balkan Peninsula): Ranunculus bertisceus is a genetically and morphologically divergent new species. Bot. J. Linn. Soc. 196, 384–406. doi: 10.1093/botlinnean/boab001

Crossref Full Text | Google Scholar

Kyogoku, D. (2020). When does reproductive interference occur? Predictions and data. Popul. Ecol. 62, 196–206. doi: 10.1002/1438-390X.12041

Crossref Full Text | Google Scholar

Landolt, E., Bäumler, B., Erhardt, A., Hegg, O., Klötzli, F., Lämmler, W., et al. (2010). Flora indicativa: Ökologische Zeigerwerte und biologische Kennzeichen zur Flora der Schweiz und der Alpen (Bern: Haupt Verlag).

Google Scholar

Marinček, P., Léveillé-Bourret, É., Heiduk, F., Leong, J., Bailleul, S. M., Volf, M., et al. (2024). Challenge accepted: Evolutionary lineages versus taxonomic classification of North American shrub willows (Salix). Am. J. Bot., e16361. doi: 10.1002/ajb2.16361

PubMed Abstract | Crossref Full Text | Google Scholar

Marinček, P., Pittet, L., Wagner, N. D., Hörandl, E. (2023). Evolution of a hybrid zone of two willow species (Salix L.) in the European Alps analyzed by RAD-seq and morphometrics. Ecol. Evol. 13, e9700. doi: 10.1002/ece3.9700

PubMed Abstract | Crossref Full Text | Google Scholar

Martin, B. T., Chafin, T. K., Douglas, M. R., Douglas, M. E. (2021). ClineHelpR: an R package for genomic cline outlier detection and visualization. BMC Bioinf. 22, 1–10. doi: 10.1186/s12859-021-04423-x

PubMed Abstract | Crossref Full Text | Google Scholar

Martini, F., Paiero, P. (1988). I Salici d’Italia (Trient: Edizione Lint).

Google Scholar

McDade, L. A. (1992). Hybrids and phylogenetic systematics II. The impact of hybrids on cladistic analysis. Evolution 46, 1329–1346. doi: 10.2307/2409940

PubMed Abstract | Crossref Full Text | Google Scholar

Merxmüller, H. (1952). Untersuchungen zur Sippengliederung und Arealbildung in den Alpen. I. Jahrbuch des Vereins zum Schutze der Alpenpflanzen und Tiere. 17, 96–133.

Google Scholar

Micevski, K. (1995). The Flora of The Republic of Macedonia (Skopje: Macedonian Academy of Sciences and Arts).

Google Scholar

Milanesi, M., Capomaccio, S., Vajana, E., Bomba, L., Garcia, J. F., Ajmone-Marsan, P., et al. (2017). BITE: an R package for biodiversity analyses. BioRxiv, 181610. doi: 10.1101/181610

Crossref Full Text | Google Scholar

Moore, A. J., Kadereit, J. W. (2013). The evolution of substrate differentiation in Minuartia series Laricifoliae (Caryophyllaceae) in the European Alps: In situ origin or repeated colonization? Am. J. Bot. 100, 2412–2425. doi: 10.3732/ajb.1300225

PubMed Abstract | Crossref Full Text | Google Scholar

Mráz, P., Gaudeul, M., Rioux, D., Gielly, L., Choler, P., Taberlet, P., et al. (2007). Genetic structure of hypochaeris uniflora (Asteraceae) suggests vicariance in the carpathians and rapid post-glacial colonization of the alps from an eastern alpine refugium. J. Biogeogr. 34, 2100–2114. doi: 10.1111/j.1365-2699.2007.01765.x

Crossref Full Text | Google Scholar

Nazareno, A. G., Bemmels, J. B., Dick, C. W., Lohmann, L. G. (2017). Minimum sample sizes for population genomics: an empirical study from an Amazonian plant species. Mol. Ecol. Resour. 17, 1136–1147. doi: 10.1111/men.2017.17.issue-6

PubMed Abstract | Crossref Full Text | Google Scholar

Neumann, A. (1981). Die Mitteleuropaeischen Salix-Arten. (Forstliche Bundesversuchsanstalt, Vienna).

Google Scholar

Neumann, A., Polatschek, A. (1972). Cytotaxonomischer beitrag zur gattung Salix. Ann. Naturhistor. Mus. Wien 76, 619–633.

Google Scholar

Ozenda, P. (1988). Die Vegetation der Alpen im europäischen Gebirgsraum. (Fischer, Stuttgart).

Google Scholar

Paris, J. R., Stevens, J. R., Catchen, J. M. (2017). Lost in parameter space: a road map for stacks. Methods Ecol. Evol. 8, 1360–1373. doi: 10.1111/mee3.2017.8.issue-10

Crossref Full Text | Google Scholar

Parisod, C. (2022). Plant speciation in the face of recurrent climate changes in the Alps. Alpine. Bot. 132, 21–28. doi: 10.1007/s00035-021-00259-6

Crossref Full Text | Google Scholar

Pawłowski, B. (1970). Remarques sur l’endemisme dans la flore des Alpes et des Carpates. Vegetatio 21, 181–243.

Google Scholar

Petit, R. J., Excoffier, L. (2009). Gene flow and species delimitation. Trends Ecol. Evol. 24, 386–393. doi: 10.1016/j.tree.2009.02.011

PubMed Abstract | Crossref Full Text | Google Scholar

Pflugbeil, G., Affenzeller, M., Tribsch, A., Comes, H. P. (2021). Primary hybrid zone formation in Tephroseris helenitis (Asteraceae), following postglacial range expansion along the central Northern Alps. Mol. Ecol. 30, 1704–1720. doi: 10.1111/mec.15832

PubMed Abstract | Crossref Full Text | Google Scholar

Pickrell, J., Pritchard, J. (2012). Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 8, e1002967. doi: 10.1371/journal.pgen.1002967

PubMed Abstract | Crossref Full Text | Google Scholar

Pickup, M., Brandvain, Y., Fraïsse, C., Yakimowski, S., Barton, N. H., Dixit, T., et al. (2019). Mating system variation in hybrid zones: facilitation, barriers and asymmetries to gene flow. New Phytol. 224, 1035–1047. doi: 10.1111/nph.v224.3

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

Rambaut, A., Drummond, A. (2012). FigTree version 1.4. 0. Available online at: http://tree.bio.ed.ac.uk/software/figtree/.

Google Scholar

Rechinger, K. H. (1964). Salix in Flora Europaea, Vol. 1. Lycopodiaceae to Platanaceae. Tutin, T. G., Heywood, V. H., Burges, N.A., Valentine, D. H., Walters, S. M., Webb, D. A. (Cambridge University Press, Cambridge) pp. 43–54.

Google Scholar

Rieseberg, L. H., Carney, S. E. (1998). Plant hybridization. New Phytol. 140, 599–624. doi: 10.1046/j.1469-8137.1998.00315.x

PubMed Abstract | Crossref Full Text | Google Scholar

Rochette, N. C., Catchen, J. M. (2017). Deriving genotypes from RAD-seq short-read data using Stacks. Nat. Protoc. 12, 2640–2659. doi: 10.1038/nprot.2017.123

PubMed Abstract | Crossref Full Text | Google Scholar

Ronikier, M., Cieślak, E., Korbecka, G. (2008). High genetic differentiation in the alpine plant Campanula alpina Jacq. (Campanulaceae): evidence for glacial survival in several Carpathian regions and long-term isolation between the Carpathians and the Alps. Mol. Ecol. 17, 1763–1775. doi: 10.1111/j.1365-294X.2008.03664.x

PubMed Abstract | Crossref Full Text | Google Scholar

Rosche, C., Schrieber, K., Lachmuth, S., Durka, W., Hirsch, H., Wagner, V., et al. (2018). Sex ratio rather than population size affects genetic diversity in Antennaria dioica. Plant Biol. 20, 789–796. doi: 10.1111/plb.2018.20.issue-4

PubMed Abstract | Crossref Full Text | Google Scholar

Schmitt, T. (2017). “Molecular biogeography of the high mountain systems of europe: an overview,” in High Mountain Conservation in a Changing World. Eds. Catalan, J., Ninot, J. M., Aniz, M. M. (Springer International Publishing, Cham), 63–74.

Google Scholar

Schneeweiss, G. M., Winkler, M., Schönswetter, P. (2017). Secondary contact after divergence in allopatry explains current lack of ecogeographical isolation in two hybridizing alpine plant species. J. Biogeogr. 44, 2575–2584. doi: 10.1111/jbi.2017.44.issue-11

Crossref Full Text | Google Scholar

Schneider, C. A., Rasband, W. S., Eliceiri, K. W. (2012). NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675. doi: 10.1038/nmeth.2089

PubMed Abstract | Crossref Full Text | Google Scholar

Schönswetter, P., Schneeweiss, G. M. (2019). Is the incidence of survival in interior Pleistocene refugia (nunataks) underestimated? Phylogeography of the high mountain plant Androsace alpina (Primulaceae) in the European Alps revisited. Ecol. Evol. 9, 4078–4086. doi: 10.1002/ece3.2019.9.issue-7

PubMed Abstract | Crossref Full Text | Google Scholar

Schönswetter, P., Stehlik, I., Holderegger, R., Tribsch, A. (2005). Molecular evidence for glacial refugia of mountain plants in the European Alps. Mol. Ecol. 14, 3547–3555. doi: 10.1111/j.1365-294X.2005.02683.x

PubMed Abstract | Crossref Full Text | Google Scholar

Schönswetter, P., Tribsch, A., Niklfeld, H. (2004). Amplified Fragment Length Polymorphism (AFLP) reveals no genetic divergence of the Eastern Alpine endemic Oxytropis campestris subsp. tiroliensis (Fabaceae) from widespread subsp. campestris. Plant Syst. Evol. 244, 245–255. doi: 10.1007/s00606-003-0096-9

Crossref Full Text | Google Scholar

Šlenker, M., Koutecký, P., Marhold, K. (2022). MorphoTools2: an R package for multivariate morphometric analysis. Bioinformatics 38, 2954–2955. doi: 10.1093/bioinformatics/btac173

PubMed Abstract | Crossref Full Text | Google Scholar

Smycka, J., Roquet, C., Boleda, M., Alberti, A., Boyer, F., Douzet, R., et al. (2022). Tempo and drivers of plant diversification in the European mountain system. Nat. Commun. 13, 2750. doi: 10.1038/s41467-022-30394-5

PubMed Abstract | Crossref Full Text | Google Scholar

Spaniel, S., Resetnik, I. (2022). Plant phylogeography of the Balkan Peninsula: spatiotemporal patterns and processes. Plant Syst. Evol. 308, 38. doi: 10.1007/s00606-022-01831-1

Crossref Full Text | Google Scholar

Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. doi: 10.1093/bioinformatics/btu033

PubMed Abstract | Crossref Full Text | Google Scholar

Stehlik, I. (2000). Nunataks and peripheral refugia for alpine plants during quaternary glaciation in the middle part of the Alps. Botanica Helv. 110, 25–30.

Google Scholar

Suarez-Gonzalez, A., Lexer, C., Cronk, Q. C. B. (2018). Adaptive introgression: a plant perspective. Biol. Lett. 14, 20170688. doi: 10.1098/rsbl.2017.0688

PubMed Abstract | Crossref Full Text | Google Scholar

Surina, B., Pfanzelt, S., Einzmann, H. J. R., Albach, D. C. (2014). Bridging the Alps and the Middle East: Evolution, phylogeny and systematics of the genus Wulfenia (Plantaginaceae). TAXON 63, 843–858. doi: 10.12705/634.18

Crossref Full Text | Google Scholar

Taberlet, P., Fumagalli, L., Wust-Saucy, A. G., Cosson, J. F. (1998). Comparative phylogeography and postglacial colonization routes in Europe. Mol. Ecol. 7, 453–464. doi: 10.1046/j.1365-294x.1998.00289.x

PubMed Abstract | Crossref Full Text | Google Scholar

Thiel-Egenter, C., Alvarez, N., Holderegger, R., Tribsch, A., Englisch, T., Wohlgemuth, T., et al. (2011). Break zones in the distributions of alleles and species in alpine plants. J. Biogeogr. 38, 772–782. doi: 10.1111/j.1365-2699.2010.02441.x

Crossref Full Text | Google Scholar

Thuiller, W. G. ,. D., Gueguen, M., Engler, R., Breiner, F., Lafourcade, B., Patin, R. (2023). biomod2: Ensemble Platform for Species Distribution Modeling. R package version 4.2-1. Available online at: https://biomodhub.github.io/biomod2/.

Google Scholar

Tribsch, A., Schönswetter, P. (2003). Patterns of endemism and comparative phylogeography confirm palaeoenvironmental evidence for pleistocene refugia in the eastern alps. TAXON 52, 477–497. doi: 10.2307/3647447

Crossref Full Text | Google Scholar

Viricel, A., Pante, E., Dabin, W., Simon-Bouhet, B. (2014). Applicability of RAD-tag genotyping for interfamilial comparisons: empirical data from two cetaceans. Mol. Ecol. Resour. 14, 597–605. doi: 10.1111/men.2014.14.issue-3

PubMed Abstract | Crossref Full Text | Google Scholar

Volf, M., Leong, J. V., de Lima Ferreira, P., Volfová, T., Kozel, P., Matos-Maraví, P., et al. (2023). Contrasting levels of β-diversity and underlying phylogenetic trends indicate different paths to chemical diversity in highland and lowland willow species. Ecol. Lett. 26, 1559–1571. doi: 10.1111/ele.14273

PubMed Abstract | Crossref Full Text | Google Scholar

Wagner, N. D., Gramlich, S., Hörandl, E. (2018). RAD sequencing resolved phylogenetic relationships in European shrub willows (Salix L. subg. Chamaetia and subg. Vetrix) and revealed multiple evolution of dwarf shrubs. Ecol. Evol. 8, 8243–8255. doi: 10.1002/ece3.2018.8.issue-16

PubMed Abstract | Crossref Full Text | Google Scholar

Wagner, N. D., He, L., Hörandl, E. (2021). The evolutionary history, diversity, and ecology of willows (Salix L.) in the european alps. Diversity 13, 146. doi: 10.3390/d13040146

Crossref Full Text | Google Scholar

Wagner, N. D., Marincek, P., Pittet, L., Hörandl, E. (2023). Insights into the taxonomically challenging hexaploid alpine shrub willows of salix sections phylicifoliae and nigricantes (Salicaceae). Plants-Basel 12, 1144. doi: 10.3390/plants12051144

PubMed Abstract | Crossref Full Text | Google Scholar