This two-year study was conducted from May to July across three populations each of D. majalis (KA, 2014-2015; SKI, 2015-2016; and SKII, 2016-2017), D. incarnata var. incarnata (ZB, 2014-2015; RO and MR, 2015-2016), and D. fuchsii (BR, CM, and GR, 2015-2016), located in northeastern Poland ( Supplementary Table S1 ; Figure 1 ). Dactylorhiza majalis grows in wet meadows, cohabitating with abundant entomophilous and rewarding plants, comprising approximately 20% of the herb layer. The abundance of D. majalis varied, with around 120–200 (to ca. 1,000) flowering individuals ( Supplementary Table S1 ). The three D. incarnata var. incarnata populations varied in sizes ranging from 35 to 200 flowering plants ( Supplementary Table S1 ). These populations inhabited sedge communities in the Biebrza Valley and Rospuda Valley, with a rewarding plant species cover of about 10%. Dactylorhiza fuchsii was found in open hornbeam forests in the Białowieża Primeval Forest and nearby areas and one population in the Biebrza Valley (84–193 flowering plants), all with a low density of rewarding plants.

Over a two-year period (2014 to 2017), open and hand pollination experiments were conducted on three Dactylorhiza taxa during peak flowering, when most flowers were open. Open pollination treatments (n = 5 inflorescences per population/year, 90–340 flowers per population/year, Supplementary Table S1 ) and four hand pollination experiments were carried out on bagged inflorescences, during the optimum of flowering (almost all freshly flowers on the inflorescences were opened). For the hand pollination treatments, we employed three treatments on nylon mesh bag-enclosed inflorescences: (1) spontaneous self-pollination to assess autogamy (n = 2–5 inflorescences per population/year, 20–154 flowers per population/year, Supplementary Table S1 ), (2) induced self-pollination to determine self-compatibility (using pollen from the same flower) (n = 2–5 inflorescences per population/year, 37–203 flowers per population/year, Supplementary Table S1 ), and (3) geitonogamy (using pollen from another flower on the same plant) (n = 1–5 inflorescences per population/year, 37–203 flowers per population/year, Supplementary Table S1 ). Additionally, (4) induced xenogamy (cross-pollination with pollen from different plants within the same population, approximately 3 m away) was performed (n = 5 inflorescences per population/year, 76–177 flowers per population/year, Supplementary Table S1 ). All flowers for each treated plant were hand pollinated. After each treatment, the inflorescences were quickly rebagged. In total, 518 inflorescences and 9994 flowers were analyzed over two years across all treatments and taxa ( Supplementary Table S1 ).

Pollination treatments were conducted in two stages: initially without division of the inflorescence, followed by targeted treatments on lower, middle, and upper sections of the inflorescence in each population to investigate the effect of flower position on fitness traits (fruit set, seed quantity and quality, germination). Flower allocation to each level was based on near-equal numbers. Statistical analyses of seed variability and the likelihood of in vitro asymbiotic germination to protocorm stage were performed for all hand pollination treatments, both divided and non-divided of inflorescence.

Upon maturity, fruit set, total seed number per fruit, mean viable seed length (mm), seed categories based on embryo shape, and viable seeds germinating into protocorms on various asymbiotic media were measured. Seeds were immediately sown post-collection and evenly spread in Petri dishes marked with a grid (6 × 6 for D. fuchsii and D. majalis, 10 × 11 for D. incarnata var. incarnata). A random selection of 11 or 20 squares (confidence level of 0.02) was used for analysis, with a pseudorandom number generator implemented in the R environment (R 3.6/4.0 software, R Core Team, 2021). Categories included seeds with deformed embryos and well-developed globular embryos (Lemus et al., 2021). Sixty well-formed seeds per capsule were randomly chosen and measured using an OLYMPUS SZX2-ILLT stereomicroscope and MultiScan v.14.02 software.

For in vitro germination, seeds were stratified at 4°C for 12 weeks, surface sterilized in 5% NaOCl with Tween 80, and air-dried before sowing on Malmgren (1996) (D. fuchsii) and FAST (Lindén, 1980) (D. majalis and D. incarnata var. incarnata) media. Germination and protocorm development were assessed under a microscope after 40 days of incubation in the dark at 19–20°C. Total seed germination and protocorm development were assessed under a light microscope.

