Interspecific Hybridization and Introgression Influence Biodiversity—Based on Genetic Diversity of Central European Viola epipsila-V. palustris Complex

Abstract

The Viola epipsila-V. palustris complex is a highly taxonomically complicated group of species in its entire circumboreal range of distribution. Habitat loss, forest flooding, and hybridization could lead to the extinction of V. epipsila. A hybrid index and principal component analysis (PCA) were used to select qualitative and quantitative morphological features to distinguish parent species and hybrids, inter simple sequence repeat (ISSR) markers to determine the genetic diversity of the populations, flow cytometry to estimate the genome size (GS), and non-coding chloroplast DNA (cpDNA) regions to indicate the directions of crosses. All taxa are very morphologically variable, and their features can change within a season. The most stable feature is the distance of the bracts on the pedicel from the rhizome. The genetic diversity of all taxa populations is low and highly influenced by selfing and vegetative propagation. The population structure is differentiated: populations of V. epipsila or V. palustris, mixed populations with both parent species, F1 hybrids and populations with introgressive forms occur in different regions. The interspecific GS variation corresponds to the ploidy level (4x = 2.52 pg, 8x = 4.26 pg, 6x = 3.42 pg). Viola epipsila is the mother plant of the hybrids. Research has shown the risk of V. epipsila extinction in Central Europe and the importance of local populations in studying the role of hybridization in reducing/maintaining/increasing biodiversity.

1. Introduction

Genetic diversity of species is constantly changing under the influence of various factors. The consequences of global climate changes are: the future species distribution, the borders of species range and fate of a species or ecological communities, and biodiversity [1]. The persistence of existing levels of genetic diversity in plant populations under climate change determines species life history, availability of habitats or migration routes, breeding system, and current range sizes. Endemic species and ecosystems with high levels of spatial isolation facing the highest risk [2]. The quaternary ice ages are considered as the most important factor determining present species ranges. It is predicted that in the future such adverse factors will lead to drastic reduction of circumboreal species range of occurrence [3]. The genetic diversity is also influenced by hybridization, leading to loss of biodiversity (in disturbed habitats species might be displaced by interspecific hybrids) or to biodiversity enhancement (origin of new ecotypes and species). Contemporary (recent) hybridization might lead to genetic diversity of populations by introducing new hybrid genotypes resulting from the introgression of novel alleles; transgressive segregation; and finally, to local adaptation [4,5,6,7,8]. To detect hybridization and introgression, different methods and approaches have been used, including studying patterns of morphological, cytological, chemical characters, and molecular markers [9]. Nowadays, DNA sequence data and other sophisticated molecular methods allow us to distinguish types of hybrids, estimate the frequency of hybridization, and detect recent and ancient hybridization [10,11].

Viola epipsila Ledeb. (2n = 4x = 24, dwarf marsh violet, sec. Plagiostigma Godr., subsect. Stolonosae Kupffer, Violaceae) is a circumboreal species with two subspecies which ranges overlap in a small part of Western Siberia. According to older floras, the type subspecies (V. epipsila subsp. epipsila) was distributed mainly in Northern Europe with scattered localities in Central Europe extending from Germany to the Ural Mountains and to Western Siberia in the east [12,13,14,15,16,17,18]. The subspecies repens (Turkz. ex Trautv. and Meyer) W. Becker prolongs the range of V. epipsila eastwards from Western Siberia to the Chukotka Peninsula and then to NW North America from Alaska, the Yukon Territory, central Canada to the Hudson Bay [14,19,20,21,22]. Within the whole European range V. epipsila subsp. epipsila may occur sympatrically with closely related species V. palustris L. (marsh violet, 2n = 8x = 48) that facilitates hybridization [14,18,23,24,25].

Both species occur on marches, wetlands, which allows to predict how peatlands species range could shift under climate change. Loss of habitats and formation of hybrids better adapted to changing environments might lead to the extinction of a rare V. epipsila in European populations [14,18,23,24,25]. In some Polish, Finnish, and Swedish populations, the frequency of V. epipsila was drastically reduced with simultaneously increased interspecific hybrids (V. × ruprechtiana Borb., syn. V. × fennica F. Nyl.) [14,25] and newly described species Viola pubifolia (Kuta) G. H. Loos [26], probably of hybrid origin.

The research aimed to estimate the frequency of V. epipsila in selected Central European populations by: (1) evaluating genetic diversity of V. epipsila-V. palustris complex; (2) confirming the role of hybridization in V. epipsila extinction by identifying interspecific hybrids (V. epipsila × V. palustris) and crosses direction; (3) establishing current western European range border of V. epipsila.

2. Material and Methods

2.1. Plant Material Collection

For morphological analysis Viola epipsila, V. palustris and their putative hybrids were collected from NE Poland (Działdowo and Krutyń regions in the Masurian Lake District and Białowieża National Park in the North Podlasie Lowland) and were grown in the experimental garden conditions for continuous observation of morphological features in subsequent seasons. The selection of individuals from Polish populations for morphological analysis was due to the availability of fresh material. From the German and Lithuanian populations, the material was sent in the form of silica gel-protected leaves for genetic diversity research.

For molecular analysis leaves of V. epipsila, V. palustris and their putative hybrids were harvested from plants collected in natural sites and growing in experimental garden or leaves were picked directly in the field. Material originated from 15 populations from three countries of Central Europe (Table S1, Figure 1). Two to three healthy and fully developed leaves/plant were harvested randomly from individuals growing in a distance of minimum 5 steps (to avoid clonality), preserved in sterile tubes (F.L. MEDICAL, Torreglia, Italy) filled with silica gel (F.H.U. “DOR-CHEM”, Cracow, Poland) at room temperature not exceeding 25 °C.

