1. Introduction

Mycorrhizal symbiosis is an ancient and prevalent interaction between plants and mycorrhizal fungi (Smith and Read, 2008; van der Heijden et al., 2015; Brundrett and Tedersoo, 2018). The pervasive influence of mycorrhizal fungi on plant physiology, population dynamics, soil processes and ecosystem functioning has been extensively explored and well recognized from plant-centric perspectives (van der Heijden et al., 1998; Powell and Rillig, 2018; Tedersoo et al., 2020). Mycorrhizal fungi affect the productivity, diversity, distribution and coexistence of the aboveground plant community through soil niche partitioning and nutrient exchange in belowground hyphal networks (van der Heijden et al., 2015; Tedersoo et al., 2020). Yet, how plant phylogenetic, physiological and eco-geographical factors influence the assembly and composition of mycorrhizal fungal communities remains relatively unknown.

The Orchidaceae is unarguably one of the largest plant families on Earth (Dressler and Rasmussen, 1996; Wang et al., 2024), with a cosmopolitan distribution (Pridgeon et al., 2009; Givnish et al., 2016; Wang et al., 2024), diverse ecological niches (Chomicki et al., 2015; Givnish et al., 2015; Pérez-Escobar et al., 2017), and contrasting physiological traits (Merckx, 2013; Jacquemyn and Merckx, 2019; Wang et al., 2021). Orchid species have obligate nutritional dependence on mycorrhizal fungi at least at the seed germination stage and associate with diverse fungal groups in the Basidiomycota and Ascomycota phyla (Dearnaley et al., 2012; Rasmussen et al., 2015; Shao et al., 2022). These fungi exhibit high taxonomic diversity, fulfill diverse ecological roles, and include both 'rhizoctonia' fungi represented by three families (Ceratobasidiaceae, Tulasnellaceae, and Serendipitaceae) and a variety of ectomycorrhizal and wood-/litter-decaying saprotrophic clades (Dearnaley et al., 2012; Li et al., 2021; Ogura-Tsujita et al., 2021; Wang et al., 2021). Together, these features make orchid mycorrhizas an ideal study system to explore the drivers of fungal diversity and community composition in plant roots.

The factors determining mycorrhizal community composition within orchid roots remain still largely unclear and their importance likely varies across taxonomic and geographical scales (Jacquemyn et al., 2017). Several studies have demonstrated that phylogenetic relationships may exert a significant influence on fungal community composition, with closely related orchid species frequently hosting similar fungal assemblages (Shefferson et al., 2007, 2010; Jacquemyn et al., 2011). Whether these findings hold across larger taxonomic scales and whether phylogeny plays an important role in fungal assembly at both sides of the interaction remains unclear, as only a well-resolved phylogeny of the plant side is available (Chase et al., 2015; Givnish et al., 2015; Li et al., 2019; Wang et al., 2021).

Plant ecophysiology is also expected to affect the assembly of fungal communities in orchid roots. The orchid family displays a continuum of trophic modes from autotrophy via partial to full mycoheterotrophy (Merckx, 2013; Jacquemyn and Merckx, 2019; Wang et al., 2021). Similar to arbuscular mycorrhizal interactions (Gomes et al., 2017; Zhao et al., 2021), orchid species with different trophic modes seem to have divergent preferences towards fungal partners (Ogura-Tsujita et al., 2021; Wang et al., 2021). Achlorophyllous, fully mycoheterotrophic orchids tend to associate with non-rhizoctonia ectomycorrhizal or saprotrophic fungi, whereas autotrophic species mostly associate with rhizoctonia fungi (Ogura-Tsujita et al., 2012; Selosse et al., 2004; Yagame et al., 2016; Wang et al., 2021). Partially mycoheterotrophic orchids, positioned between the two endpoints of autotrophy and fully mycoheterotrophy, can associate with both rhizoctonia and non-rhizoctonia fungi (Bidartondo et al., 2004; Jacquemyn et al., 2021; Suetsugu et al., 2021a, 2021b). Previous research has found that divergent fungal lifestyles are associated with orchid trophic modes, however, the taxonomic identity and composition of fungal groups were not accounted for (Wang et al., 2021). Thus, orchid partner preference at the fungal community level has yet to be adequately explored.