Given the non-normal distribution of the data, traditional transformations (arcsin, log + 1, square root) were ineffective for ANOVA requirements. Therefore, Kruskal–Wallis tests and Dunn’s post-hoc tests with Bonferroni correction were employed to identify differences in fruit production, seed number per fruit, seed length, well-developed seed frequency per fruit, and in vitro germination rates between control and hand pollination treatments, and among flower positions on the inflorescence for each taxon. Variables included treatment types (control, selfing, cross-pollination), inflorescence levels (lower, middle, upper), and dependent factors (fruit set, seed number, seed length, well-developed seed proportion, in vitro germination rates). R 3.6/4.0 software was used for analysis, and standard errors were provided with mean values (R Core Team, 2021).

The ID coefficient δ was calculated from five reproductive traits (fruit set, seed number, seed length, well-developed seed proportion, in vitro germination), both for entire inflorescences and division of inflorescence on three levels, using formulas by Ågren and Schemske (1993) and Lande and Schemske (1985). The formula δ = 1 − (Ws/Wo) was used when Ws < Wo, and δ = (Wo/Ws) − 1 for Ws > Wo, where Ws represents the average fitness of selfed progeny and Wo represents that of outcrossed progeny. The estimation was carried out using data where Ws denotes the average fitness of selfed progeny from autogamy and geitonogamy treatments and Wo denotes the average fitness of manually outcrossed progeny (Lande and Schemske, 1985; Charlesworth & Charlesworth, 1987). Cumulative δ values were calculated incorporating the correlation among data sets, using the approach by Husband & Schemske (1996): δ = 1 − (Ws_F/Wo_F x Ws_NS/Wo_NS x Ws_SL/Wo_SL x Ws_DS/Wo_DS x Ws_IV/Wo_IV), where F is fruit set; NS is seed number per fruit; SL is seed length, DS is frequency of properly developed seeds per fruit, and IV is proportion in vitro asymbiotic germinated seeds to protocorm. The values of δ range from −1 to 1, where 0 signifies no ID, positive values indicate that the outcrossed offspring outperformed the inbred offspring, and negative values indicate the opposite (Charlesworth & Charlesworth, 1987).

Pollen limitation (PL) was calculated using the Larson and Barrett (2000) formula: PL = 1 – (control fruit set/cross-pollination fruit set), where PL values range from 0 (no pollen limitation) to 1 (maximum pollen limitation), and negative values indicate greater pollen receipt in naturally pollinated flowers (Lázaro and Traveset, 2006, González-Varo et al., 2009).

There were significant differences in fruit set between control pollination and hand pollination treatments in the three Dactylorhiza taxa ( Figure 2 ; Supplementary Table S2 ). In D. majalis, D. incarnata var. incarnata, and D. fuchsii, fruit set was 3 to 4 times lower in control pollination (35.7–44.0%) compared to hand pollination (85.0–90.2% for cross-pollination and 88.9–95.6% for self-pollination, Figure 2 ; Supplementary Table S2 ). No significant decrease in fruit set was observed from the lower to the upper position on the inflorescence for each Dactylorhiza species in both control and hand pollination treatments ( Figure 2 ; Supplementary Table S1 ). Spontaneous autogamy was recorded in less than 1% of all studied flowers across the three Dactylorhiza taxa; hence, it was excluded from further analysis.

The seed number per fruit varied significantly between control and hand pollination treatments, and often between inflorescence levels within each taxon ( Figure 2 ). In D. majalis and D. fuchsii, seed numbers per fruit were higher in control pollination than in both selfing and crossing experiments. In contrast, D. incarnata var. incarnata exhibited a higher number of seeds in outcrossing compared to selfing and control pollination ( Figure 2 , Supplementary Table S2 ). The lower part of the inflorescence produced a higher seed count per fruit than the upper part in control pollination across the three taxa ( Supplementary Table S2 ).

Significant variations in seed length and the number of seeds with well-developed embryos were observed only between the different orchid treatments ( Table S2 ). Seeds from control pollination showed a higher frequency of well-developed embryos than those from selfing and outcrossing in all Dactylorhiza taxa ( Figure 2 ).

