Analyses of sexual reproductive traits in Dactylorhiza majalis: a case study from East Germany

The orchid species Dactylorhiza majalis is endangered by continuing habitat destruction and fragmentation. This requires more detailed information with respect to its sexual reproduction, which is especially relevant for Germany, where from 10 % to 30 % of the world-wide remaining populations grow. In the present study, we determined both the numbers of growing and flowering individuals per stand with regard to D. majalis at 12 localities of Upper Lusatia, Saxony, Germany, during the season 2014. For up to 25 plants per stand, sexual reproduction was assessed by checking over the numbers of blossoms and fruits per inflorescence and by calculating percentages of seed fertilities from embryo-viability stains. Applying pair-wise statistical analyses, we found correlations between two of the above-mentioned traits as well as among the above-cited population-specific reproduction parameters and four out of six Ellenberg’s indicator values, which have been calculated to characterize local site conditions. We furthermore recorded both very poor and enhanced seed fertilities, clustering into two groups which were associated with the Ellenberg’s indicator value thermal continentality. Lower seed fertilities were generally detected in the northern lowlands, whereas D. majalis is probably able to compensate the unpleasant environments of the southern highlands by bearing more fertile seeds. Conducting genetic inventories with three nuclear microsatellites, the sampled seed-producing mother plants of both fertility groups differed by the opposite frequency distribution of two prominent genotypes DD and EE at locus ms14. These findings indicate a genetic selection due to adaptation to climatical stresses. Based on the additionally detected aberrant megasporogenesis, we propose that mother plants of homozygous genotype EE and their germ-cells are less affected by both aneuploidy and large deletions on the remaining chromosomes, and we assume that a linkage disequilibrium exists between such advantageous karyotypes and the studied microsatellite locus. Regarding the challenges of global warming, repeated inventories are finally recommended at all 12 stands in order to validate the long-term indicative properties of the discovered findings.


Introduction
The perennial broad-leaved marsh orchid Dactylorhiza majalis (Reichenbach) P. F. Hunt and Summerhayes (1965) occurs in different types of Eurasian wetlands (Delforge 2006). Because such wetlands, e.g., coastal and forested wetlands, marshes, fens, bogs, peatlands, and floodplains, have been gradually removed from the landscape for human uses, D. majalis has become a rare species. Protection of natural refuges from drainage systems and intensive agriculture as well as the restoration of lost wetlands are remarkable efforts to save this taxon at least on a local scale. D. majalis suffers, however, from the conditions of an extensive agricultural management, when competitive flora is not regularly removed at the end of the season (Wotavová et al. 2004).
Due to ongoing habitat destruction and fragmentation (Riecken et al. 2006;Liu et al. 2016), D. majalis is present on many German regional Red Lists of endangered species as has been recently reported for the Upper Lusatia region, situated in East Saxony, and Thuringia (Bräutigam and Otto 2012;Rode and Schadwinkel 2014). To prevent this rare species from further decline, more efficient conservation actions are urgently needed, especially when considering the fact that from 10 % to 30 % of the world-wide D. majalis populations grow in Germany (Ludwig et al. 2007). Despite this importance, relatively little attention has been paid, however, to evaluate the current sexual fitness of German D. majalis stands by carrying-out an exact identification of flowering access and fruit production as proposed by Heinrich (2013).
Regarding the fact that D. majalis increasingly occurs in small relict populations, it could be threatened by a reduction of genetic variability in response to random drift processes and inbreeding. For statistics identifying this disadvantageous situation in plants, connected to homozygosity excesses in terms of Hardy-Weinberg equilibrium (HWE), a failure to adapt to shifting environments, and insufficient fertility, see Yeh (2000). The floral composition of Upper Lusatia is furthermore expected to be affected by forecast increasing temperatures in spring and dry summers due to global climate change (Flemming 2005;Kreienkamp et al. 2011). In order to get a better understanding in how far the sexual reproduction of D. majalis is influenced by the above-mentioned variety of genetic and environmental stresses, local population surveys have a high priority. The knowledge, gained from such studies, may provide new scientific information to implement site-specific conservation strategies in the future.
D. majalis represents a naturally occurring allotetraploid hybrid between diploid Dactylorhiza incarnata (Linné) Soó (1960) and diploid Dactylorhiza fuchsii (Druce) Soó (1960) that unifies 80 tiny roundish chromosomes within a single cellular nucleus. This complex genome combination, corresponding to marked morphological traits, has been originally confirmed by cytological investigations (pointed out by Heslop-Harrison 1953;Hunt and Summerhayes 1965) and was later on characterized using different types of genetic markers (Hedrén 1996;Devos et al. 2003;Devos et al. 2006;Paun et al. 2011;Balao et al. 2016). Allopolyploids are usually more vigorous in comparison to their parental species taking advantage from both gene-dosage and heterosis effects (Comai 2005). They can furthermore readily adapt to changing environments (Paun et al. 2010) and defend ecological niches in disturbed areas as documented by the fact that D. incarnata disappeared over time from Upper Lusatia, while D. majalis is still present (Bräutigam and Otto 2012). As has been reviewed by Comai (2005), allopolyploid hybrids are known to protect their unique gene pool by a strong reproductive isolation, which is indicated by numerous prezygotic barriers (e.g., floral isolation, highly specialized pollination vectors, or pollen-pistil incompatibility) and postzygotic barriers. The latter leads to a mosaic of differential fertility amongst outcrossing offspring, including inviable embryos and endosperm malfunctions, which is caused by altered chromosomal structures and numbers and/or incompatible gene functions with respect to nuclear genomes and organelle DNA variants (for reviewing basic principles of widespread hybridization events within the plant kingdom see Burke and Arnold 2001;Mallet 2005;Baack and Rieseberg 2007;Soltis and Soltis 2009).
In the present investigation, we examined the current sexual reproduction of 12 D. majalis stands in Upper Lusatia by a broad range of methods including environmental assessments, embryo-viability stains, chromosomal records, and genetic inventories of mother plants. Main aims of our study were i) to count the number of growing D. majalis individuals and to assess flowering, fruit-set, and embryo viability within stands; ii) to identify local site conditions and certain genotypes beneficial for sexual reproduction; iii) to find possible genetic loads and environmental factors of diminishing concern to the production of offspring.

