Roots of Cypripedium reginae and Fraxinus nigra were collected between 29 July and 18 August 2019, from six locations in western Newfoundland: Grenfell Campus fen, Corner Brook fen, Humber Gorge, Steady Brook, Humber Village, and Pasadena (Figure 2). The three westernmost sites are closer to Marble Mountain, with more rocky, calcareous soils, whereas the easternmost sites were in forest-surrounded fens and ditches. Samples were labeled with plant name (F. nigra or C. reginae), and presence or absence of the other plant growing within 15 m, yielding four comparison groups: (a) orchids with no ash nearby (ONA), (b) orchids with ash nearby (AO-Orch), (c) ash with no orchids nearby (ANO), and (d) ash with orchids nearby (AO-Ash). A total of 54 root samples (27 each of F. nigra and C. reginae) were collected (each approximately 10 root tips, 2–8 cm per sample, one sample per plant), rinsed with water, and immediately divided, with five roots from each sample preserved in CTAB (cetyl trimethyl ammonium bromide, for DNA extraction) (Gardes and Bruns, 1993) and the other five in FAA (formalin-acetic acid-alcohol, for microscopy) (Supplementary Table 1). All orchid and ash plants survived the sampling, confirmed in site visits the following year.

Root samples were cleared and stained for microscopy (Brundrett et al., 1994). In short, thicker orchid roots were sectioned transversely and longitudinally, and stained using trypan blue (0.01% w/v) in lactophenol for 1 min, followed by mounting using lactoglycerol (1:2:1 lactic acid, glycerol and water). Whole ash roots were cleared with 10% KOH (w/v) at 80°C overnight (~12 h), rinsed, and stained as above. Samples were examined and photographed with a ZeissZ1 microscope with brightfield illumination (Biotron, University of Western Ontario).

Approximately half of the CTAB root samples were randomly assigned to a surface sterilization protocol to denature the DNA of organisms on the outsides of the roots to compare endophytic fungi to non-surface-sterilized samples that included root-associated rhizosphere fungi (Supplementary Table 1). Briefly, 100 mg of 0.5–1 cm segments of root tips were cleaned by vortexing in 0.1 M sodium pyrophosphate (60 s), 100% ethanol (molecular grade; 5 s), 0.5% sodium hypochlorite (freshly diluted household bleach, 60 s), 70% ethanol (60 s), followed by three rinses in sterile molecular grade water. Samples that were not surface-sterilized were only rinsed with sodium pyrophosphate and sterile water. All root tips were plunged into and dabbed on malt extract agar (MEA) with chloramphenicol (100 μg/mL) prior to DNA extraction to assess the presence of culturable fungi on their surfaces. Plates were incubated in the dark at room temperature for 2 weeks to record any fungal growth. Extraction of DNA from root tips of Cypripedium reginae and Fraxinus nigra was conducted using a bead beating protocol from the Quick-DNA ™Plant/Seed Miniprep kit (Zymo Research, Irvine, California, United States). Extracts were quantified with a Thermo Scientific Nanodrop2000 Spectrophotometer and stored at −20°C until PCR amplification.

Each DNA sample was PCR-amplified using primers 5.8S-Fun and ITS4-Fun that amplify the ITS2 region of rDNA from most fungi (Taylor et al., 2016). The variability due to evolutionary divergence found at this locus is, in many fungi, adequate to enable identification of nucleotide sequences to the level of species (Taylor et al., 2016). Both primers were modified for Illumina MiSeq sequencing by including a forward or reverse Illumina adapter, a 4 base pair linker (NNNN), and an 8-nucleotide index barcode that allows sequences to be assigned to sample origin after multiplexing. The PCR mix included Toughmix (Quanta Biosciences) polymerase master mix with 50 × loading dye, primers, and 2 μL of each sample DNA and were run in a hot-start (2 min at 94°C), 30-cycle program of 94°C for 30 s, 55°C for 30 s, and 72°C for 30 s, yielding amplicons of approximately 380 bp in length, including the adapters, barcodes, and linkers. One negative control was used (sterile molecular grade water) and two positive controls of DNA extracts from Saccharomyces cerevisiae (Ascomycota) and Agaricus bisporus (Basidiomycota). Sequences were obtained using a 2 × 300 kit on an Illumina MiSeq sequencer at Robarts Research Institute, London, ON.

