Floral Trait and Mycorrhizal Similarity between an Endangered Orchid and Its Natural Hybrid

: Hybridization can often lead to the formation of novel taxa which can have traits that resemble either or both parental species. Determining the similarity of hybrid traits to parental taxa is particularly important in plant conservation, as hybrids that form between rare and common taxa may more closely resemble a rare parental species, thereby putting the rare parental taxon at further risk of extinction via increased backcrossing and introgression. We investigated the ﬂoral (morphological and chemical) traits and orchid mycorrhizal (OrM) fungal associations of the endangered orchid Orchis patens , its more common sister species O. provincialis , and their natural hybrid O. × fallax in natural sympatric populations. We found that both morphological and chemical ﬂoral traits of O. × fallax are shared by the parents but are more similar to O. patens than O. provincialis . OrM fungi were shared among all three taxa, indicating that the availability of OrM fungi should not represent a barrier to establishment of individuals of any of these taxa. These results suggest that O. × fallax may be able to expand its distribution within a similar niche to O. patens . This highlights the importance of quantifying differences between hybrids and parental taxon in species conservation planning.


Introduction
Hybridization among plant species is common in the wild and is thought to be a driving force in evolution [1,2]. Indeed, plant speciation through hybridization (such as allopolyploidy or homoploid hybrid speciation) is widespread in the plant kingdom [3][4][5]. However, hybridization can also lead to introgression and the loss of unique genetic identities within populations [6,7]. This is particularly true when rare species co-occur and hybridize with more common and widespread related taxa [8,9]. In these cases, hybridization may drive extinction of rarer taxa [10][11][12]. The maintenance of species boundaries depends on reproductive barriers to isolate sexually compatible species which grow sympatrically [13,14]. These barriers can be divided in pre-zygotic, which prevent individuals of different species from mating (e.g., by emitting different scents that attract different pollinators, different flowering phenology, habitat separation or geographic distribution) and post-zygotic barriers, such as hybrid sterility, fruit abortion, production

Plant Material
Flowers of Orchis patens (Figure 1a), O. × fallax (which was determined based on morphology; Figure 1b) and O. provincialis (Figure 1c) have been sampled from two populations, Breccanecca (44.3152 • N-9.35614 • E 255 m asl) and Pieve Ligure (44.37890 • N-9.08749 • E 400 m asl; Genova, Italy), where these taxa grow sympatrically. Breccanecca is characterized by olive groves (hosting both rich grasslands vegetation and semi-dry grassland vegetation), and mesophilous broad-leaved mixed woods on hillslopes. Pieve Ligure is a dry grassland and garrigue habitat and thermophilous Mediterranean woods in which artificial pine forests are also present. Untreated epidermal peels were observed by both transmitted light and epifluorescence using a Leica DM2000 Microscope coupled to a computer-driven DFC 320 camera (Leica Microsystems, Wetzlar, Germany). To detect the presence and autofluorescence of secondary metabolites, samples were examined under a UV filter (340-380 nm) according to Talamond et al. [40]. In addition, fresh tissues were stained with Sudan III to detect total lipids or with 0.001% Neutral Red to verify the presence and to measure secretory activity of the osmophores, structures involved in scent production, as performed by Stern et al. [41].
For SEM analyses, flowers were fixed in the ethanol-based working solution FineFIX (Milestone, Bergamo, Italy), and incubated overnight at 4 • C; fixed samples were then dehydrated through an ascending series of ethanol (70-80-90-100%, 1 h for each passage), following Chieco et al. [42]. Samples were processed with a Critical Point Drier (K850CPD 2M Strumenti S.r.l., Roma, Italy), and mounted on aluminium stubs where they were sputter-coated with a 10 nm gold particle layer. Plant material was finally observed with a VEGA3-Tescan-type LMU microscope (Tescan Orsay Holding, a.s., Brno, Czech Republic), at an accelerating voltage of 20 kV.
Flower parts descriptive terminology was adopted according to Bateman et al. [43] while we followed Akbulut et al. [44] for the description of papillae.

