Structure and specialization of mycorrhizal networks in phylogenetically diverse tropical communities

Background The root mycobiome plays a fundamental role in plant nutrition and protection against biotic and abiotic stresses. In temperate forests or meadows dominated by angiosperms, the numerous fungi involved in root symbioses are often shared between neighboring plants, thus forming complex plant-fungus interaction networks of weak specialization. Whether this weak specialization also holds in rich tropical communities with more phylogenetically diverse sets of plant lineages remains unknown. We collected roots of 30 plant species in semi-natural tropical communities including angiosperms, ferns, and lycophytes, in three different habitat types on La Réunion island: a recent lava flow, a wet thicket, and an ericoid shrubland. We identified root-inhabiting fungi by sequencing both the 18S rRNA and the ITS2 variable regions. We assessed the diversity of mycorrhizal fungal taxa according to plant species and lineages, as well as the structure and specialization of the resulting plant-fungus networks. Results The 18S and ITS2 datasets are highly complementary at revealing the root mycobiota. According to 18S, Glomeromycotina colonize all plant groups in all habitats forming the least specialized interactions, resulting in nested network structures, while Mucoromycotina (Endogonales) are more abundant in the wetland and show higher specialization and modularity compared to the former. According to ITS2, mycorrhizal fungi of Ericaceae and Orchidaceae, namely Helotiales, Sebacinales, and Cantharellales, also colonize the roots of most plant lineages, confirming that they are frequent endophytes. While Helotiales and Sebacinales present intermediate levels of specialization, Cantharellales are more specialized and more sporadic in their interactions with plants, resulting in highly modular networks. Conclusions This study of the root mycobiome in tropical environments reinforces the idea that mycorrhizal fungal taxa are locally shared between co-occurring plants, including phylogenetically distant plants (e.g. lycophytes and angiosperms), where they may form functional mycorrhizae or establish endophytic colonization. Yet, we demonstrate that, irrespectively of the environmental variations, the level of specialization significantly varies according to the fungal lineages, probably reflecting the different evolutionary origins of these plant-fungus symbioses. Frequent fungal sharing between plants questions the roles of the different fungi in community functioning and highlights the importance of considering networks of interactions rather than isolated hosts. Supplementary Information The online version contains supplementary material available at 10.1186/s40793-022-00434-0.


Supplementary Methods: Supplementary Methods 1: Molecular work
We amplified two nuclear regions of the rDNA operon, the 18S rRNA gene and the ITS2 region, using two sets of tagged primers: AMADf-AMDGr pair (Berruti et al., 2017) and ITS86F-ITS4 pair (White et al., 1990;Turenne et al., 1999), that respectively amplify a fragment of 380 and 280 base pairs on average. PCRs were performed using the AmpliTaq Gold™ 360 kit (Thermo Fisher Scientific) following the manufacturer's instructions and ran for a total of 35 cycles at an annealing temperature of 50°C for the 18S rRNA marker and 56.5°C for the ITS2 one. To ensure better characterization of the mycobiota, the PCR step was replicated at least 2 times. Each sample was characterized by a unique combination of tagged primers (following Taberlet et al., 2018;Petrolli et al., 2021), and some combinations of primer pairs were left empty to evaluate the amount of tagged index hopping between samples (Zinger et al., 2019). In addition, negative PCR controls were also conducted (i) to identify the possible external contaminants and (ii) to evaluate the baseline cross-contamination between samples during the library preparation.
For both libraries (constituted of either the 18S rRNA or the ITS2 amplicons), PCR products were purified using magnetic bead-based clean-up (NucleoMag™ NGS Clean-up and Size Select, Macherey-Nagel), quantified using Qubit (Qubit dsDNA High Sensitivity Assay Kit, Invitrogen)), and (up to) 3 ng of DNA amplicons per sample (including controls) were sequenced using Illumina 2x250 bp MiSeq technology (v3 chemistry, Fasteris, Geneva). 3