The composition and specificity of root fungal communities also vary between orchid species with distinct biogeographic distributions (Jacquemyn et al., 2017; Xing et al., 2020). For instance, studies have shown that fungal specificity is higher in most Australian orchids, e.g., Caladenia (Swarts et al., 2010; Phillips et al., 2016; Reiter et al., 2020), Drakaea (Phillips et al., 2011), Pheladenia (Davis et al., 2015), than in Eurasian and North American terrestrial orchids. One potential explanation for this finding is that the relatively old landscapes of Australia may allow for orchid specialization on a small range of mycorrhizal fungi that are well-adapted to local habitats (Phillips et al., 2011). Accordingly, considerable effort has been dedicated to uncovering the effect of habitat conditions, i.e., orchid ecological niche, on their fungal communities (Li et al., 2021). Previous findings at local and regional scales have suggested that fungi associated with orchids are potentially limited by abiotic habitat filtering, suggesting that a wide range of climatic, edaphic, and vegetative variables may influence fungal diversity and community dynamics (McCormick et al., 2009, 2012; Martos et al., 2012; Jacquemyn et al., 2016b; Duffy et al., 2019; Bell et al., 2020; Kaur et al., 2021).

Although previous studies have investigated the effect of orchid phylogenetic relatedness, trophic mode, biogeographic distribution and habitat condition on fungal community composition in orchid roots, their impacts have only been examined at a small phylogenetic and ecological range of orchid species. Thus, it remains to be explored how orchid mycorrhizal communities are structured at a large phylogenetic and biogeographic scale. In this study, we aim to investigate patterns of fungal community composition across the entire Orchidaceae family, using a comprehensive orchid-fungal dataset previously compiled from ecological studies on orchid-fungal associations across the globe (Wang et al., 2021). Using phylogenetic and ecological methods, we evaluate the effect of orchids' phylogenetic relatedness, trophic mode, biogeographic region, ecological biome and corresponding environmental factors on the phylogenetic diversity and community structure of orchid mycorrhizal fungi.

2. Materials and methods 2.1. Orchid-fungus associations and phylogenetic relationships