Isolation of Volatile Fraction
Samples of flowers of O. patens (6.15 g), O. provincialis (5.56 g), and O. × fallax (2.91 g) to which octyl octanoate (98%, Sigma-Aldrich, Inc., St. Louis, MI, USA) was added as internal standard, were steam distilled with odour-free water for 3 h. Due to the low yield of essential oil, it is difficult if not impossible to stratify on the aqueous phase obtained from steam distillation. Therefore, solvent extraction is necessary to isolate the essential oil. The distillate was extracted with methylene chloride (3 × 100 mL) (Merck, Darmstadt, Germany), dried over anhydrous sodium sulfate (Sigma-Aldrich, Inc., St. Louis, MI, USA), and concentrated at first with a rotary evaporator and subsequently using a gentle stream of N 2 for successive GC/FID and GC/MS analyses [45].

Fractionation and Alkylthiolation of Alkenes
Each sample was subjected to selective purification process [46] and alkylthiolation reaction [47]. The dimethyl disulfide adducts were identified, and the positions of the methyl sulfide substituents were deduced from the fragmentation pattern.

GC-FID Analysis
The analyses were carried out using a Hewlett Packard model 5980 GC, equipped with Elite-5MS (5% phenyl methyl polysiloxane) capillary column of 30 m × 0.32 mm i.d. and film 0.32 µm thick. It was used as a carrier gas at a flow of 1 mL/min. One µL aliquots of essential oil were manually injected in splitless mode. The oven temperature program included an initial isotherm of 40 • C for 5 min, followed by a temperature ramp to 260 • C at 4 • C/min, and a final isotherm at this temperature for 10 min. Injector and detector temperatures were set at 250 and 280 • C, respectively. The relative amount of each component was calculated based on the corresponding FID peak area without response factor correction. Analyses were performed in triplicate.

GC-MS Analysis
The analyses were carried out using a GC Model 6890 N, coupled to a benchtop MS Agilent 5973 Network, equipped with the same capillary column and following the same chromatographic conditions used for the GC/FID analyses. The carrier gas was He at a constant flow of 1.0 mL/min. The essential oils were diluted before analysis and 1.0 µL was manually injected into the GC system in splitless mode. The ion source temperature was set at 200 • C, while the transfer line was at 300 • C. The acquisition range was 40-500 amu in electron-impact (EI) positive ionization mode using an ionization voltage of 70 eV.

Identification of the Components of the Volatile Fractions
The identification of the volatile oil components was performed by their retention indices (RI) and their mass spectra [48], and by comparison with a NIST database mass spectral library, as well as with data from the available literature [49]. Retention indices were calculated by Elite-5MS capillary columns using an n-alkane series (C 6 -C 35 ) (Sigma-Aldrich, Inc., St. Louis, MI, USA) under the same GC conditions as for samples. The relative amount of each component of the oil was expressed as percent peak area relative to total peak area from GC/FID analyses of the whole extracts. The quantitative data were obtained from GC/FID analyses by an internal standard method and assuming an equal response factor for all detected compounds.

DNA Extraction and Amplification
To detect potential orchid mycorrhizal (OrM) fungi of the three taxa, we extracted genomic DNA from roots with the DNeasy Plant Mini Kit by Qiagen (Hilden, Germany) following manufacturer's instructions. DNA quantity and quality were assessed by spectrophotometry (ND-1000 Spectrophotometer NanoDrop; Thermo Scientific, Wilmington, Germany).
DNA amplification was performed as described in Calevo et al. [36] by means of a seminested PCR. Briefly, the entire ITS region was amplified with the primers ITS1-OFa, ITS1-OFb and ITS4-OF primers (OF), and primers ITS1 and ITS4tul (Tul) specifically designed for OM fungi [50]. The ITS2 region was amplified, from the previously amplified ITS products, with the ITS3mod and ITS4 [51,52] tagged primers. Each sample was amplified in three replicates. PCR products were checked on 1% agarose gel, and replicates were pooled together and purified using the Wizard SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA), following manufacturer's instructions. Purified amplicons were quantified with Qubit 2.0 (Thermo Fisher Scientific, Waltham, MA, USA), by means of the Qubit™ dsDNA BR Assay Kit. Sequencing libraries prepared with equimolar amounts of purified PCR products were paired-end sequenced using the Illumina MiSeq (2 × 250 bp) by IGA Technology Services Srl (Udine, Italy).