For the nuclear DNA content fresh leaves of V. epipsila, V. palustris and their putative hybrids from five populations were used and samples from six populations were selected for cpDNA sequencing (Table S1).

Law-protected Viola epipsila and plants growing in nature reserves and national parks were collected under a permission in accordance with the relevant institutions.

Plants were identified in the field on the basis of pre-selected morphological characters [25]. Voucher specimens are deposited in the Herbarium of the Institute of Botany of the Jagiellonian University in Cracow, Poland (KRA, accession numbers: 0552067-0552072).

2.2. Morphological Analyses

A total of 101 individuals of V. epipsila, V. palustris, and putative hybrids were analyzed using qualitative and quantitative traits during flowering with chasmogamous flowers (CH). A data matrix for the hybrid index was created based on the following seven qualitative features and scored according to Anderson [27]: (1) pubescence of lower leaf surface (1, abundantly pubescent; 2, pubescent; 3, glabrous), (2) shape of leaf apex (1, acute; 2, intermediate; 3, rounded), (3) type of leaf margin (1, dentate; 3, crenate), (4) number of leaves rising from rhizome (1, typically two; 3, three or more), (5) dark red veins on lower leaf surface (1, absent; 2, in part of leaf surface; 3, on the entire leaf surface), (6) length of CH spur (1, approximately three times longer than calyx appendages; 2, two times longer than calyx appendages; 3, as long as appendages), and (7) pubescence of CH flower lateral petals (1, glabrous or slightly pubescent; 3, abundantly pubescent).

Eight following quantitative traits were analyzed: (1) length of CH sepals with calycine appendages, (2) length of CH calycine appendages, (3) length of CH spur, (4) length of CH lower petal with spur, (5) ratio of leaf length to width, (6) ratio of distance of bracts from rhizome to length of pedicel, (7) length of CH flower, and (8) width of CH flower.

Principal Component Analysis (PCA) for quantitative traits was performed using function prcomp in R v. 3.6.3 package stats (R Core Team, 2020) [28].

Plants were also observed at flowering of cleistogamous (CL) flowers, capsule development and seed set to compare the variability of traits (with the exception of CH flower morphology) throughout the season. Reproductive system of both parent species and their hybrids is mixed, they represent seasonal cleistogamy. Bud-like, obligate self-pollinated CL flowering does not overlap with colorful, open, cross-pollinated CH flowering [24]. At spring (April/May), CH flowers are produced following by CL flowers which period is longer, starting from May till late autumn (October/November).

2.3. DNA Extraction and ISSR Analysis

For isolation of total DNA from dried leaves of V. epipsila, V. palustris, and their putative hybrids, Plant & Fungi DNA Purification Kit and Lyse CT buffer (EurX Sp. z o.o., Gdansk, Poland) were used after testing several methods. For assessment of extracted DNA quality electrophoresis was performed using 1% agarose gel. Concentration of DNA was measured by NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

Ten primers of ISSR markers giving the highest number of polymorphic products were selected [29,30] (Table S2). 15 μL of reaction mixture used for amplification contained: a 1.5 μL 10 × DreamTaq Green Buffer, 1.2 U of DreamTaq DNA Polymerase (Thermo Fisher Scientific, Waltham, MA, USA), 0.6 μL of a primer (10 μM), 0.3 μL of a dNTPs (10 mM; Thermo Fisher Scientific, Waltham, MA, USA), and 15 ng of a template DNA. Amplifications were done in Bio-Rad T100 Thermal Cycler (Bio-Rad, Hercules, CA, USA) according to the program: initial denaturation at 94 °C for 5 min., 42 cycles: denaturation 94 °C, 59 s.; annealing temperature depended on the primer (Table S2), 59 s.; elongation 72 °C, 59 s.; final extending 72 °C, 7 min. The PCR products were separated in 1% agarose gel with 1 × TBE and SimplySafe (EURx Sp. z o.o., Gdansk, Poland) for about 100 min. at 120 V. Band patterns of 217 samples were observed and captured with a MultiDoc-It™ Imaging System with VisionWorks^® LS Analysis Software (UVP, Upland, CA, USA). Obtained results reproducibility was tested by repeating the PCR reaction cycle on selected samples from all populations.

The analyzed plant material was a priori divided into groups that comprised V. epipsila, V. palustris, and their putative hybrids, based on morphological characters previously selected [25]. ISSR polymorphisms and genetic diversity within and among subdivided groups were analyzed using POPGENE v. 1.32 [31] and FAMD v. 1.31 [32]. To evaluate the relationships between studied individuals, species and populations, we constructed a split phylogenetic network (NeighborNet) in SplitsTree v. 4.6 [33] and performed a principal coordinates analysis (PCoA), both based on Dice coefficient. Bootstrap for NeighborNet was calculated on 2000 replicates. Identification of admixture and inference of population structure was done using clustering method on unlinked dominant markers [34] in STRUCTURE v. 2.3.4 [35], which has been widely used for inference of population structure in polyploids using dominant markers [36,37,38,39,40] also for species of the genus Viola [41,42]. STRUCTURE analysis, implemented in R v. 3.6.3 package ParallelStructure [43], assumed admixture between populations, and correlated allele frequencies between clusters. Five independent runs were performed with burn-in of 2 × 10⁵ and 2 × 10⁶ Markov chain Monte Carlo replicates after burn-in. Clustering results were summed up in CLUMPAK [44] with LargeKGreedy search method and 2000 random input order repeats. In order to detect various degrees of admixture in putative hybrid populations, we chose K = 2, corresponding to the number of admixing genomes. To evaluate statistical significance of clustering, a hierarchical analysis of molecular variance (AMOVA) was performed using Arlequin v. 3.5.2.2 [45]. To detect putative introgressants, we calculated posterior probability of three hybrid categories (F1 hybrids, V. palustris and V. epipsila backcrosses) using NewHybrids ver. 1.0 [46] with 1 × 10⁶ Markov chain Monte Carlo sweeps and burn-in of 1 × 10⁵.

