Orchid– mycorrhizal fungi interactions reveal a duality in their network structure in two European regions differing in climate

Network analysis is an effective tool to describe and quantify the ecological interactions between plants and root- associated fungi. Mycoheterotrophic plants, such as orchids, critically rely on mycorrhizal fungi for nutrients to survive, so investigating the structure of those intimate interactions brings new insights into the plant community assembly and coexistence. So far, there is little consensus on the structure of those interactions, described either as nested (generalist interactions), modular (highly specific interactions) or of both topologies. Biotic factors (e.g., mycorrhizal specificity) were shown to influence the network structure, while there is less evidence of abiotic factor effects. By using next- generation sequencing of the orchid mycorrhizal fungal (OMF) community associated to with plant individuals belonging to 17 orchid species, we assessed the structure of four orchid– OMF networks in two European regions under contrasting climatic conditions (Mediterranean vs. Continental). Each network contained four to 12 co- occurring orchid species, including six species shared among the regions. All four networks were both nested and modular, and fungal communities were different between co- occurring orchid species, despite multiple sharing of fungi across some orchids. Co- occurring orchid species growing in Mediterranean climate were associated with more dissimilar fungal communities, consistent with a more modular network structure compared to the Continental ones. OMF diversity was comparable among orchid species since most orchids were associated with multiple rarer fungi and with only a few highly dominant ones in the roots. Our results provide useful highlights into potential factors involved in structuring plant– mycorrhizal fungus interactions in different climatic conditions.