Study area and fitness assessment of D. majalis
Twelve orchid stands in the Upper Lusatia region were investigated during the growing season of the year 2014. Note that five stands represent the northern lowland and the remaining seven climb from the southern lowland to the mountains (Fig. 1). The absolute number of all growing D. majalis plants and their flowering access (absolute number of flowering individuals) were counted one-by-one in case of small populations. Extent and flowering intensity of large populations were estimated. Absolute numbers of blossoms and seed-producing fruits were determined per inflorescence and averaged over the total number of investigated sexually reproducing mother plants per stand (up to 25, collecting more was not allowed by the local environmental protection authority). Sampling of clonal groups of mother plants, derived from the same individual by asexual reproduction, was avoided by selecting unnested individuals even within small test populations. Leaves and seeds of the examined mother plants were further analyzed as described below (Tables 1 and 2).

Chromosomal smears
When D. majalis appeared in the field, roughly about a month before flowering, megasporocytes were fixed in a fluid of absolute methyl alcohol and glacial acetic acid (ratio of 3:1 by volume) for at least 1 h. They were then macerated and stained by acetocarmine at 95°C for 10 min. A stock solution was previously prepared by dissolving 1 g carmine in 100 ml of boiling 45 % (v/v) acetic acid, which was then filtered and stored at room temperature until use (Sass 1958). Treated megasporocytes were finally transferred in a fresh drop of acetocarmine and smeared by gently pressing a cover glass against the slide. The obtained chromosomal plates were subjected for microscopic investigation and photographical documentation after removing excess stain by a paper tissue.

Genetic inventories
We selected three target-specific PCR primer pairs from the literature (Nordström and Hedrén 2007), which are known to amplify nuclear-encoded microsatellites (ms) either in diploid D. incarnata (ms3) or diploid D. fuchsii (ms13 and ms14), respectively. Because these microsatellites are unique markers for only one parental genome when below-testing the allotetraploid hybrid D. majalis, easily interpretable codominant PCR product patterns were found which refer to a single Mendelian locus. This allowed the distinction of heterozygous genotypes (i.e., two different gel electrophoresis bands) from homozygous genotypes (i.e., only one band). Individual leaves, harvested per stand from the abovementioned mother plants and dried in silica gels at 6°C, were further treated for DNA extraction (Doyle and Doyle 1990) employing the NucleoSpin Plant II Kit (Macherey-Nagel, Germany). A total of 15 ng of genomic DNA was amplified in a final reaction volume of 25 μl per plant, containing 0.26 μM forward primer (5′ end-labelled with BMN-5 or BMN-6 for fluorescent detection by multiplex analysis), 0.26 μM reverse primer (not labelled), 1.6 mM MgCl 2 , 125 μM of each dNTP, 2.5 μl 10 x PCR buffer, and 0.5 units Taq DNA polymerase (VWR). All primers were provided by biomers. net (Germany). Forty PCR cycles were conducted with the following parameters: denaturation at 94°C for 60 s; annealing at 54°C for 60 s (ms3 and ms14) or alternatively at 55°C for 60 s (ms13); and extension at 72°C for 90 s.
PCR products were finally analyzed by use of the DNA separation capillary array 33 -75B, running under control of a GeXP Genetic Analyzer (both instruments from Beckman Coulter, United States), to detect fluorescent bands by electrophoresis. The operating conditions were in accordance with the executable tools implemented by the manufacturer in the program software Frag-3. Exact genotypes for each studied microsatellite sample were identified during the electrophoretic run with the aid of the fluorescence detector and 22 internal molecular-weight standards, which have been previously l a b e l l e d w i t h We l l R E D ™ f l u o r e s c e n t d y e b y GenomeLAB GeXP Consumables Chemistries (Beckman Coulter) before distributing.