2.4. Genome Size Assessment

Flow cytometry (FCM) was used to estimate the nuclear DNA content in 22 plants of V. epipsila, V. palustris, and their putative hybrids; two to three fresh leaves per plant, depending on leaf size, were harvested (Table S1). Samples of leaves were prepared as previously described [47]; nuclei isolation buffer (200 mM Tris-Cl, pH 7.5, 4 mM MgCl₂ × 6 H₂O, 0.5% v/v Triton X-100) supplemented with 2% (w/v) polyvinylpyrrolidone (PVP-10), propidium iodide (PI; 50 μg/mL), and ribonuclease A (50 μg/mL). Solanum lycopersicum cv. Stupicke (2C = 1.96 pg) [48] served as an internal standard. For each sample, PI fluorescence was measured in at least 5000 nuclei, using a CyFlow SL Green (Partec GmbH, Münster, Germany) flow cytometer. The coefficient of variation (CV) of the G₀/G₁ peak of Viola ranged between 3.32% and 5.58%. Nuclear DNA content was calculated using the linear relationship between the ratio of the 2C peak positions of Viola/Solanum on a histogram of fluorescence intensities.

2.5. Sequencing of cpDNA

PCR reactions of non-coding regions of cpDNA including: trnH-psbA, trnS-trnG were performed on 13 selected plants of V. epipsila, V. palustris and their putative hybrids from six populations (Table S1) using trnH (GUG), psbA primers [41,49] and trnS, trnG primers [50,51]. A total of 15 μL of reaction mixture used for amplification trnH-psbA contained: a 1.5 μL 10 × DreamTaq Green Buffer, 1.125 U of DreamTaq DNA Polymerase (Thermo Scientific), 0.375 μL of each primer (10 μM), 0.375 μL of a dNTPs (10 mM; Thermo Fisher Scientific, Waltham, MA, USA) and 15 ng of a template DNA. Amplifications were done in Bio-Rad T100 Thermal Cycler (Bio-Rad, Hercules, CA, USA) according to the program: initial denaturation at 94 °C for 5 min., 35 cycles: denaturation 92 °C, 45 s.; annealing 60 °C, 45 s.; elongation 72 °C, 2 min., final extending 72 °C, 10 min. 20 μL of reaction mixture used for amplification trnS-trnG contained: a 2 μL 10 × DreamTaq Green Buffer, 1 U of DreamTaq DNA Polymerase (Thermo Fisher Scientific, Waltham, MA, USA), 1 μL of each primer (10 μM), 0.8 μL of a dNTPs (10 mM; Thermo Fisher Scientific, Waltham, MA, USA), and 10 ng of a template DNA. Amplifications were done in Bio-Rad T100 Thermal Cycler (Bio-Rad, Hercules, CA, USA) according to the program: initial denaturation at 94 °C for 5 min., 35 cycles: denaturation 94 °C, 30 s.; annealing 62 °C, 1 min.; elongation 72 °C, 1 min., final extending 72 °C, 10 min. The PCR products were separated in 1% agarose gel with 1 × TBE and SimplySafe (EURx Sp. z o.o., Gdansk, Poland) for about 90 min. at 120 V. Bands were observed and captured with a MultiDoc-It™ Imaging System with VisionWorks^® LS Analysis Software (UVP, Upland, CA, USA). Sequencing was performed by Genomed S.A., Warsaw, Poland. Sequences of trnH-psbA and trnS-trnG intergenic spacers were manually checked for quality and separately aligned using MAFFT v. 7.310 [52] L-INS-i algorithm. In trnH-psbA, a hypervariable microsatellite region between 127 and 141 bp was excluded from further analysis. Phylogenetic analysis was performed using Bayesian method in MrBayes v. 3.2.6 [53] with separate partitions set for each spacer and indels. F81 substitution model was chosen for trnH-psbA spacer and HKY for trnS-trnG, based on Bayesian information criterion (BIC) selection done in jModelTest v. 2.1.10 [54]. Indels were coded by simple coding [55] in SeqState [56] and treated as restriction sites with variable coding. Two independent runs with two Markov chains each, were performed for 3 × 10⁶ generations with burn-in of 2 × 10⁵. Sequences were deposited in GenBank with accession numbers MT450443-MT450452, MT450454-MT450466, and MT450468-MT450470.

3. Results

3.1. Variability of Morphological Features

Viola epipsila, V. palustris, and their putative hybrids are cleistogamous, developing open, colorful, cross-pollinated CH flowers and closed, bud-like, inconspicuous, greenish, obligate self-pollinated CL flowers in the same season (Figure 2). Cleistogamy is seasonal, and CL flowers develop after CH flowers. The morphological features of the flowers and leaves varied intraspecifically as well as within the same season based on observations of specimens growing in garden conditions for one or two seasons.

