Genetic Consequences of Hybridization in Relict Isolated Trees Pinus sylvestris and the Pinus mugo Complex

: Scots pine ( Pinus sylvestris L.) and the taxa from the P. mugo complex can hybridize in the contact zones and produce fertile hybrids. A unique example of an early Holocene relict population of P. sylvestris and P. uliginosa (a taxon from the P. mugo complex) growing on the tops of Jurassic sandstone rocks is located in Bł˛edne Skały (Sudetes). Phenotypically, there are trees resembling P. sylvestris , P. uliginosa and intermediate forms between them. We expected that some of P. sylvestris and / or P. uliginosa -like trees could be in fact cryptic hybrids resembling one of the parental phenotypes. To address this question, we examined randomly sampled individuals, using a set of plastid (cpDNA), nuclear (nDNA) and mitochondrial (mtDNA) markers as well as biometric characteristics of needles and cones. The results were compared to the same measurements of allopatric reference populations of the P. sylvestris and the P. mugo complex ( Pinus mugo s.s , P. uncinata and P. uliginosa ). We detected cpDNA barcodes of the P. mugo complex in most individuals with the P. sylvestris phenotype, while we did not detect cpDNA diagnostic of P. sylvestris within P. uliginosa -like trees. These results indicate the presence of cryptic hybrids of the P. sylvestris phenotype. We found only three typical P. sylvestris individuals that were clustered with the species reference populations based on needle and cone characteristics. Most trees showed intermediate characteristics between P. sylvestris and P. uliginosa -like trees, indicating intensive and probably long-lasting hybridization of the taxa at this area and subsequent gene erosion of parental species.


Introduction
Hybridization and introgression are important evolutionary factors that increase species variation and may lead to speciation [1][2][3][4][5][6]. These two processes concern mostly related species, which do not express isolation mechanisms such as phenological barriers [7] and post-pollination barriers to fertilization [8]. The interbreeding between different species of Pinaceae was found to be quite as compared to some earlier works (iii) to what extent the level of genetic variation of this small and isolated population is comparable to variation in much bigger reference populations of pure species. To address those questions, we sampled individuals across the phenotypic forms present in this population and performed genetic and biometric analyses to verify taxonomic relationships between the samples and to compare morphology of trees from the mixed population to the reference populations.