Embryo-viability stains
Fruits, harvested from the above-mentioned mother plants after open pollination, were dried in silica gels at 6°C. To estimate fertility, 400 seeds per plant were assessed as recommended by the International Seed Testing Administration on the basis of a specific protocol (Van Waes and Debergh 1986), optimized for European orchids as follows: seeds were treated for 6 h in a solution, consisting of 5 % (w/v) Ca(OCl) 2 and 2% (v/v) Tween-80, and were then soaked in sterile distilled water for 24 h. They were subsequently placed into 2, 3, 5triphenyltetrazolium chloride solution (1 g substance dissolved in 100 ml sterile distilled water, pH 7.0) and stored in darkness for 24 h at 30°C. Seeds were afterwards rinsed in three changes of sterile distilled water at 5-min intervals, and examined using a stereoscopic microscope. If the embryo showed pink stain it was judged as viable (i.e., fertile seed = 1) with biochemically intact mitochondria, whereas it was judged as unviable (i.e., sterile seed = 0) when it was uncolored or showed reddishyellow colored tissues. Plant-specific seed-fertility rates, lying between 0 and 1, were obtained after dividing the sum of viable seeds by four hundred. Arithmetic means of seed-fertility rates were calculated using all tested mother plants of a given stand.
It is worthy of note that absolute seed fertilities can only be determined by germination experiments, because the relation of the percentage of germinable seeds to the corresponding percentage of pink-colored seeds is generally unknown. Despite this limitation, embryo-viability stains were successfully employed during the past to identify relative seedfertility differences among Dactylorhiza species and their hybrids (De hert et al. 2012).

Environmental assessment
Plant species, forming the habitat communities, were identified by utilizing a current version of Rothmaler's morphologic key (Jäger 2011). For this purpose, a restricted area of 225 m 2 was inspected per location. Each species is known to meet specific environmental conditions within its range of distribution, i.e., soil moisture (F), thermal continentality (K), light availability (L), nitrogen content of the soil (N), pH-reaction of the soil (R), and ambient temperatures (T) as has been documented by Ellenberg's indicator values (Ellenberg et al. 1992).
We summarized all identified plant species at a given stand and calculated unweighted arithmetic means for each of the above-mentioned indicator value. Finally, the degree of correlation was tested between the obtained means of site-specific indicator values and population-specific fitness parameters in the ensuing exploratory statistical data analysis (see below).