The period of CH flowering in all taxa studied was relatively short (2–3 weeks) and began, depending on the season, from mid-April to mid-May in the following order: V. epipsila first, followed by or at the same time as the putative hybrid, and V. palustris last. The flowering periods of the studied taxa overlapped, which allowed a comparison of the CH flower features. In contrast, the flowering period of the CL flowers is long and lasts until late autumn (October/November).

To select the features that potentially distinguish these taxa in the field, the qualitative and quantitative features of CH flowers and leaves of both parental species and hybrids were analyzed based on published floras [15,18] and previous data [25]. A field assessment of CH flower size and color, lateral petal pubescence, spur length and color (Figure 2A,B,F,G,J,K), leaf size, shape and pubescence of the lower surface appears to distinguish both parental species. Hybrids are difficult to recognize and they are often identify as V. epipsila or V. palustris. The results of the current analysis showed the weakness of some of these features as distinguishing taxa because of the high intraspecific variability.

The hybrid index based on the analysis of seven qualitative features of leaves and CH flowers clearly distinguished the two parental species and indicated high intraspecific morphological variation (Figure 3). The range of variability of third separated group consisting of hybrids (13–19), covered V. epipsila in small extent (9–14) and fell into V. palustris variability (15–20). The PCA for quantitative traits separated parental species and located hybrids mostly between both parents although they are closer to the V. palustris range of variability (Figure 4). Significant traits for parental species and hybrids delimitation were length of sepals with calycine appendages, length of lower petal with spur, ratio of distance of bracts from rhizome to length of pedicel, and length and width of flower (Figure 2A,B,F,G,J,K,O–R,

Table 1. Results of morphological analysis of V. epipsila, V. palustris and their putative hybrids during chasmogamous flowering using quantitative traits.

Taxon	N *	Mean of Length of Sepals with Calycine Appendages (mm) (±SD)	Mean of Length of Calycine Appendages (mm) (±SD)	Mean of Length of Spur (mm) (±SD)	Mean of Length of Lower Petal with Spur (mm) (±SD)	Mean of Ratio of Leaf Length to Width (±SD)	Mean of Ratio of Distance of Bracts from Rhizome to Length of Pedicel (±SD)	Mean of Length of Flower (mm) (±SD)	Mean of Width of Flower (mm) (±SD)
V. epipsila	40	6.68 a (±0.79)	1.41 a (±0.42)	9.70 a (±0.88)	17.85 a (±1.73)	0.79 a (±0.14)	0.69 a (±0.08)	21.55 a (±2.45)	24.88 a (±2.84)
putative hybrids	43	6.19 b (±0.69)	1.19 a (±0.29)	6.62 b (±0.80)	13.59 b (±1.25)	0.82 a (±0.12)	0.49 b (±0.07)	18.42 b (±1.40)	20.05 b (±2.11)
V. palustris	18	5.33 c (±0.69)	0.89 b (±0.27)	6.28 b (±0.57)	11.67 c (±1.19)	0.79 a (±0.15)	0.40 c (±0.08)	13.89 c (±2.49)	15.44 c (±2.83)

Average values followed by the same letter do not differ significantly at p ≤ 0.05 in Kruskal–Wallis test or one-way ANOVA and a Tukey HSD test post-hoc. * Number of samples depended on availability of plant material.

).

Capsules filled with seeds of CH flowers with visible straight style of pistil ended with stigma were developed in V. epipsila and V. palustris (Figure 2C,L) while the hybrid’s CH flowers were sterile; they formed small and deformed capsules (Figure 2H). Both parental species produced CL flowers (Figure 2D,M) and capsules with curved pistil style (Figure 2E,N). Hybrid CL flowers are arrested in development (Figure 2I).

The heterosis of hybrids was manifested both in the CH and CL flowering phases (Figure 2S).

At the mid and late of season, at the stage without CL flowers or CL capsules, both species and hybrids were almost indistinguishable (Figure 2T). Plants became very lush, leaf morphology was very diverse in all taxa (Figure 2Q), and some features, such as dark red veins on lower leaf surface, disappeared. The most stable feature of all taxa during the whole season was the ratio of distance of bracts from rhizome to length of pedicel.

3.2. Genetic Diversity of V. epipsila-V. palustis Complex Increased by Hybrids and Introgressants

