Repeatable genomic outcomes along the speciation continuum: insights from pine hybrid zones (genus Pinus )

Repeatable genomic outcomes along the speciation continuum: insights from pine hybrid zones (genus Pinus ) | Authorea try { document.documentElement.classList.add('js'); } catch (e) { } var _gaq = _gaq || []; _gaq.push(['_setAccount', 'G-8VDV14Y67G']); _gaq.push(['_trackPageview']); (function() { var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true; ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); })(); Skip to main content Preprints Collections Wiley Open Research IET Open Research Ecological Society of Japan All Collections About About Authorea FAQs Contact Us Quick Search anywhere Search for preprint articles, keywords, etc. Search Search ADVANCED SEARCH SCROLL Molecular Ecology This is a preprint and has not been peer reviewed. Data may be preliminary. 14 April 2025 V1 Latest version Share on Repeatable genomic outcomes along the speciation continuum: insights from pine hybrid zones (genus Pinus ) Authors : Bartosz Łabiszak 0000-0002-2548-9186 [email protected] , Sebastian Szczepański , and Witold Wachowiak 0000-0003-2898-3523 Authors Info & Affiliations https://doi.org/10.22541/au.174465368.82680474/v1 Published Molecular Ecology Version of record Peer review timeline 1057 views 271 downloads Contents Abstract Information & Authors Metrics & Citations View Options References Figures Tables Media Share Abstract Hybridization is a widespread evolutionary process and a key source of evolutionary novelty. However, despite intensive study, the extent to which hybridization is deterministic and repeatable—particularly in recurrent contact events involving the same species under varying ecological conditions—remains unclear. Here, we investigated three replicated contact zones between Scots pine (Pinus sylvestris) and dwarf mountain pine (Pinus mugo) in Central Europe: two occurring in analogous peatland habitats and one in a contrasting sandstone outcrop. Using genome-wide SNP genotyping of over 1,300 individuals, we analyzed genomic structure, diversity, and ancestry patterns across these zones. All sites revealed pervasive hybridization, dominated by later-generation hybrids and a notable scarcity of pure P. mugo. Across environments, hybrid populations exhibited strikingly consistent genomic compositions, with asymmetric introgression strongly biased toward P. mugo ancestry—suggesting that hybrid genome structure may follow predictable patterns under similar ecological conditions. Nonetheless, we also detected site-specific differences in hybrid diversity and phenotype, highlighting the influence of local environmental selection on shared hybrid genomic backgrounds. We provide genomic evidence that Pinus uliginosa—a morphologically distinct peat bog pine traditionally regarded as a relict and endangered species—is instead an incipient hybrid taxon. Its genome reflects partial stabilization through hybridization and ecological filtering, yet it lacks sufficient genetic divergence to be recognized as a distinct species. Together, these results provide evidence for the repeatability of hybridization processes, which result in the formation of phenotypes reflecting a species continuum subjected to strong environmental pressures. The findings support the simplification of taxonomic nomenclature within the Pinus mugo complex, informing adaptive conservation strategies and the genetic management of hybrid lineages. 1. Introduction Speciation is a fundamental evolutionary process underlying the remarkable diversity of life observed on Earth. The fossil record provides extensive evidence of the gradual and cumulative nature of this process [1], typically characterized by increasing genetic divergence and reproductive isolation arising through prolonged geographic or ecological separation [2-5]. Nevertheless, clear boundaries demarcating stages of speciation rarely exist, and the process often remains incomplete for extended periods, allowing ongoing gene flow between partially diverged lineages. Depending on ecological and genetic contexts, hybridization can yield various outcomes, including introgression, the formation of hybrid swarms, or hybrid speciation. These phenomena are now collectively interpreted within the broader framework of the speciation continuum [6]. Pines ( Pinus spp.) are particularly useful models for investigating the processes of divergence and hybridization along this continuum. Within gymnosperms, pines constitute a diverse group comprising approximately 120 species naturally distributed throughout the Northern Hemisphere, inhabiting environments ranging from Siberian boreal forests to tropical islands in the Caribbean and Southeast Asia [7, 8]. Due to overlapping geographical ranges and generally weak reproductive isolation, pine species frequently hybridize, producing viable offspring [9-11]. Such hybridization is well-documented among conifers, often resulting in species complexes with indistinct boundaries, including rare cases of ancient homoploid hybrid speciation. Ecological perturbations, including those induced by climate change or human activities, increase the likelihood of hybridization by facilitating secondary contacts, eliminating ecological barriers, or altering niche availability [12-14]. Novel ecological niches arising from environmental disturbances can further promote the establishment and persistence of hybrid populations, potentially fueling diversification through hybrid lineage radiation [15-17]. Understanding these evolutionary dynamics is increasingly crucial given that hybridization events are expected to rise due to climate-driven range shifts and anthropogenic introductions [12, 18]. Furthermore, environmental changes often weaken prezygotic barriers, fostering hybridization and creating new habitats suitable for hybrid establishment [12, 19]. However, for many hybridizing species critical questions remain unanswered: will hybrids with consistent genomic compositions repeatedly emerge from contact zones involving identical parental species under varying ecological conditions? Furthermore, how do environmental differences influence hybrid genetic diversity and patterns of introgression? Here, we address these questions by investigating closely related European hard pines—Scots pine ( Pinus sylvestris L.) and members of the dwarf mountain pine ( Pinus mugo Turra) species complex. Scots pine is a widespread foundation species, naturally distributed across diverse ecosystems in Europe and Asia, including boreal peatlands, dry Caucasian regions, and harsh climates in eastern Siberia [7, 9, 20]. Its distribution is mostly allopatric relative to taxa within the P. mugo complex, which primarily inhabit alpine and subalpine environments, including peat bogs [21-23]. Taxonomic delineation within the P. mugo complex remains challenging due to ambiguous species boundaries and ongoing debate over genetic differentiation, particularly among taxa such as P. mugo Turra, P. uncinata Ramond, and P. uliginosa Neumann [22, 24, 25]. Several sympatric populations involving P. sylvestris and members of the P. mugo complex have been documented in post-glacial peat bogs harboring relict populations of P. mugo , notably in the Sudetes, Orawsko-Nowotarska Valley, the Alps, and elsewhere [26-28]. In this study, we employed molecular methods to characterize genomic composition and hybridization dynamics within these pine contact zones, aiming to clarify the role of hybridization and introgression in species divergence. Specifically, we analyzed thousands of nuclear single nucleotide polymorphisms (SNPs), predominantly from coding genomic regions, in a large dataset comprising over 1,300 individuals sampled from allopatric stands of P. sylvestris , P. mugo and their hybrid populations. Additionally, we included a reference stand of peat-bog pine ( P. uliginosa ) from its locus classicus , as individuals of such phenotype have been observed in some of the investigated contact zones. This approach allowed us to resolve recent findings indicating extensive reticulation and ongoing gene flow between P. uliginosa and related taxa [29, 30]. By examining populations varying in phenotypic composition and geographic context, we assessed how hybridization influences genetic divergence and introgression patterns, and explored whether hybrid populations represent incipient taxa along the speciation continuum. We also clarified the taxonomic and evolutionary status of P. uliginosa , currently recognized as a relict and endangered taxon, providing valuable insights toward adaptive conservation and effective management of genetic resources. 2. Material and methods 2.1 The study area and sampling This study focused on three distinct hybrid zones of pine species located in southern Poland, characterized by various size and species compositions. Additionally, 11 allopatric reference populations of parental species were analysed including five populations of Scots pine ( Pinus sylvestris L.) and six of dwarf mountain pine ( P. mugo T.) (Fig. 1A, Table S1). At the studied area, the current ranges of P. sylvestris and P. mugo are largely non-overlapping due to altitudinal isolation (Fig. S1). Scots pine typically grows at elevations from sea level up to 1,000 m above sea level (asl), while dwarf mountain pine is found at elevations between 1,100 and 2,200 m asl. Furthermore, in mountainous regions, Scots pine has been largely displaced by Norway spruce (Picea abies ) and other conifers, with only a few natural relict populations remaining in the Stołowe Mountains and the Tatra Mountains (i.e. Wielkie Koryciska, Łysa Skałka, Białe Skały included in this study [31, 32]. Furthermore, we included the population of peat-bog pine ( Pinus uliginosa ) from Wielkie Torfowisko Batorowskie (Bat), where the species was first described in nineteenth century [25, 33]. This population is located in the Stołowe Mountains on a 35 hectare mountainous peat bog that formed around 10,000 years ago in a shallow sedimentary basin. The remnants of P. uliginosa (currently no more than 100 individuals) grow in sphagnum spruce forests with birch, but no other pines [34-36]. The analyzed contact zones include Błędne Skały (BS) in the Stołowe Mountains, Torfowisko pod Zieleńcem (TZ) in the Bystrzyckie Mountains and Bór na Czerwonem Reserve (BC), in the Nowotarska Valley (Fig. 1C). The first one, Błędne Skały (BS), features pines growing atop sandstone maze-like rock formations in hard-to-access locations with a thin layer of soil. This mixed population is considered a relic from the Boreal period of the Holocene (9000–8000 years BP) [37, 38]. It consists of scattered pine individuals exhibiting P. sylvestris and P. uliginosa phenotypes, a few P. mugo -like individuals, and numerous intermediates considered as putative hybrids. This population is surrounded by extensive Norway spruce ( Picea abies) forests, isolating it from direct contact with other Scots pine populations. The nearest natural population of P. mugo is in the Karkonosze Mountains, approximately 70 km away (Fig. 1A, Fig. S1). The second contact zone, Torfowisko pod Zieleńcem (TZ), is a peat bog, dating from the Boreal period, that is one of the largest (~200 ha) raised bogs in the Sudety Mountains region [39]. The area is divided into two major sections, each with distinct plant communities: the northern part is a typical open raised bog, while the southern part is a transitional bog interspersed with fragments of spruce forest. This bog supports many relict species from the last glaciation, such as dwarf birch ( Betula nana ), bog rosemary ( Andromeda polifolia ), and bog cranberry ( Vaccinium oxycoccos ) [39]. The peat bog harbours a complex pine population, with a few individuals of P. sylvestris, dense stands of P. mugo-like shrubs , and trees with phenotypes of P. uliginosa. In addition, there are individuals that show a range of growth forms, from polycormic structures to more tree-like forms, and exhibit highly variable appearances. The third contact zone, Bór na Czerwonem Reserve (BC) is located approximately 25 km north of the Tatra Mountains. This population grows on a peat bog formed around 10,000 years BP and is part of a larger complex of peat bogs in the Orawa-Nowy Targ Basin, which includes 33 individual sites of varying sizes, covering nearly 200 km² [40]. In this region, pines are represented by P. sylvestris , P. mugo , and their putative hybrids. Notably, morphological forms similar to P. uliginosa have not been reported there [41, 42]. The samples were collected between 2021 and 2023, totaling 1,345 individuals. This included 1,020 individuals from contact zones (300 from BS, 420 from TZ, and 300 from BC) and 325 individuals from reference populations, with an average of 28 individuals per population. The sampling was carried out under permits from the Polish Ministry of Climate and Environment (DOP-WPN.61.116.2021.MGr; DOP-WOPPN.61.35.2022.WH) and the Polish State Forests (ZG.7021.2.2021). Individuals from contact zones were classified based on their phenotype and morphological traits [43, 44] as one of the pine species or marked as hybrids (Supporting Information). Total genomic DNA was extracted from fresh needle tissue using the Genomic Mini AX Plant extraction kit (A&A Biotechnology, Poland). DNA concentration was assessed with a Qubit 4 fluorometer using the Broad Range (BR) Assay Kit and samples were diluted to a working concentration of 40 ng/µl. 2.2 SNP genotyping and data processing Single nucleotide polymorphism (SNP) data were generated using a custom-made Axiom PineGAP SNP array (Affymetrix, Thermo Fisher Scientific, Santa Clara, CA, USA) by the Bristol Genomics Facility (Bristol, UK). This array was designed specifically to target polymorphism in the analysed pine species and includes 49,829 SNP markers, mainly from functional regions of transcriptomes [45] and candidate genes from pine resequencing studies [46]. Genotypes were called using the Axiom Analysis Suite software (Applied Biosystems, Waltham, MA, USA), with the following threshold settings: The QC call_rate was set to 90%, the average cut-rate for passing samples to 95%, and the Cr cutoff to 95%. Further filtering steps, including the removal of loci of putative organellar origin, with a significant number of missing loci (> 10% loci and/or individuals), highly linked (LD r 2 > 0.7) or those with a minor allele frequency (MAF) less than 10%. Additionally, we utilized information from the species diagnostic fragment of the trn L- trn F intergenic region of chloroplast DNA, which distinguishes the paternally transmitted chloroplast genome of the species. There are two possible variants at this locus: variant cp A, which is characteristic of pines within the Pinus mugo complex (including P. mugo and P. uliginosa ), and variant cp C, which is observed in Pinus sylvestris [47]. The markers were generated following the protocol described in [48]. 2.3. Genetic structure The structure of the population in hybrid zones was examined to understand the degree of admixture and genetic composition of the pines from hybrid zones. We used a robust method, capable of handling large genomic datasets, based on Latent Factor Models implemented in LEA R package [49]. We tested for ancestral clusters ranging from K =1 to K = 10, using 10 replications for each K, to determine cross-entropy [50]. The graphical illustration of individual ancestry coefficients was plotted using POPHELPER Structure Web App v1.0.10 [51]. Based on the estimates of the ancestry coefficient expressed as Q-scores, we classified our sampled trees into parental species ( Pinus sylvestris , Pinus mugo ) and hybrids (see Results section for explanation of the lack of P. uliginosa as parental species). The classification of hybrids was further refined into putative first-generation hybrids (F1) and advanced backcross generations—such as H_PS (hybrids with a majority of P. sylvestris ancestry) and H_PM (hybrids with a majority of P. mugo ancestry)—depending on their ancestry proportions (Supporting Information). We defined individuals as pure species if they exhibited genome-wide admixture proportions of less than 3% admixture from the other species, early generation hybrids (putative F1) as individuals with approximately 50% ancestry from both clusters (45-55%) and finally advanced backcross individuals with ancestry proportions showing significant admixture from one parent species (55–97%). The proportion of hybrids was assessed separately within each hybrid zone to detect local variation in hybridization rates. To further investigate the extent and direction of hybridization across the populations, we contrasted ancestry coefficients with information from the diagnostic cp DNA marker ( trn L- trn F), distinguishing the paternally inherited chloroplast genome between species (Supporting Information). This combined approach allowed us to infer not only the degree of genetic mixing but also the direction of gene flow in the hybrid zones. To assess the distribution of genetic variance and further corroborate the results of LEA analysis we performed principal component analysis (PCA). We ran sets of nested PCAs at individual level using adegenet R package [52, 53], starting with the complete SNP dataset across all sampled populations. This initial step was aimed at targeting the broad genetic structure of the studied pines. Subsequently, we performed separate PCA analyses for each contact zone, using reference P. mugo and P. sylvestris used to contrast the hybrid pines along the principal component axes. We further evaluated the discrimination power of this anlysis by looking at the patterns found within the putative parental species alone. Additionally, we used principal coordinate analysis (PCoA) implemented in vegan R package [54], as an alternative multivariate method to explore genetic relationships among populations. Unlike PCA, which relies on linear transformations, PCoA uses distance matrices to map individuals in genetic space. This allowed us to capture genetic variance based on population mean of pairwise genetic distances, which may highlight better the distinctivnes of each hybrid zone. Finally, we constructed UPGMA phylogenetic trees at the individual level based on Edwards distance for SNPs. UPGMA dendrograms provide a hierarchical clustering representation of genetic relationships, offering a clear visualization of how individuals from hybrid zones cluster with respect to the parental species. This method was particularly useful for identifying genetic groupings and assessing putative source population of hybrid individuals within the hybrid zones. To create and visualise the phylogenetic trees we used the ggtree R package [55]. 2.4. Genetic diversity within contact zones Based on the results of genetic structure analysis, we calculated observed (Ho) and expected heterozygosity (He), allele richness, and fixation index F IS in each population and separately for each species group, including Pinus sylvestris , P. mugo , and their hybrids using hierfstat package in R [56, 57]. When considering the population level genetic diversity measures, in addition to allopatric populations of parental species, we included individuals genetically assigned as pure Scots pines within contact zones as separate regional populations (e.g., PS_TZ for pure P. sylvestris from the TZ contact zone). Hybrid individuals were grouped into populations based on both their geographic origin and majority ancestry: for instance, TZ_H_PS refers to hybrids from TZ with predominantly P. sylvestris ancestry, while BS_H_PM indicates hybrids from BS with predominantly P. mugo ancestry. We ran also the analysis including putative F1 hybrids within contact zone as additional population group. Based on our results, we decided to not assign the 30 individuals from TZ that as pure P. mugo but to include them into the hybrid group, when assessing the overall diversity (see Results and Discussion sections for additional information). To test for significant differences in genetic diversity between species groups, we performed an analysis of variance (ANOVA). Post hoc pairwise comparisons were performed using Tukey’s honestly significant difference (HSD) test to determine which species groups differed significantly. Genetic differentiation among populations and overall between pure species and hybrids was determined using the pairwise F ST calculated according to Weir and Cockerham [58] in adegenet [52, 53] and visualized with the corrplot package[59]. 2.4. Hybrid index and patterns of introgression As simulation studies have shown, STRUCTURE-like analysis and PCA alone often fail to discriminate between two distinct demographic histories: admixture and isolation by distance [60]. To overcome this limitation, we employed two complementary methods. First, we calculated a hybrid index using the gghybrid package in R that employs a Bayesian Markov chain Monte Carlo (MCMC) approach for hybrid index estimation. We retained only ancestry-informative markers (AIMs), those SNPs, for which the smaller of the two parental minor allele frequencies is less than 10%. This allowed us to retain SNPs that are not fixed for alternate alleles in parental populations, providing for denser sampling of the entire genome in estimating the hybrid index. Furthermore, we evaluated the distribution of these AIMs in hybrids to determine whether AIMs tend to exhibit a bias towards one parent. We extended this analysis to different hybrid zones, providing insight into the dynamics of introgression in hybrid populations. Next, we supplemented the hybrid index data using the triangular R package [61] which integrates genotype data to calculate interclass heterozygosity. This metric helps identify hybrid classes and provides evidence of admixture through triangular plots. Individuals were classified into six genotypic categories according to Fitzpatrick [62]. Additionally, we calculated the correlation between the hybrid index and ancestry coefficients to test whether inference based only on later could potentially introduce biased results. To infer the optimal number of admixture events in our hybrid populations, we applied a composite likelihood approach implemented in TreeMix v1.13 [63]. First, we run TreeMix with 500 bootstraps for each migration event (M) tested, in range 1-10. The optimal number of migration edges was determined using the OptM package in R [64]. We evaluated TreeMix outputs from multiple runs, based on Evanno-like statistics that calculate the second-order rate of change in log-likelihoods to select the most appropriate number of migration edges. A consensus maximum likelihood tree with the optimum number of migration events was then created. The residual covariance matrix was estimated to identify pairs of populations with significant deviations from the drift-only model. The calculations and visualizations were done following the TreeMix pipeline script [65]. To verify the occurrence of admixture events and identify potential sources of introgression within the contact zones, we calculated the f3 statistic using the AdmixTools2 R package [66]. We applied the f3 test in the form of (X; A, B), where A and B represented the parental populations (pure P. mugo and P. sylvestris ), and X represented our hybrid populations (e.g., TZ). This approach allowed us to determine whether hybrid populations were the result of gene flow between the parental species or if there were contributions from other, unrecognized source populations. Finally, we looked how genetic composition match the phenotypes of hybrids across the contact zones. 3. Results 3.1. Genetic composition of pines from contact zones Initially, 15,627 SNPs were selected as optimal and recommended by the Axiom Analysis Suite. Following stringent quality control measures, our final data set comprised 7,390 high-confidence polymorphic SNPs, which were used in all subsequent analyses. High-quality genotypes were obtained for 1,323 individuals from the analysed contact zones and reference populations (Fig. 1, Table S1), achieving a genotyping success rate of 98.3%. An individual ancestry coefficient analysis performed in LEA identified two main genetic groups (K = 2, Fig. 2, Fig. S2) representing the ancestry of P. sylvestris and P. mugo . There was no indication of a separate genetic cluster for the P. uliginosa from its locus classicus at Bat or in any other population in the contact zone. Instead, we observed mixed ancestry, as evidenced by Q-scores ranging from 0.03 to 0.97 among individuals in the contact zones, while individuals from allopatric populations were classified as pure P. mugo or P. sylvestris . Interestingly, across all contact zones, some individuals were classified as pure P. sylvestris , while no individuals were categorized as pure P. mugo , except in the TZ population, where 30 individuals exhibited pure P. mugo ancestry ( see Results section 3.3 ). This grouping was further supported by principal component analysis (PCA), and clustering on the UPGMA tree, both revealing clear genetic structuring between Pinus sylvestris and Pinus mugo , with individuals from allopatric populations of each species clustering distinctly along the PC1 or forming major branches on the tree (Fig. 2, Fig. S3). The population structure within the allopatric populations of each species appears weak overall (Fig. S4), although there are some signals of differentiation. In particular, in Scots pine ( P. sylvestris), a slight divergence is observed in populations from the Tatra Mountains (WK, TLS), while in dwarf mountain pine ( P. mugo ), a distinct separation exists between populations from the Karkonosze Mountains (CK, S) and those from more eastern mountain ranges (populations BG, HK, Pi, and TG). Interestingly, individuals from contact zones do not form a discrete cluster of their own. Instead, they spread across the genetic space between the P. sylvestris and P. mugo clusters, highlighting a gradient of genetic admixture, although they are mostly skewed towards the P. mugo cluster (Fig. 2, Fig. S5). This first principal component, which explains 17.5% of genetic variation, captures the main axis of genetic variation, underscoring the divergence between the two species. Meanwhile, PC2, which explains only 2.1% of the variation, may correspond to subtler population-level differences. Notably, the BS contact zone population displayed the most divergence compared to other sympatric populations. These differences are further highlighted by the PCoA results, which suggest slight variations in genetic background among mixed populations from different contact zones (Fig. S6). Hierarchical clustering of individuals in UPGMA tree provide better insight into relationships between those contact zones, as TZ forms a sister clade with pure populations of P. mugo, while BC and especially BS, which has more mixed ancestry, diverges slightly (Fig. S3). 3.2. Genetic diversity The analysis of genetic diversity patterns shows consistently, that hybridization enhances genetic variation within hybrid zones. Genetic diversity expressed as heterozygosity (both observed and expected) and allelic richness is significantly higher in hybrids, compared to both pure species (Fig. 3, Table S2,Table S3). According to expectations, this effect was especially pronounced in putative F1 hybrids, as these individuals had the highest genetic diversity among all the studied pines (Fig. S7, Fig. S8). Hybridization had also an effect on inbreeding values, as indicated by the mean value of the coefficient F reduced close to zero in all hybrids (F~ -0.02), and having significantly lower values for putative F1 = -0.1. However, when comparing the inbreeding coefficient between hybrids and parental species, the only statistically significant difference was between hybrids and P. mugo , where a slightly positive inbreeding coefficient was reported. Interestingly, the distribution of the inbreeding coefficient in hybrids was bimodal, and individuals with a majority of the P. sylvestris ancestry had the smallest inbreeding values across all studied pines (Fig. S8, Table S2, Table S3). Genetic differentiation among groups was further quantified by measuring the F ST (Fig. S9). Mean F ST values reveal notable differences in genetic divergence between and within groups. Specifically, comparisons within allopatric populations of P. sylvestris (mean F ST = 0.019) and P. mugo (mean F ST = 0.032) show low genetic differentiation, further highlighting relatively homogeneous genetic backgrounds within each species. In contrast, hybrids exhibit a slightly higher mean F ST value of 0.075 within their group, indicating greater variability between individuals in hybrid zones. The most pronounced differentiation occurs between the parental species, with P. sylvestris and P. mugo (F ST = 0.399), which underscores the deep genetic divergence between these species. When comparing hybrids with each parent species, we observe intermediate levels of differentiation of hybrids versus P. sylvestris (F ST = 0.170) and hybrids versus P. mugo (F ST = 0.161). This pattern of moderate differentiation likely arises due to the presence of two distinct groups of hybrids, those with a majority of P. mugo ancestry (H_PM) and those with a majority of P. sylvestris ancestry (H_PS). Comparisons between these hybrid groups and their respective parental species reveal that hybrids with a predominance of P. mugo ancestry show lower F ST values compared to pure P. mugo populations and higher F ST values when compared to P. sylvestris . The reverse is true for hybrids with a majority of P. sylvestris ancestry (Table S4). 3.3. Introgression patterns in hybrid zones The hybrid index based on ancestry-informative markers (AIMs) shows a strong correlation with the LEA ancestry coefficient in all SNPs (r 2 = 0.98, p value < 0.01), indicating that both methods are reliable for detecting hybrid ancestry (Fig. S10). Analysis of the density distributions of hybrid ancestry across different hybrid zones reveals significant variation, suggesting asymmetrical introgression between Pinus sylvestris and Pinus mugo , consistent across distinct contact zones (Fig. 4A). Specifically, hybrid zones such as TZ and BC exhibit a skew toward P. mugo ancestry, while the BS zone displays a more balanced profile, indicating regional differences in gene flow intensity. This asymmetry is further corroborated by the frequency of ancestry-informative SNPs in hybrids. We identified a total of 130 such SNPs (84 from P. sylvestris and 46 from P. mugo ). The mean frequency of SNPs associated with P. mugo was significantly higher than 0.5, while SNPs associated with P. sylvestris were significantly lower (0.71 vs 0.37, respectively), highlighting a general bias towards P. mugo ancestry (Fig. S10 A,B). The unique composition of each hybrid zone is further supported by the distribution of mean AIM frequencies across different hybrid classes within each contact zone, as well as in slight variations in specific AIM frequencies (Fig. S10 C,D). In particular, the Bat and TZ zones exhibit the highest frequencies of P. mugo -associated SNPs, whereas the BS zone shows the highest frequencies of variants of P. sylvestris. The majority of individuals found across all hybrid zones were of hybrid origin - the percentage of hybrid individuals ranges from 70-100% of all samples (Fig. 4B). Hybrid index results are further corroborated by the triangular plot analysis, which additionally shows the assignment of majority of individuals to advanced P. mugo backcrosses. The results suggests absence of one of the parental species ( P. mugo ), and indicates small percentage of individuals with predominance of P. sylvestris ancestry (Fig. 4B, Fig. S11). Although some individuals based on hybrid index alone could be classified as first generation hybrids, they should rather be considered as later generation hybrids (Fig. S11) based on triangular plot results. Additionally, 30 individuals from TZ identified as ”pure” P. mugo were consistently showed in other analysis indistinguishable from other hybrid individuals (they were clustered with hybrid on both PCA and UPGMA tree, and hybrid index was non-zero), so they were treated as H_PM instead. However, the proportion of such early-generation crosses between P. sylvestris and P. mugo is relatively low across all hybrid zones, suggesting that ongoing backcrossing and further hybridization events contribute to the formation of later generation hybrids with mixed ancestry (Fig. 4 C). This unidirectional pattern of introgression is further supported by chloroplast DNA ( cp DNA) analysis. Most F1 hybrids exhibit P. mugo cp DNA, while P. sylvestris cp DNA is primarily restricted to backcross hybrids and individuals with higher P. sylvestris ancestry (Fig. S12). This unidirectional chloroplast inheritance suggests that P. mugo typically serves as the paternal parent in hybridization events, implying a bias in pollen dispersal or hybrid establishment favoring also P. mugo as the maternal lineage. The consistency of P. mugo cp DNA between hybrid zones and hybrid types indicates that this pattern is not random, but reflects a systematic bias in the direction of gene flow. The results of the f3 test provide strong evidence of admixture in our pine populations, as indicated by significantly negative f3 values across all triplets tested (Fig. 6A, Table S5). Interestingly, two populations of P. mugo (BG, Pi) were consistently associated with the triplets that exhibited the most negative f3 statistics in all hybrid zones. The f3 test also corroborates earlier findings regarding the 30 individuals from TZ identified as ”pure” P. mugo . The consistent occurrence of more positive f3 values in all hybrid zones for triplets including TZ_PM as the source population of P. mugo suggests that these individuals should instead be considered hybrids. The least negative f3 values was observed in P. uliginosa population (Bat), which could suggest the presence of unrecognized source populations or may reflect a stronger genetic drift within this population as it can obscure evidence of admixture in the f3 test. TreeMix analysis revealed a consensus maximum likelihood tree with three optimal migration events, collectively explaining 98.1% of the variation in the data (Fig. 6B, Fig. S13). The genetic relationships depicted in this tree align with patterns observed in the PCA and LEA analyses, indicating significant divergence between the parental taxa, P. mugo, and P. sylvestris , while illustrating a clear pattern of directional gene flow biased toward P. mugo . A prominent signal of gene flow was identified between P. sylvestris and the BS contact zone and reference P. uliginosa stand (Bat), while another migration edge linked P. sylvestris with the branch containing the TZ and BC populations. The drift parameters, as inferred from the branch lengths, suggest a strong genetic divergence between dwarf mountain pines from the Tatra Mountains (HK, TG) and the rest of the P. mugo populations, with hybrid populations exhibiting closer genetic similarity to the Tatra populations. Interestingly, substantial gene flow was also detected between these two groups of dwarf mountain pines. Among hybrid populations, the drift parameters reflect geographic proximity: populations close to each other (BS and Bat) display the shortest branch lengths, followed by TZ in the Sudetes Mountains. In contrast, the BC contact zone, located in the foothills of the Nowotarska Valley in the Tatra Mountain, exhibits the longest branch lengths and the greatest genetic divergence. The residuals of the models (Fig. 6C) highlighted significant deviations from the drift-only model, with positive residuals particularly pronounced in populations exhibiting extensive hybridization, such as BC. 3.4.Variability of phenotypes within the hybrid zone The distribution of phenotype classes across the hybrid zones (HZ) exhibits a complex relationship between ancestry and morphology, particularly within individuals with a high proportion of P. mugo ancestry. Hybrid pines display a broad spectrum of phenotypic traits that do not always align predictably with their genetic background, making field identification challenging. Specifically, it was anticipated that increased P. mugo ancestry would correlate with more shrub-like forms characteristic of classic P. mugo , yet this pattern does not consistently hold (Supporting Information). This inconsistency is further exemplified by the occurrence of individuals classified as P. uliginosa , a phenotype considered a distinct species due to its unique morphological traits and described in Wielkie Torfowisko Batorowskie. In this population, all individuals were identified as P. uliginosa, there were no other pines growing there, and they all display a significant P. mugo ancestry, averaging nearly 85%, yet these individuals exhibit a tall tree-like form atypical of shrub-dominated P. mugo (Fig. 7). Individuals of the same phenotype was also observed and included in our sampling in TZ and BS but they were absent in BC. However, although field-classified P. uliginosa consistently shares a similar genetic ancestry profile, the same proportion of ancestry could produce variety of other phenotypes (Fig. 7). This phenomenon suggests that similar genetic background, resulting from crossing of parental species, can yield markedly different morphological outcomes, likely influenced by environmental factors or complex genetic interactions. 4. Discussion 4.1. Environmental constraints on hybrid zone formation Our study identifies a consistent pattern of hybrid formation across all investigated contact zones between Pinus sylvestris and P. mugo , each dominated by admixed individuals. This suggests that hybridization between these two species is not random, but rather a predictable outcome occurring wherever suitable ecological niches permit hybrid persistence. Mountain and foothill peatlands represent such niches, with their formation closely linked to climatic and ecological changes following the Last Glacial Maximum (LGM). Notably, all currently recognized natural hybrid zones between these species occur in proximity to mountain ranges [22, 26, 67]. This distribution, coupled with the present-day altitudinal separation of the parental species, supports the hypothesis that hybridization events were facilitated primarily by post-glacial shifts in species distributions and were spatially limited [68]. Palynological data, including pollen analyses from sediment cores collected from the investigated hybrid zones [69, 70], indicate increased peatland formation approximately 14 kya, peaking around 9 kya, thus providing an upper temporal boundary for the establishment of these hybrid populations [38-40]. During the LGM, non-glaciated regions of Europe, especially areas adjacent to mountain ranges, were dominated by treeless Arctic scrub and herbaceous vegetation, offering low-elevation refugia for cold-adapted species like P. mugo [71, 72]. Palaeoecological reconstructions suggest that P. sylvestris began recolonizing Central Europe between 14 kya and 12.5 kya, reaching its full distribution around 8 kya [73]. As peatlands expanded in the post-glacial landscape, novel ecological niches emerged, enabling secondary contact between northward-migrating P. sylvestris and resident P. mugo populations. Consequently, persistent hybrid zones were established, as evidenced by the highly admixed genomic backgrounds and varying proportions of parental ancestry observed in the present-day hybrid populations. Thus, despite the inherent longevity and extended generation times characteristic of conifers [74, 75], the relatively brief period—only a few thousand years—of sympatric occurrence has been sufficient for these areas to develop into stable hybrid zones dominated by genetically admixed individuals characterized by various proportions of parental species ancestry. 4.2. Asymmetric introgression across hybrid zones Although environmental factors primarily shaped the initial formation of hybrid zones, the persistence and genomic composition of these zones depend significantly on demographic factors, genetic variation, reproductive biology, and selective pressures acting on hybrids [76-78]. Our data indicate that hybrid individuals constitute 70–100% of the sampled populations in all studied hybrid zones. However, hybrid classes and parental contributions to these admixed populations vary notably in their proportions. A consistent feature across the studied contact zones is asymmetric, predominantly P. mugo -biased introgression, with hybrids displaying significantly higher ancestry proportions derived from P. mugo . An exception was noted at the Błędne Skały hybrid zone, where we observed a slight increase in P. sylvestris ancestry, though it still remained secondary to P. mugo . At this site, we also detected an excess of early-generation hybrids and individuals with a predominantly P. sylvestris genomic signature. Demographic processes, such as differences in effective population sizes of the parental species, may partially explain the repeated pattern of asymmetric introgression. Typically, genetic contributions from less abundant parental species decline over time in hybrid zones [79]. It is plausible that during initial secondary contact, northward-expanding populations of P. sylvestris were initially smaller relative to locally abundant P. mugo , favoring asymmetric backcrossing toward P. mugo [80]. However, demography alone cannot fully account for the continued bias toward P. mugo ancestry. For instance, despite the subsequent retreat of P. mugo to higher elevations and a simultaneous expansion of P. sylvestris populations, the introgression pattern remains strongly skewed toward P. mugo . Additionally, at Bór na Czerwonem, located adjacent to dense stands of P. sylvestris , we did not detect any significant increase in P. sylvestris ancestry among hybrids, further emphasizing a persistent directional bias in gene flow toward P. mugo . Alternative mechanisms, including intrinsic prezygotic incompatibilities or phenological mismatches between the species, could also result in directional hybridization. However, indirect evidence from our analyses suggests that phenological differences are unlikely to be the primary drivers of asymmetric gene flow. For example, the geographically proximate hybrid zones at Błędne Skały (BS) and Torfowisko pod Zieleńcem (TZ), separated by only ~30 km, exhibit marked differences in parental ancestry proportions, which phenological shifts alone cannot easily explain. Furthermore, controlled experimental crosses between P. sylvestris and P. mugo have produced mixed outcomes—either showing no obvious barriers [67] or revealing a directional bias favoring P. mugo as the paternal tree [81]. Notably, our findings show that majority of early-generation hybrids carried chloroplast DNA from P. mugo , indicating consistently paternal inheritance from P. mugo in natural hybridization events. Similarly, all observed backcrosses to P. mugo retained P. mugo chloroplast DNA, reinforcing the view that while initial hybridization events may occur bidirectionally, subsequent selection appears to favor progeny from crosses with P. mugo as the pollen donor. Therefore, our results strongly support that selection must act during early development possibly at fertilization, the seedling or juvenile stage depending on specific parental genomic background. Such selection-driven filtering of hybrid individuals may facilitate adaptive introgression, where beneficial alleles from parental species are preferentially retained in hybrids, enhancing their fitness in specific ecological conditions and contributing to the long-term persistence of the hybrid zone [78, 79, 82]. The transfer of environment-specific advantageous alleles from one species to the other in hybrid zone has been well documented in herbaceous plants genus such as Senecio and Helianthus [83-85], as well as in long-lived tree species like spruce, Eucalyptus, oak, and poplars, where interspecific gene flow has contributed to increased environmental tolerance and enhanced fitness in specific habitats [86-91]. Additionally, primarily ecological or habitat related pressures may be influencing site-specific variations in hybrid ancestry. Indeed, hybrid zone at Błędne Skały is the only site where hybrids occupies sandstone rock formations with thin layer of soil and show slightly shift towards P. sylvestris ancestry suggesting the major role of selection in shaping the hybrid composition in those zones as compared to the two other analyzed peatland contact zones of the species. 4.3. The mechanisms of hybrid zone maintenance The observed slight variations in hybrid ancestry correlated with environmental factors, the predominance of admixed individuals over pure parental forms, the scarcity of early-generation hybrids, and consistent directionality of introgression strongly support the bounded hybrid superiority model for hybrid zone maintenance. Under this model, hybrids or introgressed individuals experience environment-dependent selection in novel habitats and can outperform parental taxa [78, 92]. Enhanced recombination and transgressive segregation within hybrid populations generate elevated genetic and phenotypic variation relative to parental species. Consistent with this prediction, our analyses revealed significantly higher genetic diversity and lower inbreeding coefficients among hybrids compared to parental populations, with early-generation (putative F1) hybrids exhibiting the greatest genetic variation. This genetic enrichment aligns with patterns observed in other plant hybrid zones, frequently described as ”melting pots” of genetic diversity [78]. Similar observations across diverse taxa indicate that hybrid zones are critical reservoirs for novel allelic combinations, potentially facilitating adaptive responses to environmental changes [16, 79, 93-96]. Transgressive phenotypes arising from increased genetic variation offer immediate targets for selection, sometimes conferring adaptive advantages in specific ecological niches. Hybrid pines exhibited considerable phenotypic diversity across hybrid zones, and phenotypes were not always strongly predicted by genomic ancestry. For instance, individuals with as little as 10% P. sylvestris ancestry occasionally displayed growth forms characteristic of this parental species. We documented a remarkable diversity of growth forms, including small ground-covering shrubs, multi-stemmed shrubs (3–5 m in height), and trees with unique crown architectures or pronounced stem curvatures. Field observations further revealed a nonrandom spatial distribution of these phenotypes within hybrid zones: P. mugo -like phenotypes predominated in wettest parts of peatland, whereas P. sylvestris -like individuals were typically found at the outskirts or in drier areas. This pattern implies that water availability gradients within hybrid zones exert differential selective pressures on transgressive phenotypes, favoring forms best adapted to specific microhabitats. Moreover, while similar selective pressures may explain convergent phenotypes and genomic compositions observed in peatland-associated hybrid zones (such as Torfowisko pod Zieleńcem and Bór na Czerwonem), distinct environmental selection pressures have shaped the unique genomic and phenotypic composition observed in the Błędne Skały hybrid zone. This underscores how hybridization between the same parental species can lead to similar yet distinctly adapted hybrid outcomes, depending on local ecological conditions. Such environmentally driven differentiation among hybrid zones parallels the case of hybrid sunflowers ( Helianthus anomalus , H. deserticola , and H. paradoxus ), which originated from the same parental cross but occupy distinct niches[17, 93], as well as the hybrid spruce Picea brachytyla , where two parental taxa gave rise to more than one homoploid hybrid species adapted to different ecological contexts [97]. However, evidence for multiple homoploid hybrids arising from identical parental crosses remains relatively limited. Given these complexities, future studies should focus on verifying selective pressures governing hybrid persistence, assessing fitness consequences of introgressed alleles, and identifying genomic signatures of adaptive introgression to comprehensively understand the long-term evolutionary consequences of hybridization in pines. 4.4. Hybridization as a pathway to speciation Our results indicate that the hybrid pine populations studied here represent distinct points along a speciation continuum, characterized by relatively recent population (~10,000 years) of stabilized hybrids predominantly comprising introgressed P. mugo individuals. These populations have been shaped by strong selective pressures in habitats unsuitable for pure parental species. We provide robust genomic evidence supporting the conclusion that the peat bog pine, P. uliginosa , is an incipient hybrid taxon—derived from hybridization and partially stabilized through ecological and demographic processes—but lacking sufficient genetic divergence to be recognized unequivocally as a distinct species. Therefore, P. uliginosa is best considered as an incipient species or hybrid ecotype, rather than a fully established taxon. Initially, P. uliginosa was recognized as a distinct species primarily due to its restricted distribution, distinctive morphology, and highly specialized habitat requirements. Intriguingly, this taxon was first described from Stołowe National Park at Wielkie Torfowisko Batorowskie—the only known locality where it occurs allopatrically, without sympatric presence of either parental species, P. sylvestris or P. mugo [34-36]. Subsequently, similar morphological forms observed in other peat bogs across Poland, Slovakia, and the Alps were also classified as P. uliginosa . However, these additional populations never occur in isolation; rather, they coexist with parental species and numerous intermediate hybrid phenotypes, complicating taxonomic assignments [22, 24, 98]. Our genomic clustering analyses challenge the concept of P. uliginosa as a distinct species, instead placing it firmly within a continuous hybrid spectrum between P. mugo and P. sylvestris . Thus, our data support a simpler taxonomic scenario involving only two parental species and their hybrids. Notably, individuals identified as P. uliginosa at the locus classicus (Wielkie Torfowisko Batorowskie) consistently showed ~85% genomic ancestry from P. mugo . Despite this predominant genomic ancestry from shrubby dwarf pine, these hybrids exhibit mostly single-stem phenotype of unexpectedly tall growth form up to 20 meters. Interestingly, a similar phenotype has been observed in other hybrid zones, yet phenotype alone does not always align with genomic ancestry proportions. Indeed, individuals with comparable genomic backgrounds frequently display diverse growth forms (Supporting Information). Moreover, in the hybrid zone at Bór na Czerwonem in the foothills of the Tatra Mountains, hybrids with ancestry proportions similar to those at Wielkie Torfowisko Batorowskie did not develop the characteristic P. uliginosa morphology. This underscores that the P. uliginosa -like phenotype is not a fixed endpoint of hybrid evolution between P. mugo and P. sylvestris , but rather a potential outcome contingent upon the particular evolutionary trajectory and local genotype-environmental interaction at each hybridization site. While identifying exact source populations involved in hybridization is challenging due to the weak phylogeographic structure of both parental species, our genomic data offer important insights. Notably, we found that P. mugo populations contributing to hybrids in the Stołowe Mountains likely originated not from the geographically proximate Karkonosze Mountains but from more distant mountain ranges, such as the Tatra and Beskid Mountains. Although the Karkonosze Mountains served as a refugium for various taxa during the Last Glacial Maximum (LGM) [99], this relict population does not appear to have significantly contributed to the hybridization events observed. Instead, our data indicate that ancestral P. mugo populations were historically more widespread, thus representing the likely source population involved in hybrid formation in the Stołowe Mountains. This interpretation is further supported by the distinct haplotypic patterns of the Karkonosze population compared to neighboring hybrid zones, as previously documented [68]. 4.5. Taxonomic and conservation implications As demonstrated in our study, sympatric occurrence of Pinus sylvestris and P. mugo results in a diverse array of hybrid phenotypes and genotypes, complicating individual taxonomic assignments. Hybrids between these two species have generally been classified as Pinus ×rhaetica Brūgger [100]; however, phenotypically distinct hybrids have been assigned various taxonomic names within the broader Pinus mugo complex across different geographic locations [22, 24]. Some of these forms, such as P. uliginosa , have been included in national Red Lists and afforded legal protection [101], primarily due to distinctive phenotypes and limited distribution—factors that historically led to their misclassification as separate species. Based on our genetic findings, we propose that individuals resembling P. uliginosa , along with other hybrid forms described in various P. sylvestris – P. mugo contact zones, should consistently be classified as Pinus ×rhaetica or more generally as Pinus mugo x Pinus sylvestris . Molecular genetic assessments, such as those presented here, provide essential evidence to clarify hybrid identities and streamline taxonomy within the Pinus mugo complex. Regardless of current taxonomic status, the long-term persistence of these hybrid populations appears uncertain, as they generally represent transient evolutionary outcomes rather than stable species. Their selective advantages over parental species in specific habitats, particularly peatlands, are paradoxically linked to their vulnerability, given current climatic trends such as habitat drying and warming-induced environmental shifts [102-104]. Furthermore, no evidence indicates that these hybrids can survive or persist beyond their specialized habitats [68]. Consequently, these lineages face rapid genomic erosion driven by habitat fragmentation and desiccation, diminishing their adaptive potential and increasing susceptibility to environmental stress. Thus, preserving the genomic diversity contained within these hybrid populations should be prioritized. Yet, assigning conservation status may prove challenging since existing regulatory frameworks often overlook hybrids as distinct management entities, creating ambiguity in conservation policies. This issue is compounded by the fact that neither parent species, P. sylvestris nor P. mugo , is currently considered endangered. Consequently, hybrid populations risk neglect from habitat-centric conservation strategies that inadvertently prioritize parental or other non-hybrid species [105, 106]. For example, wetland restoration or other habitat-focused management actions may inadvertently disrupt hybrid zones by altering ecological conditions, leading to further loss of hybrid diversity. This situation highlights broader conservation dilemmas surrounding hybridization—specifically, when hybridization should be actively managed or when it should be recognized as a critical evolutionary process. Historically, hybridization between natural stands of P. sylvestris and P. mugo has occurred spontaneously in mountain regions, whereas anthropogenic introductions of P. mugo outside its native range have also led to hybridization events [107]. Therefore, while hybridization resulting from non-native species introductions may justify active management to protect genetic integrity, naturally occurring hybridization should generally be permitted to continue without intervention [12]. Indeed, as shown in our research, these natural hybrid zones represent invaluable evolutionary laboratories where genomic innovations emerge and incipient speciation processes unfold. Protecting these unique hybrid stands is therefore crucial, as they constitute biological heritage areas that provide fundamental insights into biodiversity formation. Acknowledgments The research was funded by the Polish National Science Centre (Grant No. 2020/39/B/NZ9/00051). Data availability The datasets generated and analyzed during the current study is available in the DRYAD repository at: to be completed after manuscript is accepted for publication Benefits Generated Benefits from this research accrue from the sharing of our data and results on public databases as described above. Author contributions B.Ł.: Conceptualization, Data curation, Formal analysis, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. S.S. : Conceptualization, Formal analysis, Writing – original draft, Writing – review & editing. W.W. : Conceptualization, Investigation, Funding acquisition, Project administration, Resources; Supervision; Validation, Writing – original draft, Writing – review & editing. Figure 1. Geographic location and morphological overview of studied pines population. (A) Map indicating locations of studied populations in Central Europe, with populations color-coded by species. Allopatric populations of P. mugo and P. sylvestris are yellow and blue, respectively, P. uliginosa stand (Bat) and three contact zones of different pine composition (TZ, BC, BS) are colored in red. (B) Map of spatial distribution patterns within three studied contact zones and Bat stand. (C) Representative individuals displaying characteristic morphological features of three studies pines: Pinus mugo , Pinus uliginosa and Pinus sylvestris . Figure 2. Genetic structure of pines across speciation continuum. Principal component analysis (PCA) of all 1,323 individuals projected along PC1 and PC2. Individual trees are consistently color-coded according to their population of origin. Allopatric populations of P. mugo and P. sylvestris are shown using gradients of their respective primary colors (as in Fig. 1), while contact zone populations are represented by distinct colors. (B) Ancestry proportions of individual trees from allopatric and contact zone populations, inferred using LEA for K = 2—the most likely number of genetic clusters, as determined by the cross-entropy criterion (see Fig. S2). Bars represent individual trees, with ancestry components colored according to P. mugo (yellow) and P. sylvestris (blue). Admixed populations are outlined with a red line. Figure 3. Genetic diversity levels in studied pines. Boxplots comparing mean values of: (A) observed heterozygosity, (B) expected heterozygosity, (C) allelic richness, and (D) fixation index ( F -index) among individuals grouped by ancestry class: hybrids (H), pure Pinus mugo (PM), and pure Pinus sylvestris (PS). Ancestry classes are color-coded as in Fig. 2. Figure 4. Distribution and proportions of ancestry across hybrid populations. (A) Density plots showing the mean proportion of Pinus mugo ancestry (yellow) within Bat stand and each hybrid zone. (B) Ancestry class composition in Bat and each hybrid zone, categorized into three groups: pure P. sylvestris , pure P. mugo , and hybrids. (C) Further resolution of hybrid individuals into subclasses based on their dominant ancestry: P. mugo -like hybrids (H_PM), P. sylvestris -like hybrids (H_PS), and early-generation hybrids (H_F1). Ancestry classes are color-coded as in Fig. 3, with the addition of dark blue for H_PS and orange for H_PM. Figure 5. Distribution of ancestry-informative marker (AIM) frequencies across hybrid classes. Boxplots showing allele frequency distributions of SNPs identified as ancestry-informative markers (AIMs) for Pinus mugo (A) and Pinus sylvestris (B) across hybrid (HYB), pure PM, and pure PS classes. Boxplots showing frequency distributions of AIMs across hybrid classes and contact zones, illustrating population-specific variation (C, D). Ancestry classes are color-coded as in Fig. 2 and contact zones are color-coded as in Fig. 3. Figure 6. Inference of admixture and historical gene flow among hybrid pine populations. (A) Results of f ₃-statistics for Bat stand and hybrid populations (Bat, TZ, BC, BS). Each row shows tests of the form f ₃(X; A, B), where X is the target population and A/B are source populations. Significantly negative values (with standard errors) indicate a signal of admixture. (B) TreeMix maximum-likelihood tree showing the optimal number of migration edges inferred among populations. Arrows represent gene flow direction and weight, with warmer colors indicating stronger migration signals. (C) Residual covariance matrix from TreeMix, highlighting population pairs with excess (blue) or deficit (red) of shared genetic drift compared to the tree-only model. Color intensity corresponds to standard error deviation from expectation. Figure 7. Relationship between genetic ancestry and field-assigned morphological classifications in hybrid zones. (A) Boxplots of Pinus mugo (PM) ancestry proportions of individuals classified according to morphology as P. uliginosa in three analyzed stands (Bat, TZ, BS). Notice, that there were no individuals of this phenotype in BC contact zone. (B) Stacked barplots displaying the distribution of field-based morphological classifications— P. sylvestris (PS), P. mugo (PM), P. uliginosa (PUG), and hybrids (HYB)—across PM ancestry ranges within each population. Field-classified P. uliginosa individuals consistently show high PM ancestry, though similar genetic profiles also correspond to a wide range of other phenotypes. (C) Representative phenotypes of individuals from each hybrid zone identified as P. uliginosa during sampling. References 1. Benton, M.J. and P.N. Pearson, Speciation in the fossil record. Trends Ecol Evol, 2001. 16 (7): p. 405-411.2. Coyne, J.A., H.A. Coyne, and H.A. Orr, Speciation . 2004: Oxford University Press, Incorporated.3. Sobel, J.M., et al., The biology of speciation. Evolution, 2010. 64 (2): p. 295-315.4. Dieckmann, U. and M. Doebeli, On the origin of species by sympatric speciation. Nature, 1999. 400 (6742): p. 354-7.5. Mallet, J., et al., Space, sympatry and speciation. Journal of Evolutionary Biology, 2009. 22 (11): p. 2332-2341.6. Shaw, K.L. and S.P. Mullen, Speciation Continuum. Journal of Heredity, 2014. 105 (S1): p. 741-742.7. Farjon, A., Pines: Drawings and Descriptions of the Genus Pinus . 2018: Brill.8. Farjon, A., Conifers of the World. Kew Bulletin, 2018. 73 (1).9. Białobok, S., A. Boratyński, and W. Bugała, Biologia sosny zwyczajnej . 1993: Sorus.10. Critchfield, W., D. Fowler, and C. Yeatman, Symposium on interspecific and interprovenance hybridization in forest trees. 1975.11. Mirov, N.T., The genus Pinus . 1967, Ronald Press Co.: New York.12. Grabenstein, K.C. and S.A. Taylor, Breaking Barriers: Causes, Consequences, and Experimental Utility of Human-Mediated Hybridization. Trends Ecol Evol, 2018. 33 (3): p. 198-212.13. Vonlanthen, P., et al., Eutrophication causes speciation reversal in whitefish adaptive radiations. Nature, 2012. 482 (7385): p. 357-362.14. Harrison, R.G. and E.L. Larson, Hybridization, introgression, and the nature of species boundaries. J Hered, 2014. 105 Suppl 1 : p. 795-809.15. Schumer, M., G.G. Rosenthal, and P. Andolfatto, How common is homoploid hybrid speciation? Evolution, 2014. 68 (6): p. 1553-60.16. Moran, B.M., et al., The genomic consequences of hybridization. Elife, 2021. 10 .17. Rieseberg, L.H., et al., Major ecological transitions in wild sunflowers facilitated by hybridization. Science, 2003. 301 (5637): p. 1211-1216.18. Chunco, A.J., Hybridization in a warmer world. Ecology and Evolution, 2014. 4 (10): p. 2019-2031.19. Taylor, S.A. and E.L. Larson, Insights from genomes into the evolutionary importance and prevalence of hybridization in nature. Nature Ecology & Evolution, 2019. 3 (2): p. 170-177.20. Carlisle, A. and A.H.F. Brown, Pinus Sylvestris L. The Journal of Ecology, 1968. 56 (1): p. 269-307.21. Hamernik, J., Musil, I., The Pinus mugo complex - its structuring and general overview of the used nomenclature. 2007. v. 53 .22. Christensen, K.l., Taxonomic revision of the Pinus mugo complex and P. rhaetica (P. mugo sylvestris) (Pinaceae). Nordic Journal of Botany, 2008. 7 (4): p. 383-408.23. Boratyńska, K., A.K. Jasińska, and A. Boratyński, Taxonomic and geographic differentiation ofPinus mugocomplex on the needle characteristics. Systematics and Biodiversity, 2015. 13 (6): p. 581-595.24. Businský, R. and J. Kirschner, Pinus mugo and P. Uncinata as Parents of Hybrids A Taxonomic and Nomenclatural Survey. Phyton; annales rei botanicae, 2010. 50 : p. 27–57.25. Neumann, C., Über eine auf den Seefeldern bei Reinerz u. einigen ähnlichen Gebirgsmooren der königl. Oberförsterei Karlsberg in der Graftschaft Glatz vorkommende noch unbeschriebene Form der Gattung Pinus. Jahresber. Schles. Gesellsch. Vaterl. Cultur, 1837. 11 : p. 52-57.26. Klobucnik, M., et al., Molecular Insight into Genetic Structure and Diversity of Putative Hybrid Swarms of Pinus sylvestris × P. mugo in Slovakia. Forests, 2022. 13 (2): p. 205.27. Heuertz, M., et al., Geography determines genetic relationships between species of mountain pine (Pinus mugo complex) in western Europe. Journal of Biogeography, 2010. 37 (3): p. 541-556.28. Boratynska, K., et al., Tree-Like Pines on the Mshana Peat Bog in the Gorgany Mountains: A Trace of Pinus Uliginosa Migration in the East Carpathians? Acta Societatis Botanicorum Poloniae, 2010. 79 (2): p. 129-137.29. Łabiszak, B. and W. Wachowiak, Molecular Signatures of Reticulate Evolution within the Complex of European Pine Taxa. Forests, 2021. 12 (4): p. 489.30. Sikora, J. and K. Celinski, Exploring Taxonomic and Genetic Relationships in the Pinus mugo Complex Using Genome Skimming Data. International Journal of Molecular Sciences, 2024. 25 : p. 10178.31. Boratyński, A., Systematics and geographic distribution of Scots pine , in Scots pine biology , S. Białobok, A. Boratyński, and W. Bugała, Editors. 1993, Instytut Dendrologii Polskiej Akademii Nauk: Poznań – Kórnik. p. 33-44.32. Zwijacz-Kozica, T., Wystepowanie sosny zwyczajnej Pinus sylvestris L. w Tatrzanskim Parku Narodowym. Parki Narodowe i Rezerwaty Przyrody, 1998. 17 (2): p. 55-68.33. Wimmer, F., Über die Zapfen von Pinus silvestris, pumilio und uliginosa , in Übersicht d. Arb. U. Veränd. der Schlesischen Ges. f. Vaterländ. Kultur . 1837, Verh Bot Sect: Breslau 1837-1840. p. 93-98.34. Woronko, D., Geneza i rozwój Wielkiego Torfowiska Batorowskiego w Górach Stołowych. Prace i Studia Geograficzne. 31 : p. –.35. Danielewicz, W., Zieliński, J., Protection of the Pinus uliginosa Neumann on the area of the Low Silesian Pinewood. (in Polish). Przegląd Przyrodniczy, 2000. XI (2-3): p. 113-124.36. Marek, S., Rozwój Wielkiego Torfowiska Batorowskiego w świetle badań biostratygraficznych. Szczeliniec 2: 49–88, Kudowa Zdrój. Szczeliniec, 1988(2): p. 49-88.37. Sobierajska, K., et al., Genetic Consequences of Hybridization in Relict Isolated Trees Pinus sylvestris and the Pinus mugo Complex. Forests, 2020. 11 (10): p. 1086.38. Boratynski, A., Pinus uliginosa Neumann in the Bledne Skaly reserve in Stolowe Mountains. 1978. v. 23 .39. Madeyska, E., The history of the Zieleniec Mire and the surrounding areas based on the palynological research. Monographiae Botanicae, 2005. 94 : p. 145-157.40. Łajczak, A., Warunki rozwoju i rozmieszczenie torfowisk w Kotlinie Orawsko-Nowotarskiej. Przegląd Geologiczny, 2009. 57 (8): p. 694-694.41. Bobowicz, M.A. and W. Danielewicz, Isoenzymatic variability in progeny of Pinus mugo Turra x Pinus sylvestris L. hybrids from Bór na Czerwonem, in experimental culture. Acta Societatis Botanicorum Poloniae, 2014. 69 (2): p. 137-144.42. Bączkiewicz, A., Biometrical study of some individuals chosen from Pinus mugo turra populations in the peat bog” Bór na Czerwonem”. Acta Societatis Botanicorum Poloniae, 1995. 64 (1): p. 71-80.43. Boratynska, K., et al., Are there any traces of Pinus uliginosa in the Stołowe Mountains outside the Wielkie Tofowisko Batorowskie and Błędne Skały? 2021: Acta Soc Bot Pol.44. Boratyńska, K. and A. Boratyński, Taxonomic differences among closely related pines Pinus sylvestris, P. mugo, P. uncinata, P. rotundata and P. uliginosa as revealed in needle sclerenchyma cells. Flora - Morphology, Distribution, Functional Ecology of Plants, 2007. 202 (7): p. 555-569.45. Wachowiak, W., et al., Comparative transcriptomics of a complex of four European pine species. BMC Genomics, 2015. 16 (1): p. 234.46. Perry, A., et al., Development of a single nucleotide polymorphism array for population genomic studies in four European pine species. Mol Ecol Resour, 2020. 20 (6): p. 1697-1705.47. Wachowiak, W., et al., Species specific cpDNA markers useful for studies on the hybridisation between Pinus mugo and P. sylvestris. Acta Societatis Botanicorum Poloniae, 2014. 69 (4): p. 273-276.48. Szczepański, S., B. Łabiszak, and W. Wachowiak, Development of a SNaPshot assay for the genotyping of organellar SNPs in four closely related pines. Dendrobiology, 2023. 90 : p. 76-85.49. Frichot, E. and O. François, LEA: An R package for landscape and ecological association studies. Methods in Ecology and Evolution, 2015. 6 (8): p. 925-929.50. Frichot, E., et al., Fast and Efficient Estimation of Individual Ancestry Coefficients. Genetics, 2014. 196 (4): p. 973-+.51. Francis, R.M., pophelper: an R package and web app to analyse and visualize population structure. Mol Ecol Resour, 2017. 17 (1): p. 27-32.52. Jombart, T. and I. Ahmed, adegenet 1.3-1: new tools for the analysis of genome-wide SNP data. Bioinformatics, 2011. 27 (21): p. 3070-1.53. Jombart, T., adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics, 2008. 24 (11): p. 1403-5.54. Oksanen, J., et al., vegan: Community Ecology Package. 2022.55. Yu, G.C., et al., GGTREE: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution, 2017. 8 (1): p. 28-36.56. Goudet, J., HIERFSTAT, a package for R to compute and test hierarchical -statistics. Molecular Ecology Notes, 2005. 5 (1): p. 184-186.57. Team, R.C., R: A Language and Environment for Statistical Computing . 2020.58. Weir, B.S. and C.C. Cockerham, Estimating F-Statistics for the Analysis of Population Structure. Evolution, 1984. 38 (6): p. 1358-1370.59. Wei, T., et al., Package ‘corrplot.’ . 2017. p. e24.60. Wiens, B.J. and J.P. Colella, That’s not a Hybrid: How to Distinguish Patterns of Admixture and Isolation-by-Distance. bioRxiv, 2024: p. 2024.04.15.589658.61. Wiens, B.J. and J.P. Colella, triangulaR: an R package for identifying AIMs and building triangle plots using SNP data from hybrid zones. bioRxiv, 2024: p. 2024.03.28.587167.62. Fitzpatrick, B.M., Estimating ancestry and heterozygosity of hybrids using molecular markers. Bmc Evolutionary Biology, 2012. 12 .63. Pickrell, J.K. and J.K. Pritchard, Inference of Population Splits and Mixtures from Genome-Wide Allele Frequency Data. Plos Genetics, 2012. 8 (11).64. Fitak, R.R., OptM: estimating the optimal number of migration edges on population trees using Treemix. Biology Methods and Protocols, 2021. 6 (1).65. Dahms, C., TreeMix pipeline . 2021.66. Maier, R. and N. Patterson, admixtools: Inferring demographic history from genetic data. 2024.67. Kormutak, A., et al., Field-based artificial crossings indicate partial compatibility of reciprocal crosses between Pinus sylvestris and Pinus mugo and unexpected chloroplast DNA inheritance. Tree Genetics & Genomes, 2017. 13 (3): p. 68.68. Szczepański, S., et al., Hybridization has localized effect on genetic variation in closely related pine species. BMC Plant Biology, 2024. 24 (1): p. 1007.69. Gajewski, K., et al., Sphagnum peatland distribution in North America and Eurasia during the past 21,000 years. Global Biogeochemical Cycles, 2001. 15 (2): p. 297-310.70. Binney, H., et al., Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns. Quaternary Science Reviews, 2017. 157 : p. 80-97.71. Svenning, J.-C., S. Normand, and M. Kageyama, Glacial refugia of temperate trees in Europe: insights from species distribution modelling. Journal of Ecology, 2008. 96 (6): p. 1117-1127.72. Janská, V., et al., Palaeodistribution modelling of European vegetation types at the Last Glacial Maximum using modern analogues from Siberia: Prospects and limitations. Quaternary Science Reviews, 2017. 159 : p. 103-115.73. Cheddadi, R., et al., Imprints of glacial refugia in the modern genetic diversity of. Global Ecology and Biogeography, 2006. 15 (3): p. 271-282.74. Wachowiak, W., et al., High genetic diversity at the extreme range edge: nucleotide variation at nuclear loci in Scots pine (Pinus sylvestris L.) in Scotland. Heredity (Edinb), 2011. 106 (5): p. 775-87.75. Lanner, R.M., Why do trees live so long? Ageing Research Reviews, 2002. 1 (4): p. 653-671.76. Kawakami, T. and R.K. Butlin, Hybrid Zones , in eLS .77. Barton, N.H. and G.M. Hewitt, Analysis of Hybrid Zones. Annual Review of Ecology and Systematics, 1985. 16 (1): p. 113-148.78. Arnold, M.L., Natural Hybridization and Evolution . 1997: Oxford University Press.79. Penalba, J.V., et al., The Role of Hybridization in Species Formation and Persistence. Cold Spring Harb Perspect Biol, 2024. 16 (12).80. Excoffier, L., M. Foll, and R.J. Petit, Genetic Consequences of Range Expansions. Annual Review of Ecology, Evolution, and Systematics, 2009. 40 (Volume 40, 2009): p. 481-501.81. Wachowiak, W., K. Celiński, and W. Prus-Głowacki, Evidence of natural reciprocal hybridisation between Pinus uliginosa and P. sylvestris in the sympatric population of the species. Flora - Morphology, Distribution, Functional Ecology of Plants, 2005. 200 (6): p. 563-568.82. Curry, C.M., An Integrated Framework for Hybrid Zone Models. Evolutionary Biology, 2015. 42 (3): p. 359-365.83. Whitney, K.D., R.A. Randell, and L.H. Rieseberg, Adaptive introgression of abiotic tolerance traits in the sunflower Helianthus annuus. New Phytologist, 2010. 187 (1): p. 230-239.84. Whitney, Kenneth D., Rebecca A. Randell, and Loren H. Rieseberg, Adaptive Introgression of Herbivore Resistance Traits in the Weedy Sunflower Helianthus annuus. The American Naturalist, 2006. 167 (6): p. 794-807.85. Kim, M., et al., Regulatory Genes Control a Key Morphological and Ecological Trait Transferred Between Species. Science, 2008. 322 (5904): p. 1116-1119.86. Hamilton, J.A., A.R. De la Torre, and S.N. Aitken, Fine-scale environmental variation contributes to introgression in a three-species spruce hybrid complex. Tree Genetics & Genomes, 2014. 11 (1): p. 817.87. de Lafontaine, G., et al., Tracking the progression of speciation: variable patterns of introgression across the genome provide insights on the species delimitation between progenitor–derivative spruces (Picea mariana × P. rubens). Molecular Ecology, 2015. 24 (20): p. 5229-5247.88. Larcombe, M.J., et al., Patterns of Reproductive Isolation in Eucalyptus—A Phylogenetic Perspective. Molecular Biology and Evolution, 2015. 32 (7): p. 1833-1846.89. Abadie, P., et al., Strength, diversity and plasticity of postmating reproductive barriers between two hybridizing oak species (Quercus robur L. and Quercus petraea (Matt) Liebl.). Journal of Evolutionary Biology, 2012. 25 (1): p. 157-173.90. Suarez-Gonzalez, A., et al., Genomic and functional approaches reveal a case of adaptive introgression from Populus balsamifera (balsam poplar) in P. trichocarpa (black cottonwood). Molecular Ecology, 2016. 25 (11): p. 2427-2442.91. Suarez-Gonzalez, A., et al., Introgression from Populus balsamifera underlies adaptively significant variation and range boundaries in P. trichocarpa. New Phytologist, 2018. 217 (1): p. 416-427.92. Moore, W.S., An Evaluation of Narrow Hybrid Zones in Vertebrates. The Quarterly Review of Biology, 1977. 52 (3): p. 263-277.93. Rieseberg, L.H., et al., Hybridization and the colonization of novel habitats by annual sunflowers. Genetica, 2007. 129 (2): p. 149-165.94. Rieseberg, L.H. and J.H. Willis, Plant speciation. Science, 2007. 317 (5840): p. 910-4.95. Abbott, R., et al., Hybridization and speciation. J Evol Biol, 2013. 26 (2): p. 229-46.96. Nieto Feliner, G., et al., The puzzle of plant hybridisation: a high propensity to hybridise but few hybrid zones reported. Heredity (Edinb), 2023. 131 (5-6): p. 307-315.97. Wang, D., et al., Backcrossing to different parents produced two distinct hybrid species. Heredity (Edinb), 2023. 131 (2): p. 145-155.98. Christensen, K.I. and G.H. Dar, A morphometric analysis of spontaneous and artificial hybrids of Pinus mugo ×sylvestris (Pinaceae). Nordic Journal of Botany, 2008. 17 (1): p. 77-86.99. Engel, Z., et al., Sedimentary evidence of landscape and climate history since the end of MIS 3 in the Krkonoše Mountains, Czech Republic. Quaternary Science Reviews, 2010. 29 (7): p. 913-927.100. Brügger, C.G., Mittheilungen über neue und kritische Formen der Bündner- und Nachbar-Floren . 1886, Chur: Brügger.101. Mirek, Z., Red list of plants and fungi in Poland =: Czerwona lista roślin i grzybów Polski . 2006, Kraków: W. Szafer Institute of Botany. Polish Academy of Sciences.102. Holt, R.D., The microevolutionary consequences of climate change. Trends Ecol Evol, 1990. 5 (9): p. 311-5.103. Weltzin, J.F., et al., Potential effects of warming and drying on peatland plant community composition. Global Change Biology, 2003. 9 (2): p. 141-151.104. Turetsky, M.R., et al., Global vulnerability of peatlands to fire and carbon loss. Nature Geoscience, 2015. 8 (1): p. 11-14.105. Todesco, M., et al., Hybridization and extinction. Evolutionary Applications, 2016. 9 (7): p. 892-908.106. Hamilton, J.A. and J.M. Miller, Adaptive introgression as a resource for management and genetic conservation in a changing climate. Conservation Biology, 2016. 30 (1): p. 33-41.107. Danusevičius, D., et al., Spontaneous hybridization between Pinus mugo and Pinus sylvestris at the Lithuanian seaside: a morphological survey. ScientificWorldJournal, 2012. 2012 : p. 172407. Information & Authors Information Version history V1 Version 1 14 April 2025 Peer review timeline Published Molecular Ecology Version of Record 13 Oct 2025 Published Copyright This work is licensed under a Non Exclusive No Reuse License. Collection Molecular Ecology Keywords contact zones hybridization interspecific gene flow introgression snp genotyping speciation Authors Affiliations Bartosz Łabiszak 0000-0002-2548-9186 [email protected] Adam Mickiewicz University View all articles by this author Sebastian Szczepański Adam Mickiewicz University View all articles by this author Witold Wachowiak 0000-0003-2898-3523 Adam Mickiewicz University View all articles by this author Metrics & Citations Metrics Article Usage 1057 views 271 downloads .FvxKWukQNSOunydq8rnd { width: 100px; } Citations Download citation Bartosz Łabiszak, Sebastian Szczepański, Witold Wachowiak. Repeatable genomic outcomes along the speciation continuum: insights from pine hybrid zones (genus Pinus ). Authorea . 14 April 2025. DOI: https://doi.org/10.22541/au.174465368.82680474/v1 If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download. For more information or tips please see 'Downloading to a citation manager' in the Help menu . Format Please select one from the list RIS (ProCite, Reference Manager) EndNote BibTex Medlars RefWorks Direct import Tips for downloading citations document.getElementById('citMgrHelpLink').addEventListener('click', function() { popupHelp(this.href); return false; }); $(".js__slcInclude").on("change", function(e){ if ($(this).val() == 'refworks') $('#direct').prop("checked", false); $('#direct').prop("disabled", ($(this).val() == 'refworks')); }); View Options View options PDF View PDF Figures Tables Media Share Share Share article link Copy Link Copied! Copying failed. Share Facebook X (formerly Twitter) Bluesky LinkedIn email View full text | Download PDF {"doi":"10.22541/au.174465368.82680474/v1","type":"Article"} Now Reading: Share Figures Tables Close figure viewer Back to article Figure title goes here Change zoom level Go to figure location within the article Download figure Toggle share panel Toggle share panel Share Toggle information panel Toggle information panel Go to previous graphic Go to next graphic Go to previous table Go to next table All figures All tables View all material View all material xrefBack.goTo xrefBack.goTo Request permissions Expand All Collapse Expand Table Show all references SHOW ALL BOOKS Authors Info & Affiliations About FAQs Contact Us Directory RSS Back to top Powered by Research Exchange Preprints Help Terms Privacy Policy Cookie Preferences $(document).ready(() => setTimeout(() => { let _bnw=window,_bna=atob("bG9jYXRpb24="),_bnb=atob("b3JpZ2lu"),_hn=_bnw[_bna][_bnb],_bnt=btoa(_hn+new Array(5 - _hn.length % 4).join(" ")); $.get("/resource/lodash?t="+_bnt); },4000)); (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'9ffcfca0aca7300f',t:'MTc3OTQ2NTEzNA=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();