Supplementary Methods 2: Bioinformatics
The 18S rRNA reads and the ITS2 reads were separately processed using a pipeline based on VSEARCH (Rognes et al., 2016) available on GitHub (https://github.com/BPerezLamarque/Scripts/). Pair reads were assembled, quality checked (removing reads that had on average more than 2 errors), demultiplexed (using cutadapt (Martin, 2011)), and clustered into operational taxonomic units (OTUs) using two different methods. We first used Swarm v2 (Mahé et al., 2015) with the fastidious option, a clustering approach that does not rely on a global threshold of similarity but instead uses local thresholds and amplicon abundances. Secondly, we performed a classical 97% OTU clustering using VSEARCH. We removed the chimera in both sets of OTUs using uchime_denovo in VSEARCH and assigned taxonomy to each OTU using usearch_global (BLAST algorithm). For the latter step, we used the Silva database r138 (Quast et al., 2013) for the taxonomic assignations of the 18S rRNA OTUs and the UNITE database v8.2 (Nilsson et al., 2019) for those of the ITS2 OTUs. Because the Silva database contains only a few Agaricomycetes sequences, we also assigned taxonomy at the order level to these Agaricomycetes 18S rRNA OTUs using the Fungal 18S RefSeq Targeted Loci Project (NCBI BioProject -PRJNA39195). We finally built an OTU table for each marker and removed OTUs present in less than 10 reads.
Glomeromycotina characterized using the 18S rRNA marker are often assigned to a virtual taxa (VT) based on the MaarjAM database (Öpik et al., 2010), however, the AMADf-AMDGr primer pair used in this study led to 18S rRNA amplicons that only partially overlap with virtual taxa and prevented us to do the assignations (Davison et al., 2015).
Next, the decontam pipeline in R (R Core Team, 2022) was applied to filter out the contaminants of our OTU tables, using the read abundances within the negative controls as well as the abundance profile of each OTU (Davis et al., 2018).
Finally, we reconstructed the phylogenetic tree of both fungi and plants. For the tree of the fungal OTUs, for each main fungal lineage, we aligned the OTU sequences 4 using MAFFT (Katoh & Standley, 2013), trimmed them with trimAl (Capella-Gutierrez et al., 2009), selected the best substitution model using ModelFinder (Kalyaanamoorthy et al., 2017), and reconstrued the maximum-likelihood tree using IQ-TREE (Nguyen et al., 2015) with 1,000 SH-aLRT and ultrafast bootstraps. For each main fungal lineage, a set of 3 OTUs were chosen as an outgroup (from a closely related fungal lineage) in order to root the phylogenetic trees.
For the plants, the plant mega-phylogeny from (Zanne et al., 2014) was pruned using Phylomatic (http://phylodiversity.net/phylomatic/) to obtain the phylogenetic tree of the plant species sampled in our three sampled communities. Plant taxa that were not identified at the species levels were added as polytomies at the origin of the clade.
Given that both Swarm and 97% OTU clustering gave qualitatively similar results, only the results obtained with Swarm are reported here.

Supplementary Methods 3: Evaluating cross-contaminations
First, thanks to our range of negative controls, we evaluated the amount of index hopping (using tagged primers left empty in the entire design) and crosscontaminations (using PCR negative controls). For most of the empty combinations of tagged primers, index hopping was negligible and below 5 reads (Supplementary Fig. 1a). For most of the PCR negative controls, contaminations were below 100 reads, but some samples had abundant contaminant concentrations ( Supplementary Fig. 1 b).
The contaminants mainly corresponded to human-skin Malassezia OTUs or other saprotrophs. All these negative OTUs were removed using the decontam pipelines (Davis et al., 2018). Less than 1% of the reads found in the negative controls were fungi that may form mycorrhizal interactions, suggesting that cross-contamination was low.
Thus, we considered that discarding the OTUs present in less than 5 reads in a sample represented a good balance to avoid the consideration of spurious interactions likely coming from cross-contaminations or index hopping.

5
Second, we verified that samples that were physically close in a plate (the organization of the plates was kept identical from the DNA extraction to the final amplicon pooling) did not tend to present a more similar composition than expected by chance because of cross-contaminations. To do so, we evaluated the similarities of the composition using Bray-Curtis dissimilarities computed using the vegdist function of the vegan R-package (Oksanen et al., 2016) and tested whether they were correlated with the Euclidian physical distances in each PCR plate using Mantel tests. We found no significant correlations (Mantel test: p-value>0.05) suggesting that crosscontaminations were overall negligible and did not significantly influence our results.
The significance of the checkerboard scores was evaluated using both quasiswap (a) and shuffle-sample (b) null models and using abundance, incidence, or binary networks. We considered the network to have a significantly higher checkerboard score (respectively lower checkerboard score) than expected by chance if p<2.5%, where p is the fraction of the null model values that are greater than or equal to (respectively lower than or equal to) the observed value.
(a) quasiswap null models:  The significance of the modularity values was evaluated using both quasiswap (a) and shuffle-sample (b) null models and using abundance, incidence, or binary networks.
We considered the network to be more modular (respectively less modular) than expected by chance if p<2.5%, where p is the fraction of the null model values that are greater than or equal to (respectively lower than or equal to) the observed value.
(a) quasiswap null models:  Boxplots present the median surrounded by the first and third quartiles, and whiskers extend to the extreme values but no further than 1.5 of the inter-quartile range.    using the Fruchterman-Reingold layout algorithm (Fruchterman & Reingold, 1991) from the igraph R-package. Fungal lineages (in rows) are ordered according to their network structures: networks that tend to be nested are at the top, whereas networks that tend to be modular are at the bottom. For each network, the total read abundances (L), the connectance (C), and the ratio of interactions within modules (Q) are indicated.
Q is computed from the most modular structure according to Beckett's algorithm for abundance networks, Q close to 0 indicates that most interactions are between modules (i.e. low modularity), while Q close to 1 indicates that most interactions are within modules (i.e. high modularity).  Plaine-des-Palmistes, or Dimitile), motif frequencies were computed (for all motifs containing from 2 to 5 species) and plotted in the following graph for Swarm OTUs.

Glomeromycotina
Motifs are numbered from 1 to 17 (adapted from Figure 1 of Simmons et al., 2019).    The photos at the top represent the main lycopod species in each community.