| INTRODUC TI ON
Ecological network analyses constitute a powerful tool to understand how community assemblies can form complex interactions in nature based on processes defining species coexistence (i.e., niche partitioning) and the diversity of species interactions (e.g., antagonistic or mutualistic interactions) (Bascompte, 2010;Olesen et al., 2007;Thompson, 2005). The structure of ecological networks observed across diverse communities is commonly described as nested or modular. The core of a nested network is composed of highly connected generalists that share a subset of their interacting partners with specialized species (Bascompte et al., 2003;Ings et al., 2009). In contrast, a modular network comprises distinct compartments (i.e., modules) of frequently interacting species (Bascompte & Stouffer, 2009;Guimerà & Amaral, 2005). The structural organization of networks can change due to diverse biotic and abiotic factors. Thus, the nature of species interactions (i.e., antagonistic vs. mutualistic) can influence the network topology. For example, antagonistic networks (e.g., host-parasitoid and plant-herbivore interactions) tend to be modular (Cagnolo et al., 2011;Thébault & Fontaine, 2010), probably due to more frequent specialized interactions caused by ecological and evolutionary pressures between partners (Lewinsohn et al., 2006). In comparison, mutualistic networks (e.g., plant-pollinator, plant-animal and plant-fungus interactions) tend to be highly nested (Bascompte et al., 2003;Bastolla et al., 2009;Jacquemyn et al., 2011). However, some studies have reported modular structures also in mutualistic interactions, especially in plant-mycorrhizal fungus networks (Põlme et al., 2018;Sepp et al., 2019) and in some pollination systems (Carstensen et al., 2016;Dicks et al., 2002;Joffard et al., 2019).
Plant-mycorrhizal fungi networks represent an interesting model for understanding the mechanisms underlying the coexistence of plant species. In particular, this system allows us to understand the importance of fungi in facilitating niche partitioning between plant species to reduce competition for nutrients (Jacquemyn et al., 2014;Jacquemyn et al., 2015;Van Der Heijden et al., 2008;Waterman et al., 2011). The intensity of host specificity varies greatly among mycorrhiza types (Van Der Heijden et al., 2015). Specialized interactions are commonly observed between plants and ericoid mycorrhizal fungi or orchid mycorrhizal fungi (OMF) (Li et al., 2021;Põlme et al., 2018), which lead to modular networks. Plant-ectomycorrhizal fungi (ECM) networks tend to be nested or modular (Bahram et al., 2014), while arbuscular mycorrhizal (AM) networks are commonly nested (Van Geel et al., 2018). Previous studies have revealed a duality in the network structure (with both nested and modular topologies) for plant-OMF (Xing et al., 2019) and plant-AM networks (Chagnon et al., 2012;Montesinos-Navarro et al., 2012). Thus, although host specificity varies across plant-mycorrhizal fungus interactions, both nested and modular networks have been documented, suggesting the importance of other potential drivers of network architecture, such as environmental factors (e.g., climatic and edaphic).
So far, studies demonstrating the effects of environmental factors on plant-mycorrhizal fungus network structure are still scarce.
The orchid-OMF system represents the strongest host-specific interactions among all mycorrhiza types (Põlme et al., 2018;Van Der Heijden et al., 2015). Orchid-OMF networks differ from the above-mentioned plant-mycorrhizal fungus networks by their high modularity and variation in nestedness values (Põlme et al., 2018).
These attributes may reflect the role of local abiotic conditions in shaping OMF communities in the soil and corresponding orchid-OMF networks Li et al., 2021). Orchid species critically rely on mycorrhizal fungi to acquire the nutrients essential for seed germination and subsequent life cycles (Smith & Read, 2008). Thus, orchids might be strongly conditioned by the presence and ecology of their symbionts (Jacquemyn et al., 2014;McCormick et al., 2016;Waud et al., 2016). Most temperate orchids are associated with OMF belonging to three Basidiomycete families: Tulasnellaceae, Ceratobasidiaceae (Cantharellales) and Serendipitaceae (Sebacinales) (Rasmussen, 2002;Taylor et al., 2002;Weiß et al., 2016), that is the former polyphyletic rhizoctonia group (Dearnaley et al., 2012). Across multiple plant populations and climatic regions, the orchid-OMF networks were described either with a nested topology (Jacquemyn et al., 2010;Jacquemyn et al., 2011) or a modular one Martos et al., 2012). On a local scale, coexisting orchids and OMF displayed a modular network in a Mediterranean grassland (Jacquemyn et al., 2015), while a nested (Xing et al., 2020) or both a nested and modular network (Xing et al., 2019) were found in a tropical habitat. The resulting discrepancies may be due to differences in sampling (i.e., plot size, number of orchid species, number of individuals), variations in orchid species specificity towards the fungi or changes in local abiotic conditions. Moreover, all the preceding cited studies used exclusively binary data, namely the presence or absence of OMF operational taxonomic units (OTUs), to build and interpret the ecological networks. This binary approach ignores the frequency of interactions and gives equal importance to all OTUs. In contrast, although unbalanced amplification of OTUs may arise through molecular biases, considering weighted data quantifies the frequency of a given interaction across different samples, providing more ecologically relevant information on the strength of individual interactions. Thus, the weight of rare OTUs is reduced, and more weight is given to the abundant ones, resulting in a more realistic topological structure of the network (Barrat et al., 2004;Ma et al., 2018). Hitherto, no study has yet analysed the structure of orchid-OMF networks using quantitative data.
Regarding the influence of climate on the structure of ecological networks, Kottke et al. (2013) suggested that climatic factors might influence the signal of nestedness in orchid-OMF networks. A recent meta-analysis on the effects of mean annual temperature and precipitation on network metrics found a negative correlation between the nestedness values and mean annual precipitation for the sum of all mycorrhiza types, but this correlation was absent in orchid mycorrhiza type, and no relationship was found for modularity values (Põlme et al., 2018). We may expect different signals for orchid mycorrhiza type, given their higher specialization for their fungal partners which may enable better resource partitioning (Novotná et al., 2023;Nurfadilah et al., 2013), especially in a drier environment where water resources are limited.
By using next-generation sequencing of OMF communities in 17 orchid species in total and a weighted data approach we describe the structure of four orchid-OMF networks based on the same-size plots situated in two different climatic regions (Mediterranean vs. Continental). We hypothesized that: (i) within each site, co-occurring orchid species associate with different fungal communities, which may facilitate the partitioning of limited nutrients; and (ii) the architecture of orchid-OMF networks differs between the two studied regions since the conditions of the Mediterranean climate (lower summer precipitation, higher annual temperature) may cause higher levels of network modularity than those of the Continental climate.