The ISSR analysis resulted in 167 clearly resolved bands, 110 (65.87%) of all bands were polymorphic. Within V. epipsila group, number of polymorphic bands was the highest, ranging from 6 (3.59%) in BNP to 23 (13.77%) in L. Within V. palustris group number of polymorphic bands was slightly lower, ranging from 0 in BNP to 23 (13.77%) in L. The lowest share of polymorphic bands was observed in a group of putative hybrids, ranging from 0 in BNP to 6 (3.59%) in KOZ. The highest number of private markers was found in V. epipsila group (21), the lowest in V. palustris group (5). Discriminating markers were present in each group, ranging from 2 in V. palustris to 11 in a group of hybrids (

Table 2. Parameters of genetic diversity in V. epipsila, V. palustris and their putative interspecific hybrids. N—number of specimens used in genetic analyses; P—number of polymorphic markers; %_poly—proportion of polymorphic markers; N_prt —no. of private markers present in a given population but absent in other populations; N_d—no. of discriminating markers present in all individuals of a given population, H_j—Nei’s (1973) gene diversity; H_T—total gene diversity; H_S—mean gene diversity within populations; G_ST—Nei’s (1973) gene diversity between populations.

Taxon/Pop *	N	P(%poly)	Nprt	${\mathit{N}}_{\mathit{d}}$	Hj	HT	HS	GST
V. epipsila						0.095 ($\pm$0.027)	0.035 ($\pm$0.004)	0.636
BNP	13	6 (3.59)	7	4	0.015 ($\pm$0.080)
SZ	4	21 (12.57)	7	1	0.052 ($\pm$0.138)
L	10	23 (13.77)	7	1	0.036 ($\pm$0.098)
Total	27	55 (32.93)	21	6	0.083 ($\pm$0.155)
V. palustris						0.083 ($\pm$0.024)	0.012 ($\pm$0.001)	0.850
BNP	8	0 (0.00)	1	1	0.000 ($\pm$0.000)
SZ	5	13 (7.78)	1	0	0.035 ($\pm$0.123)
N_NE	9	0 (0.00)	0	0	0.000 ($\pm$0.000)
N_B	14	4 (2.40)	0	0	0.011 ($\pm$0.070)
N_SE	5	0 (0.00)	0	0	0.000 ($\pm$0.000)
BAL	10	0 (0.00)	1	1	0.000 ($\pm$0.000)
KR	21	1 (0.60)	0	0	0.002 ($\pm$0.024)
L	10	23 (13.77)	2	0	0.051 ($\pm$0.132)
Total	82	50 (29.94)	5	2	0.074 ($\pm$0.150)
hybrids						0.101 ($\pm$0.030)	0.003 ($\pm$0.000)	0.974
BNP	73	0 (0.00)	5	5	0.000 ($\pm$0.000)
SZ	6	0 (0.00)	0	0	0.000 ($\pm$0.000)
KOZ	15	6 (3.59)	1	0	0.008 ($\pm$0.042)
OST	8	2 (1.20)	0	0	0.006 ($\pm$0.051)
L	4	0 (0.00)	6	6	0.000 ($\pm$0.000)
Total	106	45 (26.95)	12	11	0.061 ($\pm$0.130)

* Populations with sample size less than two were not included in analysis (V. palustris—OST, hybrids—N_NE). Origin of samples: Germany (N), Poland (BAL, BNP, KR, KOZ, OST, and SZ) and Lithuania (L). Detailed information in Table S1.

). Nei’s gene diversity (Hj) was slightly higher for parental species (0.015–0.052 for V. epipsila and 0.000–0.051 for V. palustris) than for hybrids (0.000–0.008). Total gene diversity (H_T) was at similar level in V. epipsila (H_T = 0.095) and hybrids (H_T = 0.101), whereas in V. palustris was lowest (H_T = 0.083). Mean gene diversity within populations (H_S) was very low ranging from 0.003 in hybrids to 0.035 in V. epipsila while the gene diversity (G_ST) between populations was very high (from 0.636 for V. epipsila to 0.974 for hybrids) (Table 2).

The PCoA presented three groups clearly separated along major component corresponding to V. palustris, hybrids and V. epipsila. The lowest overall variance along both components was measured for hybrids (Figure 5). NeighborNet analysis based on Dice coefficient genetic distance matrix also shown three diverged groups subdivided by a strongly supported splits (bootstrap separating V. epipsila: 100%, bootstrap separating V. palustris: 99.4%). Samples from different populations within groups were also separated by strongly supported splits (Figure 6).

STRUCTURE Bayesian analysis based on species clustering assumed two groups, K = 2 (Figure 7A). The first group (green color) corresponded to V. palustris, the second (red color) to V. epipsila. There were four populations consisting of only V. palustris individuals (BAL, KR, N_B, N_SE), five populations with one or both parental species and hybrids (BNP, L, N_NE, OST, SZ) and one population with only hybrids (KOZ). Hybrids from BNP consisted equally of two genetic components, while hybrids from populations KOZ, OST, SZ, and N_NE indicated slight predominance of first genetic group over the second group and vice versa for hybrids from populations L. Some individuals from three populations (OST, SZ, and L) had a significant degree of admixture (Figure 7A). NewHybrids analysis indicated presence of hybrids of three categories. The first group (black color) corresponded to first generation hybrids, the second (dark green color) to V. palustris backcross and the third group (pink color) to V. epipsila backcross. First generation hybrids were identified in BNP, KOZ, N_NE, OST, and SZ populations. Introgressive form at varying degrees of introgression were found in L, OST and SZ populations. (Figure 7B).

The AMOVA showed high percentage of variation between V. epipsila, V. palustris and hybrids (59.92%) and higher variation among populations within groups (36.53%) than within populations (3.55%). The highest variation among populations was in hybrid group (98.73%) and in V. palustris group (94.76%). Variation index between groups corresponding to two parental species was higher (F_CT = 0.79) than between three groups containing also hybrids (F_CT = 0.60). The variation within the taxon was the highest in hybrids (F_ST = 0.99) and in V. palustris (F_ST = 0.95), lower in V. epipsila (F_ST = 0.73;

Table 3. Results of the three-level AMOVA, based on ISSR profiles scored from total of 10 populations of V. epipsila, V. palustris, and their putative hybrids.

Source of Variation	df	Sum of Squares	Variance Components	Percentage of Variation	F, p-Value
V. epipsila vs. V. palustris
among groups	1	1044.50	26.30	79.18	FCT = 0.79 P = 0.012
among populations within groups	6	366.19	5.89	17.72	FSG = 0.85 a
within populations	78	80.36	1.03	3.10	FST = 0.97 a
total	85	1491.05	33.22	-
V. epipsila vs. V. palustris vs. hybrids
among groups	2	1985.09	12.80	59.92	FCT = 0.60 a
among populations within groups	13	1089.63	7.80	36.53	FSG = 0.91 a
within populations	199	150.84	0.76	3.55	FST = 0.96 a
total	214	3225.56	21.36	-
V. epipsila
among populations	2	123.81	7.20	72.90	FST = 0.73 a
within populations	24	64.27	2.68	27.10
total	26	188.07	9.88	-
V. palustris
among populations	4	242.38	5.39	94.76	FST = 0.95 a
within populations	54	16.10	0.30	5.24
total	58	258.48	5.69	-
hybrids
among populations	5	545.03	10.08	98.73	FST = 0.99 a
within populations	101	13.08	0.130	1.27
total	106	558.11	10.21	-

The analysis is based on ISSR phenotypes containing 167 band states. Significance levels are based on 1023 iteration steps. df – degrees of freedom. ^a p < 0.001.

).