No significant differences were observed in in vitro asymbiotic seed germination between control and hand pollination treatments in D. majalis and D. incarnata var. incarnata ( Figure 2 ; Supplementary Table S2 ). However, a significant difference was noted in the proportion of in vitro germinated seeds among all treatments in D. fuchsii ( Figure 2 ). Notably, in vitro asymbiotic seed germination varied significantly between inflorescence levels in control pollination for D. fuchsii and in cross-pollination for D. majalis ( Supplementary Table S2 ). The frequency of in vitro seed germination from the control fruit set was comparable to that of in vitro germinated seeds from selfing in D. fuchsii ( Figure 2 ).

The three Dactylorhiza taxa exhibited varying levels of ID. The cumulative ID at the taxa level was negative and lowest in D. majalis (δ = −1.000), nearly close to zero (δ = 0.065) in D. incarnata var. incarnata, and moderately positive in D. fuchsii (δ = 0.366). Differences in ID were observed between species for all fitness parameters and across inflorescence levels ( Figure 3 ). Dactylorhiza majalis showed negative ID values in three out of four reproductive features studied (seed length, frequency of well-developed embryos in seeds, and proportion of in vitro seed germination), both with and without division of inflorescence into three levels ( Figure 3 ). The lowest ID values were consistently noted at the upper level of the inflorescence in D. majalis.

Conversely, in D. fuchsii, the ID estimate increased up to 0.400 in in vitro asymbiotic seed germination without division of inflorescence and across the three inflorescence levels ( Figure 3 ). Other ID values in D. fuchsii were near zero ( Figure 3 ).

Moderate ID was recorded for seed length, seed number per fruit, and the proportion of seeds with well-developed embryos in D. incarnata var. incarnata ( Figure 3 ). PL was high, reaching 0.481 in D. fuchsii, 0.643 in D. incarnata var. incarnata, and the highest at 0.663 in D. majalis.

The link between mating systems and inbreeding depression (ID) has been well-established in many studies. Smithson (2006) discovered that nectarless orchids tend to exhibit high levels of ID, which may be indicative of reduced self-fertilization (autogamy and/or geitonogamy). In orchids employing food-deceptive strategies, particularly those that are self-compatible, high ID is often observed. This is likely due to the fact that deceived pollinators typically visit only a few flowers among plants within populations, thereby promoting cross-pollination and reducing the likelihood of inbreeding (Ackerman, 1986; Dafni, 1987; Gill, 1989; Nilsson, 1992; Neiland & Wilcock, 1998; Jersáková et al., 2006; Juillet et al., 2006; Kropf & Renner, 2008; Jacquemyn & Brys, 2020). However, it appears that the impact of ID varies across different fitness traits. Smithson (2006) and Sletvold et al. (2012) noted significant ID in food-deceptive Orchidaceae, particularly in the early stages of life history, ranging from fruit set to protocorm development. Our study of Dactylorhiza orchids revealed that ID varied considerably among reproductive traits. In early life stages, such as seed viability, distinct patterns of ID were observed, suggesting that deceptive strategies in these orchids cannot be attributed to a single mechanism for reproductive success. In D. fuchsii, significant ID was evident, with outcrossed progeny demonstrating superior fitness compared to selfed progeny. This shows a strong outcrossing mating system, as described by Husband & Schemske (1996). However, the frequency of in vitro germinated seeds from selfed and control-pollinated flowers was similar, yet lower than that from outcross pollination. It provided only partial support for the hypothesis that inbred progeny harboring deleterious alleles are eliminated during over-fertilization or early seed development.