| Study sites and sampling
The study was conducted in four orchid species-rich grasslands located in two different climatic regions, according to Köppen and Geiger (1930): a warm-summer Continental climate (Czech Republic, referred to as CZ1 and CZ2) and a hot-summer Mediterranean climate (Southern France, referred to as FR1 and FR2). The climatic conditions for each site (monthly temperature and precipitation) were inferred from the nearest weather stations (<40 km away, at a similar altitude, years 2008-2018) ( Table 1). The edaphic conditions were analysed from a mixture of 10 random soil cores from 10 cm of soil depth at each study site. Each sample was air dried before analyses for total nitrogen and total carbon (dry combustion on an elemental analyser), total phosphorus (perchloric acid digestion followed by a colorimetric assay), plant-available phosphorus (Mehlich III method), organic matter content (loss on ignition at 550°C), soil texture (silt and clay; by the pipette method), soil water capacity and soil pH (H2O) . The analyses were conducted in the certified laboratories of the Research Institute for Soil and Water Conservation in Prague and the Department of Biology of Ecosystems at the Faculty of Science USB in České Budějovice (Czech Republic) (Table S1).
In May-June 2018, we established a 45 × 20-m plot at each site composed of four to 12 co-occurring orchid species, including eight and six species shared in the Mediterranean region and between the two regions, respectively (Table S1). In total, we sampled 17 orchid species inside the plots, from the genera Anacamptis, Dactylorhiza, Gymnadenia, Neottia, Neotinea, Orchis, Ophrys and Platanthera. Due to the observed hybridization between Platanthera bifolia and P. chlorantha at both French sites, we sampled both species together as Platanthera spp. Within each plot, three to five root pieces were carefully sampled from a minimum of four and up to 14 plants per orchid species (Table S1). The number of orchid species was higher in the French sites than in the Czech ones. Nevertheless, the Czech site CZ2 hosted an additional two orchid species located outside the sampling plot (Anacamptis pyramidalis and Platanthera bifolia) and four species in the plot with fewer than four individuals (Gymnadenia conopsea, tetraploid G. conopsea subsp. conopsea, Neotinea ustulata and Ophrys apifera). Initially, the data set contained 309 orchid individuals, which were reduced to 238 individuals for the analyses after removing samples (i) with a low number of reads (i.e., containing fewer than 70 reads), (ii) located outside the sampling plot, or (iii) with fewer than four plants per orchid species per plot. Finally, the orchid species were sampled across the whole plot area whenever possible to reduce the effect of spatial autocorrelation.

| Molecular assessment of mycorrhizal fungi
Roots were surface sterilized for 30 s in 4.7% sodium hypochlorite, washed in sterile distilled water and inspected for mycorrhizal colonization by light microscopy. Up to 12 well-colonized 3-mm-thick cross-sections were pooled per plant and stored in cetyltrimethylammonium bromide (CTAB) for molecular identification of mycorrhizal fungi. DNA was extracted by the CTAB method, according to Doyle and Doyle (1987). The Internal Transcribed Spacer 2 (ITS2) of TA B L E 1 Geographical, climatic and biological parameters of the four sites (France: FR1 and FR2; Czech Republic: CZ1 and CZ2).
Sequences were demultiplexed using MID tags and primers, and then those presenting any mismatch to the MID tags of primers were removed using mothur 1.40 (Schloss et al., 2009). An ITS2 region was extracted from the sequences by itsx 1.2 (Bengtsson-Palme et al., 2013). Putative chimeras were identified using de novo detection followed by comparison with a reference database, uchime version 7.2 (Nilsson et al., 2015), and finally removed with vsearch. The remaining sequences were clustered into OTUs at a 95% similarity threshold using cd-hit (Li & Godzik, 2006) and then taxonomically identified by blastn queries against the UNITE database version 8.0 (Nilsson et al., 2019;UNITE Community, 2019). All singletons and doubletons were removed, as they primarily represent artificial OTUs (Brown et al., 2015). Misaligned sequences were re-aligned using codoncode aligner version 4.02 (CodonCode Corp.) and merged at a 95% similarity threshold. OTUs present only in the negative controls were removed, while the read count of OTUs present at lower magnitude in the negative controls but abundant in the samples was reduced in the samples according to their abundance in the negative controls, as suggested by Nguyen et al. (2015). Finally, we manually checked all erroneously assigned taxa, and we removed OTUs with a low abundance of reads from samples to avoid overestimation of diversity by using a 0.05% cut-off of the proportion of reads for each sample (as in Duffy et al., 2019). Blasted sequences with an e-value <e−50 and a length of ≥100 bp were considered reliable to assign OTUs in the fungal kingdom. OTUs were regarded as OMF taxa when belonging to the families Ceratobasidiaceae, Tulasnellaceae and Serendipitaceae with a similarity threshold of at least 85%.
Subsequently, other fungal taxa such as ectomycorrhizal and saprobic fungi, known to form orchid mycorrhiza (Dearnaley et al., 2012), were screened manually and kept for further analysis. One representative sequence of each OTU was deposited in GenBank at NCBI (no. SUB12222961; see Table S2).

