Network hubs in root-associated fungal metacommunities.

BACKGROUND: Although a number of recent studies have uncovered remarkable diversity of microbes associated with plants, understanding and managing dynamics of plant microbiomes remain major scientific challenges. In this respect, network analytical methods have provided a basis for exploring "hub" microbial species, which potentially organize community-scale processes of plant-microbe interactions.
METHODS: By compiling Illumina sequencing data of root-associated fungi in eight forest ecosystems across the Japanese Archipelago, we explored hubs within "metacommunity-scale" networks of plant-fungus associations. In total, the metadata included 8080 fungal operational taxonomic units (OTUs) detected from 227 local populations of 150 plant species/taxa.
RESULTS: Few fungal OTUs were common across all the eight forests. However, in each of the metacommunity-scale networks representing northern four localities or southern four localities, diverse mycorrhizal, endophytic, and pathogenic fungi were classified as "metacommunity hubs," which were detected from diverse host plant taxa throughout a climatic region. Specifically, Mortierella (Mortierellales), Cladophialophora (Chaetothyriales), Ilyonectria (Hypocreales), Pezicula (Helotiales), and Cadophora (incertae sedis) had broad geographic and host ranges across the northern (cool-temperate) region, while Saitozyma/Cryptococcus (Tremellales/Trichosporonales) and Mortierella as well as some arbuscular mycorrhizal fungi were placed at the central positions of the metacommunity-scale network representing warm-temperate and subtropical forests in southern Japan.
CONCLUSIONS: The network theoretical framework presented in this study will help us explore prospective fungi and bacteria, which have high potentials for agricultural application to diverse plant species within each climatic region. As some of those fungal taxa with broad geographic and host ranges have been known to promote the survival and growth of host plants, further studies elucidating their functional roles are awaited.

Background

Below-ground fungi in the endosphere and rhizosphere are key drivers of terrestrial ecosystem processes [1-4]. Mycorrhizal fungi, for example, are important partners of most land plant species, enhancing nutritional conditions and pathogen resistance of host plants [5-7]. In reward for the essential physiological services, they receive ca. 20% of net photosynthetic products from plants [8, 9]. Recent studies have also indicated that diverse taxonomic groups of endophytic fungi (e.g., endophytic fungi in the ascomycete orders Helotiales and Chaetothyriales) commonly interact with plant roots [10] and that those fungi play pivotal roles by providing soil nitrogen/phosphorous to their hosts [11-15], converting organic nitrogen into inorganic forms in the rhizosphere [16], and increasing plants’ resistance to environmental stresses [17-19]. Because of their fundamental roles, below-ground fungi have been considered as prospective sources of ecosystem-level functioning in forest management, agriculture, and ecosystem restoration [18-21]. However, due to the exceptional diversity of below-ground fungi [22-24] and the extraordinary complexity of below-ground plant–fungus interactions [25-27], we are still at an early stage of managing and manipulating plant-associated microbiomes [28-30].

In disentangling complex webs of below-ground plant–fungus associations, network analyses, which have been originally applied to human relations and the World-Wide Web [31, 32] and subsequently to biological systems [33, 34], provide crucial insights. By using network analytical tools, we can infer how plant species in a forest, grassland, or farmland are associated with diverse taxonomic and functional groups of fungi [25, 35–37]. Such information of network structure (topology) can be used to identify “hub” species, which are linked with many other species within networks depicting multispecies host–symbiont associations [38] (cf. [37, 39, 40]). Those hubs with broad host/symbiont ranges are expected to play key roles by mediating otherwise discrete ecological processes within a community [20, 25]. For example, although arbuscular mycorrhizal and ectomycorrhizal symbioses have been considered to involve distinct sets of plant and fungal lineages [41] (but see [42, 43]), hub endophytic fungi with broad host ranges may mediate indirect interactions between arbuscular mycorrhizal and ectomycorrhizal plant species through below-ground mycelial connections. As information of plant-associated fungal communities is now easily available with high-throughput DNA sequencing technologies [1, 22, 23], finding hub microbial species out of hundreds or thousands of species within a network has become an important step towards the understanding of ecosystem-scale phenomena.

Nonetheless, given that fungi can disperse long distances with spores, conidia, propagules, and animal vectors [44-48], information of local-scale networks alone does not provide thorough insights into below-ground plant–fungus interactions in the wild. In other words, no forests, grasslands, and farmlands are free from perturbations caused by fungi immigrating from other localities [49-53]. Therefore, to consider how local ecosystem processes are interlinked by dispersal of fungi, we need to take into account “metacommunity-scale” networks of plant–fungus associations [38]. Within a dataset representing multiple local communities (e.g., [26]), fungal species that occur in multiple localities may interlink local networks of plant–fungus associations. Among them, some species that not only have broad geographic ranges but also are associated with diverse host plant species would be placed at the core positions of a metacommunity-scale network [38]. Such “metacommunity hub” fungi would be major drivers of the synchronization and restructuring of local ecosystem processes (sensu [54]), and hence, their functional roles need to be investigated on a priority basis [38]. Moreover, in the screening of mycorrhizal and endophytic fungi that can be used in agriculture and ecosystem restoration programs [18, 21, 55], analytical pipelines for identifying metacommunity hubs will help us explore species that are potentially applied (inoculated) to diverse plant species over broad geographic ranges of farmlands, forests, or grasslands. Nonetheless, despite the potential importance of metacommunity hubs in both basic and applied microbiology, few studies have examined metacommunity-scale networks of plant–symbiont associations.