Data analysis
The following procedures were conducted in statistics and population genetics: I, pair-wise correlation analyses among fitness variables assessing sexual reproduction and with genetic diversity measures (see below); II, pair-wise exploratory statistical analyses among fitness data and environmental indicator values by estimating Kendall's correlation coefficient (Kendall 1955); III, analyzing the frequency distribution of the present embryo-viability classes over all stands by a histogram; IV, Ward's agglomerative hierarchical clustering method (Ward 1963) and non-hierarchical k-means clustering in conjunction with a scree test to dissect the investigated populations according to their seed fertilities; V, nonparametric H tests (Kruskal and Wallis 1952) comparing non-Gaussian seed-fertility distributions over all 12 stands as a whole and over all stands within a given identified seedfertility cluster; VI, chi-squared independence tests to evaluate associations between seed-fertility clusters and genotypes of each investigated marker, taking together minor genotypes of individual frequencies ≤ 5% (Sham and Curtis 1995); in the case of a significant association, multiple Fisher's exact tests, as discussed by McDonald (2014), were used to identify the microsatellite genotypes responsible for the association; VII, exact Mann-Whitney U (MWU) tests (Mann and Whitney 1947) in the case of two clusters, or H tests in the case of more than two clusters, to assess differences in seed fertilities and indicator values between the genetically distinct fertility clusters previously identified. All above-mentioned analyses were run by the software package R (R Core Team 2017). Confidence intervals for Kendall's were computed by the R package NSM3 (Hollander et al. 2015). Calculating exact MWU tests with corresponding confidence intervals was done by the R package coin (Hothorn et al. 2006).  Using POPGENE version 1.31 (F. C. Yeh, R. Yang and T. Boyle, University of Alberta and Centre for International Forestry Research, Canada, 1999), allele frequencies as well as both single-and two-locus genotype frequencies were computed for all three microsatellite markers and all investigated populations to estimate number of alleles (N A ), number of genotypes (N G ), observed heterozygosity (H o ), expected heterozygosity (H e ), F-statistics (F ST , F IS , and F IT ), genetic distances, and disequilibria measures according to Wright (1978), Nei (1972Nei ( , 1973, and Smouse et al. (1983), respectively. Departures from HWE were tested for each marker and each population by the R package genetics (Warnes 2013). To calculate the diversity for alleles (V A ), as defined by Gregorius (1978Gregorius ( , 1987, we employed the GSED computer software version 3.0 (E. M. Gillet, University of Goettingen, Germany, 2010).
All p values ≤0.05 were considered statistically significant. P values of the statistical procedures given in method II and V were not adjusted for multiple testing due to the exploratory character of the analyses. P values of the statistical operations described by methods I, VI, VII, disequilibria measures, and HWE tests were adjusted for multiple examinations by Bonferroni-Holm correction (Holm 1979).

Fitness variables assessing sexual reproduction
We registered small populations, sparsely settled with several dozens of D. majalis plants, as well as more abundant populations where hundreds and in four cases even more than thousand individuals of this taxon appeared (Tab. 1). Population no. 3 is an exception with 8000 D. majalis plants. Over all stands, flower formation was observed in May for 16.4 % to 89.8 % of the growing individuals, which provided on average 13 to 32 blossoms per inflorescence (Tab. 1). On average, 31.1% through to 69.3 % of the enumerated flowers were found to produce fruits later-on (Tab. 1). It should be noted, however, that for unknown reasons a certain number of the originally detected flowering plants or their seeds disappeared in the nature before the assessment of fruit-set took place. This phenomenon is of major concern to small populations (compare the fourth column vs. the fifth column in Tab. 1). Moreover, extremely high rates of seed sterilities were recorded from all 12 study populations, because only a minority of a total sum of 85,600 analyzed embryos were able to demonstrate viability (ranging from 4.24 % to 16.41 % on average per stand as shown in Tab. 1). We found no significant correlation among the above-mentioned fitness variables, except for the strong correlation (r = 0.965; p = 0.00000364) identified between the absolute numbers of growing and flowering individuals.

Cytological studies illustrating aberrant meiosis
The intimate structure of chromosomes in developing megasporocytes was recorded under the microscope for six D. majalis plants of population no. 9, which were collected in April 2017 and 2018. As a result, some erratic arrangements (i.e., large ring-shaped conjugations of more than two chromosomes, a fragile equatorial zone in germinal cells, a sortingout of single chromosomes, and microcytes possessing the small number of remaining chromosomes) were observed at different time points during both meiotic divisions (Fig. 2). Such irregularities can cause aneuploid chromosome numbers with respect to the produced female gametes and can therefore likely explain the high amount of sterile seeds previously determined during the growing season in 2014 (see above).

Exploratory statistical analysis aggregating fitness data and environments
We compared the above-mentioned fitness variables from 12 D. majalis stands to six environmental indicators. As can be seen by Kendall's correlation coefficient (Tab. 2), positively true (i.e., direct) associations were detected between the percentages of flowering D. majalis individuals and pH-reaction of the soil (R), between the percentages of their seed production and ambient temperature (T), between the percentages of their seed fertility and thermal continentality (K) as well as between the percentages of their seed fertility and nitrogen content of the soil (N). Negatively true (i.e., reciprocal) associations were discovered between the number of flowering D. majalis individuals and thermal continentality (K) as well as between the average number of blossoms per inflorescence and the ambient temperature (T).

Agglomerative cluster analysis detecting two seed-fertility groups in differing environments
Inspecting all plant-specific highly variable embryo viabilities, which displayed a L-shaped non-Gaussian frequency distribution (Fig. 3), the H test revealed significant differences among the stands (chi-squared = 40.516; p = 0.000029). Treating this data by Ward's agglomerative hierarchical clustering method, two distinct fertility groups were identified ( Fig. 4; upper part). Cluster I unified the samples of population nos. 1, 6, 8, 9, 10 and 11 due to their superior seed fertilities (12.88 % -16.41 % on average). Cluster II harbored the seed samples of population nos. 2, 3, 4, 5, 7 and 12, which shared very poor fertilities ranging from 4.24 % to 9.04 %. The marked bipartite classification was strongly supported by the present knee in the k-means graph ( Fig. 4; lower part), which indicated that the calculated residual sum of squares did not substantially decrease when dissecting seed fertility data in greater numbers of clusters. Performing the H test procedure, no significant fertility differences were detectable whether within cluster I (chi-squared = 1.3035; p = 0.9346) nor within cluster II (chi-squared = 8.0358; p = 0.1543). This underpins the similarity of fertility data within each identified cluster. However, significant differences (p < 0.0001) were apparent between both clusters when the seed-fertility data were subjected to the MWU test. Further, the results of additional MWU tests assessing differences in each of the six Ellenberg's indicator values between the orchid stands in seed-fertility cluster I and II can be found in Table 3. As illustrated by the calculated Hodges-Lehmann estimators, the two fertility clusters differed considerably in their thermal continentality, but not in any other assessed environmental indicator. The marked difference in thermal continentality between the two clusters is expressed by the fact that with reference to cluster I enhanced K values of 3.5, 3.6, 3.7, 3.7, 3.6 and 3.5 were registered for the thermal environments of population nos. 1, 6, 8, 9, 10 and 11, respectively. Concerning cluster II, reduced K values of 3.5, 3.0, 3.3, 3.4, 3.4 and 3.2 were on the other hand provided by population nos. 2, 3, 4, 5, 7 and 12, respectively.

DNA marker data demonstrating genetic variability in association with both fertility clusters
Conducting genetic inventories by three nuclear-encoded microsatellites, we recorded a total of 23 alleles and 39 genotypes out of 271 examined mother plants (Tab. 4). Locus ms14 was the most polymorphic, showing a diversity for alleles spanning over all stands from 1.88 to 3.39 with a number of genotypes reaching from 2 through to 9. Values were a little bit lower for markers ms13 and ms3 with a diversity for alleles spanning from 1.88 to 2.90 and from 1.52 to 2.76, respectively, and with a number of genotypes ranging from 3 to 8. No tested allelic or genotypic diversity parameter (N A , N G , V A , H o , and H e ) was found to be in linear correlation with the absolute numbers of growing individuals or the relative seed-fertility rates provided by the populations (data not shown). With respect to the F-statistics (F IS and F IT, Tab. 4), heterozygote deficits reaching from 7% to 13% in the case of loci ms3 and ms14 were detected, whereas locus ms13 revealed an opposite excess of 11% to 18% heterozygotes over all stands. Yet, statistically significant deviations from HWE were not observed except for the homozygosity excesses measured for marker ms14 in population nos. 5 (p = 0.0238) and 9 Embryo viability was dissected in six classes, as shown by the rectangles, ranging from 0.0 (i.e., 0% of the embryos were alive) to 0.6 (i.e., 60 % of the embryos were alive). Heights of the rectangles are equal to the number of mother plants representing each class. Data were obtained from 400 seeds per mother plant using arithmetic means (p = 0.0140), and both the homozygosity and heterozygosity excesses (p = 0.0036) associated with marker ms13 in population no. 5 (Tab. 4). A moderate 6 % among-population genetic differentiation was inferred from all three locus-specific F ST values (Tab. 4). Pair-wise comparisons of Nei's genetic distances, calculated for each marker locus, did not very well reflect the geographic relationships among all 12 studied locations. Specifically, geographically nearby populations (e.g. populations nos. 9 and 10, lying 5 km apart from each other) mostly did not cluster together based on genetic distances, whereas distant populations (e.g. populations no. 3 and 6, lying 57 km apart from each other) often clustered together based on genetic distances. There was, however, an informative allele frequency variation for alleles D and E of marker locus ms14 in a clinal fashion when comparing all populations (i.e., allele D dominated among the majority of mother plant populations bearing reduced embryo viabilities, whereas allele E controlled the more fertile test populations). Resulting from this diagonal line of direction, a remarkable non-random association concerning the Fig. 4 Dendrogram (upper part) displaying the results of a hierarchical agglomerative cluster analysis on the basis of Ward's algorithm. The received seed fertilities of mother plants of 12 D. majalis populations (their identification numbers are given in brackets with reference to Fig.  1; S.dorf = Seifhennersdorf) were merged into two different fertility clusters as indicated by horizontal lines and the Manhattan distance metric attributed to the ycoordinate. The two-cluster solution was underpinned by a scree test (lower part of the figure) as calculated from the residual sum of squares of a nonhierarchical k-means cluster analysis and a variable number of clusters ranging from 1 to 11. A corresponding scree test leaving out the 1-cluster solution was inconclusive. The two clusters of the hierarchical cluster analysis were the same as in the k-means analysis Table 3 Results of exact Mann-Whitney U tests evaluating differences in local site conditions between the two previously identified fertility clusters, based on embryo-viability stains. Ellenberg's indicator values Bonferroni-Holm corrected p values are shown in the first data line. The sizes of the detected differences were assessed by the Hodges-Lehmann estimators (Hodges and Lehmann 1963) given in the second data line, with bootstrapped 95 % confidence intervals (95%-CI) printed in brackets when performing 10,000 replicates. Significant differences of site conditions between fertility groups I and II are marked by a star symbol opposite frequency distribution of the prominent genotypes EE and DD was identified for marker ms14, which exactly corresponds to both seed-fertility clusters mentioned above (for detailed information see Fig. 5 and its legend). Testing multilocus departures from panmictic equilibrium, a non-random union of gametes was exclusively found in the case of population no. 7 with respect to two pairs of loci, namely between locus ms13 and locus ms14 (p = 0.0024) as well as between locus ms3 and locus ms13 (p = 0.0414), respectively.

Differential seed fertilities and the detected meiotic irregularities
In this study, we found that the orchid hybrid D. majalis exhibits a mosaic of differential fertility among the seeds, including a high number of inviable embryos, which corresponds well to the data previously released by other investigators (De hert et al. 2012). This probably resulted from the DD (4)  irregular chromosome pairing, which we detected by examining a few megasporocytes. Considering the large rings of multivalent chromosome pairs recorded (normally only rodshaped bivalents should be observed), frequent translocations among them are argued to remove centromere-linked regions from single chromosomes by unequal crossover. As a consequence of such deletions, the affected chromosomes cannot properly make contact with the spindle machinery during both meiotic divisions. Hence, they will be sorted out as furthermore discovered by our analyses, resulting later-on in aneuploid germ-cells and contributing to a reduced fertility among heterozygous offspring as reviewed by Futuyma (1986). Differences in the heterochromatin content, resulting from repetitive non-coding sequences, with pronounced effects on meiotic recombination have been already assumed to control chromosomal stability during tissue cultivation and plant regeneration (Karp 1994).
In fertile allotetraploid hybrids, two different diploid genomes of exactly the same chromosomal number are combined as outlined by Comai (2005). At meiosis, the chromosomes of both contributors segregate equally but independently into the four germ-cells since homogeneic pairing is enforced between homologous chromosomes of the same origin. This prevents intergenomic recombination and maintains the benefits of heterosis. By contrast, if there is a single unpaired chromosome (aneuploidy), moving to one or to the other cell pole, and/or heterogeneic pairing takes place between both species constituents (this is likely the case regarding the multivalent conjugations of more than two chromosomes seen in Fig. 2), subsequent generations are affected by a reduced fertility. Extraordinary heterogeneic recombination between structurally related chromosomes of both contributors is on the other hand acknowledged to be an Fig. 5 Association plot illustrating the chi-squared independence tests of both fertility clusters (named I and II) with microsatellite marker genotypes ms3, ms13, and ms14. Genotypes with frequencies ≤ 0.05 were pooled together in genotype class "NN". The contingency table with corresponding p value for the chi-squared test is given below each association plot. In the association plot, each cell is represented by a rectangle that has signed height proportional to the signed contribution to the chisquared test statistic, i.e., the difference between observed and expected counts divided by the square root of the expected counts, and width proportional to the square root of the expected counts, so that the area of the box is proportional to the difference in observed and expected frequencies (see Friendly 1992;Meyer et al. 2003). The rectangles in each row are positioned relative to a baseline indicating independence (no difference between observed and expected counts). If the observed genotype frequency of a cell is greater than the expected one, the box rises above the baseline, and falls below otherwise. In addition, signed contributions to the chi-squared test statistic (Pearson residuals) are colored according to the legend given for each subfigure (darkgrey: ≥ 2). For marker ms14, multiple comparisons using exact Fisher's tests are furthermore shown (upper panel, right). Figures in bold indicate p values ≤ 0.05. The plot was created using the R package vcd (Meyer et al. 2006(Meyer et al. , 2017Zeileis et al. 2007) advantageous evolutionary force possibly enabling D. majalis to adapt to changing environments (Hedrén 1996).

Continental climate stress, sexual reproduction and local adaptation
We discovered for D. majalis marked bipartite patterns of seed fertility in correlation with different climatic environments relying on the Ellenberg's thermal continentality indicator value (Ellenberg et al. 1992). In Upper Lusatia, winters are a little bit longer and colder, and summers are a little bit shorter and warmer compared to western Germany. This is due to a substantial continental climate impact in East Saxony (for general information see Weischet and Endlicher 2000;Kondratyev et al. 2002), which is reflected by relatively high thermal continentality indicator values spanning from 3.2 to 3.7 for 11 out of 12 studied populations. Only population no. 3 met a thermal continentality indicator value of 3, which suggests an optimal environment for D. majalis (Ellenberg et al. 1992). In line with the original observation by Ellenberg and coworkers (Ellenberg et al. 1992) that central European stands of D. majalis do not perform well when reproducing under harsh continental conditions, we detected a reciprocal statistical association between the absolute numbers of flowering individuals and increasing thermal continentality. This was supported by the fact that rising ambient temperatures lowered the average number of blossoms per inflorescence. Both phenomena, conflicting with sexual reproduction, were surprisingly compensated by a higher seed fertility in unpleasant environments to avoid local extinction. For comparison, recent investigations previously demonstrated that deceptive orchid species in Central Europe, depending on naive pollinators because not offering nectar, can increase their reproductive success by higher seed numbers in the face of a pollinator limitation . For an exact taxonomic identification of D. majalis pollinators and the floral chemical compounds attracting such insects, see Ostrowiecka et al. (2019) and , respectively, and the literature cited therein.
Other fitness variables provided by our investigation, e.g., high relative fruit-building rates, were in accordance with the results released by Sonkoly et al. (2015), but were opposed to those by Claessens and Kleynen (2011). Since the mean fruitset of deceptive orchids was found to be only 27.7 % (Neiland and Wilcock 1998), whereas we registered in mean from 31.1 % to 69.3 %, it is unlikely that our studied D. majalis populations suffered from limited services of the pollinator insects. Moreover, we found no significant positive correlation between the numbers of flowers and fruits as has been described by Vojtko et al. (2015) for D. majalis. The reason might be the fact that several flowering plants or their seeds disappeared in the nature before the fruit-assessment took placea phenomenon which can likely be explained by plant poaching.
Note that within the superior seed-fertility cluster I, differing from cluster II by enhanced continental climate conditions, four out of six localities represented stands of the southern highlands and one additionally elevated site (Schlegel) was localized in closed neighborhood to them. Compared to the northern lowlands, higher amounts of incoming solar radiation and both larger daily and annual temperature differences between night and day as well as between summer and winter, respectively, have been recently recorded for the southern highlands (Flemming 2005). It has been furthermore assumed that the landscape of the northern lowlands, characterized by numerous lakes and ponds, is able to balance out more effectively ambient temperature differences (Flemming 2005). We did not investigate local site conditions using physical instruments. Yet, the instead applied Ellenberg's indicator values have been previously shown to deliver reliable ratings of the environmental conditions under which plant communities grow up (e.g., Seidling and Fischer 2008;Guarino et al. 2012). Our finding of viability selection among the tested embryos in association with differing thermal continentalities (see below) definitely underlines that climatic changes over space and time affected the study plots. Significant temperature effects on both the male and female organs of flowers have been already reported for dozens of plant species before pollination took place and during the post-pollination stage, but the underlying mechanisms of gene expression are currently only fragmentaryly understood as has been reviewed by Zinn et al. (2010) and Hedhly (2011).

Genetic selection in response to continental climate stress
We were able to detect a genetic selection due to the current climate stress, such that the environmentally controlled bipartite fertility pattern of seed samples exactly corresponded to the opposite frequency distribution of two prominent genotypes at locus ms14 when comparing mother plants. Physiological differentiations among growing individuals are known to result from many kinds of abiotic stresses, revealing the existence of both adapted and unfit offspring, which are able to resist or perish, respectively (several case studies have been compiled by Schubert and Müller-Starck 2006). Apart from the observed clinal-fashioned genotype distribution seen along the identified climate gradient, gene-expression profiles could be useful to demonstrate that marker locus ms14 itself, or a closely linked not yet detected gene, may provide different protein variants responsible for fertility/sterility. Unfortunately, we cannot predict any possible gene function for locus ms14 containing a CTT-repeat motif. The available 287 base-pairs showed no significant similarity when performing a nucleotide sequence alignment in publicly accessible databases (current releases of GenBank, GenomeNet, and EMBL nucleotide data libraries were repeatedly scanned via WWW using the procedures described by Claverie 1997). Differential fertilities among the offspring of plant hybrids are on the other hand frequently documented as a result of both altered chromosomal structures and numbers rather than being a product of disturbed gene expression (for detailed information see Futuyma 1986).
Reflecting the higher frequency of homozygous marker genotype EE within cluster I of beneficial seed fertility and the erratic meiotic processes recorded in megasporocytes, we argue that the germ-cells of mother plants bearing more fertile seeds are probably less affected by aneuploid chromosomal numbers and/or severe mutations of the chromosomal structure (e.g., deletions due to non-reciprocal translocations). The identified homozygote advantage perfectly fits the occurrence of aneuploid chromosomal numbers, since they are expected to have a strong negative impact on fertility in heterozygous constellation as issued by Futuyma (1986). Moreover, there is evidence that imbalance of gene dosage in deletion heterozygotes alters the phenotypic expression dramatically and can be injurious to the organism (Birchler et al. 2005). Aneuploidy or mutations of the chromosomal structure in linkage disequilibrium with the present genotype at microsatellite marker ms14, lying on a separate chromosome, thus appears as a plausible hypothesis to explain our current data. Besides a natural selection among adapted and non-adapted genotypes, random variation in small populations, population subdivision, genetic hitchhiking, and gene flow are the possible causes for linkage disequilibriuma phenomenon characterized by allele frequencies which are not independently distributed at two or more loci. Such a linkage disequilibrium between two nonallelic genes on different chromosomes will only be maintained for the future if both gene products interact in metabolism by epistasis to make final protein products (for detailed information see Manly 1985). By contrast, epigenetic alterations of DNA methylation pattern in flowers have been recently shown to drive environmental adaptation in three closely related allopolyploid orchid species including D. majalis (Paun et al. 2010).
Keeping furthermore in mind that locus ms14 is present within the genome of D. fuchsii and absent within the genome of D. incarnata (Nordström and Hedrén 2007), previously undertaken crosses already suggested that D. fuchsii could control the seed fertility expressed in the natural hybrid D. majalis by an unknown maternal mechanism. Despite a high standard deviation, such crosses offered somewhat more viable seeds when D. fuchsii served as pollen recipient (mother plant) and D. incarnata was employed as pollen donor in comparison to a reverse order of operation (De hert et al. 2012).
D. majalis plants, owing disadvantageous genetic constellations to carry-out successful sexual reproduction under harsh continental climate conditions (e.g., DD or the remaining genotypes detected at locus ms14 with the exception of EE), are nevertheless expected to be retained in the nature for longer times by vegetative propagation via auxiliary bulbs (Feldmann and Finke 2013). Asexual reproduction is, however, an evolutionary dead end and can only postpone extinction of a species to a later time (Luijten et al. 2000). A knowledge of species and genotypes conferring thermal stress tolerance is currently of upmost urgency, as has been pointed out by Hedhly (2011), in order to increase our understanding of how plants cope with fluctuating environmental conditions in the face of global climate change.
What about inbreeding and gene flow?
The few significant excesses of homozygosity detected in terms of HWE for markers ms13 and ms14 do not necessarily imply inbreeding because of the following reasons. I) Affected population nos. 5 and 9 were never concerned at all three studied marker loci. II) One locus, namely ms14, was recognized to sign a non-neutral genomic variation. III) The reduced embryo viabilitya trait likely attributable to possible inbreedingis also known to be an intrinsic genetic burden of allotetraploids hybrids under the condition of heterogeneic chromosomal recombination as pointed out at the start of the discussion. IV) Selfing, an extreme form of inbreeding, is unlikely due to the specific blossom morphology of D. majalis (Schoenichen 1940, experimentally shown by Vojtko et al. 2015) and pollinaria bending (discussed in Ostrowiecka et al. 2019). V) An elevated level of gene flow, counteracting possible inbreeding, will be expected with regard to the fact that the embryo is embedded into a surrounding scaffold of collapsed cork-like cells. This enables tiny seeds, reaching from 0.000002 to 0.000008 g in mass, to sail with the wind over larger distances than the pollinator insects usually move (Schoenichen 1940). In this context, it has been shown that bumblebees, being the main pollinators of Dactylorhiza species, rarely fly distances longer than 2 km between the nest and the flower heads (Walther-Hellwig and Frankl 2000). Analyzing 27 populations of D. majalis ssp. lapponica by 15 plastid and eight nuclear marker loci, Hedrén and co-workers (Hedrén et al. 2018) concluded that long-distance seed dispersal was the driving force to colonize the Baltic island of Gotland after the Ice Age, but this type of gene flow only has a minor impact on the contemporary genetic differentiation.
Considering the fact that diploid D. fuchsii exclusively occured side-by-side with allotetraploid D. majalis at stand no. 3, a possible pollination between both species cannot be excluded and could be in part responsible for the striking low embryo viabilities found there and at stand no. 4, which is situated a few hundred meters away. The offspring of such an interspecies cross is expected to display significantly lower rates of embryo viabilities in comparison to the offspring of pure species pollinations due to postzygotic barriers as has been demonstrated by De hert et al. (2012) with regard to three Dactylorhiza species including D. fuchsii and D. incarnata. Although D. majalis has a life-span of up to 30 -40 years (Harper 1977), it is unlikely that we sampled mother plants resulting from a recent putative pollination between D. fuchsii and D. incarnata. The latter is absent from Upper Lusatia since 1950 (Hardtke and Ihl 2000), whereas D. majalis has been highly abundant in the study area for over a century (Barber 1901). Consequently, it is moreover not plausible that our data are compromised by low fertility seeds representing hidden back-crosses between D. incarnata and D. majalis. Concerning other possible mating partners, no further Dactylorhiza species grew at any of the 12 study sites.
Due to the fact that the recorded genotypic data were in part affected by selection and furthermore indicate imbalances between genetic drift and gene flow, as confirmed by the moderate among-population genetic differentiation, we refrained from the calculation of indirect gene-flow estimates (Slatkin and Barton 1989) on the basis of individual F ST values for all three investigated markers. Such an estimation is generally not applicable under the aforementioned conditions (see Hamrick and Nason 2000). The 6 % among-population differentiation identified for our study area is below the recently calculated mean for Orchidaceae (14.6 % according to Phillips et al. 2012), but is consistent with the data reported for some other orchid species, which are comparably distributed by pollination vectors and minute seeds (Hamrick and Godt 1996;Forrest et al. 2004;Chung et al. 2005). However, it should be mentioned that we only sampled mother plants at a limited regional scale in order to study their sexual reproduction. It is of note that inbreeding shapes the whole genome in contrast to selection (adaptation), which affects single genomic sites. Thus, at least 10 unlinked neutral markers (e.g., allozymes) are in sum necessary to exclude inbreeding in our study populations.

