Molecular and isotopic data were collected from seven orchid species that co-occurred in two coastal dune slacks along the Belgian coast (Ter Yde: 51°8′18″ N – 2°42′14″ E; Westhoek: 51°5′16″ N – 2°33′48″ E). Both sampling sites were devoid of trees and had similar, high light conditions. The sampled orchid species included two Epipactis species (E. neerlandica and E. palustris), three Dactylorhiza species (D. fuchsii, D. incarnata, and D. praetermissa), Herminium monorchis, and Liparis loeselii.

To assess differences in orchid mycorrhizal communities between the seven species, young roots of five replicate individuals were collected for each species in the most central part of each of the study populations in summer 2015. After roots were carefully excavated from the soil, they were transported to the laboratory and immediately surface sterilized (30 s submergence in 1% sodium hypochlorite, followed by three 30 s rinse steps in sterile distilled water). Mycorrhizal association was confirmed by detection of peloton structures using thin cross sections of orchid roots. These sections were gently cleared in hot alkali, stained using Trypan blue, and viewed with interference contrast microscopy (Brundrett et al., 1990). Subsequently, DNA was extracted from 0.5 g mycorrhizal root fragments using the UltraClean Plant DNA Isolation Kit as described by the manufacturer (Mo Bio Laboratories, Inc., Solana Beach, CA, United States) and 10 times diluted afterward. Amplicon libraries were created using the broad-spectrum internal transcribed spacer (ITS) primers ITS3 (5′-GCATCGATGAAGAACGCAG-3′) and ITS4OF (5′-TTACTAGGGGAATCCTTGTT-3′) (White et al., 1990; Taylor and McCormick, 2008). This primer pair has been shown to amplify substantial diversity of sequences and to detect a variety of orchid-associating mycorrhizal families (Waud et al., 2014). All samples (65 in total) were assigned unique MID (Multiplex Identifier) barcode sequences according to the guidelines for 454 GS-FLX Lib-L sequencing.

Polymerase chain reaction (PCR) amplification was performed in duplicate in a 25 μl reaction volume containing 0.15 mM of each dNTP, 0.5 μM of each primer, 1 U Titanium Taq DNA polymerase, 1X Titanium Taq PCR buffer (Clontech Laboratories, Palo Alto, CA, United States), and 1 μl of a 10-times diluted DNA extract. PCR conditions were as follows: initial denaturation of 2 min at 94°C followed by 30 cycles of 45 s at 94°C, 45 s at 59°C, and 45 s at 72°C. After resolving the amplicons by agarose gel electrophoresis, amplicons within the appropriate size range (∼250–500 bp) were cut from the gel and purified using the Qiaquick gel extraction kit (Qiagen, Hamburg, Germany). Purified dsDNA amplicons were quantified using the Qubit fluorometer (Invitrogen) and pooled in equimolar quantities of 1.00E+10 molecules per sample. This yielded two amplicon libraries, each representing one of the two PCR replicates. The quality of the amplicon libraries was assessed using an Agilent Bioanalyzer 2100 and high sensitivity DNA chip (Agilent Technologies, Waldbronn, Germany). Each amplicon library was loaded onto 1/8th of a 454 Pico Titer Plate (PTP). Pyrosequencing was performed using the Roche GS FLX instrument and Titanium chemistry according to the manufacturer’s instructions (Roche Applied Science, Mannheim, Germany).