| Data analysis
All analyses were performed using r software version 3.5.2 (Core Team R, 2018). We used the log-transformed number of reads to mitigate the disparities of OTU amplifications by balancing the weights between overestimated abundant OTUs and underestimated rare ones. Fungal relative frequencies were then expressed as the number of reads of an OTU within a sample divided by the total number of reads found in this sample. We explored the fungal OTU diversity associated with an orchid species calculated based on the Shannon-Wiener diversity index (referred to as the Shannon index).
We investigated if fungal diversity differs between co-occurring orchid species between and within sites using a one-way analysis of variance (ANOVA) followed by the Tukey post hoc test. Further, the sampling effort was evaluated with species accumulation curves showing the level of saturation of OMF OTU richness in sampled orchid species, calculated for each orchid species using the "specaccum" function in the vegan package.
To investigate the first hypothesis, we compared the OMF community composition between co-occurring orchid species across all sites by computing Bray-Curtis distances using the "vegdist" function in the vegan package. Precisely, we calculated a Bray-Curtis dissimilarity index of fungal communities: in the first step, we averaged distances between each individual and other heterospecific individuals at a particular site, and in a second step, the mean individual distances were averaged per each orchid species (hereafter mentioned as the interspecific dissimilarity index, IDI). Thus, the resulting IDI measured the fungal dissimilarity index across orchid species in each sampled site. Then, it was compared among sites using a one-way ANOVA followed by the Tukey post hoc test. The OMF community composition associated with orchid species at the four sites was visualized by nonmetric multidimensional scaling (NMDS) using the Bray-Curtis distances in the vegan package (Oksanen et al., 2015).
PERMANOVA tests were carried out using the "adonis" function in the same package to test whether OMF community composition differed significantly among orchid species and sites.
To test the second hypothesis on the description of the orchid-OMF network topologies, we assembled a bipartite species-level matrix at each site based on the relative frequencies of OMF OTUs associated with co-occurring orchid species. The four networks were visualized based on the algorithm Atlas 2 using the gephi software 0.9.2. (Bastian et al., 2009). The nestedness and modularity metrics (Fortuna et al., 2010;Thébault & Fontaine, 2010) were calculated for each network using the bipartite package in R (Dormann et al., 2017;Dormann et al., 2009). As most network indices in this package cannot be applied directly to the relative frequencies, we converted them to integers by multiplying them by 10,000, following Pauw and Stanway (2015). The degree of nestedness was calculated using the weighted version of the nestedness metric (wNODF), which is less influenced by the matrix size (Almeida-Neto & Ulrich, 2011). This index varies between 0 (nested) and 100 (non-nested). The degree of modularity was calculated by the Q modularity index for the weighted matrix using the DIRTLPAwb+ algorithm (Beckett, 2016). The observed values of nestedness and modularity were tested against a null distribution using the "vaznull" model (Vázquez et al., 2007)  we calculated the network specialization H2′ quantitative index in the bipartite package (Blüthgen et al., 2006), which represents the number and dominance of OMF OTUs to which an orchid species is connected. This index ranges between 0 (generalization) and 1 (specialization). Finally, we compared our general statement with the assessment of the topology in networks by calculating the nestedness and modularity indices (i) in two of the studied networks (one per region) by exclusively selecting the orchid species that occurred in both regions (see Appendix S1); and (ii) in all networks with a binary matrix (i.e., the presence or absence of an interaction between an orchid species and an OMF) (see Appendix S2).
We then analysed in detail the role of fungal OTUs and the strength of their interactions with orchid species, to complete the general network analysis. Following Guimerà and Amaral (2005), different universal roles of fungal OTUs in a network were defined based on their location in the z-P parameter space, calculated as within-module degree zi and among-module connectivity Pi using the function "czvalues" (Guimerà & Amaral, 2005;Olesen et al., 2007).
Thus, we classified the fungal OTUs into four categories of network topologies: "peripherals" representing the specialists, "module hubs" and "module connectors" as the generalists, and "network hubs" for the supergeneralist OTUs (see Appendix S3 for more details).
Evaluation of the strength and statistical significance of the interaction between the relative frequencies of OMF OTUs and orchid species was determined using the multilevel pattern analysis in the package indicspecies (function "multipatt") (De Cáceres et al., 2010) with 999 permutations. Interaction strength was calculated by an indicator value index using the association function "IndVal.g" with the "group-equalized" option to avoid the potential bias of unbalanced sampling (De Cáceres & Legendre, 2009).