Material
The mixed pine population in Błędne Skały (BS) occurs on the upper parts of sandstone rocks and is located at the summit of the north-western part of the Skalniak ridge in the Stołowe Mountains, at an elevation of 850 m a.s.l. It consists of dispersed pine individuals of P. sylvestris and P. uliginosa phenotype, with a few individuals resembling P. mugo s.s. and a number of individuals of intermediate phenotype. The pine population is surrounded by extensive Picea abies forests, which isolate it from direct contact with other pine populations.
Material for the study was taken randomly from 64 individuals in an area of about 2 ha. The sampled individuals were about 20 m distance from each other or growing on separate rocks. During sampling, part of the individuals was preliminary identified in the field based on growth form, bark color on the upper part of stem, needle color and shape and setting angle of the one-year-old conelet [15,16]. The sampling included 15 individuals of P. sylvestris-like phenotype (BS_S), 17 individuals resembling the taxa from the P. mugo complex namely P. uliginosa (BS_M), 7 individuals of intermediate character between P. sylvestris and P. uliginosa classified as hybrids (BS_H) and another 25 unclassified individuals (BS_N). As the sampled trees could potentially represent a genetic background of different pine species, we included a set of reference populations taken from across current natural ranges of P. sylvestris and the taxa of the P. mugo complex (P. mugo s.s., P. uliginosa and P. uncinata Ramond (Table 1; Figure 1). These populations have been used in previous morphometric [54][55][56][57][58][59][60][61] and genetic [30,[62][63][64][65] studies. Particular analyses were conducted using all taxa, but not always the same reference populations.
Forests 2020, 11, x FOR PEER REVIEW 3 of 18 population and performed genetic and biometric analyses to verify taxonomic relationships between the samples and to compare morphology of trees from the mixed population to the reference populations.

Material
The mixed pine population in Błędne Skały (BS) occurs on the upper parts of sandstone rocks and is located at the summit of the north-western part of the Skalniak ridge in the Stołowe Mountains, at an elevation of 850 m a.s.l. It consists of dispersed pine individuals of P. sylvestris and P. uliginosa phenotype, with a few individuals resembling P. mugo s.s. and a number of individuals of intermediate phenotype. The pine population is surrounded by extensive Picea abies forests, which isolate it from direct contact with other pine populations.  Table 1.
Material for the study was taken randomly from 64 individuals in an area of about 2 ha. The sampled individuals were about 20 m distance from each other or growing on separate rocks. During sampling, part of the individuals was preliminary identified in the field based on growth  Table 1.

DNA Extraction, Amplification and Marker Genotyping
Genomic DNA was isolated from dry needles, using a DNeasy Plant Mini Kit (Qiagen, Hilden, Germany). DNA quality was evaluated using a BioPhotometer (Eppendorf, Hamburg, Germany). The samples were analyzed with a set of organellar and nuclear markers that show different modes of inheritance. Polymorphisms at 14 mitochondrial regions, inherited in pines in the maternal line and distributed at short geographical distances by seeds, were investigated including 13 markers [66], and the insertion-deletion (IN-DEL) locus nad1 intron B/C [67]. The standardized 15 µL PCR reaction mix contained 1 × PCR buffer, 1X bovine serum albumin (BSA), 1.  elongation at 72 • C, with the final extension for 5 min at 72 • C. Variants were genotyped by polymerase chain reaction restriction fragment length polymorphism (PCR-RFLP) and digested in 10 µL reaction volume, containing 3 µL of the PCR product, 2X buffer and 2.5 U of the relevant restriction enzyme [68]. Samples were incubated for at least 10 h, following the manufacturer´s recommendations for respective enzymes (Thermo Fisher Scientific, Waltham, MA, USA), and electrophoretically separated on 1.5%-2% agarose gel, stained with GelRed in UV. The insertion-deletion polymorphisms in intron B/C of the nad1 locus was amplified and scored according to the methods described in Soranzo et al. [67]. The region contains shorter variant a that appears to be mostly fixed in P. mugo s.s. populations and it segregates with different frequency across the geographical distribution of P. sylvestris. The polymorphisms at nad1 B/C region was scored by 2% agarose gel electrophoresis of the corresponding PCR-amplified fragment.
Variation at the plastid genome, which is paternally inherited in these taxa and transmitted by pollen, was studied at a set of 12 microsatellite markers. Two multiplex reactions of the chloroplast simple sequence repeats (cpSSR) loci were run, including regions PCP1289, PCP26106, PCP30277, PCP36567, PCP41131, PCP45071, PCP87314, PCP102652 (multiplex 1) and Pt15169, Pt26081, Pt30204, Pt71936 (multiplex 2) [69]. Additionally, a species-diagnostic plastid (cpDNA) intergenic trnF-trnL region that discriminates P. sylvestris from the taxa of the P. mugo complex was used to determine the origin of the plastid genome in samples from the focal population in Błędne Skały. This PCR-RFLP marker has previously been developed based on a single-nucleotide polymorphism that leads to an undigested PCR product for P. sylvestris and a digested PCR product for P. mugo when treated with DraI enzyme [70]. The PCR was run in a total volume of 15 µL, containing 15 ng of template DNA, 10 µM of each of dNTP, 0.2 µM each of forward and reverse primer, 0.15 U Taq DNA polymerase, 1X BSA, 1.5 µM of MgCl 2 and 1X PCR buffer (Novazym, Poznań, Poland). Standard amplification conditions were applied, including 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, with a final 5 min extension at 72 • C.

Diversity and Differentiation Measures
Alleles identified within each mitochondrial (mtDNA) marker were concatenated, and connections among the obtained haplotypes were analyzed by a Median Joining network constructed in PopART 1.7 [72,73]. The number of segregating sites (S), number of private variants (S P ) and the mean within group distance (d) were calculated in DnaSP V6 [74] and MEGA 7 [75] for groups of samples defined as BS and for the reference populations. Similarly, the numbers of haplotypes (N h ) and private haplotypes (N p ) present in populations and species were calculated. To better assess the scale of genetic variation at mtDNA segregating in the taxa, the haplotype frequencies and population gene diversities (H d ) were estimated in DnaSP. The relationships among samples were evaluated based on the mean distances between populations (d xy ), using the Unweighted Pair Group Method with Arithmetic Mean (UPGMA) as implemented in MEGA. The distribution of genetic variation among populations, and among genetic structures detected in the population groups were checked by F st and G st calculated in DnaSP. The differences between populations were verified by hierarchical analysis of molecular variance (AMOVA), using Arlequin 3.5.22 [76].
To describe the variation of cpDNA SSR loci, for each defined sample group we calculated mean number of haplotypes (A h ), mean effective number of haplotypes (A e ), number of private haplotypes (A P ), haplotype richness (A R ), haplotype diversity (H d ) and mean genetic distance of individuals within populations (D 2 sh) using GenAlEx 6 [77]. Variation at the diagnostic trnF-trnL cpDNA region was scored as a presence or absence of a DraI restriction site.
Genetic variation at nSSR of the mixed population and the reference populations was analyzed based on mean number of alleles per locus (N a ), mean effective number of alleles per locus (Ne), number of private alleles (A p ), mean observed heterozygosity (H o ), mean expected heterozygosity (H e ), mean fixation index (F). The analysis of molecular variance (AMOVA) on the R st distances, as implemented in Arlequin 3.5.22, was used to evaluate the differentiation among defined groups of populations from BS contact zone and reference populations at nSSR loci. Afterward, the principal coordinates analysis (PCoA) on Nei's genetic distance was applied to visualize the genetic structure of the populations and BS group. The Bayesian clustering method implemented in STRUCTURE 2.3.4 [78][79][80] was used to assign individuals and populations into genetically distinct groups without a priori assumption of their origin. The STRUCTURE software uses the algorithm of Markov chain Monte Carlo (MCMC) to assign particular individuals to a number of genetic clusters (K) whose members share similar patterns of variation. It does not consider samples' origins but assumes that each cluster is in optimal Hardy-Weinberg and linkage equilibria. The parameter sets, that assumed correlated allele frequencies and the used admixture model, together allowed attribution of mixed recent ancestry of individuals and assignment of the proportion of inferred clusters' origin within each individual genome. The "recessive alleles" option was used, with 20 runs for each K, from K = 1 to 10, and with burn-in lengths of 500,000 and 100,000 iterations. The data probability distributions (Ln P [D]) and the ∆K values [81] were calculated using STRUCTURE HARVESTER Web v0.6.94 application and were used to determine optimal K for given dataset [82].

Biometrics
We characterized two-year-old needles using a set of traits (Table S5), following the procedures described by Boratyńska et al. [56] and Boratyńska and Boratyński [23]. Every individual was characterized using the average measurements of five needles, each from a different dwarf shoot (i.e., the morphological structure in Conifers from which needles in fascicles of 2-5 grow together). The anatomical characteristics were measured and/or evaluated on the preparations from the central part of the needle [21,23,55,56,83].
Cones were characterized using 11 measured characteristics and five proportions (Table S6) describing cone and cone scale shape [17,54,84]. Measurement procedures were adopted from Marcysiak [54] and Marcysiak and Boratyński [57]. Length, width, diameter and circumference measurements were taken from closed cones after soaking, while the number of cone scales was determined for open cones after drying. The scale characteristics were measured on one scale taken from the convex side of each cone at the maximum diameter. All measurements were done manually, using an electronic caliper.

Statistical Analyses of Biometric Data
Prior to multivariate comparisons, the data distribution and homoscedasticity of the variances were inspected, applying the Shapiro-Wilk's and Brown-Forsythe tests, respectively [85,86]. The arcsine transformation of percentage data was applied prior to statistical analyses. The data were standardized using the STATISTICA 9 software (StatSoftPolska, Kraków, Poland) to avoid a possible influence of their different types on the results of the multivariate analyses.
We determined the minimum and maximum values of the characteristics, arithmetic means, and variation coefficients. These parameters were calculated and analyzed for every individual in the BS population, and for reference populations of P. sylvestris, P. mugo, P. uliginosa and P. uncinata. Principal component analysis (PCA) was used to detect grouping of samples from the BS mixed population and their position among referenced populations of P. sylvestris, P. mugo and P. uncinata. Subsequently, discrimination analysis was conducted on these characteristics, except for the ones which were strongly biased (SRC, SRF, VB and RC -see Tables S5 and S6 for traits description) [86,87]. To assess the significance of differences for specific measured traits between groups preliminary determined in the BS population during sampling, Tukey's HSD test for traits with normal distribution and the Kruskal-Wallis test for traits with biased distribution were applied [86,88], using STATISTICA. The taxonomic position of each individual from the BS mixed population, assessed via macro-morphological traits during sampling, was verified using characteristics of the needles and cones. Finally, the results from the biometric study were compared with the results of the genetic analyses to distinguish individuals of P. sylvestris from the taxa from the P. mugo complex and their hybrids.

Genetic Diversity and Differentiation
In total, 39 mitotypes were identified across the four pine species (Figure 2, Table S1). Species-specific combination of haplotypes was observed only in P. uncinata, which was distinct when compared to the other species (Figure 2, Figure S1). Different frequencies of mitotypes were observed in individual reference populations. The highest frequency of shared haplotypes was found in individuals from the mixed BS population (Table S2). In the BS pine mixed population, we observed mitotypes of P. sylvestris, P. mugo s.s. and P. uliginosa reference populations, but not of P. uncinata. Unique mtDNA haplotypes of pines from BS population were found for five individuals from hybrid or unclassified sample groups. The remaining pines shared numerous mitotypes present in individuals from BS and reference populations. The most frequent mitotypes at the BS population were found at high frequency in the P. uliginosa population from Węgliniec Nature Reserve (Table S1). Some unique mitotypes were found at low frequency in every examined species and BS population. The samples including four reference species and a BS population formed each separate group in the UPGMA tree based on Nei's genetic distance at mtDNA regions ( Figure S2). In total, 39 mitotypes were identified across the four pine species (Figure 2, Table S1). Species-specific combination of haplotypes was observed only in P. uncinata, which was distinct when compared to the other species (Figure 2, Figure S1). Different frequencies of mitotypes were observed in individual reference populations. The highest frequency of shared haplotypes was found in individuals from the mixed BS population (Table S2). In the BS pine mixed population, we observed mitotypes of P. sylvestris, P. mugo s.s. and P. uliginosa reference populations, but not of P. uncinata. Unique mtDNA haplotypes of pines from BS population were found for five individuals from hybrid or unclassified sample groups. The remaining pines shared numerous mitotypes present in individuals from BS and reference populations. The most frequent mitotypes at the BS population were found at high frequency in the P. uliginosa population from Węgliniec Nature Reserve (Table S1). Some unique mitotypes were found at low frequency in every examined species and BS population. The samples including four reference species and a BS population formed each separate group in the UPGMA tree based on Nei's genetic distance at mtDNA regions ( Figure S2).  Table 1).
Only three individuals from the BS population had diagnostic a cpDNA marker characteristic for P. sylvestris. Those three individuals were indicated as pure P. sylvestris based on STRUCTURE analysis of cpSSR loci. The remaining individuals from BS had plastid DNA specific to the P. mugo complex (i.e., P. uliginosa), independently of the phenotype.   Table 1).
Only three individuals from the BS population had diagnostic a cpDNA marker characteristic for P. sylvestris. Those three individuals were indicated as pure P. sylvestris based on STRUCTURE analysis of cpSSR loci. The remaining individuals from BS had plastid DNA specific to the P. mugo complex (i.e., P. uliginosa), independently of the phenotype.
Individuals from BS showed reduced genetic variation at cpSSRs, with a considerably lower number of haplotypes and lower haplotype diversity as compared to the reference species (Table 2A) and individual populations (Table S3). The individuals unclassified taxonomically and those defined as P. sylvestris from the BS mixed pine population showed the highest mean genetic distance between individuals within populations (D 2 sh range 32.1-51.9) as compared to the other groups in this area (7.3-9.1) and the reference populations (4.4-13.6). Haplotype and gene diversity at mtDNA markers were also lowest in the BS population (Table 2B). However, in contrast to organelle markers, variation at biparentally inherited nSSR loci was similar in the BS population as in the reference samples (Table 2C,  Table S4). There was no evidence of reduced heterozygosity among samples in the BS population (Table 2C). According to the PCoA based on genetic distance at nSSR loci, all defined groups of individuals from BS showed close genetic similarity to each other, to P. uliginosa from the Węgliniec Nature Reserve and to P. uncinata samples (Figure 3). This group of samples was distinct from the reference populations of P. mugo s.s. and P. sylvestris. Among two outlier populations detected appeared one P. uliginosa from the Large Batorów Peatland, with high genetic similarity to group of three P. mugo s.s. populations and the fourth P. mugo stand localized in Sudetes (M8), which was distinct as compared to all other samples analyzed (Figure 3). from BS showed close genetic similarity to each other, to P. uliginosa from the Węgliniec Nature Reserve and to P. uncinata samples (Figure 3). This group of samples was distinct from the reference populations of P. mugo s.s. and P. sylvestris. Among two outlier populations detected appeared one P. uliginosa from the Large Batorów Peatland, with high genetic similarity to group of three P. mugo s.s. populations and the fourth P. mugo stand localized in Sudetes (M8), which was distinct as compared to all other samples analyzed (Figure 3).  Table 1).
Structure analysis at nSSR loci indicated the presence of five genetically distinct groups corresponding approximately to the four reference species and one sample group from BS (Figure 4,  Table 1).
Structure analysis at nSSR loci indicated the presence of five genetically distinct groups corresponding approximately to the four reference species and one sample group from BS (Figure 4, Figure S3). However, the pattern was not homogenous, and the samples showed varying levels of admixture among different clusters. Three genetic clusters, including P. sylvestris, P. mugo s.s. and a mixture of samples from populations of other taxa, were found at cpSSR loci ( Figure S4).
Forests 2020, 11, x FOR PEER REVIEW 9 of 18 Figure S3). However, the pattern was not homogenous, and the samples showed varying levels of admixture among different clusters. Three genetic clusters, including P. sylvestris, P. mugo s.s. and a mixture of samples from populations of other taxa, were found at cpSSR loci ( Figure S4).  Table 1).