Data Analysis and Bioinformatics
Data were analysed as described in Calevo et al. [36]; summarizing, paired-end reads from each library were merged using PEAR v.0.9.2 [53], setting the quality score threshold at 28 and the minimum length of reads after trimming at 200 bp. The assembled reads were then processed with the Quantitative Insights into Microbial Ecology (QIIME) v.1.8 software package [54]. Operational Taxonomic Units (OTUs) were defined with an open referencebased clustering strategy, by means of the USEARCH61 method [55], at 98% similarity; OTUs encompassing less than 10 sequences were removed [52]. Taxonomy assignment and OTU picking were performed by using the full "UNITE + INSD"; dataset database version 7.2 for QIIME as a reference [56,57]; the BLAST algorithm [58] was used for taxonomy assignment as described in Calevo et al. [36]. The OMF OTU representative sequences were submitted to GenBank and registered under the following string of accession numbers: OK335158-OK335175.
Further analyses were focused on tulasnelloid, ceratobasidioid, serendipitoid, and sebacinoid fungi, which are common OrM symbionts in the genus Orchis [34,36]. Their OTU representative sequences were used to create a phylogeny with a maximum likelihood (ML) approach. Sequences were aligned using the program MUSCLE [59] with default conditions for gap opening and gap extension penalty in MEGA v.7.0.26 [60]. ML estimation was performed with RAxML v.8 [61] through 1000 bootstrap replicates [62] using the GTR + GAMMA algorithm on CIPRES Science Gataway [63]. Bootstrap values were mapped on the best tree using the -f option of RAxML and -x 12345 as a random seed. Nodes with a bootstrap value ≥70% were considered as well supported.
All further analyses were performed with R v4.1.1 [64] and R studio v1.4.1103 [65]. OTU table was rarefied to 2359 reads in order to make samples comparable with the rrarefy function from the 'vegan' R package v2.5-7. We calculated OTU frequency as: (number of reads per plant taxa)/(total number of reads), and then we plotted the OTUs on a ternary plot using 'ggtern' v3.3.5 R package [66]. Ternary plots (similarly to Venn-diagrams) assign elements to one or more sets, considering their relative frequency. We used a random forest [67] approach to classify plant taxa basing on their mycorrhizal community, this approach is useful to understand the most influencing (the relative importance being calculated as mean decrease Gini score) OTUs shaping a microbial community. We classified orchids mycorrhizal communities using the randomForest function with the 'randomForest' v4.6-14 R package [68] with a Breiman's algorithm. We used orchid taxa as variables and rarefied OTU table as matrix of predictors.