| Diversity and composition of the OMF community across orchid species and sites
The saturation curves of orchid species showed that the total fungal diversity is likely to increase with sampling effort across the sites, as no plateau was reached for most orchid species curves (except for Anacamptis pyramidalis and Orchis militaris in the FR2 site, Figure S1). Most orchid species were associated with multiple OMF OTUs across all sites in which they occurred; for example,  Table S3a). Regarding fungal diversity among co-occurring orchid species, the Shannon index showed a significant difference at the FR2 site (ANOVA, F = 4.42, df = 9, p = .0001; Figure S2, Table S3b), while no significant differences were observed at other sites.
There were significant effects of both site and orchid species identity on OMF community composition (PERMANOVA test,  Table 2; Table S4, Figure S4). Although most orchid species within the site did not share fungal OTUs, some species-pairs shared more OTUs than others, for example cooccurring Or. militaris and Or. purpurea (IDI <0.4 at both French sites) ( Figure S4). Other species-pairs displayed a low IDI only at one of the sites of co-occurrence (e.g., A. morio with G. conopsea at CZ1 and Neotinea ustulata with Ophrys passionis at FR1) ( Figure S4).

| The orchid-OMF BIPARTITE networks
The four orchid-OMF BIPARTITE networks were both significantly nested (p < .05) and modular (p < .001) ( Table 3)  were slightly more nested than those in France, as indicated by their higher Z-scores (Table 3). In contrast, the networks at the French sites (mostly FR1) were more modular (higher Z-scores) than those at the Czech sites ( Table 3). The number of identified modules varied from four at the Czech sites to seven and nine at the French sites ( Figure 2). The degree of network specialization shown by the H2′ index was lowest at the CZ1 site (H2′ = 0.57) and the highest at the CZ2 and FR2 sites (H2′ = 0.81; Table 3), indicating that co-occurring orchid species at the CZ1 site shared most of their OTU associates compared to the other sites. Moreover, both French networks were characterized by OTUs significantly associated with one particular orchid species compared to the Czech networks (see Table S5 Figure 2).
The most common and generalist fungal OTUs across both regions were T3, T4, C14 and Serendipitaceae S1, forming a high number of binary links (14, 15, 10 and 10 links, respectively, Figure 3). In certain orchid species, these generalist OTUs were highly dominant (Figures 2 and 3), such as T4 in G. conopsea (relative frequency in the roots of 52% in CZ1), T3 in A. pyramidalis (46% in FR1; 74% in FR2) and S1 in Neottia ovata (77% in CZ1; 30% in CZ2). These dominant fungi served as module connectors within the French networks ( Figure 2) but not in the Czech ones. In detail, the FR1 network was composed of eight generalist module connectors and one module hub, while the FR2 network was composed of three generalist module connectors (Figure 2). In all networks, most OTUs were specialists, referred to as peripherals, and formed all or most of their links within their modules. In the Czech networks, the peripherals represented 100% of all OTUs, while in the French networks, they represented 74.3% and 90% of OTUs (at FR1 and FR2, respectively). The largest modules were composed of 11 and 10 OTUs in the French sites, associated with Neottia ovata (FR1) and A. morio (FR2), followed by a module in CZ1 composed of nine OTUs associated with A. morio and G. conopsea (Figure 2). The smallest modules were composed of highly specific and dominant OTUs associated with particular orchids, such as Orchis mascula, G. conopsea and Platanthera spp.
at the FR1 site (Figure 2; see Table S5 for details).
The sampling sites in the Mediterranean region were richer in the number of orchid species than the Czech ones (varying from four to 12 species). We compared one network per region by using orchid species that occurred in both regions (see Appendix S1; Table S6). The reduced French network was still significantly nested and modular, revealing higher modularity and lower nestedness values than the Czech network. In addition, orchids growing in the reduced network in the Mediterranean grassland still shared fewer fungal partners than those growing in the Continental region (results not shown). When comparing the reduced and unreduced FR1 networks, the level of modularity (observed Q modularity and Note: Nestedness was analysed with the weighted NODF index (wNODF), modularity with the Q modularity index (Q Mod) and network specialization by the H2′ index.
Z-score) was lower when the orchid species number was reduced, while the level of nestedness was higher (observed wNODF and Zscore; Table S6). Finally, we observed that the same orchid species do not share many fungal OTUs (Table S7) and have a different fungal community when growing in different regions.

