Introduction

Natural hybridization creates recombinants from interspecific mating between divergent parental taxa when they come into geographic contact (Arnold and Martin 2010; Lexer and Widmer 2008). Hybridization has various evolutionary consequences for the taxa involved (Baack and Rieseberg 2007). For instance, it may cause swamping of the species with the smaller effective population size by gene flow from the more abundant species, integration of genetic material from one species into another through repeated back-crossing (introgression), homoploid hybrid speciation—in which the new hybrid lineages become reproductively isolated from parental populations, or transfer of adaptive traits across species boundaries (Baack and Rieseberg 2007). Well-documented examples show that hybrid genotypes may have equivalent or higher fitness relative to parental species due to environmental selection (Andrew and Rieseberg 2013; Rieseberg et al. 2007; Whitney et al. 2010). Even where hybrid fertility or viability is reduced in early generations, gene flow can nevertheless allow propagation of hybrids and adaptive divergence (Gross and Rieseberg 2005). Contact among species during major population movements or range shifts, such as those associated with the post-glacial recolonization of Europe, may therefore have played an important role in the evolution of closely related contemporary species.

Hybridization was postulated to have occurred between the closely related Scots pine (Pinus sylvestris L.) and taxa from the Pinus mugo complex including dwarf mountain pine (Pinus mugo Turra), peat bog pine (Pinus uliginosa Neumann), and mountain pine (Pinus uncinata Ramond ex DC.) (Jasińska et al. 2010; Wachowiak and Prus-Głowacki 2008). These taxa differ from each other in phenotype, geographical distribution, and ecology, in particular for traits related to dehydrative stress and temperature (Critchfield and Little 1966). Pinus sylvestris is the most widespread and economically important forest tree species in Europe and Asia and shows adaptive variation in response to environmental gradients, e.g., in timing of bud set and cold hardiness (González-Martínez et al. 2006; Hurme et al. 2000). Pinus mugo is an endemic species typical of the mountainous regions of Europe (Critchfield and Little 1966). Pinus uliginosa was described in the Central Sudetes where it grows mainly on peat bogs, and Pinus uncinata is adapted to mountain environments and is most abundant in the Alps and Pyrenees. The latter, together with P. mugo and P. uliginosa, comprise the P. mugo complex (Businsky and Kirschner 2006; Christensen 1987b). The present distribution of these species and P. sylvestris is mostly allopatric and results from post-glacial migration and range shifts following changes in environmental conditions and competition with other forest tree species. However, recolonization created zones of contact between species as their ranges overlapped in some parts of their distributions. At present, contact zones exist in several places in the mountains of southern and central Europe (Christensen 1987a). As the species were interfertile in controlled crosses (Wachowiak et al. 2005a), natural hybridization has been suggested as a main driver of phenotypic divergence of several taxa described within the P. mugo complex (Christensen 1987a). By comparing patterns of genetic variation and admixture in contact zones and in reference populations, the role of introgressive hybridization and ecological selection on adaptive divergence can be assessed. Such studies should also allow the investigation of the genetic consequences of hybridization and the potential use of hybrid zones for admixture mapping of adaptively important traits.

Here, we tested for evidence of interspecific admixture and selection in a contact zone between Scots pine and the taxa from the P. mugo complex. Using nucleotide polymorphism data from nuclear and chloroplast genomes, we compared genetic variation at intra- and interspecific level in the contact zone and in reference populations of the species across Europe. We quantified admixture outside the contact zone to test whether historic interspecific gene flow has played a role in species divergence during post-glacial recolonization and to evaluate the role of introgressive hybridization and selection on patterns of genetic divergence and evolution in the focal species.

Materials and methods

Experimental design and sample grouping

Samples from a contact zone and 30 reference populations of P. sylvestris and the taxa from the P. mugo complex from their European distribution were used in the study (Fig. 1, Table 1). Previous biometric studies suggested the existence of relict P. uncinata populations in the Silesian Lowlands. Those locations are geographically distant from the contemporary range of the species (Marcysiak and Boratyński 2007) but close to the contact zone. Therefore, we included P. uncinata in our analysis.

Fig 1
figure 1