Final conclusions for site-specific conservation management
Counting the available numbers of flowering individuals per stand as well as checking over the numbers of blossoms and fruits per inflorescence as has been proposed by Heinrich (2013), did not deliver reliable information concerning the reproductive success of D. majalis. In the same line, no significant correlation existed between flowering, fruit-set, and seed-fertility. After manufacturing laborious embryo-viability stains, we surprisingly discovered that the outstanding large lowland population no. 3, which seemed to be attached to an optimal environment during the past, had the smallest seedfertility rate among all 12 investigated stands. The consequences of this result are not yet forseeable, but ongoing global warming is expected to provide additional stresses likely interfering with reproduction (for reviewing both local and world-wide effects of global warming see Flemming 2005 as well as Alcamo and Olesen 2012; resulting challenges in persuit of the welfare of orchids are outlined by Gale et al. 2018). Since the extent and effect of potential backcrossing of tetraploid D. majalis and diploid D. fuchsii at stand no. 3 and probably also no. 4 might in part be responsible for the very low embryo viabilities, this topic should thoroughly be investigated in the future. Further, populations no. 8 and no. 10 of the southern highlands, although characterized by superior embryo viabilities, could also be endangered due to their little population sizes (44 and 35 growing individuals, respectively). Fifty flowering individuals are at least necessary to produce a sizeable number of seedlings as has been estimated for the food deceptive orchid species Orchis purpurea by Jacquemyn et al. (2007). Focussing on the remarkable twolocus non-random linkage disequilibria exclusively identified in population no. 7, the data could reflect the relatively small population size, or indicate genetic selection. Even allele frequency differences alone are sufficient to generate linkage disequilibrium (Smouse et al. 1983). Finally, taking into account the alarming facts described above, sexual reproduction of D. majalis should be regularly monitored within all 12 stands in the future, including embryo-viability stains and genetic inventories, to validate the long-term properties of our findings and to further refine monitoring strategies and conservation plans to prevent local extinction.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.