By compiling Illumina sequencing datasets of root-associated fungi [56], we herein inferred the structure of a metacommunity-scale network of below-ground plant–fungus associations and thereby explored metacommunity hubs. Our metadata consisted of plant–fungus association data in eight forest localities across the entire range of the Japanese Archipelago, including 150 plant species/taxa and 8080 fungal operational taxonomic units (OTUs) in temperate and subtropical regions. Based on the information of local- and metacommunity-scale networks, each of the fungal OTUs was evaluated in light of its topological position. We then examined whether fungal OTUs placed at the core of local-level plant–fungus networks could play key topological roles within the metacommunity-scale network. Overall, this study uncover how diverse taxonomic groups of mycorrhizal and endophytic fungi can form metacommunity-scale networks of below-ground plant–fungus associations, providing a basis for analyzing complex spatial processes of species-rich host–microbe systems.

Methods

Terminology

While a single type of plant–fungus interactions is targeted in each of most mycological studies (e.g., arbuscular mycorrhizal symbiosis or ectomycorrhizal symbiosis), we herein analyze the metadata including multiple categories of below-ground plant–fungus associations [56]. Because arbuscular mycorrhizal, ectomycorrhizal, and endophytic fungi, for example, vary in their microscopic structure within plant tissue [41], it is impossible to develop a general criterion of mutualistic/antagonistic interactions for all those fungal functional groups. Therefore, we used the term “associations” instead of “interactions” throughout the manuscript when we discuss patterns detected based on the Illumina sequencing metadata of root-associated fungi. Consequently, our network data could represent not only mutualistic or antagonistic interactions but also neutral or commensalistic interactions [25, 57, 58]. Our aim in this study is to gain an overview of the metacommunity-scale plant–fungus associations, while the nature of respective plant–fungus associations should be evaluated in future inoculation experiments.

Data

We compiled the Illumina (MiSeq) sequencing data collected in a previous study [56], in which community-scale statistical properties of below-ground plant–fungus associations were compared among eight forest localities (four cool-temperate, one warm-temperate, and three subtropical forests) across the entire range of the Japanese Archipelago (45.042–24.407 °N; Fig. 1). In each forest, 2 cm segment of terminal roots were sampled from 3 cm below the soil surface at 1 m horizontal intervals [56]. Those root samples were collected irrespective of their morphology and mycorrhizal type; hence, the samples as a whole represented below-ground relative abundance of plant species in each forest community. Host plant species were identified based on the sequences of the genes encoding the large subunit of ribulose-1,5-bisphosphate carboxylase (rbcL) and the internal transcribed spacer 1 (ITS1) of the ribosomal RNA region, although there were plant root samples that could not be identified to species with the rbcL and ITS1 regions [56]. The sequencing data are available through DDBJ Sequence Read Archives (accession DRA006339).

Fig. 1

Study sites examined in this study. Across the entire range of the Japanese Archipelago, root samples were collected in four cool-temperate forests (sites 1–4), one warm-temperate forest (site 5), and three subtropical forests (sites 6–8)