Samples for isotopic analysis were collected from the same set of orchids at Ter Yde and Westhoek in summer 2015, at the same time that roots were collected for mycorrhizal analyses. At Ter Yde, we sampled Epipactis palustris, E. neerlandica and Dactylorhiza praetermissa, and at Westhoek, in addition to the two first previous species, we sampled Dactylorhiza fuchsii, D. incarnata, H. monorchis, and L. loeselii. At both sites, we also sampled representatives of several surrounding non-orchid autotrophic plants (see names in Figure 4). For each species, we collected five leaves, each from different plants. Samples were ground in 2 mL Eppendorf tubes in a ball mill (MM200, Retsch Gmbh, Haan, Germany) and analyzed for total N concentration, as well as ¹³C/¹²C and ¹⁵N/¹⁴N ratios using an elemental analyser (EA) coupled to a ThermoFinnigan DeltaV Advantage Continuous-Flow Isotope-ratio mass spectrometer, and expressed as δ-values (Hynson et al., 2013). Isotope values were calibrated using internal calibrated standards (EDTA and ammonium oxalate). The standard deviations of the replicated standard samples were 0.029‰ for ¹³C and 0.215‰ for ¹⁵N.

The quality-filtered pyrosequencing data set comprised 556 fungal OTUs (89900 sequences), of which 104 (64732 sequences – 72.0%) were assigned to putatively orchid mycorrhizal OTUs according to Dearnaley et al. (2012) and information from previous studies that detected mycorrhizal fungi from the roots and protocorms of these and related orchid species (Supplementary Table S1). Individual samples contained on average 996 sequences of putatively orchid mycorrhizal OTUs, although there was substantial variation between samples (minimum: 101, maximum: 5489). Representative sequences for each mycorrhizal OTU found in this study were submitted to GenBank under the Accession Numbers KY083558 – KY083678 and MF567576 – MF567600.

The majority of the detected fungi belonged to the Ceratobasidiaceae (8 OTUs, 20671 sequences) and Tulasnellaceae (10 OTUs, 9701 sequences) among rhizoctonia fungi, and ectomycorrhizal fungi from the Thelephoraceae (23 OTUs, 13297 sequences), Sebacinaceae (26 OTUs, 8420 sequences) and Inocybaceae (14 OTUs, 8261 sequences) (Figure 1). Besides, a number of other ectomycorrhizal taxa previously shown to associate with these orchid species were detected, including members of the Cortinariaceae (8 OTUs, 2709 sequences), Tuberaceae (1 OTU, 376 sequences), and Pezizaceae (4 OTUs, 97 sequences). Additionally, a number of ectomycorrhizal or possibly endophytic fungi belonging to the Psathyrellaceae, Tricholomataceae, and Russulaceae were only sporadically observed (Figure 1).

Mycorrhizal communities associating with each species were diverse (Figure 2) (average number of OTUs per species per population: 37.8 ± 7.6, range: 23–48) and all sampled individuals associated with multiple fungi simultaneously. The relative abundances of the fungal genera differed between species (Figure 2). Members of the Tulasnellaceae and Ceratobasidiaceae were retrieved from all sampled species, but only in very low abundances in E. neerlandica and L. loeselii (Figure 2). Tulasnella fungi were the dominant community members of Dactylorhiza species, whereas members of the Ceratobasidiaceae fungi dominated the mycorrhizal community associating with E. palustris (Figure 2), none of which clustered with ectomycorrhizal Ceratobasidiaceae taxa (sensu Veldre et al., 2013). Members of the Inocybaceae reached high relative abundances in L. loeselii and E. neerlandica and to a lesser extent in H. monorchis. Tuber reached high relative abundances in E. neerlandica, but was absent in all other species. Finally, members of the Sebacinaceae reached high relative abundances in E. neerlandica, whereas members of the Thelephoraceae had the highest relative abundances in L. loeselii.

Because Tulasnella strains have very derived ITS sequences (Dearnaley et al., 2012), additional analyses were performed using primers specifically designed for ITS detection of other Tulasnella species, in particular Tulasnella calospora. This strain appeared to be a major fungal associate in L. loeselii (Ouanphanivanh et al., 2008; Illyés et al., 2009) and was sporadically observed in Dactylorhiza (Jacquemyn et al., 2012, 2016b). However, while these primers successfully detected positive controls, the Liparis and Dactylorhiza root samples displayed no detectable strains of the T. calospora clade (Supplementary Figure S1).

