Genetic and morphometric variability between populations of Betula ×oycoviensis from Poland and Czechia: A revised view of the taxonomic treatment of the Ojców birch

Birches are generally known for their high genetic and morphological variability, which has resulted in the description of many species. Ojców birch was described in 1809 by Willibald Suibert Joseph Gottlieb Besser in Poland. Since then, several studies assessing its taxonomy were conducted. Today, various authors present Ojców birch at different taxonomic ranks. In Czechia, the Ojców birch is classified a critically endangered taxon and confirmed at one locality consisting of several tens of individuals. However, before a strategy for its conservation can be applied, we consider it necessary to assess the taxonomic position of the endangered Czech population and to evaluate its relationship to the original Polish population. This study aimed to evaluate the morphometric and genetic variability between populations of B. ×oycoviensis in Poland and the Czechia and their relationship to regional populations of B. pendula, one of the putative parental species of the Ojców birch. Altogether, 106 individuals were sampled, including the holotype of B. szaferi, the second putative parental species of B. ×oycoviensis, received from the herbarium of W. Szafer, which is deposited at the Institute of Botany in Kraków. Morphological analyses identified differences in leaves between B. ×oycoviensis and B. pendula. However, no significant differences were found in genome size between selected taxa/working units except for B. pendula sampled in Czechia. The identified difference of the Czech population of B. pendula is probably caused by population variability. Genetic variability between all the taxa under comparison, regardless of their origin, was also very low; only the benchmark taxa (B. nana and B. humilis) clearly differed from all samples analyzed. The results indicate minute morphological and negligible genetic variability between the Czech and Polish populations of B. ×oycoviensis. In light of our results, the classification of B. ×oycoviensis as B. pendula var. oycoviensis seems more accurate than all hitherto presented alternatives (e.g. B. ×oycoviensis as a separate species).