Geographical location of the contact zone of pine species from Zieleniec reserve and reference allopatric populations of Pinus sylvestris and the taxa from the P. mugo complex. Distribution range of P. sylvestris provided by Euforgen network is marked in dark gray

Full size image
Table 1 Location of populations from Zieleniec reserve and the four reference pine species
Full size table

The contact zone comprised sympatric populations of the pine species at the Zieleniec reserve which is the largest peat bog complex in the Sudety Mountains, southwest Poland. For taxonomic classification of trees in the contact zone, we used several phenotypic traits, i.e., growth form, bark color of the upper part of stem and main branches, color and shape of needles, and setting angle of conelet from the previous year (when available) (Christensen 1987b). Based on these traits, we identified pure bushy P. mugo individuals (PMZ), monocormic (single-stemmed) P. sylvestris individuals (PSZ), mono- and oligocormic (one-three stemmed) P. uliginosa individuals (PUZ), and a group of oligo- and polycormic (multi-stemmed) P. uliginosa-like pines of a height of over 2 m and untypical morphology that could not be classified as either of parental species (HBZ). Needles from 40 to 60 trees representing each phenotypic group were sampled. The samples were collected over about 50 ha. Trees of different phenotypes were intermixed across the area and the distance between sampled trees was from 5 to 100 m.

The sampled trees were screened with a diagnostic chloroplast DNA marker (cpDNA) from the trnF–trnL region (Wachowiak et al. 2000). This PCR-RFLP marker (using DraI as a restriction enzyme) is paternally inherited and dispersed through pollen and seeds (Wachowiak et al. 2006a). It was previously developed based on a single-nucleotide polymorphism that leads to an undigested PCR product for P. sylvestris (S marker) and a digested PCR product for P. mugo (M marker) giving two bands after electrophoresis on agarose gel. Analysis of many populations from across the distribution range of the pine species in Europe confirmed that this marker discriminates P. sylvestris from taxa of the P. mugo complex (Wachowiak et al. 2000). PCR-RFLP analysis of the samples from Zieleniec reserve indicated the presence of the M marker in all P. mugo, P. uliginosa, and the group of taxonomically unclassified polycormic (multi-stemmed) P. uliginosa-like pines of untypical morphology (HBZ). However, the M marker was also found in some individuals that had been initially classified as P. sylvestris. Therefore, we distinguished another group of samples at Zieleniec reserve, namely monocormic P. sylvestris-like individuals carrying cpDNA diagnostic for the P. mugo complex (HPSZ) (Table 1).

A subset of 16 to 20 individuals from each group was selected for molecular analyses at 26 nuclear gene loci (Supplementary Table S1). We compared levels of nucleotide variation in the contact zone versus reference populations from allopatric zones of the species. The reference populations were treated either separately or were grouped by geographical location to compare a similar number of samples relative to the putative hybrid zone. Pinus sylvestris reference populations were grouped into Northern Europe (Finland and Sweden), Central Europe (Poland and Austria), Southern Europe (Spain), and Northwest Europe (Scotland; Table 1). Pinus mugo reference populations were grouped into Central Europe (Sudetes and Alps), Carpathians (Eastern and Southern), and Balkans (Pirin and Durmitor Mts.), and Italy (Carnic Alps and Abruzzi). Seeds from at least ten trees separated at a distance of a minimum 50 m were sampled from each reference population providing a set of 384 trees (Table 2). Using the QIAGEN DNeasy Plant Mini Kit, DNA was extracted from either haploid megagametophyte of seeds (samples from reference populations) or from needles (samples from the contact zone).

Table 2 Nucleotide and haplotype diversity of the groups of pines defined at Zieleniec reserve (Table 1) and reference populations of pure taxa of Pinus mugo, P. sylvestris, P. uliginosa, and P. uncinata (average values across 26 nuclear genes are reported)
Full size table

PCR amplification and sequencing