Microscopy
The longitudinal section of the spurs revealed that the inner surface of Orchis patens and O. × fallax were covered with papillae having different length (Figure 2a  In O. patens, papillae had an expanded base, were conical in shape, and displayed straight radiated cuticular striations. Near the spur entrance, papillae were medium-sized and more or less pyramidal ( Figure 2d). These elements gradually increased in length while approaching to the spur apex, where they were more elongated ( Figure 2g); some swollen elements, defined as uniseriate trichomes by Naczk et al. [70], were occasionally present (Figure 2g, arrow). This complex of characters (type of cuticular striations, their morphology, and the increasing size of the papillae along the spur internal surface) has also been observed in O. × fallax (Figure 2b In the groove's centre of O. patens and O. × fallax, papillae were conical to conicalvilliform, while in the groove margins they were conical-elongated to villiform (Figure 3d,e). In all of the three taxa, towards the labellum leg portion along the groove margins, papillae gradually decreased their length and changed into elongated/conical to conical. In the torso region of O. provincialis, villiform papillae were diffused both in the groove centre and in the margins (Figure 3c In the torso, arms, and legs of the labellum of O. patens, papillae were more flattened to slightly conical. Radial linear to vermiform cuticular striations have also been recorded (Figure 3h,k). In O. × fallax, papillae from these regions were more frequently conicalshaped, sometimes villiform, and showed radial linear to vermiform striations (Figure 3i,l). Orchis provincialis showed more often conical to villiform papillae, with the latter ones as the most diffused type in other zones of the labellum's torso ( Figure 3c); their cuticular striations were radial linear to vermiform at the papillae tip (Figure 3j,m).
Histochemical features and the presence of cytoplasmatic contents have been recorded by LM. All of the size classes of conical papillae inside the spur of O. patens were stained by Neutral Red, indicating osmophoric activity (Figure 4a,e). Under UV excitation, all these structures emitted blue autofluorescent signals (Figure 4b-d). Autofluorescence varied from weak to intense light blue when cytoplasmatic content was present, indicating that polyphenols (occasionally hydroxycinnamic acids) were abundant in the secretion. In addition, a weak affinity to Sudan III was observed (not shown). A similar reaction to all of the treatments has been observed also for the papillae inside the spur of O. × fallax (Figure 4f-k). However, SEM analyses highlighted that in this taxon papillae from the median-apical region of the spur were significantly covered by exudated organic residues (Figure 4h). In O. provincialis, no osmophoric activity has been recorded in spur epidermal cells. Their cytoplasmatic content showed only a weak affinity to Sudan III (Figure 4i). Under UV light, only a weak blue autofluorescent signal coming from the lignified cell walls has been detected (Figure 4m,n).

Chemical Composition of Essential Oils
Steam distillation of inflorescences of O. patens, O. provincialis, and O. × fallax yielded 83.32 mg, 37.46 mg, and 18.78 mg of yellowish oils representing 0.34%, 0.16%, and 0.13% w/w yields on dry plant materials, respectively. GC-MS and GC-FID analyses led to the identification of 31 components, listed in order of their elution and reported as percentages of the total oil. Table 1 shows the results of qualitative and quantitative analyses of essential oil. The main bulk of constituents of volatile fractions were found to be saturated hydrocarbons, accounting for 67.63%, 49.40%, and 62.59% of the total essential oils of O. patens, O. provincialis, and O. × fallax, respectively. This class of compounds was mainly constituted of saturated linear chain hydrocarbons in the range C 11 -C 27 : pentadecane (20.15% in O. patens), octadecane (8.97% in O. × fallax), and pentacosane (8.61% in O. provincialis) were the most abundant. Among saturated hydrocarbons, the odd homologs were the most representative. A series of unsaturated linear chain hydrocarbons ranging from C 23 to C 27 was also identified (13.74% O. patens, 7.33% O. provincialis, 3.72% O. × fallax). The double bond position in these linear isomers was determined by GC/MS after alkylthiolation reaction [47]. The most representative compounds were found to be 9-heptacosene (

Fungal Metabarcoding from Roots
Illumina MiSeq sequencing was performed on roots from the three taxa (O. provincialis, O. × fallax and O. patens). After extensive quality filtering and chimera removal, reads were clustered into 870 OTUs. All the analyses reported in this study are limited to the 18 OTUs assigned to the OrM taxa Tulasnellaceae, Ceratobasidiaceae and Sebacinales (which includes Serendipitaceae and Sebacinaceae).
Most of the 18 mycorrhizal OTUs were present in all the three taxa, with the exception of five of them ( Figure 6A). Only two were uniquely associated with O. patens (OTU524 and OTU431), belonging to Serendipitaceae and Ceratobasidiaceae respectively. Four OTUs were only shared by O. patens and O. × fallax (the ceratobasidioid OTU380 and the tulasnelloid OTU5) or by O. patens and O. provincialis (the serendipitoid OTU565 and the ceratobasidioid OTU 310). The most abundant OTUs (in term of number of reads) were present in the roots of the three taxa and a group of them were almost equally distributed in terms of reads frequency (Figure 6A), revealing a core of OTUs representing a common fungal symbiotic microbiota across the two parental species and their hybrid. Using a random forest approach, we classified orchid taxa based on their mycorrhizal community composition, obtaining a OOB estimate of error rate = 71.43% which reveals a strong similarity among the single orchid mycorrhizal communities. The most frequent OTUs were almost the same equally distributed across taxa ( Figure 6B), as reported by the Mean Decrease Gini score, used to assess the importance of variables (in this case OTUs) in model prediction. OTU28 and OTU35 were assigned to Ceratobasidiaceae and were retrieved in all samples. Tulasnelloid fungi were strongly represented by OTU1 (Tulasnella helicospora cfr), the most abundant OTU; Serendipitaceae abundance was instead more variable among taxa (see Figure 6).

Discussion
This study investigated floral traits and OrM fungal associations between two orchid species and their natural hybrid. The two parental species showed distinct morphological traits in their labella and spur morphology. While some hybrid traits are clearly associated with only one of the parental species, other traits appear as intermediate between the two. In other cases, the hybrid shows qualitative and quantitative traits not detected in the two parental species. Formation of novel traits in hybrids have been shown to be drivers of ecological adaptation to new environmental conditions [71], but further work is required to see if this is the case in our study.
Orchis × fallax showed the presence of osmophoric papillae with polyphenolic secretions inside the spur, indicating that the character is morphologically and functionally similar and has been inherited from O. patens. Likewise, secretory uniseriate trichomes were detected in these two orchids. However, the abundant presence of secreted substances on the papillae surface was detected only in O. × fallax. As proposed by Kowalkowska et al. [72], the presence of polyphenolic secretions on cuticles may enhance pollinator visits through an intensification of the scent production. In O. provincialis, no secretory structure nor activity has been detected in the spur. As highlighted by histochemical tests, the secretory/osmophoric activity is mainly accomplished by labellar papillae. The long spur has the role of visual attraction in the deceit of pollinators by mimicking spurs or other structures typical of rewarding species [69]. As stated by Figueiredo and Pais [74], Stpiczyńska and Matusiewicz [75] and Naczk et al. [70], components being part of the floral fragrance belong to plastoglobuli which can be present in osmophores. The transportation of these products outside secretory tissues can occur through cuticular microchannels [70,[76][77][78][79][80]. A high presence of cytoplasmatic globules and secondary metabolites has been noticed in spur papillae of O. patens, but especially in O. × fallax. This could be at the basis of the abundant secretion recorded by SEM. Some of the cytoplasmatic inclusions, as also confirmed by observations under UV light, could correspond to phenolic bodies already detected in species of the Dactylorhiza maculata complex by Naczk et al. [70]. These authors and Wiśniewska et al. [81] stated that dihydroxyphenolic globules may be involved in fragrance production, and thus we can hypothesize a similar role for cytoplasmatic inclusions that we detected in the spur and labellar papillae of O. patens and O. × fallax. In the hybrid, cuticular ornamentations of both spur and labellar papillae seemed to be intermediate between the two parental species. It has been proposed that cuticular striations could be involved in light reflection, acting as a visual cue for pollinators [70,82]. The hybrid showed intermediate characters also in the case of the distribution and types of labellar papillae. It has been widely reported that osmophores appear, generally, in the shape of conical papillae, and that sometimes are not discernible from other flower features [80,83,84]. By Neutral Red and autofluorescence tests, we observed that osmophoric cells and secreted substances were more frequent in the labella of O. provincialis and in the hybrid than in O. patens.
The scent of the three taxa showed dissimilarities: although the volatile spectra presented comparable percentages of compound classes (Table 1), differences were reported at the single component level. It has been proposed that even subtle fragrance variations could be determinant for floral isolation of different species growing in sympatry and sharing one or more pollinators [71,85]. While two compounds were found shared only by O. provincialis and O. × fallax (nonanoic acid and an unidentified compound), the hybrid did not present compounds uniquely in common with O. patens. On the other hand, three compounds (tetradecane, 9-pentacosene, and 7-heptacosene) were shared uniquely by O. patens and O. provincialis. Four compounds were prerogative of O. patens only (pentadecane, cinnamaldehyde, 7-pentacosene, and 12-pentacosene). Among these, pentadecane was recognized as the predominant compound for this species (20.15%). Two compounds were detected only in O. provincialis (trans-Cinnamic acid and 7-tricosene), while the hybrid did not present unique molecules. Each orchid was also distinguishable for the relative abundance of single compounds. Among the various examples, one of the most notable is oxo-β-ionone, recorded as predominant in O. provincialis (20.66 ± 0.3 %) but present in lower amounts in O. patens and in intermediate quantities in O. × fallax. The occurrence in orchid scent of such cleavage product of β-carotene [86] is here reported for the first time. Results seems to suggest that unique molecules detected in the floral scent could contribute to floral isolation in the parental species. The abundant pentadecane in O. patens could also be involved in pollinator attraction. A further contribution by single components present in relative higher amounts in the three taxa cannot be excluded. However, the high level of hybridization and number of hybrids in the wild (Supplementary Figure S1) and a similar level of fruit set (personal observation J.C. and M.B.) suggest that pre-mating barriers are weak or absent among these taxa. It is known, indeed, that in food deceptive flowers, visual attraction is generally regarded as the primary cue to attract pollinators, while the chemical attraction plays only a secondary role in pollination [87]. While for O. provincialis few pollinators have been recorded such as the Anthophoridae Eucera hungarica and Eucera nigrescens or the Apidae Bombus humilis [88,89], no records were available for O. patens. The recording of the deposition of pollinia into O. patens stigma by Eucera sp. (Supplementary Figure S2) confirms that pre-mating species boundaries are not present for these species except for phenology. Indeed, the two taxa are supposed to have different flowering seasons, with O. provincialis flowering in late March/April and O. patens in late April/May [23]. However, climatic conditions more often cause a later flowering for O. provincialis and earlier for O. patens resulting in partially overlap and potential hybridization. It is important to note that fruit development in orchids is a direct result of pollination and is independent of later fertilization since pollen deposition causes ovary enlargement. [90], Consequently, estimates of fruit production can serve as indirect estimates of plant pollination success [4]. We might be careful, therefore, in considering the role of the hybrid in the distribution of the rare O. patens. The hybrid, indeed, always and regardless of weather conditions, overlap its flowering with both O. provincialis and O. patens (personal observation). Since there is no genetic data on the origin of O. patens fruit set, we cannot rule out the possibility that reproductive success was contributed by either the pollen of parental or hybrid plant and future investigations should verify this hypothesis through a manipulated cross experiment.
The analysis of putative OrM fungi by fungal metabarcoding showed that, as previously hypothesized [36], there seems to be no partitioning of OrM fungal barrier by the three taxa. Several dominant OTUs were abundant (in term of read number) in and shared by all taxa (Supplementary Figure S3). Only two OTUs were exclusive of Orchis patens. It is interesting to note, therefore, that contrary to what Těšitelová et al. observed in Gymnadenia with different ploidy levels [37], ploidy seems not to drive a shift in the OrM fungal assemblages among the sympatric taxa we investigated. Calevo et al., [36] already found that Tulasnella helicospora, the main OrM fungus found in O. patens but also O. provincialis and their hybrid ( Figure 6; Supplementary Figure S4), is capable of inducing germination of seeds of all these taxa. The phylogenetic tree of tulasnelloid fungi (Supplementary Figure  S4) revealed that also OTU5, only shared between O. patens and O. × fallax in fact belongs to the clade of T. helicospora. It is tempting to hypothesized that some "secondary" OrM fungi might play a role in restricting the distribution of O. patens, which is particularly rare compared to the widespread O. provincialis, and we cannot exclude a possible role of the two OTUs found only in O. patens in its distribution and persistence. The frequency of these OTUs in O. patens samples (75% and 50%), however, suggests that they do not play a crucial role in the ecology of this species. Comparing the OTUs sequences obtained from this study with the ones obtained in the previous work by Calevo et al. [36], we notice that not all the OTUs found by these authors were retrieved in the roots of O. patens (or the other orchids we investigated) (Supplementary Figures S4-S6); indeed, while the three tulasnelloid OTUs had been previously found in orchid roots (Supplementary Figure S4), only ceratobasidioid OTUs 28, 35 and 59 (the most abundant in our samples), and the serendipitoid OTUs 524 and 517, had been previously retrieved from orchid roots, suggesting that the others OTUs found by Calevo et al. (2020) [36] are putative soil saprotrophs or plant endophytes ( Supplementary Figures S5 and S6). However, in vitro experiments are required to understand their ecological roles for these orchid taxa.

Conclusions
Our results suggest that different floral micromorphology and chemical composition of essential oils do not differ among taxa in the hybrid zones of Orchis patens and Orchis provicialis and that OrM fungi do not vary among these taxa with different ploidy levels. These results are congruent with what was previously hypothesised by Calevo et al., [36] suggesting that OrM fungi do not limit the distribution of the endangered O. patens. The narrow distribution of O. patens seems, therefore, to be the result of abiotic factors and its lack of adaptability to environmental change. It is still not clear if the increasing abundance of O. × fallax relative to O. patens plays a role in limiting the distribution of O. patens. However, given that they share both floral traits and OrM fungi, it can potentially backcross with O. patens or exploit its niche and expand its range. Future studies should focus on the determination of pollinators range and observation of their behaviour on the parental species and their hybrid on fitness components among these hybridizing Mediterranean orchids to test how hybrids influence the distribution of their parental taxa (and if they constitute a threat to their long-term survival).
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d13110550/s1. Figure S1: Stacked column chart of the percentage of hybrid plants within populations in Liguria, Italy. Populations are listed geographically from East to West. Figure S2: Eucera sp. depositing pollinia on Orchis patens stigma. Figure S3: Relative abundance in term of read number of the main OTUs found in Orchis provincialis, Orchis × fallax and Orchis patens.  Acknowledgments: Authors are very thankful to Laura Negretti (University of Genova) for SEM assistance and Arianna Cassetti (CREA-OF) for microscopy technical support.

Conflicts of Interest:
The authors declare no conflict of interest.