3.3. Genome Size as a Good Marker to Distinguish Parental Species and Their F1 Hybrids

Genome size of V. epipsila from BNP (2C = 2.52 pg ±0.03) and V. palustris from BNP, OST and KR populations (2C = 4.26 pg ±0.02) differed significantly (Table S3). Fifteen selected hybrids from four populations (BNP, SZ, KOZ, OST) possessed DNA content (mean 2C = 3.42 pg ±0.06) close to the sum of holoploid (1C) genome sizes of parental species (2C = 3.39 pg), confirming their hybridity. The Cx-values (DNA contents of a monoploid genome with chromosome base number x) for parental genotypes and their hybrid were 0.63, 0.53, and 0.57, respectively Table S3).

3.4. Viola epipsila Is Maternal Species of Hybrids Based on Non-Coding cpDNA Sequences

Two haplotypes were identified, first corresponded to V. epipsila and putative hybrids, second to V. palustris (Table S4) based on trnH-psbA and trnS-trnG intergenic spacer obtained from 13 individuals from BNP, KR, L, N_NE, N_B, N_B! populations, including V. epipsila, V. palustris, and putative hybrids (Figure 8). Sequences of the trnH-psbA were from 298 to 302 bp and trnS-trnG were from 560 to 587 bp in length. The two haplotypes diverged by total of four substitutions (3 – trnH-psbA, 1 – trnS-trnG) and five indels (3 – trnH-psbA, 2 – trnS-trnG) ranging from one to 11 nucleotides (Table S4).

4. Discussion

The results of our research along with the previous data on V. epipsila and V. palustris distribution in Poland [25], indicated a threat of extinction for V. epipsila in Central Europe. The reasons are environmental changes, habitat disappearance as a consequence of climate change and hybridization causing hybrids to displace one or both parental species, depending on the region. Sympatric Polish populations of V. epipsila and V. palustris are good local models for studying the role of hybridization in the decrease in the number of parental species.

The Viola epipsila-V. palustris complex is a highly taxonomically complicated group of species in its entire circumboreal distribution range. The genetic diversity of these closely related species may also be influenced by hybridization [14,18,23,24,25,57,58]. The situation is even more complicated because one of the putative hybrid parents (V. epipsila 2n = 4x = 24) likely participated in the formation of a species at a higher ploidy level (V. palustris 2n = 8x = 48) [23,24].

4.1. Difficulties in the Identification of V. epipsila, V. palustris and their Putative Hybrids in the Field

Intraspecific morphological variations result from habitat conditions but are also associated with the life cycle stage in a season and influenced by hybridization. These factors make species identification in the field very difficult and lead to incorrect species designations. The most critical period in a season is a stage of early cleistogamous flowering. During this period, the leaf features are very variable and the capsules of the CL flowers are not yet developed, which allows the identification of fully or partially sterile interspecific hybrids. In sympatric populations, distinguishing species and hybrids is particularly difficult. Our observations of both parental species and the hybrids growing in garden conditions in two subsequent seasons allowed us to conclude that the qualitative and quantitative features indicated by the authors of floras e.g., [15,18] as typical for the species are variable and should be applied with caution. For V. epipsila, the following features were assigned: leaf lower surface abundantly pubescent, acute leaf apex, dentate leaf margin, two leaves rising from the rhizome, CH flowers lateral petals glabrous or slightly pubescent, spur approximately three times longer than the calyx appendages, and capsules of CH and CL flowers filled with seeds; for V. palustris: leaf lower surface glabrous, rounded leaf apex, crenate leaf margin, three or more leaves rising from the rhizome, lateral petals of CH flowers abundantly pubescent, spur as long as the appendages, and capsules of CH and CL flowers filled with seeds. The hybrid morphology was described as intermediate to those of the parental species or more or less similar to either species, CH and CL flower capsules deformed, dwarf and without seeds.

During the CH flowering period, pubescence on the leaf lower surface and the lateral petal, flower size, bract position on the pedicel, capsule size, and seed production are taxonomically useful traits for distinguishing species and hybrids in the field. However, finding plants with flowers or capsules is not always possible in the field. The high diversity of plants, the similarity of the hybrids to the parents, and the morphological variability across the whole season might result in hybrids being unrecognized or erroneously identified as either parental species, which makes determining the actual frequency of hybrids in floras of different regions difficult [59,60]. The phenomenon of interspecific hybrids having new characteristics that do not occur in the parents or being morphologically very similar to one of the parents, e.g., in interploidy crosses to parents with higher ploidy levels or in introgressants originating via unidirectional introgression, is well known [9,61]. Spontaneous hybridization is nonrandomly distributed among taxa; the frequency often does not depend on the size of the family or the genus [60]. The Viola genus is a phylogenetic group that is biologically predisposed towards the generation and maintenance of hybrids via homoploid and/or alloploid speciation [14,18,23,24,25,41,57,58,62,63,64,65,66,67,68,69,70,71,72,73,74].