| DISCUSS ION
In our study, we provide evidence that the structure of orchidmycorrhizal fungus networks presents both nested and modular topologies simultaneously, yet each region was marked by the prevalence of one of the topologies. By sampling four orchid speciesrich grasslands in two different climatic regions, our study is the first to compare the orchid-fungus network architecture using quantitative data, providing greater ecological information than the frequently used binary data. This perspective offers new leads for understanding the possible effects of environmental conditions on the orchid community's assembly and their interaction with mycorrhizal fungi.

| OMF diversity and variation among orchid species and sites
Most of our studied orchid species harboured a high number of mycorrhizal fungi in their roots, such as Anacamptis morio (similarly reported in Bailarote et al., 2012;Ercole et al., 2015), Gymnadenia conopsea (Těšitelová et al., 2013;Vogt-Schilb et al., 2020), Neotinea ustulata and Neottia ovata (Djordjević et al., 2016;Kotilínek et al., 2015). Despite the occurrence of multiple fungal associations per species, we demonstrated that only a single or a few fungi were constantly more abundant in orchid roots than the others based on their relative frequencies. We hypothesize that dominant OTUs might represent fungi essential for seed germination of these orchids (see Těšitelová et al., 2022;J. Jersáková et al., unpublished data).
Most of the mycorrhizal fungi present in all studied orchid species belonged to the basidiomycetous family Tulasnellaceae, as previously reported by Jacquemyn et al. (2017). The second most common basidiomycetous family, Ceratobasidiaceae, occurred at various frequencies in nearly all our studied species,

F I G U R E 2
Orchid-OMF modular matrices and species roles. Red rectangles represent the modules (aggregated sets of highly interacting species). Coloured squares show the presence of the interaction, and the blue colour intensity refers to the OMF relative frequencies associated with the orchid species, classified into six classes of abundance (from <2% to 80%-100%). Coloured circles below the matrix represent generalist OTUs connecting the different modules within a network: green for the module connectors and yellow for a module hub. The absence of a circle represents specialist OTUs called peripherals that form all or most links within its module (see Appendix S2 for details). while Serendipitaceae OTUs were confined mainly to Neottia ovata Vogt-Schilb et al., 2020). Some studied orchid species harboured ECM, which are usually detected in mycoheterotrophic orchid species growing in shaded habitats (Selosse & Roy, 2009), but also found in orchids in open stands (Vogt-Schilb et al., 2020). However, their role in the nutrition of grassland orchids has not been unequivocally demonstrated (Rasmussen, 2002;Selosse et al., 2022). Future studies could explore their implications since we found abundant ECM in some orchid roots.
Although orchids formed multiple fungal associations, we found that most coexisting orchid species did not share OTUs.
One exception was Orchis militaris and Or. purpurea, which shared a similar fungal community. In their study, Jacquemyn et al. (2010) showed that Or. militaris shared the most dominant OTUs with Or. purpurea, which can be explained by the close phylogenetical signal between orchid species (Rezende et al., 2007). However, in our study, Or. purpurea was associated with a completely different fungal community than the phylogenetically related Or. simia and Or. anthropophora. Surprisingly, G. densiflora was the only species that was both very specific and shared its fungus with the other orchid species in the network. Although G. densiflora was available at only one site in this study and we could have overestimated species specificity, Těšitelová et al. (2013) found associations of G. densiflora with the same fungal OTU (GenBank accession no. KC243932) at two additional sites in the same region. This may indicate specialization for a fungus that is geographically widespread, similar to Or. mascula, which seems to conserve a high specificity for the same Tulasnellaceae isolate across different