Phenotypic Variability and Differentiation
The level of the intraspecific morphological variation differed among the taxa and populations from BS in terms of needle (Table S5) and cone characteristics (Table S6). The highest variability in most of the needle characteristics showed P. uliginosa and the population from BS; reference taxa variation was at a lower level (Table S5). Regarding the cone characteristics, the most variable ones were P. sylvestris and P. uliginosa (Table S6).
Most of the analyzed needle and cone characteristics differed significantly (p ≤ 0.01) between samples from BS and reference populations. The BS population significantly differed in terms of needle characteristics from every reference species, while in terms of cone characteristics it differed significantly from P. uncinata and P. mugo (p ≤ 0.01). The groups of individuals described in the putative hybrid population during material collection based on macro-morphology as P. sylvestris-like (BS_S), P. uliginosa-like (BS_M), hypothetical hybrids intermediate between these two taxa (BS_H) and individuals not determined (BS_N), appeared close to each other in terms of needle and cone characteristics. Each of these four groups did not reveal significant pairwise differences  Table 1).

Phenotypic Variability and Differentiation
The level of the intraspecific morphological variation differed among the taxa and populations from BS in terms of needle (Table S5) and cone characteristics (Table S6). The highest variability in most of the needle characteristics showed P. uliginosa and the population from BS; reference taxa variation was at a lower level (Table S5). Regarding the cone characteristics, the most variable ones were P. sylvestris and P. uliginosa (Table S6).
Most of the analyzed needle and cone characteristics differed significantly (p ≤ 0.01) between samples from BS and reference populations. The BS population significantly differed in terms of needle characteristics from every reference species, while in terms of cone characteristics it differed significantly from P. uncinata and P. mugo (p ≤ 0.01). The groups of individuals described in the putative hybrid population during material collection based on macro-morphology as P. sylvestris-like (BS_S), P. uliginosa-like (BS_M), hypothetical hybrids intermediate between these two taxa (BS_H) and individuals not determined (BS_N), appeared close to each other in terms of needle and cone characteristics. Each of these four groups did not reveal significant pairwise differences neither in terms of needle nor in terms of cone characteristics. The individuals determined during sampling as P. uliginosa-like (BS_M) in terms of needle characteristics revealed a relatively low number of traits differentiating them from reference populations of that species, but a relatively high number of characteristics differentiating them from P. mugo s.s. and P. uncinata. Twice as many characteristics differentiated the group BS_M from P. sylvestris populations (Table S7). However, there was no general rule in terms of cone characteristics (Table S8). Surprisingly, there were no differences between each group from BS and the P. uliginosa population from the Zieleniec Nature Reserve (UL2).
The centroids of four groups from BS were close to each other on the scatterplots of the PCA and fell between reference populations of P. sylvestris and the taxa representing the P. mugo complex.
The group assigned to P. sylvestris during sampling appeared only slightly similar to these species reference populations, and only in terms of needle characteristics ( Figure 5). The needle characteristics that mostly contributed to group differentiations were NRC, SC, SF, TN/WN, VBF and VBL, while the most important cone characteristics were CD, CL, CDM, DAU, TA, CVX and WA ( Figure S5, code descriptions in Tables S5 and S6).

