1. Fully mycoheterotrophic plants are thought to obtain carbon exclusively from their root-associated fungal partners. The general enrichment of these plants in the heavy isotopes ¹³C and ¹⁵N suggests that fungi are the main nutrient source for these plants. Yet, the majority of studies have targeted mycoheterotrophic plants associated with ectomycorrhizal, orchid mycorrhizal and saprotrophic fungi, while mycoheterotrophic plants living on arbuscular mycorrhizal fungi remain understudied.
2. Here, we sampled 13 species of arbuscular mycorrhizal fully mycoheterotrophic plants from five families and co-occurring autotrophic reference plants growing in forests of tropical South America, tropical South East Asia, and temperate Australasia. We measured stable isotope natural abundances (δ¹³C, δ¹⁵N, δ²H and δ¹⁸O), determined total nitrogen concentrations, and used high-throughput DNA sequencing to characterize the arbuscular mycorrhizal fungal communities associated with the sampled mycoheterotrophic plants.
3. We observed a general enrichment in ¹³C and ¹⁵N isotopes across mycoheterotrophic plant families and geographic regions. We confirm cases where no ¹⁵N enrichment is present, but we show that in general arbuscular mycoheterotrophic plants are enriched in ¹⁵N. Moreover, we demonstrate for the first time that these plants are significantly enriched in ²H but not in ¹⁸O in relation to their autotrophic references. The fungal communities targeted by the mycoheterotrophs mainly consist of Glomeraceae and show strong association with the isotopic signatures and geographic origin of the plants.
4. Synthesis: Our finding enlarges the limited knowledge on the multi-element stable isotopic signatures of mycoheterotrophic plants living on arbuscular mycorrhizal fungi. We show that these plants are enriched in ¹³C and ²H as expected due to their mycoheterotrophic nutrition, and that in general they are also enriched in ¹⁵N, despite some exceptions. Variation in stable isotope signatures is likely influenced by plant taxonomy, geography, and fungal community composition.

Sampling

In total, 17 plots were sampled in French Guiana, Malaysia, Australia, and New Zealand. Within these plots, a total of 13 arbuscular mycorrhizal fully mycoheterotrophic plant species were collected (see details in Table S1), from five different families: Burmanniaceae (Dictyostega orobanchoides, Hexapterella gentianoides, Gymnosiphon breviflorus), Gentianaceae (Voyria aphylla, Voyriella parviflora), Polygalaceae (Epirixanthes cylindrica), Thismiaceae (Thismia clavarioides, T. hillii, T. rodwayi, and an undescribed Thismia species), and Triuridaceae (Sciaphila densiflora, Soridium spruceanum, and Triuris hyalina). Within each plot, samples of the target mycoheterotrophic plant species were taken in five replicates, if available, (resembling five 1 m² subplots). Within each of these subplots at least three reference autotrophic plant species were sampled as proposed by Gebauer & Meyer (2003). For sites in Tasmania, autotrophic reference plants were identified by local expert Mark Wapstra. For all other sites, leaf samples were identified to family level or below by sequencing matK or trnL following the methods described in Gomes et al. (2017). Reference plants that belong to Fabaceae or rely on C4 or CAM photosynthesis were excluded from the analyses, since these taxa have distinct stable isotope signatures of N and C respectively (Hynson et al., 2013). A total of 265 leaf samples obtained from autotrophic understorey species served as reference plants (Table S1). For DNA analysis of arbuscular mycorrhizal fungi, roots of the sampled mycoheterotrophs were rinsed with water and stored in CTAB buffer (cetyltrimethylammonium bromide) at −18°C until further processing.

Analysis of stable isotope abundance

Aboveground samples of the 13 mycoheterotrophic species (n = 78) and autotrophic references (n = 265) were dried at 105°C and ground to a fine powder in a ball mill (Retsch Schwingmühle MM2, Haan, Germany). Relative C and N isotope natural abundances were measured in a dual element analysis mode with an elemental analyser (Carlo Erba Instruments 1108, Milano, Italy) coupled to a continuous flow isotope ratio mass spectrometer (delta S; Finnigan MAT, Bremen, Germany) via a ConFlo III open-split interface (Thermo Fisher Scientific, Bremen, Germany) as described in Bidartondo, Burghardt, Gebauer, Bruns, & Read (2004). Relative H and O isotope natural abundances of the samples were obtained only for a selection of representative species within each family because these measurements require the preparation of five samples in total, and for many species not enough material was available. Sufficient material was obtained for E. cylindrica (Polygalaceae), G. breviflorus (Burmanniaceae) and V. parviflora (Gentianaceae), and reference autotrophic plants (n = 31). Representatives of the families Thismiaceae and Triuridaceae were not measured for H and O because insufficient material was available. The H and O isotope natural abundances were measured with thermal conversion through pyrolysis (HTO; HEKAtech, Wegberg, Germany) coupled to a continuous flow isotope ratio mass spectrometer (delta V advantage; Thermo Fisher Scientific) via a ConFlo IV open-split interface (Thermo Fisher Scientific) as described in Gebauer et al. (2016).