The observed cumulative δ value in D. fuchsii (0.366) aligns roughly with previous estimates for other species, such as D. sambucina (δ = 0.42, Nilsson, 1980), but is lower than reported by Jersáková et al. (2006), (δ = 0.63) and Juillet et al. (2006), (δ = 0.47). This value was similar to Smithson (2006) findings for nectarless orchids (δ = 0.32 ± 0.05) and lower than that reported by Husband & Schemske (1996) for allogamous angiosperms (δ = 0.53). The remarkably low ID in D. majalis and the absence of ID in D. incarnata var. incarnata raise intriguing questions. In D. majalis, strong outbreeding depression, evidenced by differences in seed length, well-developed embryos, and in vitro seed germination between outcross and self-progeny, was observed. However, no significant differences were noted across three inflorescence levels. Dactylorhiza incarnata var. incarnata exhibited outbreeding depression in seed germination, while hand-selfed flowers produced seeds with similar viability to naturally pollinated flowers. These results suggest that in D. majalis and D. incarnata var. incarnata, pollinator limitation is likely, and ID is negligible, with biparental inbreeding being more common. The negligible ID in D. incarnata var. incarnata (δ = 0.065) may reflect a historical prevalence of selfing, leading to a reduced genetic load and decreased inbreeding depression (Husband & Schemske, 1996). Hedrén & Nordström (2009) proposed that populations of D. incarnata var. incarnata might consist of several inbred lines. This hypothesis is supported by studies from Naczk et al. (2018) and Vallius et al. (2004). The extremely low ID in the allopolyploid D. majalis, potentially exacerbated by within-population relative pollination, suggests reduced fitness due to detrimental alleles at a local geographic scale (Husband et al., 2008). There is evidence suggesting a lack of selection against selfed progeny in this polyploid. It is plausible that inbred progeny harboring deleterious alleles are not eliminated during over-fertilization or early germination phases (Ortiz-Barney & Ackerman, 1999). Dactylorhiza majalis, an allotetrapolyploid relative to its diploid progenitors D. incarnata and D. fuchsii (Pillon et al., 2007), represents a complex taxonomic group that evolved through multiple and independent hybridization events. The genus’s taxonomic complexity is likely attributable to its migration history during glaciations and interglacial periods, as well as episodes of polyploidization. Otto & Whitton (2000) noted that taxa belonging to polyploidy, which are at or near mutation-selection equilibrium, are expected to harbor greater genetic loads than comparable diploids, potentially leading to higher levels of ID in the progenitors (Ronfort, 1999; Pannell et al., 2004).

A mixed mating system was observed in all studied Dactylorhiza taxa (Hedrén & Nordström, 2009; Ostrowiecka et al., 2019). Ostrowiecka et al. (2019) found that pollinator behavior in D. majalis likely promotes geitonogamy, explaining the development of selfed seeds in fruits at different levels of the inflorescence with germination potential similar to that of outcrosses within populations. Our results, particularly concerning the level of fruit set, corroborate previous studies on Dactylorhiza (Neiland & Wilcock, 1998; Vallius, 2000; Claessens and Kleynen, 2011). We did not confirm spontaneous autogamy in the three Dactylorhiza taxa, and our findings affirm that these taxa are highly dependent on pollinators, with a notable presence of pollen limitation, especially in D. majalis. This observation aligns with the common occurrence of pollen limitation in deceptive orchids (Tremblay, 2006; Henneresse et al., 2017). The foraging behavior of deceived pollinators likely results in the transfer of low amounts of pollen per visit, increasing the risk of pollen limitation and reducing the likelihood of receiving pollen from different individuals. Ostrowiecka et al. (2019) reported that mean visitation frequency by Apis mellifera was extremely low and time of visits of A. mellifera on inflorescences lasted from 11-40 s in studied D. majalis. Additionally, our data suggest that pollinator foraging behavior remained consistent across the inflorescence during the flowering season, indicating that flower position on the inflorescence did not significantly impact fitness components such as seed length, well-developed seeds, and in vitro asymbiotic seed germination. However, the observed decrease in seed number per fruit from the lower to the upper level of the inflorescence is typically attributed to the resource hypothesis, suggesting increased competition for resources among pollinated flowers (Devlin, 1989; Berry & Calvo, 1991; Obeso, 1993; Diggle, 1995). This decrease in seed set across the inflorescence has also been noted in Dactylorhiza maculata and Epipactis helleborine, potentially due to limited resources and morphological differences in the upper flowers (Light & MacConaill, 1998; Vallius, 2000). In summary, our study as the first revealed that fitness traits in early life stage in three Dactylorhiza taxa depend significantly on type of mating system, but little influence on flower position on inflorescence. These hypotheses warrant further investigation in the context of the three taxa studied.