Phenotypic Variability and Differentiation
The level of the intraspecific morphological variation differed among the taxa and populations from BS in terms of needle (Table S5) and cone characteristics (Table S6). The highest variability in most of the needle characteristics showed P. uliginosa and the population from BS; reference taxa variation was at a lower level (Table S5). Regarding the cone characteristics, the most variable ones were P. sylvestris and P. uliginosa (Table S6).
Most of the analyzed needle and cone characteristics differed significantly (p ≤ 0.01) between samples from BS and reference populations. The BS population significantly differed in terms of needle characteristics from every reference species, while in terms of cone characteristics it differed significantly from P. uncinata and P. mugo (p ≤ 0.01). The groups of individuals described in the putative hybrid population during material collection based on macro-morphology as P. sylvestris-like (BS_S), P. uliginosa-like (BS_M), hypothetical hybrids intermediate between these two taxa (BS_H) and individuals not determined (BS_N), appeared close to each other in terms of needle and cone characteristics. Each of these four groups did not reveal significant pairwise differences neither in terms of needle nor in terms of cone characteristics. The individuals determined during sampling as P. uliginosa-like (BS_M) in terms of needle characteristics revealed a relatively low number of traits differentiating them from reference populations of that species, but a relatively high number of characteristics differentiating them from P. mugo s.s. and P. uncinata. Twice as many characteristics differentiated the group BS_M from P. sylvestris populations (Table S7). However, there was no general rule in terms of cone characteristics (Table S8). Surprisingly, there were no differences between each group from BS and the P. uliginosa population from the Zieleniec Nature Reserve (UL2).
The centroids of four groups from BS were close to each other on the scatterplots of the PCA and fell between reference populations of P. sylvestris and the taxa representing the P. mugo complex. The needle characteristics on the dispersion diagram between the two first discrimination variables, responsible for 80% of the total variation, indicate the intermediate position of BS samples between P. uncinata and individuals of P. uliginosa and P. sylvestris. Some of BS individuals intermixed with individuals of P. uncinata, P. uliginosa and P. mugo, and entered the 95% confidence intervals of these species. Interestingly, BS individuals did not intermix with P. sylvetris specimens ( Figure S6A). In the cone characteristics, pines from BS partly intermixed with reference P. uncinata, P. uliginosa and P. sylvestris, but did not enter P. mugo s.s. at a 95% confidence interval ( Figure S6B).
The dispersion of BS samples described by the two first PCA variables also indicates their intermediate character among the compared taxa. In the needle characteristics, a great part of BS individuals, determined during sampling as P. uliginosa-like, was located between P. mugo s.s., P. uliginosa and P. uncinata, while the remaining ones were located between P. uliginosa and P. sylvestris. Individuals determined as BS_S were located between reference P. sylvestris and taxa of the P. mugo complex, while those determined as hybrids were dispersed among all the others ( Figure 6A). In the cone characteristics, BS_M individuals mostly revealed close connection to reference P. uliginosa, P. uncinata and P. mugo s.s., although some of them were located between P. sylvestris and taxa from the P. mugo complex, with high affinity to P. uncinata. P. mugo complex, with high affinity to P. uncinata.
In PCA analysis, trees from BS assigned to P. sylvestris were located between reference P. sylvestris and P. uncinata populations, while the putative hybrids were more dispersed ( Figure 6B). The three individuals identified as pure P. sylvestris in genetic analysis were close to this species´ reference populations on the PCA, both based on needle and cone characteristics (Figures 6) and in the discrimination analysis ( Figure S6).  In PCA analysis, trees from BS assigned to P. sylvestris were located between reference P. sylvestris and P. uncinata populations, while the putative hybrids were more dispersed ( Figure 6B). The three individuals identified as pure P. sylvestris in genetic analysis were close to this species´reference populations on the PCA, both based on needle and cone characteristics ( Figure 6) and in the discrimination analysis ( Figure S6).