The Illumina sequencing reads of the fungal ITS1 region were processed using the program Claident [59, 60] as detailed in the data-source study [56]: the UNIX scripts used are available as Additional file 1. After filtering and denoising processes, operational taxonomic units (OTUs) representing less than 10 sequencing reads were discarded. The primers used were designed to target not only Ascomycota and Basidiomycota but also diverse non-Dikarya (e.g., Glomeromycota) taxa [61]. In most studies analyzing community structure of Ascomycota and Basidiomycota fungi, OTUs of the ITS region are defined with a cut-off sequence similarity of 97% [23, 62, 63] (see also [64]). Meanwhile, Glomeromycota fungi generally have much higher intraspecific ITS-sequence variation than other taxonomic groups of fungi [65]. Consequently, we used 97% and 94% cut-off sequence similarities for defining non-Glomeromycota and Glomeromycota fungal OTUs, respectively [56]. The OTUs were then subjected to reference database search with the query-centric auto-k-nearest-neighbor algorithm [59, 60] and subsequent taxonomic assignment with the lowest common ancestor algorithm [66]. For the molecular identification process, the nt database ver. 2015-11-11 was downloaded from the NCBI FTP server (https://ftp.ncbi.nlm.nih.gov/blast/db/) and sequences lacking genus-level information were removed. Based on the inferred taxonomy, the functional group of each fungal OTU was inferred using the program FUNGuild 1.0 [67].

After a series of bioinformatics and rarefaction procedures, 1000 fungal ITS reads were obtained from each of the 240 samples collected in each forest locality (i.e., 1000 reads × 240 samples × 8 sites). In the process, plant–fungus associations whose read counts represented less than 0.1% of the total read count of each sample were removed to minimize the effects of PCR/sequencing errors [68]. A sample (row) × fungal OTU (column) data matrix, in which a cell entry depicted the number of sequencing reads of an OTU in a sample, was obtained for each local forest (“sample-level” matrix) (Additional file 2: Data S2). Each local sample-level matrix was then converted into a “species-level” matrix, in which a cell entry represented the number of root samples from which associations of a plant species/taxa (row) and a fungal OTU (column) was observed: the binary (presence/absence) information was used in the conversion from sample-level matrices to species-level matrices. In the species-level matrix of each forest, 17–55 plant species/taxa and 1149–1797 fungal OTUs were detected from the local species-level matrices (Additional file 3: Data S3). In total, the matrices included 150 plant species/taxa and 8080 fungal OTUs (Additional file 4: Data S4).

Local networks

Among the eight forest localities, variation in the order-level taxonomic compositions was examined with the permutational analysis of variance (PERMANOVA; [69]) and the permutational analysis for the multivariate homogeneity of dispersions (PERMDISP; [70]) with the “adonis” and “betadisper” functions of the vegan 2.4–3 package [71] of R 3.4.1 [72], respectively. The β-diversity values used in the PERMANOVA and PERMDISP analyses were calculated with the “Bray-Curtis” metric based on the sample-level matrices (Additional file 2: Data S2). Note that the “Raup-Crick” β-diversity metric [73], which controls α-diversity in community data but requires computationally intensive randomization, was not applicable to our large metadata. Geographic variation in the compositions of fungal functional groups was also evaluated by PERMANOVA and PERMDISP analyses. The R scripts for the PERMANOVA and PERMDISP analyses are available as Additional file 5.

For each of the eight local forests, the network structure of below-ground plant–fungus associations was visualized based on the species-level matrix (Additional file 3: Data S3) using the program GePhi 0.9.1 [74] with the “ForceAtlas2” layout algorithm [75]. Within the networks, the order-level taxonomy of fungal OTUs was highlighted. Although the dataset of each local forest had the information of plant–fungus association frequency (i.e., the number of root samples from which respective plant–fungus associations were observed) (Additional file 3: Data S3), all the links were shown equally in the network visualization because varying line (link) thickness could result in considerable overlaps of links, reducing the visibility of the network figure.

To evaluate host ranges of each fungal OTU in each local forest, we first calculated the d’ metric of interaction specificity [76]. However, estimates of the d’ metric varied considerably among fungal OTUs observed from small numbers of root samples (Additional file 6: Figure S1) presumably due to overestimation or underestimation of host preferences for those rare OTUs. Therefore, we scored each fungal OTU based on their topological positions within each local network by calculating network centrality indices (degree, closeness, betweenness, and eigenvector centralities metrics of network centrality; [32]). Among the centrality metrics, betweenness centrality, which measures the extent to which a given nodes (species) is located within the shortest paths connecting pairs of other nodes in a network [77], is often used to explore organisms with broad host or partner ranges [38]. Thus, in each local network, fungal OTUs were ranked based on their betweenness centrality scores (local betweenness). Note that there was clear correlation between degree and betweenness centrality scores in each of the eight forests studied (see below; Pearson’s correlation r, 0.853–0.950; P < 0.0001 for all the eight sites). While binary (presence/absence) link information is used in analyses with the original betweenness metric [77], quantitative (frequency) link information (i.e., the number of root samples from which respective plant–fungus associations were observed) can be taken into account by using a newly developed “weighted” betweenness metric [78]. However, as far as we know, few biological interpretations have been made on the use of weighted betweenness in analyses of plant–microbe associations. Therefore, we used the original betweenness metric for binary data in this study. Note that binary and quantitative betweenness scores (Additional file 4: Data S4) were highly correlated with each other in each local network (Pearson’s correlation r, 0.481–0.826; P < 0.0001 for all the eight sites).

Metacommunity-scale network

By compiling the species-level matrices of the eight local forests, the topology of the metacommunity-scale network of plant–fungus associations was inferred. In general, species interaction (association) networks of local communities can be interconnected by species that appear in two or more local networks, thereby merged into a metacommunity-scale network [38]. In our data across the eight local forests, 2109 OTUs out of the 8080 fungal OTUs appeared in two or more localities. Therefore, we could infer the topology of a metacommunity-scale network, in which the eight local networks were combined by the 2109 fungal OTUs. In the metacommunity-scale network, plant species/taxa observed in different localities were treated as different network nodes because our purpose in this study was to explore fungi that potentially play key roles in synchronizing local ecosystem processes [38]. In total, 227 plant nodes representing local populations of 150 plant species/taxa were included in the metacommunity-scale network.

We then screened for fungal OTUs with broad geographic and host ranges based on the betweenness centrality scores of respective fungal OTUs within the metacommunity network (metacommunity betweenness, Bmeta). In general, species with the highest metacommunity betweenness scores not only occur in local communities over broad biotic/abiotic environmental conditions but also are associated with broad ranges of host/partner species [38]. Possible relationship between local- and metacommunity-scale topological roles was then examined by plotting local and metacommunity betweenness scores (Blocal and Bmeta) of each fungal OTUs on a two-dimensional surface. To make the betweenness scores vary from 0 to 1, betweenness centrality of a fungal OTU i was standardized in the local- and metacommunity-scale networks, respectively, as follows:

where Blocal, and Bmeta, were raw estimates of local- and metacommunity-scale betweenness of a fungal OTU i, and min() and max() indicated minimum and maximum values, respectively. For local betweenness of each OTU, a mean value across local networks was subsequently calculated (): the local communities from which a target OTU was absent was omitted in the calculation of mean local betweenness. On the two-dimensional surface, the OTUs were then classified into four categories: metacommunity hubs having high betweenness in both local- and metacommunity-scale networks ( ≥ 0.5; B'meta,  ≥ 0.5), metacommunity connectors that had broad geographic ranges but displayed low local betweenness ( <  0.5; B'meta,  ≥ 0.5), local hubs that had high betweenness in local networks but not in the metacommunity-scale network ( ≥ 0.5; B'meta,  <  0.5), and peripherals with low betweenness at both local and metacommunity levels ( <  0.5; B'meta,  <  0.5) [38]. Approximately, 1–2% of fungal OTUs can show betweenness scores higher than 0.5 in each local or metacommunity network, while the threshold value may be changed depending on the purpose of each study [38].

In addition to metacommunity hubs within the metacommunity-scale network representing all the eight localities, those within the metacommunity-scale network representing northern (sites 1–4) or southern (sites 5–8) four localities were also explored. This additional analysis allowed us to screen for fungal OTUs that potentially adapted to broad ranges of biotic and abiotic environments within northern (cool-temperate) or southern (warm-temperate or subtropical) part of Japan.

Results

Among the eight forest localities, order-level taxonomic compositions of fungi varied significantly (PERMANOVA; Fmodel = 35.7, R2 = 0.116, P < 0.001), while the differentiation of community structure was attributed, at least partly, to geographic variation in among-sample dispersion (PERMDISP; F = 13.2, P < 0.001) (Fig. 2a). Compositions of fungal functional groups were also differentiated among the eight localities (PERMANOVA; Fmodel = 34.9, R2 = 0.113, P < 0.001), while within-site dispersion significantly varied geographically (PERMDISP; F = 9.2, P < 0.001) (Fig. 2b). The proportion of ectomycorrhizal fungal orders (e.g., Russulales, Thelephorales, and Sebacinales) was higher in temperate forests than in subtropical forests, while that of arbuscular mycorrhizal fungi increased in subtropical localities (Fig. 2). The proportion of the ascomycete order Helotiales, which has been known to include not only ectomycorrhizal but also endophytic, saprotrophic, and ericoid mycorrhizal fungi [79], was higher in northern localities. In contrast, Diaporthales, which has been considered as predominantly plant pathogenic taxon [80] (but see [81]), was common in subtropical forests but not in others.

Fig. 2

Compositions of fungal taxa and functional groups in each forest. a Order-level taxonomic composition of fungal OTUs in each locality. The number of fungal OTUs detected is shown in a parenthesis for each forest. b Functional-group composition. The fungal functional groups were inferred by the program FUNGuild [67]

In each of the eight local networks depicting plant–fungus associations, some fungal OTUs were located at the central positions of the network, while others are distributed at peripheral positions (Additional file 7: Figure S2). Specifically, fungal OTUs belonging to the ascomycete orders Chaetothyriales (e.g., Cladophialophora and Exophiala) and Helotiales (e.g., Rhizodermea, Pezicula, Rhizoscyphus, and Leptodontidium) as well as some Mortierella OTUs had high betweenness centrality scores in each of the cool-temperate forests (Fig. 3a, b). In contrast, arbuscular mycorrhizal fungi (Glomeromycota) were common among OTUs with the highest betweenness scores in subtropical forests (Fig. 3a–c). Some fungi in the ascomycete order Hypocreales (e.g., Trichoderma, Ilyonectria, Simplicillium, and Calonectria) also had high betweenness scores in some temperate and subtropical forests (Fig. 3b).

Fig. 3

Fungal OTUs with highest local betweenness. a Order-level taxonomic composition of top 20 OTUs with highest local betweenness in each forest. (Additional file 4; Data S4) for betweenness scores of all fungal OTUs in respective local forests. b Genus-level taxonomic composition of top 20 OTUs with highest local betweenness. c Functional-group composition of top 20 OTUs with highest local betweenness

In the metacommunity-scale network representing the connections among the eight local networks, not only arbuscular mycorrhizal but also saprotrophic/endophytic fungi were placed at the central topological positions (Fig. 4; Additional file 8: Figure S3). Among non-Glomeromycota OTUs, Mortierella (Mortierellales), Cryptococcus [Trichosporonales; the Blast top-hit fungus in the NCBI database was recently moved to Saitozyma (Tremellales); [82]], Malassezia (Malasseziales), Oidiodendron (incertae sedis), Trichoderma (Hypocreales), and a fungus distantly allied to Melanconiella (Diaporthales) displayed the highest metacommunity betweenness (Table 1). Among the OTUs with high metacommunity betweenness, only a Mortierella OTU was designated as a metacommunity hub (i.e.,  ≥ 0.5; B'meta,  ≥ 0.5) and others had low betweenness scores at the local community level ( <  0.5; Fig. 5a).

Fig. 4

Metacommunity-scale network including all the eight local forests. The size of circles roughly represents relative scores of betweenness centrality. The functional groups of fungi inferred with the program FUNGuild [67] were organized into six categories, i.e., arbuscular mycorrhizal (bue), ectomycorrhizal (red), ericoid mycorrhizal (skyblue), saprotrophic/endophytic (yellow), plant pathogenic (purple), and other/unknown fungi (gray) (Additional file 4: Data S4). For plant species/taxa (green), the geographic information of source populations is indicated in (Additional file 8: Figure S3)

Table 1

Top 10 list of non-Glomeromycota OTUs with the highest betweenness within the metacommunity networks

OTU	Score	Phylum	Class	Order	Family	Genus	Category	NCBI Blast top hit	Accession	Cover	Identity
Full (8 sites)
F_0042*	1.000	–	–	Mortierellales	Mortierellaceae	Mortierella	Saprotroph/endophyte	Mortierella humilis	KP714537	100%	100%
F_0381	0.844	Basidiomycota	Tremellomycetes	Trichosporonales	Trichosporonaceae	Cryptococcus	Others_unknown	Saitozyma podzolica†	KY320605	92%	99%
F_0079	0.679	Ascomycota	Sordariomycetes	Hypocreales	Nectriaceae	–	Saprotroph/endophyte	Ilyonectria protearum	NR_152890	99%	100%
F_0489	0.576	–	–	Mortierellales	Mortierellaceae	Mortierella	Saprotroph/endophyte	Mortierella sp.	KM113754	100%	100%
F_0010	0.411	Ascomycota	Leotiomycetes	–	Myxotrichaceae	Oidiodendron	Ericoid_mycorrhizal	Oidiodendron maius	LC206669	100%	100%
F_0368	0.403	Basidiomycota	Malasseziomycetes	Malasseziales	Malasseziaceae	Malassezia	Others_unknown	Malassezia restricta	KT809059	100%	100%
F_0623	0.382	–	–	Mortierellales	Mortierellaceae	Mortierella	Saprotroph/endophyte	Mortierella gamsii	KY305027	100%	100%
F_1188	0.334	Basidiomycota	Tremellomycetes	Trichosporonales	Trichosporonaceae	Cryptococcus	Others_unknown	Saitozyma podzolica†	KY320605	92%	99%
F_0007	0.333	Ascomycota	Sordariomycetes	Diaporthales	Melanconidaceae	Melanconiella	Saprotroph/endophyte	Melanconiella elegans	KJ173701	100%	85%
F_0485	0.288	Ascomycota	Sordariomycetes	Hypocreales	Hypocreaceae	Trichoderma	Saprotroph/endophyte	Trichoderma sp.	HG008760	100%	100%
Northern 4 sites (cool-temperate)
F_0042*	1.000	–	–	Mortierellales	Mortierellaceae	Mortierella	Saprotroph/endophyte	Mortierella humilis	KP714537	100%	100%
F_0034*	0.795	Ascomycota	Eurotiomycetes	Chaetothyriales	Herpotrichiellaceae	Cladophialophora	Saprotroph/endophyte	Cladophialophora chaetospira	KF359558	100%	99%
F_0079*	0.767	Ascomycota	Sordariomycetes	Hypocreales	Nectriaceae	–	Saprotroph/endophyte	Ilyonectria protearum	NR_152890	99%	100%
F_0015*	0.763	Ascomycota	–	–	–	–	Others_unknown	Cadophora orchidicola	KX611558	100%	99%
F_0202*	0.749	Ascomycota	Eurotiomycetes	Chaetothyriales	Herpotrichiellaceae	Cladophialophora	Saprotroph/endophyte	Cladophialophora chaetospira	HQ871875	100%	99%
F_0195*	0.611	Ascomycota	Eurotiomycetes	Chaetothyriales	Herpotrichiellaceae	Cladophialophora	Saprotroph/endophyte	Cladophialophora chaetospira	EU035405	100%	100%
F_0181*	0.580	Ascomycota	Leotiomycetes	Helotiales	Dermateaceae	Pezicula	Endophyte	Pezicula melanigena	LC206665	100%	99%
F_0010	0.545	Ascomycota	Leotiomycetes	–	Myxotrichaceae	Oidiodendron	Ericoid_mycorrhizal	Oidiodendron maius	LC206669	100%	100%
F_0103*	0.538	Ascomycota	Eurotiomycetes	Chaetothyriales	Herpotrichiellaceae	Cladophialophora	Saprotroph/endophyte	Cladophialophora chaetospira	EU035403	100%	97%
F_0489*	0.534	–	–	Mortierellales	Mortierellaceae	Mortierella	Saprotroph/endophyte	Mortierella sp.	KM113754	100%	100%
Southern 4 sites (warm-temperate and subtropical)
F_0381*	1.000	Basidiomycota	Tremellomycetes	Trichosporonales	Trichosporonaceae	Cryptococcus	Others_unknown	Saitozyma podzolica†	KY320605	92%	99%
F_0042*	0.569	–	–	Mortierellales	Mortierellaceae	Mortierella	Saprotroph/endophyte	Mortierella humilis	KP714537	100%	100%
F_0610*	0.460	Ascomycota	Sordariomycetes	Hypocreales	Hypocreaceae	Trichoderma	Saprotroph/endophyte	Trichoderma spirale	KU948158	100%	100%
F_1188*	0.443	Basidiomycota	Tremellomycetes	Trichosporonales	Trichosporonaceae	Cryptococcus	Others_unknown	Saitozyma podzolica†	KY320605	92%	99%
F_0029	0.414	Ascomycota	Eurotiomycetes	Chaetothyriales	Herpotrichiellaceae	–	Others_unknown	Cladophialophora sp.	LC189029	100%	99%
F_0017	0.319	Ascomycota	–	–	–	–	Others_unknown	Scleropezicula sp.	KT809119	100%	98%
F_0007	0.314	Ascomycota	Sordariomycetes	Diaporthales	Melanconidaceae	Melanconiella	Saprotroph/endophyte	Melanconiella elegans	KJ173701	100%	85%
F_0485	0.303	Ascomycota	Sordariomycetes	Hypocreales	Hypocreaceae	Trichoderma	Saprotroph/endophyte	Trichoderma sp.	HG008760	100%	100%
F_0112	0.266	Basidiomycota	Agaricomycetes	Thelephorales	Thelephoraceae	Tomentella	Ectomycorrhizal	Tomentella stuposa	KR019860	100%	98%
F_0073	0.258	Ascomycota	Sordariomycetes	–	–	–	Others_unknown	Rhexodenticula acaciae	KY173442	94%	95%

In each of the three metacommunity-scale networks examined (full, cool-temperate, and warm-temperate/subtropical), fungal OTUs were ranked based on their betweenness centrality scores. As taxonomic information of Glomeromycota OTUs with high betweenness scores was redundant (e.g., Glomus spp. or Glomeraceae spp.), the top 10 list of non-Glomeromycota OTUs is shown. Taxonomy information of each OTU was inferred based on the query-centric auto-k-nearest-neighbor algorithm of reference database search [59, 60] and subsequent taxonomic assignment with the lowest common ancestor algorithm [66]. The results of the NCBI nucleotide Blast are also shown. For simplicity, the functional groups of fungi inferred with the program FUNGuild [67] were organized into several categories. See Additional file 4; Data S4 for details of the categories and for full results including Glomeromycota and other fungal OTUs

*Fungal OTUs classified as metacommunity hubs (mean local betweenness > 0.5; metacommunity betweenness > 0.5)

†Synonym, Cryptcoccus podzolica

Fig. 5

Relationship between local- and metacommunity-level betweenness. a Full meatcommunity. On the horizontal axis, the mean values of betweenness centrality scores across all the eight local forests are shown for respective fungal OTUs. On the vertical axis, the betweenness scores within the metacommunity-scale network consisting of the eight localities (Fig. 4) are shown for respective OTUs. b Metacommunity of cool-temperate forests. For the sub-dataset consisting of the four cool-temperate forests (Additional file 9: Figure S4), mean local betweenness and metacommunity betweenness are shown on the horizontal and vertical axes, respectively. c Metacommunity of warm-temperate and subtropical forests. For the sub-dataset consisting of the warm-temperate forest and the three subtropical forests (Additional file 10: Figure S5), mean local betweenness and metacommunity betweenness are shown on the horizontal and vertical axes, respectively

In the metacommunity-scale network representing the four cool-temperate forests (sites 1–4), many saprotrophic/endophytic fungal OTUs were associated with diverse plant species/taxa, located at the central topological positions within the network topology (Additional file 9: Figure S4; Fig. 5b). The list of these fungi with high metacommunity betweenness involved OTUs in the genera Mortierella, Cladophialophora (Chaetothyriales), Pezicula (Helotiales), and Oidiodendron as well as OTUs allied to Ilyonectria protearum (Nectriales) and Cadophora orchidicola (Helotiales) (Table 1). Most of these fungal OTUs also had high metacommunity betweenness designated as metacommunity hubs (Fig. 5b).

In the metacommunity-scale network consisting of the warm-temperate and subtropical forests (sites 5–8), arbuscular mycorrhizal and saprotrophic/endophytic fungi had the highest betweenness scores (Additional file 10: Figure S5; Fig. 5c). The list of non-Glomeromycota OTUs with the highest metacommunity betweenness included Saitozyma (Cryptococcus), Mortierella, Trichoderma, and Tomentella as well as OTUs allied to Cladophialophora, Scleropezicula (Helotiales), Melanconiella (Diaporthales), and Rhexodenticula (incertae sedis) (Table 1). Among the taxa, Saitozyma and Mortierella included OTUs classified as metacommunity hubs (Fig. 5c; Table 1). In an additional analysis of a metacommunity-scale network including only the three subtropical forests (sites 6–8), similar sets of fungal taxa were highlighted (Additional file 11: Table S1). Results also showed that similar sets of fungal taxa were highlighted in the analyses with binary and weighted betweenness metrics (Additional file 12: Table S2). The detailed information of the network index scores examined in this study is provided in Additional file 4; Data S4 .

Discussion

Based on the metadata of root-associated fungi across the Japanese Archipelago, we herein inferred the structure of a network representing metacommunity-scale associations of 150 plant species/taxa and 8080 fungal OTUs. Our analysis targeted diverse functional groups of fungi such as arbuscular mycorrhizal, ectomycorrhizal, ericoid-mycorrhizal, saprotrophic/endophytic, and pathogenic fungi, which have been analyzed separately in most previous studies on plant–fungus networks. The comprehensive analysis of below-ground plant–fungus associations allowed us to explore metacommunity hub fungi, which not only occurred over broad geographic ranges but also had broad host ranges in respective local communities. Consequently, this study highlights several taxonomic groups of fungi potentially playing key roles in synchronizing metacommunity-scale processes of temperate and/or subtropical forests.

In the metacommunity-scale network representing all the eight local forests (Fig. 4), fungi in several saprotrophic or endophytic taxa showed higher betweenness centrality scores than other fungi (Table 1). Mortierella is generally considered as a saprotrophic lineage [83], but it also includes fungi contributing to the growth and pathogen resistance of plants [84-86]. A phosphate solubilizing strain of Mortierella, for example, increases shoot and root growth of host plants under salt stress, especially when co-inoculated with an arbuscular mycorrhizal fungus [84]. In addition, polyunsaturated fatty acids produced by some Mortierella species are known to increase resistance of plants against phytopathogens [85, 86]. Fungi in the genus Trichoderma are commonly detected and isolated from the rhizosphere [83, 87]. Many of them inhibit the growth of other fungi, often used in the biological control of phytopathogens [88-90]. Some of them are also reported to suppress root-knot nematodes [91] or to promote root growth [92]. The analysis also highlighted basidiomycete yeasts in the genus Saitozyma or Cryptococcus (teleomorph = Filobasidiella), which are often isolated from soil [23, 93] and above-/below-ground parts of plants [94-97].

Along with possibly saprotrophic or endophytic taxa, ericoid mycorrhizal and phytopathogenic taxa of fungi displayed relatively high betweenness scores within the metacommunity-scale network representing all the eight local forests (Table 1). Specifically, Oidiodendron (teleomorph = Myxotrichum) is a taxon represented by possibly ericoid mycorrhizal species (O. maius and O. griseum) [98, 99], although fungi in the genus have been found also from roots of non-ericaceous plants and soil [100]. On the other hand, fungi in the family Nectriaceae are known to cause black foot disease [101], often having serious damage on economically important woody plants [102, 103]. Although we collected seemingly benign roots in the eight forests studied, some samples may be damaged by those pathogens. Alternatively, some fungi in the family Nectriaceae may be associated with plant hosts non-symptomatically, having adverse effects context-dependently.

Although the fungi mentioned above are candidates of metacommunity hubs, which are characterized by broad geographic ranges and host plant ranges, none except but a Mortierella OTU had high betweenness scores at both local and metacommunity levels (Fig. 5a). This result suggests that even if some fungi have broad geographic ranges across the Japanese Archipelago, few played important topological roles in each of the local networks representing plant–fungus associations. In other words, fungi that can adapt to biotic and abiotic environments in forest ecosystems throughout cool-temperate, warm-temperate, and subtropical regions are rare.

Therefore, we also explored fungi with broad geographic and host ranges within each metacommunity representing northern (cool-temperate) or southern (warm-temperate and subtropical) part of Japan. In the metacommunity consisting of the four cool-temperate forests (Additional file 9: Figure S4), fungal OTUs in the genera Mortierella, Cladophialophora, and Pezicula as well as those allied to Ilyonectria and Cadophora had the highest betweenness at both local and metacommunity levels, classified as metacommunity hubs (Fig. 5b; Table 1). Among them, Cladophialophora is of particular interest because it has been known as a lineage of “dark septate endophytes” [104-106] (sensu [15, 16, 107]). A species within the genus, C. chaetospira (= Heteroconium chaetospira), to which highest-betweenness OTUs in our data were closely allied, has been known not only to provide nitrogen to host plants but also to suppress pathogens [13, 17, 108]. Likewise, the Helotiales genus Pezicula (anamorph = Cryptosporiopsis) includes endophytic fungi [109-111], some of which produce secondary metabolites suppressing other microbes in the rhizosphere [112, 113]. Our finding that some of Cladophialophora and Pezicula fungi could be associated with various taxonomic groups of plants over broad geographic ranges highlights potentially important physiological and ecological roles of those endophytes at the community and metacommunity levels.

In the southern metacommunity network consisting of warm-temperate and subtropical forests (Additional file 10: Figure S5), some arbuscular mycorrhizal OTUs and Saitozyma (Cryptococcus) and Mortierella OTUs had high betweenness scores at both local and metacommunity levels, designated as metacommunity hubs (Fig. 5c; Table 1). Given the abovementioned prevalence of fungal OTUs allied to Cladophialophora chaetospira in the cool-temperate metacommunity, the contrasting list of metacommunity hubs in the southern (warm-temperate–subtropical) metacommunity implies that different taxonomic and functional groups of fungi play major metacommunity-scale roles in different climatic regions. This working hypothesis is partially supported by previous studies indicating endemism and vicariance in the biogeography of fungi and bacteria [114, 115], promoting conceptual advances beyond the classic belief that every microbe is everywhere but the environment selects microbial species/taxa colonizing respective local communities [116].

The roles of possible metacommunity hubs highlighted in this study are of particular interest from the aspect of theoretical ecology. Hub species connected to many other species in an ecosystem often integrate “energy channels” [117] within species interaction networks, having great impacts on biodiversity and productivity of the ecosystems [38]. The concept of “keystone” or “foundation” species [118, 119] can be extended to the metacommunity level, thereby promoting studies exploring species that restructure and synchronize ecological (and evolutionary) dynamics over broad geographic ranges [38]. Given that below-ground plant–fungus symbioses are key components of the terrestrial biosphere [1, 2], identifying fungal species that potentially have great impacts on the metacommunity-scale processes of such below-ground interactions will provide crucial insights into the conservation and restoration of forests and grasslands. We here showed that the list of metacommunity hubs could involve various lineages of endophytic fungi, whose ecosystem-scale functions have been underappreciated compared to those of mycorrhizal fungi. As those endophytic fungi are potentially used as inoculants when we reintroduce plant seedlings in ecosystem restoration programs [21, 55], exploring fungi with highest potentials in each climatic/biogeographic region will be a promising direction of research in conservation biology.

The finding that compositions of metacommunity hubs could vary depending on climatic regions also gives key implications for the application of endophytes in agriculture. Although a number of studies have tried to use endophytic fungi and/or bacteria as microbial inoculants in agriculture [18, 19, 120], such microbes introduced to agroecosystems are often outcompeted and replaced by indigenous (resident) microbes [121, 122]. Moreover, even if an endophytic species or strain increases plant growth in pot experiments under controlled environmental conditions, its effects in the field often vary considerably depending on biotic and abiotic contexts of local agroecosystems [18] (see also [123]). Therefore, in the screening of endophytes that can be used in broad ranges of biotic and abiotic environmental conditions, the metacommunity-scale network analysis outlined in this study will help us find promising candidates out of thousands or tens of thousands microbial species in the wild. Consequently, to find promising microbes whose inocula can persist in agroecosystems for long periods of time, exploration of metacommunity hubs needs to be performed in respective climatic or biogeographic regions.

For more advanced applications in conservation biology and agriculture, continual improvements of methods for analyzing metacommunity-scale networks are necessary [10]. First, while the fungal OTUs in our network analysis was defined based on the cut-off sequence similarities used in other studies targeting “species-level” diversity of fungi [63, 65], physiological functions can vary greatly within fungal species or species groups [15, 124]. Given that bioinformatic tools that potentially help us detect single-nucleotide-level variation are becoming available [125], the resolution of network analyses may be greatly improved in the near future. Second, although some computer programs allow us to infer functions of respective microbial OTUs within network data [67, 126], the database information of microbial functions remains scarce. To increase the coverage and accuracy of automatic annotations of microbial functions, studies describing the physiology, ecology, and genomes of microbes should be accelerated. With improved reference databases, more insights into the metacommunity-scale organization of plant–fungus associations will be obtained by reanalyzing the network data by compiling enhanced information of fungal functional groups. Third, as the diversity and compositions of plant–fungus associations included in a network can depend on how we process raw samples, special care is required in the selection of methods for washing and preparing root (or soil) samples. By sterilizing root samples with NaClO [127], for example, we may be able to exclude fungi or bacteria that are merely adhering to root surfaces. Meanwhile, some of those fungi and bacteria on root surfaces may play pivotal physiological roles in the growth and survival of plants [128]. Accordingly, it would be productive to compare network topologies of plant–microbe associations among different source materials by partitioning endosphere, rhizoplane, and rhizosphere microbial samples with a series of sample cleaning processes using ultrasonic devices [129]. Fourth, although this study targeted fungi associated with roots, our methods can be easily extended to network analyses involving other groups of microbes. By simultaneously analyzing the prokaryote 16S rRNA region [129-131] with the fungal ITS region, we can examine how bacteria, archaea, and fungi are involved in below-ground webs of symbioses. Fifth, not only plant–microbe associations but also microbe–microbe interactions can be estimated with network analytical frameworks. Various statistical pipelines have been proposed to infer how microbes interact with each other in facilitative or competitive ways within host macroorganisms [40, 132, 133]. Overall, those directions of analytical extensions will enhance our understanding of plant microbiome dynamics in nature.

Conclusions

By compiling datasets of below-ground plant–fungus associations in temperate and subtropical forest ecosystems, we explored metacommunity-hub fungi, which were characterized by broad geographic and host ranges. Such metacommunity-scale analyses are expected to provide bird’s-eye views of complex plant–microbe associations, highlighting plant-growth-promoting microbes that can be applied to diverse plant taxa in various environments. Given that endophytic fungi promoting the growth and pathogen resistance of host plants can be isolated from forest soil (e.g., Cladophialophora chaetospira [105]), the list of metacommunity-hub endophytic fungi featured in this study itself may include prospective species to be used in agriculture. By extending the targets of such network analyses to diverse types of plant-associated microbes (e.g., phyllosphere fungi and bacteria [81, 130, 134]) in various climatic/biogeographic regions, a solid basis for managing plant microbiomes will be developed.

Data S1. UNIX scripts for the bioinformatic pipeline. (ZIP 6 kb)

Data S2. Sample-level matrices of the eight forests examined. (XLSX 8388 kb)

Data S3. Species-level matrices of plant–fungus associations. (XLSX 9256 kb)

Data S4. Information of 8080 fungal OTUs analyzed. (XLSX 4319 kb)

Data S5. R scripts for the PERMANOVA and PERMDISP analyses. (ZIP 65 kb)

Figure S1. Number of sequencing reads, interaction specificity, and local betweenness. (DOCX 428 kb)

Figure S2. Structure of plant–fungus networks in each local forest. (DOCX 4270 kb)

Figure S3. Locality information within the full metacommunity-scale network. (DOCX 4949 kb)

Figure S4. Metacommunity-scale network of cool-temperate forests. (DOCX 3939 kb)

Figure S5. Metacommunity-scale network of warm-temperate and subtropical forests. (DOCX 4792 kb)

Table S1. Top 10 list of non-Glomeromycota OTUs with the highest betweenness within the subtropical metacommunity network. (DOCX 114 kb)

Table S2. Top 10 list of non-Glomeromycota OTUs with the highest weighted betweenness within the metacommunity networks. (DOCX 125 kb)