Introduction species, B. szaferi, was available to us (original acknowledged specimen-holotype) only as a herbarium item from the collections in Kraków (Poland), see Material and methods.
Besides the taxonomic perspective, such studies are also desirable from a practical standpoint. Management steps should be taken towards the effective conservation of the Czech Ojców birch population because it is growing old and its natural regeneration is poor. However, these steps should reflect the taxonomical position of the Ojców birch, so it is important to determine whether its Czech population is taxonomically identical to that in Poland. If it is, Polish individuals could theoretically be used to strengthen the Czech population.
For assessing the variability between B. ×oycoviensis and closely related B. pendula in Poland and Czechia, we employed leaf morphometry, genome size analysis and microsatellite analysis. Previously, leaf morphometry was employed for the study of relationships between birch taxa, for example by Gardiner et al. [29], Atkinson [30] or Gill and Davy [31], whose studies served as a basis for the selection of parameters used in this study. The presented microsatellite analysis was recently used for the determination of taxonomic relationships of Sample locations for additional cytometric analyses (see the Discussion) are depicted using grey circles. Locations with the acknowledged (generally accepted) occurrence of B. ×oycoviensis are marked with filled red triangles (or red dots inside of black circles, if these locations were sampled within our study). Empty red triangles depict locations of proposed (unverified) occurrence of B. ×oycoviensis according to Database of the Czech flora and vegetation PLADIAS (https://pladias.cz/en/ to date of 2020/04/07) and study by Staszkiewicz [18].

Plant material
Individuals of B. ×oycoviensis for our study were sampled between 2017 and 2019 at Dolina Kobylańska and Skielek in Poland and at Volyně u Výsluní (in the Ore Mts) in Czechia. The specimens of B. ×oycoviensis were complemented by samples taken in the Botanical Garden of the Jagiellonian University (BGJU) in Kraków (Poland) and at Chomutov ZOO (Czechia; the sampled trees were transplanted from Volyně in the past). One birch showing several traits of B. ×oycoviensis was found and sampled in the Křivoklátsko region (Czechia). Reference samples of B. pendula originate from the localities Božídarské rašeliniště, Jezerka, Novodomské rašeliniště and Ú hošť in the Ore Mts. (Czechia), and Skielek (Poland). The Czech samples of B. pendula, together with a sample of B. nana described below, were collected in 2012 and used in our previous study [32]. As benchmark species for the evaluation, we included two additional birch taxa (B. nana, also collected in 2012 at Božídarské rašeliniště in the Ore Mts., Czechia, ca. 24 km to the west of the town of Volyně, and B. humilis, which was provided by the BGJU). The exact location (using Leica GPS 1200 at Volyně and Garmin GPSMAP 62S in other cases) and assignment of each sampled individual was recorded in the field. Specific individuals showing some traits of B. szaferi were cultivated at the Research Station Truba near the town of Kostelec nad Černými lesy in Czechia ('cult. B. szaferi' for further reference). The ancestors of the 'cultivated B. szaferi' specimens from Truba were individuals of B. ×oycoviensis growing at Volyně u Výsluní, whose seeds were collected in 2017. One dry leaf of the holotype (KRAM 303846) of B. szaferi [19] was provided by the herbarium of the W. Szafer Institute of Botany of the Polish Academy of Science (PAS) in Kraków (its photo is included as S1 Fig).
The determination of the taxa was done in the field as follows: As B. ×oycoviensis we determined individuals that exhibited at least 80% of the traits described in the relevant literature [16,19], and the rest were determined as B. pendula. All other transitional individuals, which combined traits of the two taxa, were determined as B. pendula/B. ×oycoviensis (a 'working unit' sensu Kuneš et al. [32]). Besides the benchmark species B. nana and B. humilis, we also collected sample material from B. atrata and B. obscura, which are visually distinct darkbarked birches occurring at the sites of sampling (conducted in 2019), although the taxonomic significance of these dark-barked birches remains questioned. We determined these taxa on the basis of morphological traits described by Domin [10], Hejtmánek [36] and Franiel [37]. Individuals of cult. B. szaferi were classified according to traits referred by Staszkiewicz [19]. However, in contrast to the holotype specimen of B. szaferi from the herbarium in Kraków, the individuals of 'cult. B. szaferi' grown at Truba should be viewed more as a working unit than as an acknowledged species. The determination of the remaining taxa (B. nana and B. humilis) is straightforward. From most of the sampled tree individuals, we took three short branches (two for morphometric analysis and one for flow cytometry and molecular analyses).
In total, we sampled twelve taxa and working units from twelve localities ( Table 1, Fig 1). We provided a simple table with basic description of all involved taxa in this study as (S1 Table). To simplify the text, we use the term 'group' to refer to each taxon or working unit with the country of origin. The limited number of specimens of B. ×oycoviensis in the study is given by the fact that the wild Czech and Polish populations are small. Moreover, the individuals with sufficiently manifested traits of Ojców birch are scattered in stands of B. pendula and other species, and their proportion is low. Even in BGJU in Kraków, only a few trees retained traits of the Ojców birch to such an extent that it allowed us to classify them as unequivocally identifiable B. ×oycoviensis.

Morphological analyses
Whenever possible, we collected two branches from each sampled individual for analysis of leaf morphology using telescopic shears. From each branch, we randomly selected two leaves for taking measurements (in total, we assessed four leaves per each analyzed individual). We measured following traits on each leaf: blade length [ we obtained and averaged only two measurements because it was not possible to take two branches from these samples (the individuals were too young), and in the case of B. szaferi from the herbarium of PAS we obtained only two measurements because we used a scanned image of the herbarium specimen. We analyzed this specimen by taking measurements of its image using the software ImageJ 1.52 [38]. We did not measure the morphological features of the benchmark taxa (B. nana and B. humilis), because it was not possible to measure all parameters on the leaf (4 th veins were not present on their leaves). Except for B. szaferi (Herbarium), we measured all specimens using a simple ruler and protractor. We did not include some individuals in the morphological analyses, because it was not possible to take enough material under field conditions; for example, the crown was sometimes too high to allow the sampling of sufficiently developed, insolated leaves for morphometric measurements (e.g. B. atrata and B. obscura in our study). Further methodological information on the measuring of morphological parameters of birch leaves are provided by previous studies [39,40]. In total, we selected 87 individuals for morphological analyses.

Genome size analysis
The genome size of sample plants was determined by propidium iodide flow cytometry [41] with Solanum pseudocapsicum L. (2C = 2.61 pg) as the internal standard. Leaf tissue from two petioles of each sample was chopped together with about 1.5 cm 2 of S. pseudocapsicum leaf tissue in 0.5 ml of Otto I buffer [42]. The resulting suspension was filtered through a 42-μm nylon mesh and left still for ca 20 minutes at 20˚C. After that, the suspension was stained with a solution of the following composition: 1 ml of Otto II buffer [42], β-mercaptoethanol (2μl/ ml), propidium iodide and RNase IIA. Flow cytometry was performed using a Partec CyFlow flow cytometer (Partec, Germany) equipped with a green solid-state laser (Cobolt Samba, 532 nm, 100 mW). Holoploid genome size and 1Cx-values (i.e. holoploid genome size divided by the number of chromosome sets [43]) was computed from raw cytometric data using FloMax software and evaluated statistically.

Molecular analyses (SSR genotyping)
Before DNA extraction, individual samples (leaves), stored at −80˚C, were frozen in liquid nitrogen and ground by an oscillation mill Retsch MM400 (Retsch GmbH, Haan, Germany). DNA was then extracted from the powdered samples using the QIAGEN DNEasy Plant Mini Kit (QUIAGEN, Hilden, Germany) following the standard protocol. The concentration and quality of DNA was checked prior to PCR using a NanoDrop 2000 (ThermoFisher Scientific, Waltham, Massachusetts, USA) spectrophotometer. For PCR, each DNA sample was diluted to a concentration of 10 ng/μl.
The genetic diversity of the taxa under study was assessed by microsatellite analysis of nuclear DNA. Microsatellite markers were selected based on publications by Tsuda et al. [34,35] and Kulju et al. [33]. In total, 50 markers were tested, from which 12 polymorphic loci were selected (see S2 Table) and optimized for two multiplex groups, identically as in a previous study by Kuneš et al. [32].
The PCR reactions were run in a total volume of 20 μl: 15 ng of DNA, primers (0.25 μM of each primer), 200 μM of dNTP, 2.5 mM MgCl 2 and 1 μM buffer PCR multiplex mix with polymerase.
The volume of 1 μl of PCR product was added to 14 μl of a solution of formamide with Gen-eTrace500 DNA ladder (Carolina Biosystems, Prague, Czechia), prepared according to the manufacturer's protocol. The solution was denatured (5 minutes, 95˚C) and quickly cooled on ice.
For determining the length of the amplicons containing the microsatellites, the genetic sequencer Genetic Analyser 3500 (Applied Biosystems, Foster City, California, USA) was used. Raw sequence data were analyzed by GeneMapper 4.1 software (Applied Biosystems, Foster City, California, USA).

Statistical analyses and computations
Morphometric analysis. To compare foliar features between groups (taxa/working units), we performed multivariate analysis of variance (MANOVA) and canonical discriminant analysis (CDA). The results of the CDA are presented as scatterplots. We evaluated the accuracy of CDA discrimination based on ratios of successfully classified individuals for each working group. We performed all computations in R software [44] and produced scatterplots using the R package 'ggplot2' [45] and the MorphoTools function set for R [46].
We tested for differences in individual parameters between unequivocally identifiable B. pendula and B. ×oycoviensis individuals across all populations taken together using the t-test (when all assumptions were met) or the Wilcoxon rank-sum test (when they were not). We evaluated the normality of our data was using the Shapiro-Wilk normality test. When the assumption of normality was met but the datasets differed significantly in variances (tested by the F-test of equality of variances), we used Welch test instead. We performed all statistical computations in R [44], with the significance level set to α = 0.05.
Genome size analysis. We did not evaluate the benchmark species (B. nana, B. humilis) and dark-barked birches (B. atrata and B. obscura) in this analysis, as only a few samples of these taxa were available. The Betula szaferi sample obtained from the Kraków herbarium was analyzed without success, as it was impossible to perform flow cytometry on this old dry sample. Flow cytometry was also not successful for five samples of B. pendula from Czechia, one sample of B. ×oycoviensis from Poland, two samples of B. ×oycoviensis from Czechia and one sample of 'cultivated B. szaferi'. Therefore, flow cytometric analysis was accomplished in 92 samples.
Differences in 1Cx values between working groups were tested for by a Kruskal-Wallis test with multiple comparisons [47], as the assumptions of ANOVA were not met. A box-plot presenting statistically significant differences was produced in R software [44] using the package 'ggplot2' [45]. All statistical tests were performed using the significance level of α = 0.05. Molecular analyses. Genotype data were checked for errors and invalid records. Basic parameters (total allele count, number of used loci and total number of individuals) were evaluated together with basic F-statistics (F ST , F IS , F IT ) [48] between selected taxa (working groups). For each of the selected taxa (working groups), the numbers of private alleles (those which did not occur in any other group at the respective loci) across all loci were also evaluated.
To visualize relationships between individuals in the study, discriminant analysis of principal components (DAPC) [49] was performed. The structure of the dataset was also evaluated by clustering analysis in STRUCTURE 2.3.4 [50][51][52] using the following settings: K-1 to 15 (number of groups + 3), number of runs for each K-30, burn-in-100,000 repeats, MCMC repeats-100,000 repeats. An admixture model with correlated allele frequencies was used. The remaining parameters were left at their default settings. The actual number of populations in the dataset was determined by the method presented by Evanno [53] using Structure Harvester [54]. Cluster matching and permutation was performed in CLUMPP 1.1.2, utilizing the Greedy algorithm with 1,000 repeats [55]. STRUCTURE plots were produced using DIS-TRUCT [56]. Testing for differentiation between selected taxa (working groups) was done using Goudet's G-statistic Monte Carlo test [57].

Morphometric analysis
As is apparent from the CDA scatterplot (Fig 2), there was relatively low variability in morphological traits of leaves among most of the individuals except for the sample of B. szaferi (taken from the herbarium), which is clearly separated on the CDA scatterplot. When the original sample of B. szaferi was excluded from the analysis, samples of 'cultivated B. szaferi' appeared outside the point cloud. After the removal of B. szaferi, also the distinction between B. ×oycoviensis and B. pendula, and to some extent between the Czech and Polish populations of B. pendula and B. ×oycoviensis, became more visible. Significant differences between B. pendula and B. ×oycoviensis were identified in several leaf parameters by statistical analyses ( Table 2).
Trees classified as B. pendula and B. ×oycoviensis (those which showed at least 80% of traits described in the relevant literature, see the Material and methods) differed significantly (p < 0.05) from each other in 13 out of the 16 morphometric parameters analyzed (see Table 2).

Molecular analyses
Altogether, 116 alleles at twelve loci were identified in a total of 106 individuals. Private alleles were found in seven groups out of the twelve groups examined. The greatest numbers of private alleles per individual were found in the cases of B. humilis and B. nana (benchmark taxa). Among the other groups, the greatest number (proportion) of private alleles per individual was found for B. pendula from Czechia ( Table 3).
The overall DAPC analysis (with all samples included in the analysis) clearly separated the benchmark taxa (B. nana and B. humilis) whereas all other groups formed a compact cloud with relatively low variation (Fig 4).
The homogeneity between all groups (working units) except for B. nana and B. humilis was also supported by analyses in STRUCTURE and CLUMPP software. The number of clusters in the whole dataset was set to 4, based on an analysis according to Evanno et al. [53] in Structure Harvester [54]; for more information, see S2 Table. The benchmark taxa were included in a separate group (denoted by orange color in Fig 6) whereas all other groups formed quite a homogeneous cluster.
Like Staszkiewicz [19], we found, at sites of reported occurrences of B. ×oycoviensis the individuals showing the traits of B. ×oycoviensis scattered across stands of typical B. pendula or mixtures of typical B. pendula and other tree genera (e.g. Picea, Acer and Fagus). Our search for B. szaferi in the wild of Czechia and Poland was not successful.
However, we examined several individuals resembling B. szaferi that were grown from seeds of specimens classified as B. ×oycoviensis, sampled in Czechia (locality Volyně u Výsluní). Similar results were obtained by Jentys-Szaferowa [28], who described three types of progeny obtained by controlled crossings of B. ×oycoviensis: type oycoviensis, type pendula and type 'nova' (lately 'szaferi'). Both B. ×oycoviensis and B. szaferi were originally classified as species. The existence of three types of progeny could imply a 'hybrid' origin of B. ×oycoviensis but is definitely not sufficient to prove such a hypothesis. Moreover, no difference between B. szaferi and B. pendula have been observed using molecular methods (see . Therefore, B. ×oycoviensis could have originated by the crossing of visually somewhat distinct parents belonging to the population of B. pendula. This theory would also explain the occurrence of B. ×oycoviensis in countries other than Poland. The question is what traits actually differentiate species from taxa of lower taxonomical ranks, which is rather a conceptual issue [63,64], and how these traits are controlled and linked. The testing for differences in leaf morphology via MANOVA yielded significant results. Subsequent CDA analysis revealed distinct separation of the sample of B. szaferi taken from  the Krakow herbarium. The separation of this particular sample is not surprising, as the peculiarity of its leaf shape is apparent to the naked eye. When the sample of B. szaferi was removed from the CDA analysis, samples of 'cultivated B. szaferi' appeared as outliers, which could be expected too, as those samples were selected based on their distinctive leaf morphology. The overall leaf shape of the samples of B. pendula and B. ×oycoviensis (including samples with mixed traits) showed some minor differences; for example the Czech population of B. ×oycoviensis as that of B. pendula appeared to be slightly distinctive from its Polish counterparts. This could be a 'random effect' of geographical distance, but it could also be a result of minor phenotypic and genetic differences, possibly reflected in slightly different genome size of the Czech population of B. pendula compared to that of other selected working units (see below). This topic, however, surely merits a detailed study, and our data are not sufficient to draw any definite conclusion. Our testing of individual morphological traits revealed significant differences between clearly identifiable samples of B. pendula and B. ×oycoviensis in thirteen parameters out of the sixteen tested, namely all parameters tested except for blade fitting angle [˚], petiole length [mm] and basal angle [˚]. These outcomes of our morphological studies are more or less consistent with the results obtained by Baláš et al. [27], who found significant differences in seven  Each bar depicts one individual and its estimated probability of affiliation to a group denoted by its color. The number of groups in the dataset was estimated to be 4 using the method of Evanno et al. [53] using STRUCTURE HARVESTER [54]. The plot was created by DISTRUCT software [56].
https://doi.org/10.1371/journal.pone.0243310.g006 out of sixteen morphological traits, but it is important to mention that the study by Baláš et al. [27] was performed only on samples collected at Volyně u Výsluní, Czechia. Blade length and the number of leaf veins are also mentioned as important for the determination of B. ×oycoviensis by the keys to the floras of Poland [20] and Czechia [65], together with the growth of leaves in groups on brachyblasts in B. ×oycoviensis. An overview of the whole habitus is needed for the determination of B. ×oycoviensis in the field, as the leaf shape differences between the taxa in question do not have to be distinctive enough. Some habitus traits are also described in the keys to the floras mentioned above. Still, many individuals exhibit mixed traits of B. pendula and B. ×oycoviensis, so it is quite difficult to classify them [66]. That is also the main reason why we distinguished the B. pendula/B. ×oycoviensis working unit in this study.
Our genome size analysis revealed no significant differences between the defined groups except for B. pendula originating from the Czech population (the Ore Mts. and their surroundings), which displayed significantly greater values compared to the other groups. Because B. pendula from Czechia does not differ from other groups in the results of our molecular analyses, the reason behind the differences in genome size could reside not only in intraspecific ploidy-related or monoploid-level variation [67] but also in some kind of chromosomal disorder. We, however, attribute these differences to intraspecific variability because our as yet unpublished data indicate relatively extensive variation in genome size (1Cx values, respectively) between populations of B. pendula in Czechia (see Fig 7). However, the difference could easily be also a result of some level of methodological inaccuracy. Therefore, a survey focused on genome size across B. pendula populations may be desirable to check for the existence of a pattern of genome size with respect to geographic distribution.
Our molecular analyses detected relatively low inter-group variability, which suggests that all of the groups considered are genetically highly homogeneous (including the Polish and Czech B. ×oycoviensis and B. pendula). The same results concerning B. ×oycoviensis were reported in a previous paper by Kuneš et al. [32], based, however, only on a small number of B. ×oycoviensis samples from Czechia. The only clearly distinguished samples were those of the diploid benchmark species (B. humilis and B. nana), which were visibly separated from the other samples in the DAPC scatterplot (Fig 4). When these benchmark taxa were removed, only small differences could be observed, for example between samples of 'cultivated B. szaferi' and samples of B. ×oycoviensis from Poland or between samples of B. ×oycoviensis from Czechia and Poland. However, in spite of some gentle difference, the Polish and Czech populations of the Ojców birch most probably belong to the same taxon. The low-level variation is probably caused by some population variability or random effects.
Aside from the benchmark taxa, the greatest mean number (proportion) of private alleles per sample was found in case of B. pendula sampled in Czechia (0.5), which is in line with the observed difference in genome size. This fact might suggest that B. ×oycoviensis originated from B. pendula which was 'slightly different' from B. pendula sensu stricto, but as already said, further research is needed to draw solid conclusions. Nevertheless, it can be stated that B. pendula, B. ×oycoviensis and B. szaferi in our dataset do not differ genetically at the species level. Ashburner & McAllister [2] argued that the Ojców birch is a weak growing, precocious and heavily fruiting form of B. pendula. More recently, the chapter on birches in the Flora of the Czech Republic [21] classified B. ×oycoviensis as a variety of B. pendula-B. pendula var. oycoviensis. The experimental data summarized in our study support these opinions.
In Poland, B. ×oycoviensis is a species strictly protected by law (Ministry Decree: Poz. 1409 Rozporządzenie ministra środowiska z dnia z dnia 9 października 2014 r. w sprawie ochrony gatunkowej roślin). In Czechia, this taxon is included on the Red List of Vascular Plants, although it is not protected by law [68]. On the other hand, the Ojców birch is included neither on the European Red List of Trees [69] nor on the Red List of Betulaceae [70]. The exclusion of the Ojców birch from the latter two Red Lists might be attributable to its expected downgrading from the species level to a significantly lower taxonomic rank.
Nonetheless, we are convinced that the Ojców birch still deserves some level of protection and distinction in terms of conservation regardless of the taxonomic rank at which this birch is most appropriately classified. It is an interesting birch with a characteristic 'broomy' habitus whose origin is yet to be completely resolved. To answer the remaining questions or to test standing hypotheses related to the Ojców birch and its putative parent (B. szaferi), we should keep protecting this rare taxon and its habitats. At present, conservationists are reducing the emphasis on species conservation and are becoming increasingly aware of biodiversity at all levels of the hierarchy of life [64]. We suggest that this should be reflected on some reasonable level also in case of the Ojców birch, at least before its taxonomic position is resolved satisfactorily. In our opinion, this holds true, especially if the rare populations of the Ojców birch in Czechia and Poland diminish in size. A good example of the usefulness of this approach is that of the curly birch (Betula pendula Roth. var. carelica [Merklin] Hämet-Ahti), whose unique traits are highly valued even though its taxonomic value may be low. Comparison of genome size of selected populations of B. pendula from Czechia, including some hitherto unpublished data obtained by an identical procedure as described in the material and methods section. Relatively close locations within Central Bohemia were selected to illustrate the relatively large genome size variability on a small area; samples from Southern Moravia are included as an outlier population. The number of samples analyzed is indicated below each box. The box and whiskers plots are presented in standard Tukey's design: Whiskers depict the minimum and maximum excluding outliers, black dots represent outliers (less than lower quartile-1.5 times the inter-quartile range and more than upper quartile + 1.5 times the inter-quartile range, respectively). https://doi.org/10.1371/journal.pone.0243310.g007

Conclusions
Betula pendula and B. ×oycoviensis in our dataset do not differ genetically at the species level despite being distinct morphologically. On a scatterplot produced by the DAPC method, excluding diploid benchmark species (B. nana and B. humilis), the holotype of B. szaferi rests in a cloud consisting of individuals of B. pendula and B. ×oycoviensis, although it is shifted to the part of the cloud represented by specimens of Czech origin.
Our molecular analyses detected low variability between the groups under comparison (taxa and working units of the Czech and Polish provenance) after the exclusion of the benchmark species. This low variability suggests that the Polish and Czech populations of B. ×oycoviensis are genetically very close even though some small differences related to geographic origin may be traced. The classification of the Ojców birch as B. pendula var. oycoviensis seems more accurate than its treatment as a separate hybridogenous species under the name B. ×oycoviensis.