The NMDS ordination showed that the investigated orchid species associated with distinctive mycorrhizal communities (Figure 3). The mycorrhizal communities associating with E. neerlandica and L. loeselii were clearly different from the other species, which tended to cluster together in the central part of the plot. The PERMANOVA analysis confirmed that community composition of mycorrhizal fungi differed significantly between the investigated orchid species (R² = 0.34, pseudo-F = 3.81, P < 0.001). Species Indicator Analysis revealed 24 indicator species. Dactylorhiza species were significantly associated with Tulasnella spp. (members of the rhizoctonias). E. neerlandica was significantly associated with ectomycorrhizal Tuber, Inocybe and several Thelephora OTUs, whereas rhizoctonia Ceratobasidium was a significant indicator species for E. palustris. Finally, L. loeselii significantly associated with several ectomycorrhizal Thelephora and Inocybe OTUs.

At both sites, E. neerlandica was significantly more enriched in ¹³C and ¹⁵N than all other species (P < 0.001 for all comparisons), supporting a mixotrophic nutrition (Figure 4). δ¹³C values for the species collected at the two sites did not differ significantly (namely for E. neerlandica, P = 0.084; E. palustris, P = 0.227; and Rubus fruticosus agg., P = 0.342). δ¹⁵N values for E. neerlandica were significantly higher at Westhoek (6.48 ± 0.40; mean ± standard deviation) than at Ter Yde (4.57 ± 0.21; P < 0.001), whereas E. palustris showed lower values at Westhoek (-3.91 ± 0.20) compared to Ter Yde (-2.67 ± 0.22; P < 0.001), while R. fruticosus agg. had similar values at both sites (P = 0.365). Regarding autotrophic references, Prunella vulgaris at Westhoek had very low δ¹³C (-32.01 ± 0.19; Figure 4A) compared to the mean value for other autotrophic species (-28.98 ± 0.75). At Westhoek, all orchids except E. neerlandica were depleted in ¹³C relative to the autotrophic references except P. vulgaris (P < 0.05 for all comparisons). Notable exceptions were D. fuchsii and E. palustris that were not distinguishable from Oenothera biennis (P = 0.616; P = 1; Figure 4A). At Ter Yde, all orchids besides E. neerlandica were either depleted (compared to Eupatorium cannabinum and R. fruticosus agg., P < 0.05 for all comparisons) or enriched in ¹³C (compared to Senecio vulgaris, P < 0.05 for both comparisons; Figure 4B) or showed no difference in ¹³C (compared to Leontodon sp., P > 0.05 for both comparisons). As for ¹⁵N, the orchids outside E. neerlandica had higher δ¹⁵N than most of the autotrophic references (P < 0.05 for all comparisons) at both sites; the only exceptions were D. incarnata and H. monorchis at Westhoek that were not distinguishable from O. biennis and Salix repens for the former (P = 0.121; P = 0.998) and Salix sp. for the latter (P = 1).

Considering total N concentrations, the orchids fell in three distinct groups at Westhoek (Figure 5A): E. neerlandica and L. loeselii had the highest N concentrations (P < 0.001 for all comparisons with other species); D. fuchsii and E. palustris had intermediate N concentrations (P < 0.01 for all comparisons with other species); and all other orchids did not significantly differ from most autotrophic non-orchid species (P > 0.05 for all comparisons with other species except O. biennis). The situation was similar at Ter Yde with E. neerlandica having the highest N concentration (P < 0.001 for all comparisons; Figure 5B), E. palustris having intermediate N concentrations (P < 0.001 for all comparisons), whereas D. praetermissa and autotrophic references had similar N concentrations (P > 0.05 for all comparisons) lower than those of the two former species. Overall, all orchids tended to be enriched in ¹⁵N compared to surrounding autotrophic reference plants and four out of seven had higher total N concentrations, but only E. neerlandica also displayed the ¹³C enrichment characterizing mixotrophy.