Measured isotope relative abundances are denoted as δ values that were calculated according to the following equation: δ¹³C, δ¹⁵N, δ²H or δ¹⁸O = (R_sample/R_standard − 1) × 1,000 [‰], where R_sample and R_standard are the ratios of heavy to light isotope of the samples and the respective standard. Standard gases were calibrated with respect to international standards (CO₂ vs. V-PDB, N₂ vs. N₂ in air, H₂ and CO vs. V-SMOW) with the reference substances ANU sucrose and NBS19 for carbon isotopes, N1 and N2 for nitrogen isotopes, CH7, V-SMOW and SLAP for hydrogen isotopes and IAEA601 and IAEA602 for O isotopes, all provided by the International Atomic Energy Agency, Vienna, Austria. Reproducibility and accuracy of the C and N isotope abundance measurements were routinely controlled by measuring the laboratory standard acetanilide (Gebauer & Schulze, 1991). For the relative C and N isotope natural abundance analyses, acetanilide was routinely analysed with variable sample weight at least six times within each batch of 50 samples. The maximum variation in δ¹³C and δ¹⁵N both within and between batches was always below 0.2‰. For relative H and O isotope natural abundance analyses, benzoic acid was routinely analysed with variable sample weight at least six times within each batch of 40 samples. The maximum variation in δ²H and δ¹⁸O both within and between batches was always below 4‰ and 0.6‰, respectively.

Total N concentrations in mycoheterotroph and reference plant aboveground samples were calculated from sample weights and peak areas using a six-point calibration curve per sample run based on measurements of the laboratory standard acetanilide with a known N concentration of 10.36% (Gebauer & Schulze, 1991).

Identification of arbuscular mycorrhizal fungi

Fungal DNA was extracted from the CTAB preserved roots with the KingFisher Flex Magnetic Particle Processors (Thermo Fisher Scientific, Waltham, MA, USA) using the NucleoMag 96 Plant Kit (Machery-Nagel Gmbh & Co., Düren, Germany). Because morphological observations suggest that the mycoheterotrophic plants of this study are associated with fungi that belong to arbuscular mycorrhizal fungi (Merckx, 2013 and references therein), we targeted Glomeromycotina fungi. The ITS2 region (approximately 250 bp) was sequenced using the primers fITS7 (Ihrmark et al., 2012) and ITS4 (White, Bruns, Lee, & Taylor, 1990) with a Personal Genome Machine (Ion Torrent; Life Technologies, Guilford, CT, USA) with 850 flows. These primers potentially exclude mycorrhizal fungi within the Mucoromycotina (Ihrmark et al., 2012), yet there is no prior indication that these fungi associate with mycoheterotrophic plants. Raw reads were processed using the UPARSE algorithm (Edgar, 2013) incorporated in USearch v.7. Briefly, raw reads were screened for quality control and trimmed at the first base with a Phred score of Q < 20. Then, reads were dereplicated with singletons and sequences with length < 100 bp removed. The resulting reads were clustered into OTUs at 97% similarity (Blaalid et al., 2013), and taxonomy was assigned by querying against the UNITE + INSD database (released on 10.09.2014). We only retained fungal OTUs classified as Glomeromycotina in the subsequent analysis. The OTUs that were represented by less than five reads in each sample were excluded to avoid spurious OTUs (Lindahl et al., 2013). The fungal OTUs were aligned together with sequences from 40 reference Glomeromycotina taxa from Krüger, Krüger, Walker, Stockinger, & Schüßler (2012), and the MaarjAM database (Opik et al., 2010) using the MUSCLE algorithm (Edgar, 2004). The phylogenetic relationships were inferred with RAxML v.8.2.12 using the GTRCAT model of evolution (Stamatakis, 2014). Paraglomeraceae and Archaeosporaceae reference sequences were selected as outgroups and the phylogeny was transformed into an ultrametric tree with TreePL (Smith & O’Meara, 2012) fixing the root at 505 million years ago (Davison et al., 2015). Since the ITS2 region is likely variable below species level in Glomeromycotina, we used the phylogeny to estimate evolutionarily independent clusters of DNA sequences by differentiating branches within clusters, associated with neutral coalescent (merging of lineages within a population) processes, from branches between clusters, associated with speciation events (Powell, Monaghan, Öpik, & Rillig, 2011), using  a single-threshold Generalized Mixed Yule Coalescent (GMYC) approach (Pons et al., 2006) implemented in the splits R package (Ezard, Fujisawa, & Barraclough, 2009). The GMYC table is available in Table S2.