Twenty-six nuclear gene loci related to regulation of gene expression, metabolism, signal transduction, and transport were analyzed (Supplementary Table S1). PCR primers are given in Supplementary Table S1. PCRs were performed with Thermo MBS thermal cyclers and carried out in a total volume of 15 μl containing 15 ng of template DNA, 10 μM of each dNTP, 0.2 μM each of forward and reverse primers, 0.15 U Taq DNA polymerase, 1× BSA, 1.5 μM of MgCl2, and 1× PCR buffer (BioLabs). Standard amplification procedures were used with initial denaturation at 94 °C for 3 min followed by 35 cycles with 30 s denaturation at 94 °C, 30 s annealing at 60 °C and 1 min 30 s extension at 72 °C, and a final 5 min extension at 72 °C. PCR fragments were purified using Exonuclease I-Shrimp Alkaline Phosphatase enzymatic treatment. About 20 ng of PCR products were used as templates in 10 μl sequencing reactions with the BigDye Terminator v3.1. DNA Sequencing Kit (Applied Biosystems) performed by Genomed S.A. (Warsaw, Poland) on a 3130xl Genetic Analyzer (Applied Biosystems). CodonCode Aligner software ver. 3.7.1 (Codon Code Corporation, Dedham, MA, USA) was used for editing of the chromatograms, visual inspection of all polymorphic sites detected, and alignments. GenBank reference sequences from Pinus taeda were used as an outgroup in neutrality tests (Supplementary Table S1). The haplotype phase (combined multi-locus polymorphism) of samples amplified from needles was determined using PHASE software and reference haplotype sequences of the pure species derived from DNA sequencing of the haploid megagametophyte tissue. Haplotype sequence data are deposited in the NCBI repository (Supplementary Table S1).

Nucleotide variation

We looked at the patterns of nucleotide and haplotype polymorphism among groups of samples from the contact zone, and between the samples from the contact zone and reference pure species populations. Nucleotide diversity was measured as the average number of differences per site (π) between two sequences (Nei 1987). The number of haplotypes (N h) and haplotype diversity (H d) were computed for each gene using DnaSP v.5 (Librado and Rozas 2009). The number and frequency of unique and shared haplotypes in pairwise comparisons between species was calculated with Arlequin v.3.5 (Excoffier and Lischer 2010). Locus-by-locus estimates of net divergence between the groups of samples (Nei 1987) were determined using SITES 1.1 (Hey and Wakeley 1997). Genetic relationships between reference pure species populations and groups of samples from the contact zone were assessed based on pairwise net average genetic distance at all polymorphic sites detected at 26 loci using Unweighted Pair Group Method with Arithmetic Mean (UPGMA) as implemented in Mega 6.06 (Tamura et al. 2013). Standard errors of the interior branches were calculated using 1000 bootstrap procedure and the Maximum Composite Likelihood model as implemented in Mega.

Population structure and taxonomic relationships

Population structure was studied using Bayesian clustering methods and analysis of molecular variance. In addition, we assessed genetic differentiation between groups of samples from the contact zone and the reference populations of the species. Clustering of individuals and populations and estimates of admixture were explored using BAPS 6.0 (Corander and Tang 2007). Ten independent runs were conducted for each K (1–35) to estimate the number of clusters for all samples combined, for the three taxa from the P. mugo complex (P. mugo, P. uliginosa, and P. uncinata), P. sylvestris, and samples from the contact zone based on all polymorphic sites. The linear linkage model was used, the number of iterations used to estimate admixture coefficients was set to 100, the number of reference individuals was set to 100, and the number of iterations used to estimate admixture coefficients for the reference individuals was 10. Genetic differentiation in pairwise comparisons between populations was measured as Wright’s fixation index (Weir and Cockerham 1984) over all polymorphic sites and tested for significance by 1000 permutations of the samples between populations using Arlequin v.3.5 software (Excoffier and Lischer 2010). To further assess among-population differentiation, we used principal coordinate analysis (PCoA) based on the mean net genetic distance between populations at all polymorphic sites.

Tests for natural selection