4.2. Is There a Risk of Extinction of V. epipsila in Central Europe?

The range of V. epipsila in Central Europe has been shrinking. According to older floras [12,13,14,15,16,17,18,75], the western European border of this species was located in Germany, north and east of the river Elbe, with the highest concentration found in the regions around the Baltic Sea. The new Verbreitungsatlas der Farn- und Blütenpflanzen [76] presents a distribution map of V. epipsila in Germany with its historical and current locations. This shows that the majority of population sites have not been confirmed in recent years, which suggests the disappearance of this species, resulting in an eastward shift of the western Central European border. Recently, a small population of V. epipsila was found in eastern Germany near the Polish border [77]. Of the 29 specimens from Germany in our study, no sample was revealed to be V. epipsila by ISSR markers.

In Poland, V. epipsila populations are located in the north-eastern part of the country [78]. There are lack of confirmed pure V. epipsila populations in western Poland as was shown based on analysis of herbarium specimens from Polish herbaria. Majority of specimens designed as V. epipsila were verified as hybrids [25].

In our studies only 17 (10.36%) of the 164 individuals analyzed by ISSR markers from north-eastern Poland were identified as V. epipsila, 39 (23.78%) as V. palustris, 6 (3.66%) as introgressive forms and the vast majority (102, 62.20%) were hybrids (Figure 7B; Table S3). In some Polish populations (BAL, KR) in which V. epipsila was observed previously [25], the species had disappeared. In the light of these data, the western boundary of V. epipsila occurrence is shifting to north eastern Poland, and the presence of hybrids in Lithuanian populations may suggest that V. epipsila is also endangered in this region (Figure 6 and Figure 7).

4.3. Local Sympatric Populations of V. epipsila and V. palustris as a Good Model for Studying the Hybridization Effect on Population Structure

Knowledge of the existence and frequency of hybrids in nature and their characteristics is important for learning about species interactions and maintenance. The occurrence of parental taxa and hybrids at one location may not be informative about other sympatric populations and the established genetic variation in parents and hybrids is not sufficient for determining the genotypes, phenotypes, and fitness of hybrids [79]. However, detailed, well-documented empirical studies of sympatric local population structures are the basis for further research with advanced molecular techniques.

The Białowieża National Park and Masurian Lake Region violet populations at the west margin of the V. epipsila range are good models for studying hybridization because on these regions, the environmental conditions are more extreme and stressful than those in the center of the distribution and they are influenced by human activity (especially in the Masurian Lake District). These conditions lead to the weakening of ecological barriers between species and create "hybridized" (intermediate) habitats that can be occupied by hybrids [79,80,81,82,83].

Environmental monitoring conducted for several decades in the BNP indicates the disappearance of V. epipsila and V. palustris habitats (flooded forests). The increasing number of Carpinus trees has reduced the biodiversity of ground vegetation [84], and some predictions can be made about the responses of both species to environmental changes based on the existing genetic variation and the ecological environment. Moreover, an analysis of the pollen collected in peat bogs in the BNP enabled the reconstruction of the vegetation history and climate fluctuations of the Late Glacial and Holocene. The most important finding from this analysis is that the impact of human activity was weak and dated to the late Neolithic or the early Bronze Age and that the disturbances (agricultural activity) did not notably break the continuity of the forest cover in the core area of the Białowieża forest throughout the Holocene to the present day [85,86]. The vegetation history, plant communities, and species fluctuations are correlated with climatic and environmental changes [87], rather than with the impact of human activity in the BNP. Furthermore, in the forest in the BNP, the areas of occurrence of both species overlap, which creates the possibility of hybridization [25,88]. In the Masurian Lake District, mixed populations with numerous hybrids were observed during the population exploration in 2016–2019 by the authors (LM & JZ). In the BNP, the frequency of hybrids was very high, as estimated by ISSR markers.

Hybrids represented 76.68% of mixed populations (102 individuals), whereas V. epipsila represented only 12.78% (17 individuals) and V. palustris represented 10.53% (14 individuals) (Figure 5, Figure 6 and Figure 7). In some populations (KOZ) in the Masurian Lake District, only hybrids were identified, whereas in other (OST and SZ) hybrids, introgressive forms (resulting from backcrosses to one parent or the other parent) were identified.

4.4. Low Intrapopulation Genetic Diversity of the V. epipsila-V. palustris Complex and Unidirectional Crosses

Genetic variation was estimated by ISSR, rather than by next-generation sequencing (NGS) markers due to the large sample number, the simplicity of this method and the ability to distinguish species and hybrids. Highly polymorphic ISSR markers have successfully been used in research on genetic diversity, phylogeny, and genome mapping [89].

The genetic diversity between populations was very high, with very low intrapopulational variation. The highest variation among hybrid populations (G_ST = 0.974) occurred due to their low fertility and vegetative propagation. For both parental species, the within-population variation was lower than that in plants with a predominance of cross-pollination [90,91,92], but the parental species produced cross-pollinated CH flowers. This resulted from the relatively short period of CH flowering (2–3 weeks in spring), whereas obligate self-pollinated CL flowering extended until the end of the season. This was most evident in the vegetatively propagated hybrids (H_S = 0.003) but also in V. palustris (H_S = 0.012) and V. epipsila (H_S = 0.035). Surprisingly, all studied hybrids originated from one crossing direction: V. epipsila as the maternal plant and V. palustris as the pollen donor, based on the trnS-trnG and trnH-psbA sequences. Viola epipsila and hybrids have a common haplotype that is different from V. palustris. The one-way crossing can be explained by the very low frequency of V. epipsila in the mixed populations studied.