Discussion
Due to the relatively recent speciation history and the possible gene exchange during speciation [13], the taxa of the P. mugo complex and P. sylvestris show high genetic similarity at nuclear and organellar genomes [62][63][64]71]. So far, species diagnostic markers were found at cpDNA discriminating P. sylvestris from the taxa of the P. mugo complex [33,35,70]. In advances on earlier studies, we provide another diagnostic marker at mtDNA, which discriminates between P. uncinata and the remaining species from the P. mugo complex. However, the unique P. uncinata mitotype was not observed in the studied BS population, indicating that P. uncinata did not contribute, at least as a seed donor, to the development of the group of pines at the study area. Pines from the BS stand had most of the mtDNA mitotypes observed at high frequency in the small, relict population of P. uliginosa surrounded with extensive P. sylvestris stands in the WęgliniecNature Reserve. However, the discriminating power of the mtDNA is still too low to be conclusive about any further genetic relationships between the samples and between P. mugo s.s. and P. uliginosa.
The indication of a hybrid origin of the individuals comes from the analysis of the paternally inherited plastid and the biparentally inherited nuclear genome [30,62,64,71]. Interestingly, only three individuals carrying the cpDNA genome typical for P. sylvestris were identified among the 64 individuals sampled: one of them was assigned based on the macro-morphological phenotype to P. sylvestris (BS_S), while the other two remained undetermined during sampling (BS_N). All remaining individuals analyzed from BS possessed plastid DNA characteristic of the P. mugo complex, independently of phenotype.
The plastid genome derived from the P. mugo complex observed in the 14 out of 15 trees classified phenotypically as P. sylvestris means that those samples represent cryptic hybrids, which resemble pure species and cannot be distinguished in the field via macro-morphological characteristics. Furthermore, despite the vast diversity of phenotypic forms sampled from Błędne Skały and preliminarily classified into three distinct phenotypic groups, they all showed a uniform chloroplast background and high genetic similarity to each other at nuclear loci. For instance, they had the lowest genetic distance and high similarity at the allelic frequency spectra that placed them together in clustering analysis as compared to the reference species populations. These results indicate that pines from the BS contact zone represent a highly hybridized population, and the presence of genetically pure individuals of P. sylvestris seems marginal, despite a number of macro-phenotypically trees deceptively resembling this species. Overall, pines from that area showed higher genetic similarity to P. uliginosa and P. uncinata, but not to P. mugo s.s. and P. sylvestris (Figure 2 and Figures S2 and S4). The predominance of hybrid individuals and their morphological differentiation could indicate either hybridization lasting for many generations with restricted inflow of genes of P. sylvestris, or more recent hybridization without competition of P. sylvestris individuals from the P. mugo complex. However, both hypotheses shall be verified in the special study applying more genetic markers at nuclear loci to reveal the proportion of genetic background from each parental species at much finer resolution.
Some discordance of needle characteristics and cpDNA markers in pines has been reported earlier for P. uncinata. Four individuals sampled as P. sylvestris appeared to be cryptic hybrids with cpDNA typical for the taxa of P. mugo and a set of needle characteristics that placed them at a marginal position of P. sylvestris individuals [31]. In the present biometric investigations, the majority of individuals determined during sampling as P. sylvestris did not resemble this species from the reference populations, neither in terms of needle ( Figure 6A), nor in terms of cone characteristics ( Figure 6B), falling between P. sylvestris and P. uliginosa reference populations (see also Figure S5). The only three individuals carrying P. sylvestris cpDNA also appeared close to P. sylvestris in the biometrical analyses on the needle and cone characteristics. During sampling, only one of them was classified as P. sylvestris based on the main phenotype traits used in the field of taxonomic classification of trees [15,16]. However, the remaining individuals classified as this species were intermediate between P. sylvestris and P. uliginosa, indicating a lack of correlation between habit form vs. needle and cone traits. Nevertheless, this observation has not been verified using statistical methods. The lack of correlation between morphological characteristics of hybrids has been reported earlier by Bobowicz et al. [89], but on restricted and relatively young material.
Interestingly, the individuals bearing cpDNA of P. sylvestris were close to the reference individuals of the species in terms of PCA and discrimination scatter plots for needle and cone characteristics ( Figure 6A,B and Figure S6) showing that the results of the biometrical analyses are similar to those of cpDNA used for distinguishing between P. sylvestris and the P. mugo complex. In spite of that, the 64 individuals representing the BS mixed population formed one group situated between P. sylvestris and the taxa of the P. mugo complex, without a clear division between P. sylvestris and P. uliginosa, as detected by Jasińska et al. [31] for P. sylvestris and P. uncinata.
The presence of P. uliginosa × P. sylvestris hybrids bearing the phenotype of P. sylvestris and the lack of a vice versa morphological combination indicates a potential unidirectional gene flow from P. uliginosa to P. sylvestris, without phenotype expression of the former taxon (e.g., [30,31] and the literature cited herein). The putative gene flow from P. uliginosa (also from other taxa of the P. mugo complex) to P. sylvestris could be a specific evolutionary attribute, which allow the taxa of the P. mugo complex to diverge from P. sylvestris. The known invasions due to gene flow by pollen are a frequent phenomenon in nature. The gene inflow of Quercus petraea (Matt.) Liebl. into Q. robur L. by pollen is a good example of such a process [90][91][92], which also takes place in the populations of other wind-pollinated central-European trees ( [93,94] and the literature cited herein). Examples of invasion by hybridization and gene flow were described also for other plant species [95][96][97], but only certain types of hybrids have adaptive advantages over other hybrid types and parental species in diverse environments. This hypothesis, however, shall be verified in the study on the genetic basis of hybrid adaptation, which requires in-depth investigation of polymorphisms at various genomic regions [98].
Predominance of individuals bearing cpDNA barcodes of the P. mugo complex suggests long-lasting hybridization in isolation from other populations of P. sylvestris, as a result of vegetation changes during the Holocene [44,45,48]. The occurrence of pines in the Sudetes was strongly reduced during the Boreal period of the Holocene. This was mainly due to the expansion of Picea, which pushed pines to specific sites, such as rocky ridges or rock tops and, consequently, involved the spatial isolation between their populations. Thus, in Błędne Skały we observe the result of hybridization between P. sylvestris and P. uliginosa without or with highly restricted gene inflow from other populations. Isolation together with unidirectional gene flow from the P. mugo complex to P. sylvestris, could explain the considerable presence of hybrid individuals there and the restricted number of specimens with a genome typical for the latter species.
Despite its hybrid origin, we detected a relatively low level of genetic variation in the mixed BS population in terms of mitochondrial and plastid DNA markers, lower than in the reference populations of P. mugo s.s., P. sylvestris and P. uliginosa [63][64][65]68,99,100]. Some discrepancy between the level of genetic variation at the genomes of different inheritance modes in pines suggests that a limited number of paternal trees contributed to the reproductive success of the population. However, this observation and scale of possible inbreeding would need some additional investigations including genetic variation of the progeny derived from seeds sampled in the population.