We tested for the effects of selection on genes in parental species at both reference and contact zone populations. Tests were based on a comparison of allelic frequency spectra between different groups of populations, and departures from neutral expectations of polymorphisms versus divergence at the interspecific level. Deviations from the standard neutral model of evolution were assessed using the frequency spectrum test and coalescence-based approaches. Tajima’s D (Tajima 1989) was computed using the difference between two distinct estimates of the scaled mutation parameter theta for each locus, and statistical significance was evaluated by a comparison to a distribution generated by 1000 coalescent simulations using DnaSP 5.1 (Librado and Rozas 2009). Species may undergo different population histories that can influence the assumptions of standard neutrality tests. Therefore, deviations from neutrality were also tested using two compound neutrality tests that are robust to demographic processes, HEW and DHEW (Zeng et al. 2007). The HEW and DHEW tests are a compound of Fay and Wu’s H and Tajima’s D/Fay and Wu’s H with the Ewens-Watterson neutrality test, respectively (Zeng et al. 2007). Significance levels were determined by 10,000 coalescent simulations based on Watterson’s estimator of theta as implemented in the dh package (http://zeng-lab.group.shef.ac.uk/wordpress/?page_id=28). The distribution of the test statistic was investigated for each locus in all samples combined from each species. Genetic differentiation in pairwise comparisons within and between species was studied locus by locus, and significance thresholds of the F ST values were set at 99 % and estimated by 1000 permutations of the samples between populations, regional groups, and species using Arlequin v.3.5. The false discovery rate (FDR) adjustment for multiple testing (lambda = 0.15, FDR level = 0.01) was conducted using QVALUE software based on the distribution of P values for the set of F ST statistics (Storey and Tibshirani 2003).

Results

Nucleotide and haplotype diversity

Across the 26 genes, ~10 kbp were aligned providing a set of 579 polymorphic sites (Supplementary material, Appendix file). Average nucleotide diversity ranged between πtot = 0.0036 and πtot = 0.0055, and it was slightly higher in samples from the contact zone versus pure species populations (Table 2, Supplementary Table S2). Out of 789 haplotypes detected across 26 loci in all 384 samples analyzed, the majority (57 %) were unique (present only once). Exclusive haplotypes were found in all species and the contact zone groups. The within-population percentage of unique haplotypes ranged from 2.3 % (P. uliginosa-like individuals from the contact zone) to 0.2 % (P. uliginosa from Batorów reserve) (Supplementary Table S3). The average haplotype diversity was similar for each species (H d = 0.59–0.66), regional groups of populations (H d = 0.61–0.70), and groups of samples from the contact zone (H d = 0.54–0.71) (Table 2). In the contact zone, no net divergence was observed between P. mugo, P. uliginosa, and polycormic P. uliginosa-like pines (Fig. 2, Supplementary Table S4 and S5). Divergence of these three groups was lower from the reference P. uliginosa (0.004–0.013) and P. mugo populations (0.005–0.013) and slightly higher from P. uncinata (0.009–0.023) (Fig. 2). Pinus sylvestris from Zieleniec and the reference populations of the species showed low divergence (0.001–0.031), and they were grouped together with P. sylvestris-like trees with M cpDNA haplotype (Supplementary Table S4 and S5, Fig. 2). The group of P. sylvestris-like trees from Zieleniec showed a similar level of divergence (0.010–0.015) from P. mugo, P. sylvestris, and P. uliginosa from that area (Supplementary Table S4).

Fig 2
figure 2

Unweighted pair group method with arithmetic mean (UPGMA) tree based on average net genetic distances between groups of samples from a contact zone at Zieleniec reserve (asterisks) and reference populations of the species including Pinus mugo (filled square), P. sylvestris (filled circle), P. uliginosa (filled diamond), and P. uncinata (filled triangle) (see Table 1 for details). Numbers indicate branch length. Genetic distance and its standard errors for each pairwise comparisons are shown in Supplementary Table S4

Full size image

Population structure

Of the reference populations, P. mugo and two of the P. uliginosa populations formed one group in the cluster analysis, while P. sylvestris and P. uncinata formed separate groups (total K = 3; Fig. 3a). The P. uliginosa population from Węgliniec showed greater similarity to P. uncinata. Two P. sylvestris individuals from Central and Northern Europe, two P. mugo, and several P. uliginosa and P. uncinata samples were identified as potentially admixed (Fig. 3a).

Fig 3
figure 3

a Bar plot from Bayesian assignment (BAPS, Corander and Tang 2007) of the group of pines from Zieleniec reserve and the reference samples indicating three genetic clusters (K = 3) corresponding to pure Pinus mugo, P. sylvestris, and P. uncinata species. Different colors represent proportional assignment of each sample to individual clusters, the black vertical lines separate the corresponding populations, and numbers above vertical bars refer to populations given in Table 1. Evidence of admixture (P < 0.01) was found across the samples from the sympatric population of the species from Zieleniec reserve. At the intraspecific level, heterogeneous population structure was observed between P. uliginosa populations. b Assignment of samples according to trnF–trnL cpDNA marker diagnostic to the P. mugo complex (M marker—green) and P. sylvestris (S marker—red)

Full size image

In the contact zone at Zieleniec, three clusters were identified: one included P. mugo (PMZ), P. uliginosa (PUZ), and polycormic P. uliginosa-like pines (HBZ) from that area; the second contained P. sylvestris (PSZ); and the third contained P. sylvestris-like samples (HPSZ) that were admixed with P. mugo/P. uliginosa and P. uncinata (Fig. 3a). There was evidence of admixture in a group containing P. mugo, P. uliginosa, and polycormic P. uliginosa-like pines (HBZ) that appears to result from hybridization with P. sylvestris and P. uncinata. Only M and S haplotypes of the trnF–trnL cpDNA region were observed across all samples analyzed (Fig. 3b). Clear division between species was found in the PCoA analysis based on the genetic distance between populations (Fig. 4). Groups of P. mugo pines from Zieleniec reserve (PMZ and PUZ) and oligo- and polycormic P. uliginosa-like pines (HBZ) were placed between P. mugo and P. uliginosa populations. Pinus sylvestris-like hybrids (HPSZ) were placed between P. mugo/P. uliginosa groups and P. sylvestris. Pinus sylvestris from the contact zone clustered with allopatric populations of the species (Fig. 4).

Fig 4
figure 4

Plot of the first two axes of a principal coordinates analysis (PCoA) based on a genetic distance matrix for groups of samples from Zieleniec reserve and reference populations of the Pinus mugo, P. sylvestris, P. uliginosa, and P. uncinata at all polymorphic sites from 26 nuclear genes. Population acronyms as in Table 1

Full size image

Significant differentiation between each pair of taxa was found across all polymorphic sites with P. mugo and P. sylvestris from Zieleniec reserve as the most diverged (F ST = 0.40, P < 0.01) (Table 3). Populations in regional groups of P. mugo and P. sylvestris showed low but significant differentiation across all polymorphic sites except P. sylvestris populations from Central and Northern Europe. No differentiation was found between P. mugo, P. uliginosa, and the group of polycormic pines from Zieleniec reserve (Table 3). The other groups showed similar levels of differentiation as between populations of allopatric zones of P. mugo and P. sylvestris.

Table 3 Pairwise F ST between groups of Pinus samples from Zieleniec reserve and the reference populations of pure taxa at all polymorphic sites (significant values at P <0.01 are marked in bold)
Full size table

Neutrality tests

An excess of singleton mutations across genes was evident as multilocus Tajima’s D was negative in most groups (D = −0.859 to −0.039) except P. sylvestris-like samples from Zieleniec reserve and reference P. uncinata (Table 2). Significant deviations from neutrality at Tajima’s D and/or HEW and DHEW compound neutrality tests were found at eight loci in P. sylvestris and P. mugo, seven at P. uliginosa and at one locus in P. uncinata (Supplementary Table S6). In the samples from Zieleniec reserve, significant excess of singleton mutations was found at locus Pr4_5 in a group of oligo- and polycormic P. uliginosa-like pines, at Pr4_17 in P. sylvestris-like individuals, at Pr1_28, Pr2_17, Pr4_17 in P. sylvestris, and Pr2_28 in P. mugo. An excess of intermediate-frequency variants was found at Pr2_47 in P. sylvestris and Pr4_19 in polycormic pines from Zieleniec reserve (Supplementary Table S6).

Locus-specific genetic differentiation and selection

The loci studied showed high divergence between species and very low differentiation between populations within species (Supplementary Table S7). No significant differentiation at any locus was observed between samples from Zieleniec reserve defined as P. mugo/P. uliginosa and a group of taxonomically unclassified P. uliginosa-like polycormic pines. These three groups of samples were highly diverged from P. sylvestris and P. sylvestris-like samples at most of the loci (Supplementary Table S8). Divergence of populations from Zieleniec from the reference populations (Supplementary Table S9) was higher than divergence between pure species populations (Supplementary Table S7). At the SNP level, several loci were significantly differentiated between P. sylvestris and the P. mugo complex: 46 SNPs were found in 22 genes from all categories between P. sylvestris and P. mugo, 40 SNPs in 16 genes between P. sylvestris versus P. uliginosa, and 24 SNPs in 11 genes between P. sylvestris versus P. uncinata. Within the P. mugo complex, P. mugo versus P. uliginosa and P. uncinata were highly diverged at 5 and 12 SNPs, respectively, and P. uliginosa and P. uncinata at 6 SNPs (Supplementary Table S10). These polymorphic sites and genes form a valuable set of new markers for tracking hybrid genotypes. The most diverged genes were those involved in regulation of gene expression, signal transduction, transport, and cellular metabolism and appear to be of high importance in speciation and local adaptation (Supplementary Table S1 and S10).

Discussion

Here, we used nucleotide polymorphism data to compare patterns of genetic variation in pine trees from the contact zone with those from allopatric zones of four species from across Europe. We focused on the contact zone between species at a peat bog complex located in the Sudety Mts. in Poland. The aim was to evaluate the role of introgressive hybridization and selection on nucleotide diversity in the analyzed populations. To date, such studies have been limited by the lack of diagnostic biometric and biochemical characters for tracking interspecific gene flow and the identification of individual hybrid trees. For instance, many anatomical traits of needles show overlapping frequency distributions among taxa (Boratyńska et al. 2015). In our study, we used polymorphisms at nuclear genes and a cpDNA marker to identify interspecific admixture. Nuclear markers that experience high rates of intraspecific gene flow are especially relevant for species delimitation (Petit and Excoffier 2009). Our data clearly indicate that the markers applied in our study can accurately discriminate pure parental species. These markers provide a large resource of SNP information for future use in tracking interspecific gene flow and evaluating species composition in other contact zones where individuals with mixed morpho-anatomical characteristics have been described (e.g., in the Alps (Christensen 1987a), Slovakia (Kormutak et al. 2014)).

Hybrids can exhibit intermediate trait values, combine traits from both parents, and/or exhibit extreme trait values as compared to the parental species (Gross and Rieseberg 2005). In our study, high genetic similarity was observed between samples classified as P. mugo, P. uliginosa, and a group of polycormic P. uliginosa-like pines (nucleotide diversity πtot = ~0.0044–0.0049, zero net divergence, no significant differentiation in allele frequency spectra). The three groups formed a uniform genetic cluster and showed low divergence (about 1 %) from pure-species populations of the taxa from the P. mugo complex. In contrast, they differed clearly from P. sylvestris from Zieleniec reserve and from reference populations of the species. Those oligo- and polycormic P. uliginosa-like pines that had cpDNA diagnostic to the P. mugo complex should be considered as hybrids between P. mugo and P. uliginosa with no evidence of a substantial contribution from P. sylvestris. In addition, we identified a group of P. sylvestris-like trees that had cpDNA markers characteristic of the P. mugo complex. This group of monocormic pines showed clear evidence of admixture between P. mugo/P. uliginosa, P. sylvestris, and also P. uncinata. These trees represent a second group of phenotypically distinct hybrids with P. sylvestris acting as a mother tree at an early stage of hybridization, and cpDNA delivered by pollination from the polycormic pines of the P. mugo complex. All hybrids shared some of the nuclear gene haplotypes observed in the reference populations of parental species. Different levels of introgression of those individuals, as evident from variable assignment to the parental species, indicate that this group of pines includes F1 and subsequent generations of hybrids, although always fertilized by pollen carrying the cpDNA diagnostic for the P. mugo complex. This asymmetric introgression between taxa may be the result of selection, but different mechanisms of incompatibility between hybrids cannot be excluded (Currat et al. 2008).

Pinus uliginosa had a heterogeneous genetic background, with one population from Węgliniec reserve showing close genetic similarity to P. uncinata. This population showed patterns of admixture not observed in P. uliginosa from the species locus classicus at Batorów reserve (Neumann 1837) and the population from Mittelwalde. This suggests that some individuals classified as P. uliginosa from Węgliniec but also Zieleniec reserve may represent remnants of marginal populations of P. uncinata. Interestingly, the pattern of admixture observed in our study is in line with some biochemical and biometric data on cone and needle morphology (Marcysiak and Boratyński 2007 and references therein). These studies showed that P. uliginosa from Węgliniec shares some biometric traits with P. sylvestris, P. mugo, and P. uncinata and is distinct as compared to other allopatric stands of the species (Lewandowski et al. 2000; Prus-Głowacki et al. 1998). Therefore, the existence of remnant populations of P. uncinata in the Silesian Lowlands seems possible and may result from the expansion of west-European refugia across the northwestern pre-Alpine territories during the late Dryas/early Holocene (Marcysiak and Boratyński 2007). Alternatively, considering the large number of shared ancestral alleles still segregating in the species from the P. mugo complex (Wachowiak et al. 2013), interspecific gene flow could create combinations of alleles in individual hybrid trees clustered in our analysis as different taxonomic units. These results provide another dimension to the very complex demographic history of the taxa within the P. mugo complex

Our data contribute to the assessment of the genetic relationships of the taxa from the P. mugo complex showing evidence of close genetic identity of P. uliginosa as compared to P. mugo and P. uncinata (net divergence of 0.8–0.9 %, respectively). Pinus uliginosa does not harbor a distinct gene pool as compared to P. mugo (Fig. 3a), and therefore it should not be considered as a separate species but rather subspecies within the P. mugo complex. The taxon was originally described from the peat bog population at Batorów reserve in the Sudety Mts. (Neumann 1837). In a biometric revision of the complex, it was proposed as a subspecies of P. uncinata (Businsky and Kirschner 2006) or as a synonym of Pinus rotundata—a pine species forming high-altitude populations in the Alps and Northern Carpathians and small isolated populations in the Pyrenees and Massif Central (Christensen 1987b). Our molecular analysis did not aim to resolve the taxonomic position of P. uliginosa. However, the data showed closer genetic similarity of P. uliginosa to P. mugo than P. uncinata. It is clear that more molecular studies are needed to clarify the taxonomic status of several taxa described within the P. mugo complex and to evaluate the role of interspecific gene flow with P. sylvestris (Christensen 1987b).

Barriers to interspecific hybridization and a lack of evidence for bidirectional gene flow between P. sylvestris and P. mugo were suggested in previous studies that found hybrid seeds derived only from P. sylvestris-like individuals pollinated with P. mugo but not from reciprocal crossings (Wachowiak et al. 2005a). A lack of hybrids from crossings between P. mugo as a maternal and P. sylvestris as a paternal tree, but putative hybrid individuals from reverse crossing combinations (with P. mugo as a pollen donor), were found based on a joint analysis of cpDNA, isozymes, and phenotypic characteristics of trees (Wachowiak and Prus-Głowacki 2008) and at nuclear genes in a P. sylvestris and P. mugo population (Kormutak et al. 2014). Cryptic hybrids between P. sylvestris and P. uncinata were found in the sympatric populations of the species (Jasińska et al. 2010). So far, the only evidence of reciprocal hybridization was found in a sympatric population of P. sylvestris and P. uliginosa (Wachowiak et al. 2005b).

Our results suggest that hybrids express distinct phenotypic variability as compared to parental species. In previous biometric and biochemical studies, the variety of morphological forms observed in sympatric populations of P. sylvestris and the P. mugo complex was explained as either the result of intensive hybridization and introgression that changed the population into a hybrid swarm or as a mixture of mostly pure pine species from the P. mugo complex and P. sylvestris (Bobowicz 1990; Odrzykoski 2002; Wachowiak et al. 2006b). Our data indicate that hybridization takes places in contact zones of the species and leads to propagation of viable hybrid trees. They exist together with pure parental species and maintain their phenotypic identity.

Natural selection can cause the fixation of advantageous alleles (or chromosomal segments) in ecologically diverged hybrids (Beaumont and Balding 2004; Buerkle and Lexer 2008; Lexer et al. 2003). Introgressed alleles may often have a positive fitness effect in their new genetic background and traits responsible for adaptation can be transferred between species (Martin et al. 2006). Foreign alleles in different genetic or ecological backgrounds will show a range of fitness outcomes, but only those that increase the adaptive optimum in a given environment will effectively introgress. Consequently, introgression of alleles derived from other species has the potential to speed adaptation (Gompert and Buerkle 2010), which may be particularly influential in populations undergoing spatial or temporal transitions into new environments. Our study provides evidence of successful hybridization within the sympatric study population but no evidence for interspecific gene flow outside the contact zone. It is possible that hybrids have reduced fitness in environmental conditions occupied by parental species, and they may persist best in new habitats. Indeed, our results show that hybrid genotypes have succeeded in peat bogs close to mountain regions, which are environments untypical of either parental species.

Co-existence of morphologically variable taxa and hybrids together with asymmetric gene flow indicates the role of selection in maintaining certain phenotypes. Strong directional selection on loci underlying fitness-related adaptation in the ecologically diverged hybrids should increase the frequency of advantageous alleles. We expect, however, that the effect of selection should be localized in the genome and the genetic background of a species should not be affected by the spatial expansion of an advantageous allele (Currat et al. 2008). Two distinct groups of hybrids seem to have maintained their phenotypic differentiation, and we found signatures of selection at some loci as compared to background variation. Pinus sylvestris-like hybrids showed increased frequency of alleles specific to P. sylvestris and alleles specific to P. mugo at different genes. Similarly, at one polyol transporter gene, oligo- and polycormic hybrids showed differentiation from P. mugo but not from P. uliginosa. An increase in frequency of alleles unique to one of the parental species at some loci may result from selection of particular alleles in the hybrids that increase their fitness in the peat bog environment. In contrast, no evidence of differentiation at any locus was found between P. mugo and P. uliginosa, and between P. uliginosa and a group of polycormic hybrids from Zieleniec reserve. In the case of parental species, strong directional selection at some loci due to local adaptation in ecologically diverged peat bog environments should increase differentiation between the peat bog and the reference pure-species populations. In the absence of population structure in the parental species, significant differences in allele frequency spectra and/or departures from neutrality between reference and contact zone populations were found at nine loci in P. mugo and eight in P. sylvestris. These loci showed no evidence of differentiation between pure-species populations. Assuming different patterns of diversity at selectively influenced loci relative to background genetic variation, this increased population differentiation suggests selection in response to specific peat bog environments not optimal for either of the parental species (e.g., Eveno et al. 2008; Kujala and Savolainen 2012; Wachowiak et al. 2009).

In many cases, as shown in theoretical models and experimental studies of contact zones, introgression is highly asymmetric (Currat et al. 2008; Petit et al. 2004) and may go from the local to the invading species. If this scenario holds for the investigated pine taxa, then we could consider P. sylvestris and possibly P. uliginosa as an invading taxon as compared to the local P. mugo. At present, however, this contact zone is isolated from the continuous range of any of the parental taxa. In the case of invading species, the pattern of introgression at some neutral loci resulting from the range expansion of a species into an already occupied territory may mimic the effect of selection. Therefore, more studies of the contact zone dynamics (e.g., Cinget et al. 2015) at Zieleniec reserve and its neighborhood are needed. Such studies would help to test invasion models and evaluate the role of demographic processes on the patterns of genome-wide nucleotide sequence variation.

With the presence of different taxa and hybrid groups in environmental conditions not optimal for either of the putative parental species, the contact zone at Zieleniec reserve is a relevant biological system for studying the role of hybridization on adaptation to new environments at the genetic level. This approach has recently been developed in several model plant species demonstrating the role of hybridization and adaptive introgression in the evolution of irises (Iris; Arnold et al. 2004), ecological divergence of sunflowers (Helianthus; Rieseberg et al. 2007), and the signatures of divergent and balancing selection in campions (Silene; Minder and Widmer 2008) and poplar (Populus; Lexer et al. 2010). So far, a few genes involved in adaptation or speciation have been identified in plants including hybrid sterility loci (Lexer and Widmer 2008), determinants of flower color-linked pollinator shifts (Hoballah et al. 2007), and genes involved in hybrid necrosis (Bomblies and Weigel 2007). Our study shows that the investigated taxa maintain genetic and phenotypic differentiation in the presence of extensive gene flow. Considering the abundance of trees growing on peat bog in our focal populations, both P. sylvestris-like hybrids and oligo- and polycormic P. uliginosa-like pines could serve as suitable mapping populations in the search for loci underlying local adaptation and genetic and phenotypic differentiation between taxa.