| Asymmetric interactions and highly specific modules in orchid-OMF networks
Our study provided evidence that orchid-fungus networks can express both nested and modular topologies, highlighting the importance of addressing both parameters to describe the network structure (Lewinsohn et al., 2006). Both nestedness and  Table S2 for more information on OMF OTUs. modularity can provide stability to ecological networks (May, 1972;Teng & McCann, 2004) and possibly enhance the number of coexisting species (Bastolla et al., 2009;Olesen et al., 2007;Stouffer & Bascompte, 2011). Both topologies were found in some plantpollinator networks (Olesen et al., 2007), and it has been suggested that this structural duality is an outcome of trade-offs in the way densely connected communities can organize their interactions (Fortuna et al., 2010). Cagnolo et al. (2011) suggested that structural duality in plant-AM fungal networks might consist of specialized modules with internal nestedness.
High modularity in all four studied networks indicates the presence of specific subsets of orchid species and fungal OTUs interacting more frequently within a subset than with the rest of the network. Similarly, co-occurring orchid species were associated with distinct OMF communities, as shown in previous studies in the Mediterranean and other European habitats Jacquemyn et al., 2015). These results support the assumption that in species-rich orchid communities, fungal species are partitioned into subsets across co-occurring species, thus allowing niche partitioning and facilitating orchid coexistence (for a review see Põlme et al., 2018). Fungal differentiation between niches may regulate the coexistence of associated orchid species due to the ability of different mycorrhizal fungi to utilize distinct nutrient sources (Novotná et al., 2023;Nurfadilah et al., 2013).
Furthermore, we found that the four networks were also significantly nested, as expressed by asymmetric interactions between fungal OTUs and orchid species. Our results indicate that most OTUs were specialists with peripheral roles, while most orchids established multiple links with the OMF. The higher proportion of specialist OTUs contrasts with plant-AM fungus networks, in which the fungi are considered rather to be generalists (Põlme et al., 2018). A few generalist OTUs, such as Tulasnellaceae T3 and T4, Serendipitaceae S1 and Ceratobasidiaceae C14, displayed a module connector role in the networks as they interacted with different orchid species, creating high interconnectivity among orchid species within the network. They can be considered keystone species (Banerjee et al., 2018;Ings et al., 2009;Jordan, 2009), as their absence may weaken the networks and affect the dynamics and persistence of orchid populations.

| Comparison of the network architecture between the two climatic regions
When comparing the network structure across the two climatic regions, we demonstrated that the orchid-fungus networks displayed higher modularity and lower nestedness in the Mediterranean grasslands than in the Continental habitats. These results align with our hypothesis on the increase in orchid specialization towards their fungus species in the Mediterranean region.
Mean annual temperature and precipitation are known to affect both fungal diversity and community in the soil, as previously reported for ECM richness (Tedersoo et al., 2012) and AM community composition (Rasmussen et al., 2018). However, the effects of climatic factors on the structure of plant-mycorrhizal fungus networks are largely unexplored, especially in the orchid-OMF system (see Põlme et al., 2018). Here, we found differences in climatic factors between the two regions, with the Mediterranean region being a more droughtstressful environment, due to two-fold lower summer rainfall and higher mean annual temperature than the Central Europe region. Equally, the orchid-fungus networks in the Mediterranean region showed higher modularity, suggesting that greater hydric stress may be involved in this specialization process. This pattern was maintained when we compared two networks across regions composed of a similar number and assembly of orchid species, which supports our assumption.
Edaphic factors are also known to affect the orchid mycorrhizal fungal communities in the soil (Illyés et al., 2009;Mujica et al., 2016;Vogt-Schilb et al., 2020). In this study, we observed some differences in edaphic conditions between regions, namely strong soil acidity of the CZ1 site compared to slightly alkaline soils of all other sites. The CZ2 site had a low amount of available phosphorus and rather high amount of available calcium, which might be attributed to its clayish soil structure that has a great capacity to adsorb large quantities of nutrients, including calcium.
The carbon/nitrogen ratio, which is an indicator of organic matter decomposition rate, was similar between the sites, but additional unmeasured parameters such as soil temperature and moisture content also affect the microbial decomposition (Curiel Yuste et al., 2007). Due to an incomplete picture of soil characteristics and processes, we can only speculate which edaphic parameters could be responsible for the observed network patterns between the two regions. Though a higher nestedness signal was found in orchid-fungus networks in the Continental region. The association with multiple fungi may reduce the dependence of orchids on a single fungus, which can be particularly advantageous in a fluctuating environment and may occur more frequently in a nutrientpoor habitat (McCormick et al., 2006). However, this contradicts the expectation of increased nestedness in more drought-stressed conditions with limited water resources (Jacquemyn et al., 2010;Põlme et al., 2018).
We acknowledge that our results on edaphic and climatic con-