Amplicons were initially demultiplexed using a custom BASH script (https://github.com/nweerasu/primer_pull). Forward and reverse FASTQ files for each primer were again demultiplexed following a Python script (https://github.com/ggloor/miseq_bin/blob/master/demultiplex_dada2.pl) to separate samples, where primers and barcode sequences were removed. Sample FASTQs were processed through DADA2 (Callahan et al., 2016), following a combination of scripts including the DADA2 workflow for Big Data: Paired-end (1.4 or later) (https://benjjneb.github.io/dada2/bigdata_paired.html) and a modified DADA2 pipeline for paired-end sequences that provided summary statistics (https://github.com/ggloor/miseq_bin/dada2_workflow_1.4.R). Taxonomic assignments used the UNITE ITS General FASTA release (v 8.3) using singletons set as RefS (Abarenkov et al., 2021). Parameters for filtration, trimming, merging, and chimera checking steps within DADA2 are provided in Supplementary Table 2. ITS2 data were additionally filtered to remove singletons and minimize sample bleeding using both positive controls (S. cerevisiae and A. bisporus) as guides. A minimum read threshold of ≥10 reads and a minimum relative abundance threshold of ≥0.03% in each sample was used to reduce index-hopping, or sample bleeding, of amplicon sequence variants (ASVs) based on our positive and negative controls.

Taxonomy was validated using reference sequences downloaded from BLASTn and included with the ASVs in a neighbor-joining tree built with a MAFFT Online alignment (Kuraku et al., 2013; Katoh et al., 2018). Taxonomic clades were created for heatmaps and network analysis to label ASVs from the same individual. Maximum-Likelihood trees of shared orchid/ash ASVs found in SS samples were created and ASVs were grouped into clades where distances between nodes were <0.01. Non-target sequences were filtered from each ASV table prior to analysis. Nomenclature follows IndexFungorum.org.

The packages ggmap, microeco, ggplot2, ggpubr, rstatix, and vegan were used for analyses and visualization of community similarities and differences within R (v 4.1.1) (Kahle and Wickham, 2013; Wickham, 2016; Kassambara, 2020, 2021; Oksanen et al., 2020; Liu et al., 2021). Alpha diversity analyses were done by calculating all applicable diversity indices: observed richness, Shannon's index H′, Simpson diversity index, inverse Simpson diversity index, and Fisher's alpha index for each group of samples (ONA, ANO, AO-Ash, and AO-Orch), split by surface sterilization (sterilized—SS, not sterilized—NSS), and sampling location (six levels, see Methods). An ANOVA and Duncan's Multiple Range Test (DMRT; agricolae package) (De Mendiburu, 2021) determined whether observed differences of the group means were significant, and which groups were significantly different for each diversity measure.

ANOSIM (vegan) with 999 permutations was used to detect any differences in ASV composition between samples by comparing the hosts (orchid or ash; presence or absence of the other host nearby), differences in surface sterilization, sampling location, and any batch differences between DNA extraction/PCA protocols (An et al., 2019). All comparisons used a Bray-Curtis distance matrix with a Kruskal-Wallis test for group comparisons, with between-group comparisons made using Dunn's Test for multiple comparisons with Benjamini-Hotchberg adjusted p-values.

The significant orchid-ash presence and location factors were visualized using PCoA (Bray-Curtis distance) and NMDS ordinations. The most relevant ASVs are included as vectors for each plot type, with the top 10 ASVs that have the highest correlation to the sample matrix included in the PCoA plots, and the most significant 10 ASVs (p < 0.05) included in the NMDS plots.

Prior to network analysis, 66 ASVs that were found in both orchid and ash roots were clustered into 45 clades (distance <0.01) in a ML tree using MEGA X (Kumar et al., 2018; Stecher et al., 2020). Network analyses were done using a non-parametric Spearman's correlation matrix. The correlation optimization parameter (COR_optimization) in the microeco package was used, where the Random Matrix Theory is applied to select the optimal correlation cutoff for the data set (Deng et al., 2012; Liu et al., 2021). Network files were exported to visualize in Gephi (v 0.9.2). Edge widths were scaled to represent Spearman correlation size, undirected network diameters for node sizes were calculated using a betweenness centrality parameter that emphasizes nodes that are central to the network. Nodes were colored by module group calculated by the “cluster_louvain” algorithm, a community assignment for vertices (nodes or groups of nodes) that maximizes contribution to modularity (Csardi and Nepusz, 2006; Blondel et al., 2008). Modules were then queried through FUNGuild (Nguyen et al., 2016) and FungalTraits (Põlme et al., 2020) to provide functional characteristics for each module, split by the weighted relative abundance of each ASV.

Sections of both orchid and black ash roots showed non-mycorrhizal associations and probable orchid mycorrhizal and arbuscular mycorrhizal structures (Figure 3). Within orchid samples, several examples of collapsed pelotons are visible along with some thin intra- and intercellular hyphae, which are typical of orchid mycorrhizal associations (Figures 3A–D). The black ash sections show possible arbuscular mycorrhizal (Glomeromycota) vesicles and hyphae (Figure 3E) and external hyphae of root-associated fungi (Figure 3F). Most sections showing putative AMF structures were from black ash roots, whereas orchid roots primarily had collapsed peloton structures, with other intra- and intercellular hyphae present less often. No mantle or Hartig net, characteristic of ectomycorrhizae were observed in either root type, nor were clamp connections, diagnostic of Basidiomycota, present on any mounts. No growth of fungi or bacteria was recorded on MEA plates from any of the SS or NSS root tips 2 weeks after stabbing them into the agar prior to DNA extraction; the time in CTAB had rendered all root-associated microbes non-viable, whether surface-sterilized or not (data not shown).

After manual filtration and removal of non-target ASVs, there were 993 ASVs with a total of 630,513 sequences across 51 samples (not including controls or three SS orchid samples with <1,000 read depth) that were used for downstream analyses (Supplementary Table 2). Surface sterilized orchid root samples had lower average read depth (an average of 1,073 sequences/sample) than ash roots (average 22,764 sequences/sample). Sample reads were not rarified since there was biological relevance to the nearly 10-fold difference in sequencing depths of orchid and ash.

Alpha diversities were measured for orchid and ash samples split by SS and NSS roots. Ash roots had greater ASV richness (observed) and diversity (Shannon and Fisher's alpha) than orchid roots in both SS and NSS samples (Figure 4; Supplementary Table 3). Relative proportions of fungi at the order level differed between sample groups after surface-sterilization (SS). Glomerales dominated ash roots after surface-sterilization, particularly in ash with orchids nearby (Supplementary Figure 1). In non-surface-sterilized roots (NSS), most sequences (62.4%) but only 112 ASVs were found in both orchid and ash roots, whereas 37.3% of the sequences and many more (625) ASVs were found only in ash roots, and 0.6% of sequences and 37 ASVs were found only in orchid roots (Supplementary Figure 2A). After surface-sterilization, there was a 30% (232 ASV) reduction of fungal sequences, with the majority (51.2% of sequences and 449 ASVs) found only in SS ash roots, 0.7% of sequences and 21 ASVs found only in SS orchid roots, and 48.1% and 66 ASVs found in both orchid and ash roots (Supplementary Figure 2B).

Mean ASV comparison between groups showed no significant difference between SS and NSS samples (p = 0.48) or sample processing batch effects (p = 0.31), but there were significant differences between orchid and ash treatment groups (p = 1 × 10⁻⁴, ANOSIM R = 0.5153) and sampling location (p = 9 × 10⁻⁴, ANOSIM R = 0.1997). Comparing sample (Bray-Curtis) distances within treatment groups using a Kruskal-Wallis test yielded no significant differences between treatment groups (p = 0.9), but there were significant differences by location (p = 3.7 × 10⁻⁴) (Supplementary Figure 3). Location is likely linked with both orchid and ash treatment groups and ASV community structure since sometimes only one treatment type was found in a location (e.g., Grenfell Campus fen and Corner Brook fen only had ONA samples), and there were noticeable site-specific differences as mentioned in the Sampling section.

A PCoA (based on linear mapping) shows greater within sample diversity in both orchid root treatments (ONA and AO-Orch) than roots of black ash (Figure 5). In comparison to ash roots that had multiple fungal associates forming more similar communities, roots of orchids were dominated by quite different fungi (Supplementary Figure 4).

Surface-sterilized ash roots have a substantially higher read depth than orchids (all of the top 10 correlated ASVs were specific to ash samples) (Figure 5). Notable taxa in ash included the broad-host root endophyte Cadophora orchidicola (Fernando and Currah, 1996), Pyricularia, and Hymenoscyphus epiphyllus, and the AMF Glomus, Rhizophagus irregularis, R. intraradices, and Dominikia (Glomeraceae). Ordination in NMDS (based on rank distribution and non-linear mapping) revealed that community proportions between samples are more similar, with slight variations due to location and orchid/ash presence. Relative to sample differentiation, there were 10 significant (p < 0.05) ASVs, of which seven are Glomus in ash roots, plus Tetracladium maxilliforme, Leptosphaeria, and Alatospora.

Non-surface-sterilized roots of ash had substantially more Dactylonectria macrodidyma and Xylomyces aquaticus than orchids (PCoA) (Supplementary Figure 4A) and were significantly (p < 0.05) dominated by AMF (Glomus spp., Dominikia); other significant taxa without any obvious location or presence-based influences included Amphinema sp., Hymenoscyphus epiphyllus, Ilyonectria radicicola, and Thelephora alnii (NMDS) (Supplementary Figure 4B).

There were 112 ASVs shared between NSS roots of orchid and ash and 66 ASVs shared between SS roots (Supplementary Figure 2); all 66 shared SS ASVs were also present in NSS shared samples (Supplementary Table 4). The top ASV clade (by relative abundance) shared between both SS and NSS orchid and ash roots was Cadophora orchidicola C2 (Figure 6). Root pathogens Dactylonectria macrodidyma, D. pauciseptata, Ilyonectria radicicola and Pyricularia, and saprotroph Tetracladium maxilliforme are also present in high abundance in both SS and NSS roots. Twenty-six of the 66 shared ASVs are arbuscular mycorrhizal Glomus, G. macrocarpum, Rhizophagus intraradices, R. irregularis, with lower read numbers of Claroideoglomus claroideum and Dominikia (Supplementary Table 4). Roots of SS had higher relative abundances of AMF, whereas NSS roots had more root pathogens. Several ectomycorrhizal fungi were also present in SS roots of both Cypripedium and Fraxinus, including two Tomentella sp. clades (C4 and C24), one T. galzinii C4, and Sebacina incrustans C5. Both SS and NSS orchid root samples had single, dominating ASVs or clades (Tomentella C4, Hyaloscypha finlandica, Tomentella C9, Ceratobasidium, Dominikia C1, Leptosphaeria) nearing 75–100% relative abundance in single samples—i.e., orchid roots had more differential, unique taxa. In contrast, ash roots had single ASVs or clades with more equal abundances across samples, without single dominants (Supplementary Figure 5).

A network analysis was performed using all 66 shared SS ASVs, labeled with 38 distinct taxonomic clades (evolutionary distance <0.01) in a ML tree (Figure 7). Overall, the network of fungi found to be shared by orchid and ash roots was composed of nine separate modules, with several peripheral nodes that interconnected six modules, and three entirely disconnected modules. All network members were positively correlated (ρ ≥ + 0.57). The most common nodes were those of arbuscular mycorrhizal clades (Glomus and Rhizophagus), followed by Cadophora orchidicola, and a plant pathogen (Pyricularia). Peripheral nodes with high among-module connectivity (Pi > 0.4) and high within-module connectivity (Zi > 0) included Glomus C6, Alatospora, Rhizophagus C4, R. irregularis C10, Glomus C15, and T. maxilliforme C1 (Supplementary Figure 6). Highly correlated clusters of C. orchidicola (module 2), Glomus sp. 3 SUN 2011 G2 (module 9), Spirosphaera cupreorufescens G1 (module 8) suggest the presence of artifactual or genetic variants of the same organism. Neither FungalTraits nor FUNGuild had matches for Dominikia, and Cadophora orchidicola ASVs had to be searched under its synonym Leptodontidium orchidicola in FungalTraits, with FUNGuild lacking matches to either synonym, but instead matching for Helotiaceae (Supplementary Figure 7).