Our results showed that different orchid species co-occurring within dune slacks associated with diverse fungal communities. These fungi encompassed true rhizoctonia fungi (sensu Dearnaley et al., 2012), namely Tulasnellaceae and Ceratobasidiaceae (Serendipitaceae were absent), as well as various ectomycorrhizal clades (Thelephoraceae, Sebacinaceae, Inocybaceae, Tuberaceae, and Pezizaceae). The latter are likely supported by surrounding willows (Salix repens), which are ectomycorrhizal hosts to these taxa in dune ecosystems (Geml et al., 2014; Boonen et al., 2015) and which can provide them carbon. Their abundance in terms of reads and OTUs, given that they do not belong to usual contaminant fungi in soil samples, suggests that their presence has a true biological meaning. In E. neerlandica at least, the mycorrhizal status of these taxa can be supported since this species is phylogenetically close to, although ecologically isolated from E. helleborine (Jacquemyn et al., 2017). The latter commonly associates with ectomycorrhizal fungi (Bidartondo et al., 2004; Ogura-Tsujita and Yukawa, 2008; Těšitelová et al., 2012). Moreover, its isotopic values (see below) suggest that it belongs to mixotrophic species that receive organic matter from ectomycorrhizal fungi. In L. loeselii, the absence of rhizoctonia fungi also suggests that ectomycorrhizal fungi are the orchid’s mycorrhizal partners.

Surprisingly, we did not observe a clear separation between orchids harboring rhizoctonia and ectomycorrhizal fungi, but rather a ‘spectral’ continuum of mycorrhizal communities was found, challenging the idea of a strict dichotomy between these two types of orchids. Our molecular analyses showed that all investigated species associated with both rhizoctonia and ectomycorrhizal fungi, although the prevalence of each fungal guild differed among species. At one extreme, E. neerlandica was nearly devoid of rhizoctonia fungi, whereas L. loeselii was very poor in rhizoctonia fungi; the same trend is also supported by analyses of this species at other sites (Waud et al., 2017). This trend contrasts with previous results from fen habitats in Hungary (Illyés et al., 2005, 2009), perhaps because the latter were obtained after in vitro isolation of the fungus: this procedure often screens against slow-growing ectomycorrhizal fungi. At the other extreme, Dactylorhiza spp. predominantly associated with rhizoctonia fungi, which largely supports previous works (e.g., Rasmussen, 1995; Kristiansen et al., 2001; Bidartondo et al., 2004; Shefferson et al., 2008; Bailarote et al., 2012; Jacquemyn et al., 2012, 2016b). E. palustris also showed a preference for rhizoctonia fungi, mainly fungi from the Ceratobasidiaceae and to a lesser extent from the Tulasnellaceae (Bidartondo et al., 2004; Jacquemyn et al., 2016a). An intermediate situation where both fungal guilds were present was found in H. monorchis (whose mycorrhizal fungi were never studied previously to the best of our knowledge; Rasmussen, 1995).

The presence of ectomycorrhizal fungi in rhizoctonia-associated orchids has only been reported sporadically before the rise of high-throughput barcoding methods (e.g., Kristiansen et al., 2001; Shefferson et al., 2007; Jacquemyn et al., 2010, 2011; Liebel et al., 2015; Těšitelová et al., 2015) and it was suggested that their presence should not be neglected in analyses of partner diversity (Selosse et al., 2010). The rise of next-generation sequencing, which significantly enhanced the sensitivity of detection and the ability to handle samples that are colonized by multiple fungi, confirmed this trend (e.g., Jacquemyn et al., 2014, 2015; Oja et al., 2015). Still, the exact interaction of ectomycorrhizal fungi with rhizoctonia-associated orchids, in terms of root (or even cell) colonization, and their functional importance remains unclear. It can be questioned whether they are simple intercellular endophytes or whether they form intracellular pelotons similar to true mycorrhizal fungi in orchids. More detailed microscopic and physiological investigations are now needed to investigate whether ectomycorrhizal fungi in rhizoctonia-associated orchids are simple root commensalists or whether they bring nutrients, as can be expected from their roles in mixotrophic orchid species that exclusively associate with them.