Conclusions
The investigated population consists mostly of hybrids, and the presence of pure P. sylvestris individuals seems marginal. Hybridization between P. sylvestris and the taxa of the P. mugo complex is predominantly unidirectional from the latter species to P. sylvestris. The macro-morphological characteristics of the latter parent tree are frequently conserved in hybrids, but we did not find correlation of those traits with morphological and anatomical characteristics of the needles and cone traits in our study.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/10/1086/s1, Figure S1: Median Joining network of the mitotypes of pines from the mixed population in Błędne Skały and the reference four pine species populations based on 14 mtDNA regions, Figure S2: UPGMA tree based on Nei genetic distance at 14 mtDNA regions showing genetic relationships between defined groups of samples from the mixed pine population in Błędne Skały and reference pure species populations, Figure S3: ∆K and L(K) indices showing the most likely partitioning of the genetic variation at nSSR loci in Structure analysis, Figure S4: STRUCTURE results indicating three genetic groups of samples based on genetic variation at cpSSR loci (above) with optimal number of clusters detected using the method of Evanno's ∆K among compared pine populations (below), Figure S5: Influence of needle (A) and cone (B) characteristics on the dispersion of samples of four sample groups from Błędne Skały and reference pine taxa, Figure S6: Results of discrimination analysis for mixed populations from Błędne Skały (BS), P. sylvestris (PS), P. mugo (PM), P. uliginosa (PUL) and P. uncinata (PUN) based on the characteristic of needles (A) and cones (B), Table S1: mtDNA haplotypes and their frequencies in the four groups of pines from Błędne Skały, Table S2: Frequency of the 20 most frequent cpDNA haplotypes (present in > 3 individuals) out of 302 detected ones in the four groups of pines from Błędne Skały, Table S3: Genetic variation at individual populations based on cpSSR loci, Table S4: Genetic variation of samples from the mixed pine population in Błędne Skały and reference populations of the four pine species based on nDNA SSR loci, Table S5: Basic summary statistics of analysed needle traits of pines from the Błędne Skały (BS) and of the taxa compared-P. sylvestris (S), P. mugo (M), P. uliginosa (UL) and P. uncinata (UN), Table S6: Basic summary statistics of analysed cone characteristics of pines from Błędne Skały (BS) and of the taxa compared-P. sylvestris (S), P. mugo (M), P. uliginosa (UL) and P. uncinata (UN), Table S7: Needle characteristics differing at p ≤ 0.01 and p ≤ 0.05 in the Tukey's HSD and Kruskal-Wallis tests between subpopulations from Błędne Skały and compared populations of P. sylvestris, P. mugo, P. uliginosa and P. uncinata, Table S8: Cone characteristics differing at p ≤ 0.01 and p ≤ 0.05 in Tukey's HSD and Kruskal-Wallis tests between subpopulations from Błędne Skały and compared populations of P. sylvestris, P. mugo, P. uliginosa and P. uncinata.