The orchid-fungal association dataset, compiled in a previous study (Wang et al., 2021), includes 8860 fungal ITS sequences obtained from 50 countries and regions across the globe and covers 750 orchid species spanning the major phylogenetic clades within the Orchidaceae family. Here, we use an ITS-based phylogeny of orchid fungal associates reconstructed in our earlier work (Wang et al., 2022). Because the phylogenetic position of orders within Basidiomycota and Ascomycota is relatively robust (Zhao et al., 2017; Tedersoo et al., 2018; He et al., 2019), we inferred separate order-level phylogenies for fungal clades within the two phyla. Briefly, within each fungal order, we first performed OTU clustering based on 97% sequence similarity using USEARCH v.11 (Edgar, 2010), retrieved taxonomic information for OTUs by blasting against the UNITE reference database (UNITE Community, 2019), built alignments of OTUs and UNITE reference sequences with MUSCLE v.3.8.425 (Edgar, 2004), and reconstructed the Maximum likelihood (ML) tree for each order using RAxML v.8.2.12 (Stamatakis, 2014). After divergence time estimation using a penalized likelihood approach with treePL v.1.0 (Smith and O'Meara, 2012), time-calibrated trees of fungal orders were combined as a single phylogeny according to a backbone phylogeny of fungal orders extracted from available phylogenetic studies (Kohler et al., 2015; Zhao et al., 2017; Tedersoo et al., 2018; He et al., 2019; Mao and Wang, 2019). Datasets of fungal ITS sequences, OTUs, alignments, phylogenetic trees and related metadata are provided in Wang et al. (2022).

In this study, we only included fungal OTUs that were (1) identified as potential mycorrhizal fungi and (2) annotated with precise orchid host information. Although orchid roots host diverse fungal communities, many fungi detected via molecular methods can be endophytes, saprotrophs, or incidental colonizers without mycorrhizal function. To avoid overestimating the diversity of 'true' orchid mycorrhiza fungi, we adopted a conservative approach by retaining only fungal groups with experimentally confirmed mycorrhizal capability (e.g., those demonstrated to promote seed germination, support seedling development, or form peloton structures; Rasmussen et al., 2015). In accordance with previous studies (Dearnaley et al., 2012; Ogura-Tsujita et al., 2021; Wang et al., 2021), we considered the following functional groups as orchid mycorrhizal fungi: 'rhizoctonia' fungi (Ceratobasidiaceae, Tulasnellaceae, and Serendipitaceae), ectomycorrhizal fungi (e.g., Sebacinaceae, Russulaceae, Thelephoraceae, Tuberaceae, Bankeraceae), and wood-/littler-decaying saprotrophic fungi (e.g., Psathyrellaceae, Meruliaceae, Hydnodontaceae, Hyaloscyphaceae). To account for biases in sampling frequency (see Fig. S1) and sequencing methods across studied orchid species (Wang et al., 2021), we created an orchid-fungal association matrix based on presence-absence information of fungal OTUs for each orchid species. The final dataset comprises 836 fungal OTUs (spanning 11 orders) associating with 365 orchid species across all five subfamilies, representing 1481 possible interactions between orchids and fungi (0.5% of all possible links; Table S1). The dataset has been reduced primarily due to frequent gaps in host-species information from original publications and NCBI accessions. Both fungal and orchid phylogenies were incorporated: the phylogenetic relationships of fungal OTUs were extracted from the previously reconstructed fungal phylogeny (Wang et al., 2022), while the phylogenetic relationship of those orchid species were sourced from Wang et al. (2021).

2.2. Orchid trophic mode, subfamily, biogeographic region and biome

Each orchid species was assigned to one category of trophic mode, subfamily, biogeographic region and biome (Table S2). The classification of trophic mode (autotrophic, partially mycoheterotrophic, and fully mycoheterotrophic) was based on morphological traits and stable isotope signatures of orchid species (Wang et al., 2021). Achlorophyllous species lacking visible chlorophyll were labeled as fully mycoheterotrophic. Green-leaved species were classified as autotrophic unless the ¹³C and ¹⁵N isotope signatures suggested they were partially mycoheterotrophic in certain habitats (Gebauer and Meyer, 2003). Notably, our classification strategy for trophic mode is probably conservative because an increasing number of putatively autotrophic orchids were found to be partially mycoheterotrophic with the adoption of the ²H enrichment factor (Gebauer et al., 2016; Schiebold et al., 2018; Schweiger et al., 2019).

We assigned biogeographic region and biome for each orchid species based on available descriptions in eFloras (http://www.efloras.org/), World Checklist of Selected Plant Families (WCSP, http://wcsp.science.kew.org) and the global biodiversity information facility (GBIF, http://www.gbif.org). According to previous research on orchid biogeography (Givnish et al., 2016), the biogeographic region of orchid species was classified into the following: Eurasia, North America, South America (including Central America), Southeast Asia, Australia (including Pacific islands), Africa, and 'others' if an orchid species occurs in two or multiple biogeographic regions. The ecological biome of orchid species was assigned with one of four major types: "Boreal Forests/Taiga", "Mediterranean Forests, Woodlands & Scrub", "Temperate Forests and Grasslands" and "Tropical and Subtropical Forests and Grasslands" according to the RESOLVE Ecoregions dataset (http://ecoregions2017.appspot.com). Additional descriptions for the classification of biogeographic regions and biomes can be found in Appendix S1.

2.3. Orchid occurrence records and environmental niche variables

To assess the environmental niche of each orchid species, we downloaded all available occurrence data of orchid species recorded in the orchid–fungus association matrix from GBIF (GBIF.org, 2020). We discarded records of occurrence with missing coordinates in longitude and latitude or with coordinate precision beyond 10 km. After removing duplicated records for each orchid species, we discarded species occurrences with less than 10 records (Raes and Aguirre-Gutiérrez, 2018), resulting in a final dataset of 240 orchid species for niche analyses (Table S3). The occurrence data were processed with R packages "rgbif" (Chamberlain, 2019) and "CoordinateCleaner" (Zizka et al., 2019).

The environmental raster layers used in this study included variables related to climate, soil and vegetation. A full list of the variable descriptions of each raster layer can be found in Appendix S1. The 19 bioclimatic variables were retrieved from WorldClim online database (http://worldclim.org/version2) with a spatial resolution of 5 arc-minutes. Edaphic variables of top soil (0–30 cm) were extracted from the Harmonized World Soil Database (Batjes, 2009) with a spatial resolution of 1 km. Eight continuous variables of land cover representing broad vegetation types were obtained from the Global 1-km Consensus Land Cover (http://www.earthenv.org/landcover). Raster layers of top soil and land cover were resampled to the spatial resolution of bioclimatic variables using the nearest neighbor sampling method in ArcGIS Desktop 10. Subsequently, all raster layers with the same spatial resolution were stacked together, from which the niche variables were extracted for species occurrences. Because species occurrences are not equally recorded in GBIF, we calculated mean niche variables to represent the common environmental niche space for each orchid species following Hendrix and Vos (2019). The mean niche variables (Table S3) were further standardized for each species using the "dudi.pca" and "niche" functions with the "ade4" R package (Dray and Dufour, 2007).

2.4. Fungal community analyses 2.4.1. General patterns of orchid-fungal associations

To investigate the phylogenetic alpha diversity of fungal OTUs, we calculated the phylogenetic species variability (PSV) for each orchid species. The PSV index summarizes diversity, taking into account the degree of phylogenetic relatedness among fungal OTUs (Helmus et al., 2007). The PSV index is one when all fungal OTUs of one orchid species are phylogenetically unrelated (i.e., a star phylogeny), and approaches zero as fungal OTUs become more related. If not specifically mentioned, the "picante" R package (Kembel et al., 2010) was used for all community analyses.

To access the phylogenetic distribution (clustering vs. overdispersion) of fungal OTUs, we computed the net relatedness index (NRI) and the nearest taxon index (NTI) (Webb, 2000; Webb et al., 2002). Both indices were used to examine if fungal partners of orchid species were more phylogenetically related than expected by chance. NRI quantifies the phylogenetic structure of a species set based on the mean pairwise distances, whereas NTI measures the terminal structure of the species set by computing the mean phylogenetic distance to the nearest taxon of every species. Significant positive or negative NRI and NTI values indicate phylogenetic clustering or overdispersion, respectively.

To measure the community structure of fungal OTUs between groups (beta-diversity), we used the unweighted UniFrac metric (Lozupone and Knight, 2005) to include the phylogenetic relatedness between fungal OTUs. UniFrac distance is the fraction of the phylogenetic distances not shared between two samples, with larger values indicating a greater phylogenetic difference of fungal OTUs between two orchid species. The mean UniFrac distance of one species was calculated by averaging the distances with all other species.

Differences in PSV, NRI, NTI, and mean UniFrac values were assessed among all groups (orchid trophic modes, subfamilies, biogeographic regions, and biomes) by using non-parametric Kruskal–Wallis rank-sum tests and pairwise Mann–Whitney U tests. Differences in fungal community structure (UniFrac metric) were tested among groups using multivariate analysis of variance (MANOVA) with the "rrpp" R package (Collyer and Adams, 2018). The pattern of fungal community structure was visualized using the non-metric multidimensional scaling (NMDS) with the "vegan" R package (Oksanen et al., 2019). In addition, we performed variation partitioning on fungal communities using the 'varpart' function of "vegan". Variation partitioning allows us to investigate the contribution of the four groups (trophic mode, subfamily, biogeographic region, and biome) to the total variation in the dataset.

2.4.2. Phylogenetic relatedness

To test whether the phylogenetic relatedness of orchid species influences their fungal community composition, we computed Mantel tests between the orchid phylogenetic distance metric and the fungal UniFrac metric. The phylogenetic distances of orchid species were extracted from the orchid phylogeny, and the dissimilarity matrix was calculated using the function "cophenetic.phylo" using the "ape" R package (Paradis and Schliep, 2019). Mantel tests were computed using the "vegan" R package with 999 permutations. To visualize the general pattern of orchid-fungal associations, the presence of associations (links) was mapped onto the tangled phylogenies of orchid species and fungal OTUs using the function "cophyloplot" with "ape".

2.4.3. Niche variables and correlated effects with trophic modes, biogeographic regions and biomes

To explore and visualize the correlations between categorical groups (trophic modes, biogeographic regions, and biomes) and environmental niche variables, orchid species with a complete set of information were projected on a two-dimensional trait space using principal component analyses (PCA). Furthermore, we used two-dimensional kernel density estimation (Wand et al., 2015) to visualize the occurrence probability (0.2, 0.5, 0.95, and 0.99 quantiles) of orchid species. Kernel density estimation was implemented with the 'kde' function in the R package 'ks' (Duong, 2007).

The marginal effects (without considering the effects of other predictors) of niche variables on fungal community composition (the UniFrac metric) were tested using distance-based redundancy analysis (db-RDA). Niche variables with a Pearson's R correlation coefficients above 0.7 were identified using the R package "corrplot" (Wei, 2013). The non-correlated niche variables that contributed most to the first two axes of the previous PCA plot were used in the db-RDA model as explanatory variables. Niche variables that best fit the db-RDA model were selected using the "ordiR2step" function from "vegan". The significance of selected variables was assessed by the function 'anova.test' in "vegan". Furthermore, to test the partial effect of significant niche variables on fungal community composition, the biogeographic region was used as a covariate in a partial db-RDA model. The significance of the model and the first two constrained axes was assessed by the function 'anova.test' in "vegan".

3. Results 3.1. Fungal diversity and community structure in Orchidaceae

Phylogenetic alpha-diversity of fungal communities, as measured by PSV values, differed significantly between orchid trophic modes and biogeographic regions, but not between subfamilies and biomes (Table 1; Fig. 1a, S2a, S3a and S4a). Pairwise comparisons showed that partially mycoheterotrophic species on average had the highest PSV values, whereas fully mycoheterotrophic species had the lowest PSV (Fig. 1a). Significantly lower PSV values were found for orchid species distributed in "Australia" compared with those in Eurasia, Southeast Asia and Africa (Fig. S3a).

Positive NRI and NTI values representing phylogenetic clustering were found for orchid trophic modes and biogeographic regions, but not for subfamilies and biomes (Table 1; Figs. 1b, S2b, S2c, S3b, S3c, S4b, S4c and S5). Fully and partially mycoheterotrophic species had the highest and the lowest level of phylogenetic clustering, respectively (Figs. 1b and S5). Orchid species distributed in Australia had significantly higher levels of phylogenetic clustering than species in Eurasia and Africa (Fig. S3b and S3c).

Significant differences in mean UniFrac distance values were found for orchid trophic modes, biogeographic regions and biomes, but not for subfamilies (Table 1; Figs. 1c, S2d, S3d and S4d). Larger UniFrac distances were found for partially and fully mycoheterotrophic species than autotrophic species (Fig. 1c), for species distributed in Eurasia than those in Africa, Southeast Asia and Australia (Fig. S3d) and for species distributed in "Temperate Forests and Grasslands" than those in "Tropical and Subtropical Forests" (Fig. S4d). Furthermore, multivariate analyses showed significant differences in fungal community composition between trophic modes (MANOVA: F = 7.978, p-value = 0.001; Fig. 1d), subfamilies (MANOVA: F = 3.789, p-value = 0.001), biogeographic regions (MANOVA: F = 4.357, p-value = 0.001) and biomes (MANOVA: F = 1.905, p-value = 0.001). In addition, variation partitioning models showed that orchid categorical groups (trophic mode, subfamily, biogeographic region and biome) contributed to 5.1%, 4.0%, 7.6% and 3.6% of the total variation in fungal communities (Fig. 2).

3.2. Phylogenetic distribution of orchid-fungal associations

Mantel test showed no significant correlation between fungal community structure and orchid phylogenetic distance (R = −0.027, p-value = 0.921), indicating no effect of host phylogenetic relatedness on fungal community assembly. The phylogenetic mapping of orchid–fungus associations (Fig. 3) suggests that fungal associations display divergent distribution patterns among orchid trophic modes rather than subfamilies.

3.3. The effect of environmental niche on fungal community composition

The PCA showed three apparent clusters of orchid species in the two-dimensional trait space (Fig. 4). With a total of 43% of the overall variation represented by the first two components, we observed segregations between trophic modes, biogeographic regions and biomes of orchid species within the trait space (Figs. 4 and S6). In the first cluster (upper left), South America and Southeast Asia were overlaid with "Tropical and Subtropical Forests and Grasslands" and were represented by niche variables including "Evergreen Broadleaf Trees", "Annual Mean Temperature" and "Annual Precipitation". In the second cluster (upper right), Eurasia and North America were overlaid with "Boreal Forests/Taiga" and "Temperate Forests and Grasslands" and represented by niche variables including "Evergreen/Deciduous Needleleaf Trees" and "Temperature Annual Range". Australia in the third cluster (bottom) was represented by Mediterranean Forests, Woodlands & Scrub, "Temp. Mean Diurnal Range" - Mean of monthly (max temp - min temp), "Topsoil Sand Fraction" and "Topsoil Bulk Density". Autotrophic and fully mycoheterotrophic orchids were scattered over the three clusters, while partially mycoheterotrophic orchids mainly occurred in the second and third clusters.

We used phylogenetic distance-based redundancy analysis (db-RDA) to explore whether niche variables affected the fungal community structure of orchids. A set of non-correlated niche variables that contributed most to the first two axes of PCA (Figs. S7 and S8) were implemented as explanatory variables in the db-RDA model. Model selection indicated that "Temp. Mean Diurnal Range", "Evergreen Broadleaf Trees" and "Topsoil Sand Fraction" were significant variables explaining a small portion of the variation in fungal community structure (db-RDA: R² = 0.049; F = 4.059; p-value = 0.001). The partial dbRDA model (Fig. 5) using biogeographic region as a covariate showed that half of the explained variation was assigned to niche variables (R² = 0.024; F = 2.170; p-value = 0.001). The significance of the first axis was supported by permutation tests (dbRDA1: F = 4.440, p-value = 0.007).

4. Discussion

In this study, we used a global dataset of orchid-fungal associations to study how mycorrhizal fungal communities vary in relation to orchid phylogeny, trophic mode, biogeographic distribution and broad niche variables. Our findings show that the diversity and community structure of orchid mycorrhizal fungi are likely influenced by the trophic mode and biogeographic distribution, but not by the overall phylogenetic relatedness of orchid species (Figs. 1–5 and S2–S4; Table 1). These findings indicate that fungal community composition is shaped by host plant ecophysiology and biogeography on a global scale.

The relationship between orchid phylogeny and the structure of their fungal communities has been suggested for small-scale phylogenetic clades (Shefferson et al., 2007, 2010; Jacquemyn et al., 2011). However, we detected no effects of host phylogeny on fungal community structure at a large phylogenetic scale. Similar results have been found for mycorrhizal associations in the orchid subfamily Cypripedioideae (Shefferson et al., 2019), arbuscular mycorrhizal associations (Zhao et al., 2021) and other host–microbial interactions on a large scale (Malacrinò, 2022). There are likely two explanations for why studies have detected phylogenetic signals in small subclades of orchids but not in a large phylogenetic framework. First, the multiple independent evolutionary origins of mycoheterotrophy in subclades of Orchidaceae (Wang et al., 2021) may diffuse the cophylogenetic signal between orchids and mycorrhizal fungi. Consequently, in this study we see the effect of trophic mode instead of host phylogeny on fungal community assembly (Fig. 1 and Table 1). Second, orchids have been diversifying into different habitats (Chomicki et al., 2015; Givnish et al., 2015), particularly the contrasting growth conditions of terrestrial and epiphytic species (Martos et al., 2012; Xing et al., 2019; Qin et al., 2020). This habitat specialization has also likely led to the recruitment of locally available fungal partners (Li et al., 2021), which may obscure the phylogenetic signal in mycorrhizal interactions.

Our finding that fungal community composition differs between orchid trophic modes is consistent with previous studies on subclades of orchids (Ogura-Tsujita et al., 2012; Těšitelová et al., 2015; Yagame et al., 2016) and arbuscular mycorrhizal plants (Gomes et al., 2017; Perez-Lamarque and Selosse, 2020; Zhao et al., 2021). Corroborated with orchid family-wide evolutionary analyses on the associated fungal lifestyles (Wang et al., 2021), our global analyses at the fungal community level reflect orchids' consistent physiological constraints on symbiotic association regardless of their phylogenetic and ecological distribution. Orchid partner preferences in mycorrhizal interactions suggest different nutritional needs from their fungal partners (Bever et al., 2009; Kiers et al., 2011; Bever, 2015). Unlike their autotrophic relatives, fully mycoheterotrophic orchids have obligate nutritional demands on mycorrhizal fungi throughout their entire life cycle, and thus tend to recruit ectomycorrhizal fungi or wood/litter saprobes to access specialized substrates (Taylor and Bruns, 1997; Selosse and Martos, 2014; Wang et al., 2021). Recent comparative omics studies and ecophysiological experiments have confirmed that fungal guilds or families differ in the ability to take up nutrients from organic material (Nurfadilah et al., 2013; Kohler et al., 2015; Fochi et al., 2017; Nehls and Plassard, 2018; Tedersoo and Bahram, 2019).

Interestingly, we found that partially mycoheterotrophic orchids associate with a broad partner breadth that was more phylogenetically dispersed than those of autotrophic and fully mycoheterotrophic orchids (Fig. 1). Recent high-throughput sequencing studies have determined that many partially mycoheterotrophic orchids are mycorrhizal generalists associated with both rhizoctonia and non-rhizoctonia fungi (Jacquemyn et al., 2016a, 2021; May et al., 2020; Xing et al., 2020; Wang et al., 2023). A more diverse fungal community may represent an advantage (complementarity effect) for those orchids to access broad niche dimensions (Batstone et al., 2018). Moreover, dual associations and partner switching may serve as necessary intermediate steps for evolutionary transitions to later symbiotic stages (Jacquemyn and Merckx, 2019; Wang et al., 2021, 2023). It has been inferred that the ultimate loss of photosynthesis in fully mycoheterotrophic orchids was likely promoted by gradually discarding fungal taxa that had already colonized roots of their ancestors (Jacquemyn and Merckx, 2019; Wang et al., 2021) — an evolutionary trajectory vividly termed the 'waiting-room' hypothesis (Selosse et al., 2022). However, these evolutionary inferences require direct experimental validation, e.g., long-term co-culture assays with synthetic fungal communities to simulate the stepwise exclusion of fungal partners observed in nature.

The biogeographic distribution of orchid species may have an important impact on mycorrhizal fungal diversity and composition on the global scale (Fig. 2). In line with previous studies (Phillips et al., 2011; Jacquemyn et al., 2017), our global analyses show that Australian orchid species tend to specialize in a small range of phylogenetically related fungi compared with other biogeographic regions (Table 1 and Fig. S3). Because orchid mycorrhizal fungi are geographically widespread (Jacquemyn et al., 2017), the fungal specialization observed in Australian orchid species is probably not limited by fungal availability but rather related to landscape characteristics and mycorrhizal ecology (Phillips et al., 2011). Intriguingly, O'Donnell et al. (2024) demonstrated that fungal specificity can serve as an inherent trait to clarify phylogenetic relationships among Australian orchids, complementing molecular and morphological data where they fall short. This emphasizes that fungal specificity is probably a hallmark of Australian orchids. In this study, we found that different biogeographic regions were represented by different sets of environmental factors (Fig. 4) and differences in the diversity and structure of orchid mycorrhizal communities were likely related to broad climatic, edaphic and vegetative factors that influence orchid distribution ranges (Fig. 5), indicating the effect of large-scale environmental filtering on mycorrhizal associations in the orchid family.

Our multivariate model was only able to explain a small portion of variation in fungal communities (Fig. 5). This is unsurprising given the diverse environmental niches of both orchid species (Givnish et al., 2015) and orchid mycorrhizal fungi (Wang et al., 2021). Under these conditions, ecological factors may influence mycorrhizal communities differently among species and habitats. Thus, it might be difficult to find a clear, cohesive pattern of environmental filtering on mycorrhizal associations within the orchid family and across the globe. Instead, microhabitats may play a greater role in shaping orchid mycorrhizal associations, as indicated by the usually patchy distribution of orchids and fungi within local habitats (Waud et al., 2016; McCormick et al., 2018; Kaur et al., 2021). It is likely that the encounter with suitable partners and the successful establishment of orchid mycorrhizal associations rely on a complex set of biotic and abiotic factors that may act differently on plants, fungi and substrates within microsites (McCormick et al., 2012; McCormick and Jacquemyn, 2014; Rasmussen et al., 2015). Future investigations on orchid physiology, fungal ecology and their interaction with the environment may help to reveal the decisive environmental factors for establishing orchid mycorrhizal associations in certain ecological contexts.

Similar to other global datasets of mycorrhizal types (Tedersoo and Bahram, 2019), there is a sampling bias in orchid mycorrhizal data, as well as in species occurrences in the GBIF database (Fig. S1), showing sampling paucity in tropical forests harboring extremely high orchid diversity. These biases leave much of the variation in orchid mycorrhizal communities unexplained (Fig. 2). New species and new records of orchid species including mycoheterotrophs are being continuously described globally, especially in the tropics (Metusala and Supriatna, 2017; Suetsugu et al., 2018; Aung et al., 2020; Yang et al., 2021). Furthermore, some putative autotrophic species have been shown to be in fact partially mycoheterotrophic according to the newest adoption of ²H and ¹⁸O isotopes in combination with ¹³C and ¹⁵N isotopes (Gebauer et al., 2016; Schiebold et al., 2018; Schweiger et al., 2019). Thus, future studies employing standardized protocols and more complete datasets should better elucidate patterns in fungal community variation across biogeographic distributions and trophic modes. Given that half of current orchid-fungal association data is unusable due to incomplete reporting, standardized documentation of host information, fungal identity, and habitat characteristics is critically important.

Acknowledgements

We would like to acknowledge the helpful discussions with Dr. Rutger Vos from Naturalis Biodiversity Center on environmental niche analyses. We sincerely thank for the funding provided by the China Scholarship Council (Grant No. 201804910634) and the Ecology Fund of the Royal Netherlands Academy of Arts and Sciences (KNAWWF/807/19039) to Deyi Wang.

CRediT authorship contribution statement

Deyi Wang: Writing – original draft, Visualization, Funding acquisition, Formal analysis, Conceptualization. Vincent S.F.T. Merckx: Writing – review & editing, Supervision, Conceptualization. Hans Jacquemyn: Writing – review & editing, Methodology. Sofia I.F. Gomes: Writing – review & editing, Supervision, Methodology, Conceptualization.

Data availability

The data that supports the findings of this study are available in the supplementary material of this article.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.08.002.