Isotopic abundances showed that, although most orchids were enriched in ¹⁵N abundance and total N, ¹³C enrichment was only present in E. neerlandica, supporting the idea that it was the only mixotrophic orchid species based on these features. This species is the only one that grows at the edges of dune slacks, where the water table is probably lower than in dune slacks, and it was mostly colonized by ectomycorrhizal fungi of the Thelephoraceae, Tuberaceae, and Inocybaceae. Thus, E. neerlandica perfectly fits the mixotrophic syndrome described in the introduction. At the opposite, D. fuchsii, D. incarnata and D. praetermissa, E. palustris and H. monorchis, with total N content and ¹³C enrichment close to that of autotrophic non-orchid plants, fitted the syndrome of autotrophic orchids (see the possible reservations and current debates over the exact autotrophy in this syndrome in Introduction). Orchids that mostly associated with rhizoctonia fungi only tended to have slightly lower ¹³C enrichments as compared to non-orchid plants (discussed below), suggesting that ectomycorrhizal fungi did not contribute massively to their carbon budget. This is especially true for H. monorchis that, in spite of a high abundance of ectomycorrhizal fungi, turned out to be even significantly depleted in 13C as compared to autotrophic non-orchid plants.

The situation was less clear for L. loeselii, which mainly associated with members of the Thelephoraceae and Inocybaceae, with minor and infrequent occurrence of rhizoctonia fungi (see also Waud et al., 2017). Although its total N content was similarly high to that of the mixotrophic E. neerlandica, it was not enriched in ¹³C, and was in this respect similar to the orchids that mostly associated with rhizoctonia fungi. L. loeselii mycorrhizal partners are normally enriched in ¹³C as other ectomycorrhizal fungi (Kohzu et al., 1999), although we did not measure this in the investigated sites, so that no evidence for mixotrophy was found in L. loeselii. Thus, the investigated species set from dune slacks display a disconnection between the presence of ectomycorrhizal fungi and mixotrophy, as estimated from ¹³C and total N enrichment.

Our results further support previous results that have shown that E. palustris is autotrophic (Bidartondo et al., 2004; Hynson et al., 2013; Schiebold et al., 2017). The confirmation of autotrophy in E. palustris contrasts with the mixtrophy usually observed in other Epipactis spp. (Selosse et al., 2004; Bidartondo et al., 2004; Gonneau et al., 2014; see Hynson et al., 2013 for review). The genus Epipactis therefore represents another example of orchid genera that encompass both mixo- and autotrophic species, such as also observed in Cymbidium (Ogura-Tsujita et al., 2012) and Neottia (Těšitelová et al., 2015). These genera are relevant to study the shift to mixotrophy (since autotrophy is likely the ancestral state), and the possibility of reversions from mixo- to autotrophy.

Our data confirm previous reports that have shown the presence of ectomycorrhizal fungi in orchids that primarily associate with rhizoctonia-like fungi. Even if their exact role and influence remain unclear in autotrophic orchids, they may potentially indicate a first step in the evolution of mixotrophy assisted by ectomycorrhizal fungi, and suggest that the evolution to mixotrophy occurs in two consecutive steps. In such a scenario, ectomycorrhizal fungi would first gradually replace rhizoctonia fungi as mycorrhizal partners, without necessarily providing carbon, as suggested by the investigated L. loeselii. Subsequently, their indirect access to substantial nutrient resources from tree photosynthesis may allow the evolution of mixotrophy. However, whether ectomycorrhizal fungi in rhizoctonia-associated orchids represent a true predisposition remains to be assessed in other orchid lineages where mixotrophy has evolved. In this framework, the genus Epipactis represents, among others, an interesting taxon where the mix of species with different trophic modes and different symbionts, rhizoctonia and/or ectomycorrhizal fungi, offers possibilities to test this scenario thanks to phylogenetic reconstructions.