| Methodological considerations
The computation of network metrics of nestedness and modularity can be influenced by factors involving statistical approaches (i.e., data categories) and sampling design (sampling size, number of species).
Using binary data (Appendix S2;  (2015); however, they did not find a signal of nestedness. This structural duality was observed only in one previous study on an orchid-fungus network in the tropics (Xing et al., 2019). Compared to our main results based on weighted data, we observed that the nestedness Z-score values were significantly higher when using binary data, while the modularity values were lower. These trends were previously demonstrated for plant-ant networks (Miranda et al., 2019) and plant-ECM networks (Bahram et al., 2014), suggesting that the data type may influence the network structure.
The computation of nestedness and modularity is known to be sensitive to a few factors, such as the size of the network (i.e., the number of interactions), or the sampling area (Bascompte et al., 2003;Doré et al., 2020;Olesen et al., 2007;Põlme et al., 2018). While we analysed network architecture from equal sampling areas, the number and size of the module might have been influenced by the number of interacting species (Olesen et al., 2007) and the species composition (Rezende et al., 2007), as genetically related species tend to associate with similar mycorrhizal fungi (Li et al., 2021). Interestingly, our comparison of two networks composed of a subset of five orchid species occurring in both regions showed a decrease in modularity and an increase in nestedness values, which may confirm the effect of network size on both modularity and nestedness (Bascompte et al., 2003;Olesen et al., 2007). However, selecting a subset in a network composed of a larger number of interacting species does not reveal all the biases that size (i.e., number of species) can cause in the network structure. In this regard, more caution should be paid to network comparisons by considering the same sampling protocol at all sites, as suggested by Li et al. (2021), and sampling as many plant species as possible that are similar.

| CON CLUS IONS
Our study confirmed that coexisting orchid species within speciesrich sites have a tendency to associate with different fungal OTUs, although multiple fungi were shared among some orchid species. In the Mediterranean and Continental climatic regions, the four networks were both significantly nested and modular, though (i) there was no significant difference in OMF diversity and richness of fungi associated with orchids, (ii) nestedness tended to be slightly higher in the Czech networks than in the Mediterranean ones, and (iii) Mediterranean orchid species shared significantly fewer fungal partners within the network than at the Czech sites. This suggests that modular network architecture is more likely in the Mediterranean region, which differs mainly from the Continental one by its higher annual temperature and two-fold lower summer precipitation due to a pronounced summer drought.
Our results align with the hypothesis that low overlap in OTUs between plant species may decrease competition by resource partitioning (Jacquemyn et al., 2014;Waterman et al., 2011), though we did not provide a direct test of this. Further investigations on functional traits of both OMF fungi (e.g., enzymatic activities and growth rate) (see Novotná et al., 2023;Nurfadilah et al., 2013) and orchids (e.g., specificity towards fungal associates) may provide experimental evidence for the differentiation of resource uptake of plants via their fungal partners. Finally, we recommend using preferentially weighted data in future investigations to consider more biologically realistic interactions by distinguishing the presence of some highly dominant OMF in orchid species that may play a more important role in growth and germination processes than rare fungal taxa (Těšitelová et al., 2022;Vogt-Schilb et al., 2020).

ACK N O WLE D G E M ENTS
This work was supported by the Czech Science Foundation (GAČR 31-18-11378S). We thank Jitka Kocková, Linda Šternerová and Joana Pimentel for field and laboratory assistance, Jiri Košnar for technical assistance, and Marie-Pierre Dubois and Roula Zahab for assistance at the GEMEX platform of the CEFE laboratory.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflicts of interest.

O PEN R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The datasets used in this study are available in Dryad Repository https://doi.org/10.5061/dryad.0vt4b8h3m].

DATA AVA I L A B I L I T Y S TAT E M E N T
Our results were based on fungal DNA metabarcoding and obtained sequences are shared in Genbank (n° SUB12222961). The datasets used in this study are available in Dryad Repository (https://doi. org/10.5061/dryad.0vt4b 8h3m).