The results indicating the occurrence of F1 hybrids in mixed populations (BNP) as well as introgressive forms with the predominance of either one or the other parental genome (L, OST, SZ), which is clearly indicated in STRUCTURE at K = 2 and in the NewHybrids analysis, are interesting. These forms could be introgressants towards V. palustris (L, OST, SZ) with a predominance of V. palustris in their genomes or towards V. epipsila (L). The question is how introgressants can arise if F1 hybrids are sterile, they do not produce seeds neither by CH nor by CL flowers? Both male and female lines of hybrids are strongly disturbed leading to formation of cytologically unbalanced male and female gametes and further to abnormal male and female gametophytes. Nevertheless, triploid hybrids produce stainable (viable) pollen, in very low frequency, with chromosome numbers 12, 24 (also unreduced with 36 chromosomes and aneuploid numbers) which could pollinate the parental species and could lead to the formation of introgressants. The abnormal development of the female gametophyte of F1 hybrids prevents the formation of the next generation of hybrids [24].

Introgression can be considered an intermediate stage for the stabilization of hybrids and the emergence of a new species, introducing new, possibly adaptive, genetic variation into a population. This process can also eliminate rare species that are less well adapted to a changing environment than newly originated species [11,93]. Introgression may be a more common phenomenon in nature than previously expected because of the difficulty of recognizing these forms by classical methods. New molecular techniques, genetic mapping, NGS technologies, pooled barcoded amplicon sequencing, and restriction site-associated DNA sequencing (RADseq) will allow research on wild hybridization and introgression as well as the genetic basis of hybrid fitness in changing environments [10,94,95,96,97].

4.5. Lack of Intraspecific Variation in Genome Size

Estimating the nuclear DNA content by FCM is a reliable method of distinguishing parental species from their hybrids because they differ in ploidy. Interspecific GS variability can be an effect of different ploidies, different activities and amplification of transposable elements and environmental and ecological factors, such as growing season length [98,99,100]. Intraspecific GS variability did not occur in V. epipsila, V. palustris and their hybrids which is in line with the suggestion by Greilhuber et al. [101] that intraspecific GS variation data result from methodological errors (“plastic genome or plastic data”) rather than from the influence of environmental conditions.

In GS evolution of polyploids contraction/downsizing or expansion/upsizing results in a lack of correlation between GS and ploidy level [102,103]. The mean 2C value of the octoploid V. palustris (4.26 pg) established here is not exactly that predicted (2-fold the 2C value of V. epipsila, 2.52 pg), suggesting genome downsizing during polyploidization. However, in the origin of V. palustris, in addition to V. epipsila, a second unknown species having a GS that is not exactly the same as that of V. epipsila could have been involved. Although hybridization has been shown to have the potential to alter the GS (homoploid hybrid speciation, allopolyploidy), hybridization itself did not increase the GS, as demonstrated in studies on Helianthus hybrids [104,105]. The GS of the V. epipsila × V. palustris hybrids studied in the present research confirms this statement, because their 2C values were the sum of the haploid parental genomes.

5. Conclusions and Further Perspectives

Viola epipsila and V. palustris occur on marches, wetlands, and could be considered as indicators of peatlands. The range of both species, especially V. epipsila, in Central Europe drastically decreased which indicates that the habitats of these species, i.e., peatlands, are disappearing due to climate change. Hybrids likely have other ecological preferences helping them to survive in changing environment.

We assume that the presence of hybrids (V. epipsila × V. palustris) and introgressants in the absence of V. epipsila or/and V. palustris is an indicator of changing peatlands (disappearance of peat bogs and marshes).

This research provides a good basis for further studies on the role of hybridization in Viola evolution: its impacts on variation maintenance, biodiversity reduction by species extinction, and the introduction of new variation by new hybrid-derived species.

Viola epipsila is declining in Central Europe, and its western boundary of distribution has shifted from Germany to north-eastern Poland. At the western species range limit, interspecific hybrids with the closely related species V. palustris have displaced V. epipsila. Reserves should be created, especially in Germany and in the Masurian Lake District region, where a few individuals still exist, and determining the structure of the Ukrainian and Russian populations to monitor V. epipsila ssp. epipsila occurrence and its eastern range limit is advisable. Ex situ conservation, including in vitro culture, which has successfully been used for endangered violets [106,107], could also help in species preservation.

The well-documented genetic structure of local mixed populations is an important source of knowledge of rare species extinction via hybridization and/or the identification of introgressive forms that introduce new variation.

The selection of non-variable morphological features in species and the hybrid delimitation proposed in this study will allow us to estimate the real frequency of hybrids and avoid errors in species designation.

The identification of hybrids by ISSR markers with unequal parental genome admixtures suggests that introgression in both directions influenced hybrid/species variation. If so, the sterility of F1 hybrids should be re-examined [24], especially in terms of pollen viability. Introgressants found in several Polish and Lithuanian populations showed that hybridization followed by backcrosses could lead to a new hybrid-derived species origin. This should be confirmed by using advanced molecular techniques, as glucose-6-phosphate isomerase (GPI) genes successfully confirmed the origin of V. pluviae [73] and restriction site associated DNA sequencing (RAD-Seq) to resolve the hybrid origin of Viola pubifolia (Kuta) G. H. Loos [25,26].