Figures

Abstract

Diverse clades of mycorrhizal and endophytic fungi are potentially involved in competitive or facilitative interactions within host-plant roots. We investigated the potential consequences of these ecological interactions on the assembly process of root-associated fungi by examining the co-occurrence of pairs of fungi in host-plant individuals. Based on massively-parallel pyrosequencing, we analyzed the root-associated fungal community composition for each of the 249 Quercus serrata and 188 Quercus glauca seedlings sampled in a warm-temperate secondary forest in Japan. Pairs of fungi that co-occurred more or less often than expected by chance were identified based on randomization tests. The pyrosequencing analysis revealed that not only ectomycorrhizal fungi but also endophytic fungi were common in the root-associated fungal community. Intriguingly, specific pairs of these ectomycorrhizal and endophytic fungi showed spatially aggregated patterns, suggesting the existence of facilitative interactions between fungi in different functional groups. Due to the large number of fungal pairs examined, many of the observed aggregated/segregated patterns with very low P values (e.g., < 0.005) turned non-significant after the application of a multiple comparison method. However, our overall results imply that the community structures of ectomycorrhizal and endophytic fungi could influence each other through interspecific competitive/facilitative interactions in root. To test the potential of host-plants' control of fungus–fungus ecological interactions in roots, we further examined whether the aggregated/segregated patterns could vary depending on the identity of host plant species. Potentially due to the physiological properties shared between the congeneric host plant species, the sign of hosts' control was not detected in the present study. The pyrosequencing-based randomization analyses shown in this study provide a platform of the high-throughput investigation of fungus–fungus interactions in plant root systems.

Funding: This work was financially supported by the Funding Program for Next Generation World-Leading Researchers of Cabinet Office, the Japanese government (to H.T.; GS014). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

In terrestrial ecosystems, various functional groups of root-associated fungi interact with plants [1]. Ectomycorrhizal and arbuscular mycorrhizal fungi are mutualistic partners of more than 80% of terrestrial plant species, enhancing host plant growth and survival by transporting soil nutrients [2]-[5], protecting host plants from pathogens and herbivores [6]–[8], and reducing competition among co-occurring plant individuals/species [9], [10]. Diverse clades of root endophytic fungi, which do not form mycorrhizae, also interact with diverse phylogenetic groups of plants [11], [12]. Some clades of those endophytic fungi are known to enhance the nutritional conditions of host plants [13], [14], whereas many others inhabit plant roots as commensalistic symbionts or parasites [15]. Those different functional groups of root-associated fungi often co-occur in plant roots [16], [17], and understanding the assembly processes of those ecologically and phylogenetically diverse fungi in root systems is one of the major challenges in fungal ecology.

Competitive and facilitative interactions between fungal species are considered to be the major factors responsible for the community organization of fungi in roots [18]–[20]. Fungal species in roots can compete with each other by impeding the colonization of others [21]–[24] or by expelling other species from host roots [25]. These competitive interactions between fungal species in roots, importantly, are expected to result in fine-scale segregated distributions of competing species [26], [27]. Interactions between fungal species can be competitive even between fungi in different functional or phylogenetic groups. In the roots of a Eucalyptus plant, for example, ectomycorrhizal fungi prevent arbuscular mycorrhizal fungi from infecting and proliferating in the host [28], [29]. On the other hand, the presence of a fungal species in roots does not always negatively affect the colonization of others. A morphological observation of fungal hyphae in roots revealed that ectomycorrhizal and root-endophytic fungi coexisted in a single root system of a Pinus tree, presumably because the two functional groups of fungi occupied different habitats within the roots [30]. Likewise, recent molecular studies have demonstrated the coexistence of ectomycorrhizal and endophytic fungi within roots [31]–[33]. In these cases, interactions between fungal species may be neutral or even facilitative.

Although competitive or facilitative interactions between fungal species can be assessed by inoculation experiments [22], [24], [34], many of the mycorrhizal and endophytic fungi that dominate root-associated fungal communities in natural forests are unculturable. Therefore, such inoculation-based experimental studies are applicable to only a part of fungus–fungus interactions occurring in the wild. An alternative research approach for elucidating the nature of fungal interspecific interactions is to examine the pattern of the presence/absence of fungal species in host individuals and thereby examine the spatial segregation and aggregation (co-occurrence) patterns. Spatial segregation could result from competitive interactions or differences in niches. On the other hand, spatial aggregation could arise from either species sorting [35], in which a pair of fungi with shared habitat requirements come to exist in a particular root, or from facilitative interactions. Therefore, community-wide analyses of the spatial segregation and aggregation of root-associated fungi are expected to reveal the species pairs that have a role in shaping fungal community structures. Indeed, some recent studies on root-associated fungi have examined such segregation/aggregation patterns and inferred the patterns of possible fungus–fungus ecological interactions [27], [36]

While previous studies analyzed the segregation/aggregation of pairs of root-associated fungi regardless of the effects of host plant species [26], [27], [36], such spatial patterns representing fungus–fungus competitive/facilitative interactions are expected to vary depending on host plant species. For example, if plant species with lower photosynthetic rates provide less carbohydrate to their root symbiont communities, it may promote competition for limited carbon resource and cause competitive exclusion between root-associated fungi [37]. Host plants' preference for a particular microenvironment (e.g., soil moisture) may also indirectly affect the relative competitive ability of root-associated fungi that interact with each other in root systems. By comparing segregation/aggregation patterns of root-associated fungi among different host plant species, we can examine such hypothetical plant-mediated processes of fungus–fungus interactions.

By identifying segregated and aggregated distributions of pairs of root-associated fungi on two oak species, we determined the patterns of fungal interspecific interactions and examined the dependence of such fungus–fungus interactions on background plant species. In a temperate forest in Japan, we first analyzed the community composition of root-associated fungi for the seedlings of co-occurring deciduous and evergreen oak species based on the pyrosequencing of fungal internal transcribed spacer (ITS) sequences. We then conducted a randomization test to detect the spatial segregation and aggregation of pairs of root-associated fungi. Furthermore, we examined whether or not the identity of host-plant species was associated with the spatial segregation or aggregation patterns of root-associated fungi.

Materials and Methods

Study area and sampling

The seedling samples were collected in a warm-temperate secondary forest on Mt. Yoshida located in Kyoto City, Japan (35.026°N, 135.786°E). No specific permissions were required for the location/activity. We confirmed that the field study did not involve endangered or protected species. The climate of Kyoto City is characterized by humid summers and dry winters: mean temperature and precipitation over the recent 30 years are 26.0 °C and 566.5 mm in summer (from June to August), and 5.6 °C and 164.9 mm in winter (from December to February) [38]. Mt. Yoshida (alt 121 m), which is a small hill with an area of 14.3 ha, is covered mainly by the two oak species Quercus serrata and Q. glauca (Fagaceae), while Pinus densiflora (Pinaceae) and Ilex pedunculosa (Aquifoliaceae) co-occur in the canopy layer. The two dominant oak species belong to different subgenus and have different ecological properties: Q. serrata is a deciduous species that occurs in the early stages of the secondary succession of temperate forests [39], while Q. glauca is an evergreen species whose seedlings occur both sunny and shaded understory of warm-temperate secondary forests [40]. Note that our field research site is not privately owned and sampling seedling in the research area is not banned.

From 20 to 31 May 2011, seedlings of each Quercus species were sampled at a minimum interval of 1 m: the number of sampled seedlings was 261 for Q. serrata and 199 for Q. glauca. The 460 sampling positions were recorded with a GPS device (Germin, GPSMAP 62S; Fig. 1). The size of sampled seedlings was 20–30 cm in height. To sample the seedlings, we dug to a depth of ca. 25 cm, taking great care not to damage the root tips of the seedlings. The amount of the dug soil of each root system was approximately 3,000 cm3. The sampled seedlings were individually stored in sealed plastic bags in an ice chest. On the same day of the fieldwork, we randomly collected ten 2-cm fragments of terminal roots per seedling in the laboratory: note that there were seedlings with less than ten 2-cm fragments of terminal roots and hence the number of root fragments collected per seedling ranged from five to ten. The terminal roots were stored in 1.5-ml tubes with 70% ethanol at –20°C until DNA extraction.

DNA extraction, PCR, and pyrosequencing

To remove soil adhering to roots, 1-mm zirconium balls were introduced into the sample tubes and then the tubes were shaken at 18 Hz for 3 min using TissueLyser II (Qiagen, Hilden, Germany) as detailed elsewhere [33]. After the cleaning, five terminal roots per sample were subjected to fungal DNA extraction. Terminal roots were transferred to a new tube and were pulverized with 4-mm zirconium balls by shaking at 20 Hz for 3 min. DNA extraction was conducted with the cetyltrimethylammonium bromide (CTAB) method [41].

To selectively amplify the fungal ITS2 region from the extracted DNA of the plant roots, a nested polymerase chain reaction (PCR) method was applied [33]. In the first PCR step, we amplified the entire ITS region using a fungus-specific primer (ITS1-F_KYO2 [42]) and a universal primer (ITS4 [43]) with the Ampdirect Plus (Shimadzu, Kyoto, Japan) buffer system. The ITS2 region was then amplified using a fusion primer of ITS3_KYO2 [42], which included the 454-adapter-A sequence and a 8-mer multiplex identifier (MID) tag sequence designed by Hamady et al.[44] to identify the source seedling (5′-CCA TCT CAT CCC TGC GTG TCT CCG ACT CAG [adapter A]–NNNNNNNN [MID]–GAT GAA GAA CGY AGY RAA [ITS3_KYO2] -3′), and a reverse fusion primer of ITS4, which included the 454-adapter-B sequence (5′-CCT ATC CCC TGT GTG CCT TGG CAG TCT CAG [adapter B]–TCC TCC GCT TAT TGA TAT GC [ITS4]-3′). The PCR products of all the 460 seedling samples were pooled in a new tube. We then purified the pooled PCR amplicons using ExoSAP-IT (GE Healthcare, Little Chalfont, Buckinghamshire, UK) and a PCR purification kit (Qiagen, Venlo, The Netherlands) before pyrosequencing using 454 GS Junior (Roche Diagnostics, Indianapolis, IN, USA). Because we did not obtain enough data in the first pyrosequencing run, we conducted an additional emulsion PCR and a pyrosequencing run. The pyrosequencing data were deposited to a public repository (DDBJ DRA: DRA000926).

Constructing operational taxonomic units and taxonomic identification

The pyrosequencing data were processed following the method of Toju et al.[33]. The low-quality 3′-tails of the pyrosequencing reads obtained were trimmed based on a threshold sequence quality value of 27. The reads were then filtered by a minimum sequence length of 150 bp excluding forward primer sequences and MID tags. To obtain the molecular operational taxonomic units (OTUs), the remaining reads were assembled as follows. The reads were sorted by seedling samples using the sample-specific MID tags and assembled into contigs for each sample with a minimum sequence similarity cutoff of 97% using Assams v0.1.2012.05.24 [45] (see also [46] for detailed assembling process), which is a parallelized pipeline for implementing the assembler program Minimus [47]. This within-sample assembling helped to avoid overestimates of OTU richness [48]. Possible chimeric sequences were detected and removed using UCHIME v4.2.40 [49] with a minimum score to report a chimera of 0.1. The within-sample contigs that passed the chimera removal process were subjected to further assembling of among-sample contigs with a sequence similarity cutoff of 97% using Assams, and the among-sample contigs were analyzed in the following statistical analyses as fungal OTUs. Note that the downstream statistical results did not qualitatively differ from those obtained based on 93% and 95% cutoff similarity settings in the among-sample clustering (see Results).

We conducted a BLAST-search of the fugal OTUs using the NCBI nt database on 22 February 2013. We also attempted identification based on the lowest common ancestor (LCA) algorithm [50], of which the results were much more conservative than the BLAST top-hit matches. Specifically, the query-centric auto-k-nearest-neighbor (QCauto) method implemented in the program Claident v0.1.2012.05.21 [51], [52] was applied using the reference-sequence information of the “all_genus” and “all_underclass” sequence databases and the NCBI-Taxonomy information of the “all_genus” and “all_underclass” taxonomy databases (see [51] for details of those databases). The query OTU sequences are shown in Data S1.

Binary data matrices

Based on the pyrosequencing dataset, we obtained a binary matrix that depicted the presence (1) or absence (0) of fungal OTUs in each of the 261 Q. serrata and 199 Q. glauca seedling samples (Data S2). Before obtaining the binary matrix, seedling samples with less than 20 pyrosequencing reads were excluded: 10 and nine seedling samples were excluded for Q. serrata and Q. glauca, respectively. In addition, OTUs representing less than 5% of the sample-total reads were excluded from each sample to reduce among-sample variance in α-diversity that resulted from variance in sequencing effort (i.e., variance in the number of sequencing reads among samples [mean = 129.2, SD = 69.7]). In this process, singletons and rare OTUs, which were expected to contain high proportions of pyrosequencing errors in their sequences [53], were eliminated. The eleven seedling samples that included the sequences of plants other than Fagaceae or Quercus spp. were also excluded from the data set: these samples were contaminated by DNAs of Ilex, Prunus, and Ericaceae plants, which commonly occurred in the study forest. Consequently, 249 Q. serrata and 188 Q. glauca seedlings were subjected to the following analyses (Data S2).

Fungal diversity and spatial autocorrelation

Based on the presence/absence data matrix of fungal OTUs (Data S2), the diversity and spatial structure of the fungal communities on the two Quercus species were evaluated. To assess the species richness of fungi, species accumulation curves (Mao Tau curves) were drawn for each of the two host species using the function ‘specaccum’ in the Vegan package [54] of R (version 3.0.2 [55]). To evaluate spatial autocorrelation in fungal OTU composition within the study site, a Mantel correlogram analysis was applied to each host plant species. In the analysis, we calculated Mantel's correlation (r) between dissimilarity in fungal OTU composition (i.e., β-diversity) and Euclidean distance spanning sampling positions (999 permutations). For the calculation of β-diviersity, we used Raup-Crick metric [56], which could minimize statistical artifacts resulting from difference in α-diversity among samples (see [56]). In addition, to test the spatial autocorrelation of the occurrence of each fungal OTU within the study site, we conducted a Moran's I analysis [57] for each fungal OTU that occurred on 10 or more seedling samples using the R package Ape 3.0-11 [58], [59].

Comparison of the root-associated fungal community structure between Q. serrata and Q. glauca

Prior to the statistical analysis of spatial segregation and aggregation of root-associated fungi, we examined differences in fungal community structure between Q. serrata and Q. glauca by PERMANOVA [60]. In this analysis, we measured dissimilarity in fungal OTU composition between seedling samples based on Raup–Crick β-diversity, and then tested for differences in the centroid of the fungal community structure of each Quercus species in multivariate space (9,999 permutations). The difference in the homogeneity of multivariate dispersion (the variance of β-diversity) between the fungal communities of the two hosts was also examined by PERMDISP [61] as implemented in Vegan.

We also examined the presence of host-specific fungal OTUs in the data set. A test using the multinomial species classification method (CLAM [62]) was performed to classify fungal OTUs preferentially associated with either host species and OTUs commonly found on both host species. A multinomial model was used to examine the statistical significance of respective fungal OTUs' preferences for host plants with a specialization threshold value of 2/3 (“supermajority” rule [62]). Because Bonferroni correction generally returns too stringent results in CLAM analysis [62], an α value of 0.001 was used as the threshold of statistical significance. The Vegan package of R was used in this analysis.

Segregation and aggregation of pairs of fungi in roots

We used the Checkerboard and Togetherness scores (C score and T score [63], [64]) as indices of spatial segregation and aggregation of fungal OTUs in the seedling samples, respectively. The C score is defined as (Ri – S) × (Rj – S), where Ri and Rj represent the total number of occurrences of species i and j, respectively, and S is the number of co-occurrences [63]. The T score is defined as S(N + S – Ri – Rj), where N is the number of seedlings analyzed [64]. In the presence of antagonistic interspecific interactions or the differentiation of niches, the observed C score of each pair of species is expected to be greater than that obtained by randomization under a null model. However, in the presence of facilitative interspecific interactions or shared habitat requirements, the observed T score is expected to be greater than that obtained by randomization under a null model. In contrast, a lack of significance suggests that species co-occur randomly.

For the dataset of each host plant, we tested the significance of C and T scores with 100,000 randomizations using the Bipartite package 2.0-1 [65] of R (Test 1). Observed and randomized C and T scores were standardized to range from 0 (the possible lowest level of segregation in terms of C scores and the possible lowest level of aggregation in terms of T scores) to 1 (the possible highest level of segregation in terms of C scores and the possible highest level of aggregation in terms of T scores) [66]. In each test for Q. serrata or Q. glauca dataset, OTUs observed in five or more seedlings and OTU pairs whose sum of seedling-sample counts were 25 or more were used because the statistical significance of C and T scores was difficult to examine for fungal OTU pairs with fewer sample counts. In addition to the examination for each host plant, the randomization analysis was applied to the whole dataset including both Q. serrata and Q. glauca seedling samples: OTUs observed in 10 or more Q. serrata and Q. glauca seedlings and OTU pairs whose sum of seedling-sample counts were 50 or more were used (Test 2). False discovery rate (FDR; Benjamini-Hochberg method) control [67] was applied to each randomization analysis.

We further assessed whether or not each pair of fungal OTUs was associated with each other in different ways on different host plant species. To this end, for each pair of fungal OTUs, we calculated the difference of C scores between the two host plant species (i.e., Cserrata – Cglauca where Cserrata and Cglauca were standardized C scores on Q. serrata and Q. glauca, respectively). Likewise, the difference of T scores between the two host plant species (Tserrata – Tglauca, where Tserrata and Tglauca were standardized T scores on Q. serrata and Q. glauca, respectively) was calculated for each pair of fungal OTUs. The significance of the difference of C or T scores on different host plants was tested based on 100,000 randomizations (Test 3). In this analysis, we used the fungal OTU pairs that were included in the dataset of both host plants in the Test 1. The fungal OTU pairs analyzed in the Test 1 of both host plants were used. FDR control was also applied to the analysis.

Finally, to assess whether whole community of root-associated fungi in the study site show segregated/aggregated patterns, we tested the significance of C and T scores with 100,000 randomizations. We also applied the analysis to each sub-dataset including a taxonomic or functional group of fungi (Ascomycota, Basidiomycota, and ectomycorrhizal fungal sub-datasets).

Results

Assembling and identification of molecular OTUs

In total, 65,150 reads were obtained by pyrosequencing (18,667 and 46,483 reads in the first and second GS Junior runs, respectively). Mean length of those reads were 348 (SD = 69.8) bp for the first run and 362 (SD = 55.29) bp for the second run. Only 0.175% of the reads were those of plants. The total numbers of OTUs were 1869, 1079 and 785 based on sequence cutoff similarities of 97, 95 and 93%, respectively; the numbers of singletons were 940, 414 and 270, respectively. After removing seedling samples with less than 20 pyrosequencing reads, those with ITS reads of plants other than Quercus, and OTUs representing less than 5% of the sample-total reads, the binary data matrix of the fungal community included 319, 274 and 242 OTUs with 97, 95 and 93%-cutoff similarities, respectively (Data S2).

Of the 319 OTUs detected with a cutoff sequence similarity of 97%, 89.7% were identified at the phylum level, 58.9% at the order level, 50.4% at the family level, and 40.7% at the genus level (Fig. 2). From Q. serrata and Q. glauca seedlings, 94 and 86 basidiomycete, 96 and 98 ascomycete, and three and one glomeromycete OTUs were detected, respectively. Of the 319 OTUs, 34 occurred in 10 or more seedling samples. Among these 34 most common OTUs, 14 were assigned to ectomycorrhizal genera. Of the remaining 23 OTUs, three and 19 were respectively assigned to Basidiomycota and Ascomycota at the phylum level but their genera remained unidentified; the remaining one OTU could not be identified even at the phylum level by the QCauto method (Table 1). BLAST searches against the NCBI nr/nt database (Table 1) indicated that the commonly observed ascomycete OTUs were allied to genera or species that had been generally detected from living plant tissue in previous studies (i.e., possibly endophytic ascomycetes; e.g., Catenulifera[68], Pezicula [teleomorph of Cryptosporiopsis[69]], Lophodermium[70], and Cladophialophora [teleomorph of Capronia[71]), except for an OTU (OTU 121) related to soil fungi in the genus Archaeorhizomyces[72]. The list of commonly observed fungi included Mycena, Oidiodendron, and Glomeromycota, which are known as saprobes or ericoid/arbuscular mycorrhizal fungi (Table S1). We found no inconsistency between the QCauto-based and BLAST-based identification results, although results by the QCauto method were more conservative than those of BLAST (Table 1). Taxonomic diversity of OTUs were qualitatively similar among identification results with different sequence-similarity cutoffs (see Tables S1 and S2 for results at 93% and 95% cutoffs).

Fungal diversity and spatial autocorrelation

Of the 319 OTUs detected with a cutoff sequence similarity of 97%, 103 occurred on both Q. serrata and Q. glauca, while 115 and 101 OTUs occurred only on either Q. serrata or Q. glauca, respectively (Fig. 1a). Species accumulation curves did not reach a plateau for either host-plant species (Fig. 1b). The average number of OTUs per seedling was almost similar between the two host plant species but it was statistically higher on Q. serrata than on Q. glauca (3.51 and 3.20 on Q. serrata and Q. glauca, respectively; t = 2.12, df = 414.377, P < 0.05).

A Mantel correlogram analysis indicated that the fungal OTU composition of the examined seedling samples was spatially auto-correlated at a very small spatial scale within the study site: i.e., the scale of the autocorrelation was < ca. 8 m for Q. serrata seedlings and < ca. 2 m for Q. glauca seedlings (Fig. 3c, d). In a Moran's I analysis, significant spatial autocorrelation within the study site was observed for only seven of the 34 common fungal OTUs examined (Table 1).

Comparison of root-associated fungal community structures between Q. serrata and Q. glauca

Although PERMANOVA analysis indicated a significant difference in community structures between the two host species (F = 4.79, P = 0.0002), the R2 value was very low (R2 = 0.0109), suggesting that the effects of host plant identity on fungal OTU compositions were small in the present dataset. In fact, the community compositions of root-associated fungi on Q. serrata and Q. glauca were largely similar to each other (Fig. 2). The PERMDISP analysis showed significant differences in the among-sample β-diversity (Raup-Crick) of fungal communities between the two host-plant species. The β-diversity of the fungal communities was greater on Q. glauca than on Q. serrata (average distance between each fungal community and the centroid of fungal community were 0.5765 for Q. serrata and 0.6283 for Q. glauca, F = 12.61, P = 0.0004). Consistent results were obtained in the analyses based on 95% and 93% cutoff sequence similarities: for simplicity, results at 97% cutoff sequence similarity are shown in the following statistical analyses.

The CLAM test revealed the presence of 18 OTUs commonly associated with both Quercus species and an OTU exclusively associated with Q. serrata, although many OTUs were too rare to be assigned host preference (Table 1). The OTU that exclusively occurred on Q. serrata (OTU 1089) showed 99% sequence similarity to Lactarius quietus (JF273529). The list of fungi commonly associated with both plant species included fungal OTUs related to various ectomycorrhizal fungi, e.g., Cenococcum geophilum (OTU 167), Russula cerolens (OTU 193), and Thelephora terrestris (OTU 211), and those related to possibly endophytic fungi in the orders Chaetothyriales and Helotiales, e.g., Catenulifera luxurians (OTU 757) and Pezicula sp. (OTU 329).

Segregation and aggregation of pairs of fungi in roots

When the seedling samples of the two host plants were analyzed separately (Test 1), no pair of fungal OTUs showed significant C or T scores after adjusting P values based on FDR control. However, many pairs of fungal OTUs displayed low (< 0.005) P values without FDR control (Table 2). In the C score analysis, an ectomycorrhizal ascomycete in the genus Cenococcum (OTU 167) and a possibly endophytic ascomycete (OTU 757) displayed segregated distribution on Q. serrata (Table 2). In the T score analysis, eight pairs of fungal OTUs displayed aggregated patterns with low (< 0.005) P values on either of the two host plants. Among the eight pairs, six were pairs of an ectomycorrhizal basidiomycete (Russula [OTU 185, 193, 509, or 1135] or Lactarius [OTU 205]) and a possibly endophytic ascomycete fungus (OTU 115 or 331; Table 2). The remaining two pairs were those of possibly endophytic ascomycetes (Table 2).

In the analysis in which the datasets of the two host plants were combined (Test 2), three pairs of fungal OTUs displayed statistically significant aggregated patterns even after FDR control (Table 3). Among the three pairs, two were those of an ectomycorrhizal basidiomycete (Russula [OTU 193] or Lactarius [OTU 205]) and a possibly endophytic ascomycete (OTUs 115 or 331). The remaining pair was that of possibly endophytic ascomycetes (OTUs 199 and 1845). Note that such aggregation of pairs of ectomycorrhizal and endophytic fungi or those of endophytic fungi was also observed in the datasets based on 93% or 95% cutoff similarities (Table S2). In addition to the abovementioned pairs, two pairs of an ectomycorrhizal basidiomycete (Russula [OTU 1135] or Lactarius [OTU 1089]) and a possibly edophytic ascomycete (OTU 115 or 331) and a pair of possibly root-endophytic ascomycete fungi displayed aggregated patterns with low (< 0.005) P values (Table 3). Likewise, segregated patterns with low P values were observed for a pair of an ectomycorrhizal ascomycete (Cenococcum [OTU 167]) and a possibly endophytic ascomycete (OTU 757), although the patterns were non-significant after FDR control (Table 3).

When the effects of host plants on the segregation/aggregation of fungal OTUs were examined (Test 3), no pair of fungal OTUs showed significant difference of C or T scores between Q. serrata and Q. glauca after FDR control; even without FDR control, no fungal pair showed low (< 0.005) P values (Table S3). These results indicate that the identity of host plant species did not affect the segregation/aggregation patters of root-associated fungi in the present dataset.

In the community-scale analysis of C scores (Table 4), significantly segregated patterns were observed within the entire community of the observed fungi and within the basidiomycete and ectomycorrhizal fugal sub-communities, while no sign of segregation was observed within the ascomycete sub-community. As expected by the C score analysis, observed values of togetherness (T) were very low (Table 5).

Discussion

In the present study, the sign of segregated or aggregated patterns was detected for a small number of fungal OTU pairs. Although the results may reflect the rarity of competitive or facilitative interactions between root-associated fungi in the study forest, there are potential statistical issues that may have hampered the detection of significant patterns. To make clear the potential pitfalls in pyrosequencing-based high-throughput analyses of segregated/aggregated distributions of fungi, we start with discussing problems related to multiple comparisons and the power of randomization tests.

When the checkerboard (C) and togetherness (T) scores of pairs of fungal OTUs were examined on each of the two host plant species (Test 1), no fungal pair showed statistically significant segregation nor aggregation patterns after FDR control (Table 2). Given that more than 179 pairs of fungal OTUs were examined in our pyrosequencing-based analysis, the application of the multiple comparison method might make the results prone to type II errors [73]. In the studies of segregation/aggregation patterns of root-associated fungi, such a statistical issue has been generally underappreciated. For example, a previous study discussed spatial segregation or aggregation of fungal species without any multiple-comparison adjustment of P-values [26], thereby making their results prone to type I errors (false positives) rather than type II errors (false negatives). Therefore, by applying multiple comparison methods to the datasets of the previous studies, one may be able to screen for pairs of fungi with strong sign of segregated/aggregated patterns.

In addition to the possible effects of multiple comparison methods, problems related to sample size could affect the results of randomization test. Intriguingly, when the samples of the two host species were pooled in the C and T score analysis (Test 2), three pairs of fungal OTUs showed statistically significant aggregation even after FDR control (Table 3). Given that the identity of host plants was not associated with the segregation/aggregation patterns of root-associated fungi in the current dataset (Test 3; Table S3), the fact that significant segregation/aggregation was observed in the simultaneous analysis of Q. serrata and Q. glauca seedlings (Test 2) but not in the independent analysis of them (Test 1) might be attributed to the larger number of samples in Test 2 than in Test 1. Importantly, when the sample counts of examined fungi (i.e., the number of seedling samples in which respective fungi occurred) are small, C or T scores of randomized data matrices frequently take the minimum (0) or maximum (1) values and hence deviate from normal distribution (Fig. S1), possibly reducing the power of the randomization tests. Therefore, in the present study, sample counts of fungi would have been insufficient to detect segregation/aggregation patterns independently on each host plant (Test 1), despite the intensive sampling of seedlings in the forest (249 Q. serrata and 188 Q. glauca seeding samples). Overall, the comparison of the results between Test 1 and Test 2 suggests that fungal pairs with strong segregation/aggregation signs become detectable by increasing sample size (and thereby sample counts; Fig. S1).

Although the abovementioned statistical issues should be treated with caution, the present results have important implications for competitive or facilitative interactions between root-associated fungi. Albeit non-significant when FDR control was applied, a segregated pattern with a low P value (< 0.005) was observed between an ectomycorrhizal fungus in the ascomycete genus Cenococcum and a possibly endophytic ascomycete fungus (Tables 2 and 3). This result is intriguing, given that most previous studies investigated the potential competitive interactions within ectomycorrhizal or root-endophytic fungal communities [19], [20], [26], [27] (see also [30], Tables 4 and 5) and did not assume that ectomycorrhizal fungi could affect the spatial distribution of endophytic (or other types of non-ectomycorrhizal) fungi and vice versa. However, as suggested by our present results, competition for space or resources within root systems could occur not only between fungi with similar ecological or physiological properties but also between fungi in different functional groups. Similar phenomena were reported in a study of root-associated fungi on an ericaceous plant [36]. Combined with the findings of the previous study, our present results suggest that simultaneous analysis of multiple functional groups would give novel insights into the community dynamics of root-associated fungi.

In addition to spatial segregation, aggregated patterns of pairs of fungi were observed in this study (Tables 2 and 3). These pairs included those of an ectomycorrhizal fungus in the genus Russula or Lactarius and a possibly-endophytic ascomycete fungus (Tables 2 and 3). Such co-occurrence of ectomycorrhizal and possibly endophytic fungi in single root system has been inferred from other lines of evidence [30]-[33], suggesting that specific pairs of fungi in different functional groups can interact with each other, potentially in facilitative or mutualistic ways. However, aggregated patterns can be observed between fungal species with similar physiological requirements for root environments, and the observed coexistence of possibly endophytic ascomycete fungi (Tables 2 and 3) may be attributed to shared preferences for habitats. Given that we focused on the segregated/aggregated patterns of fungi at the scale of host seedling individuals, potential vertical heterogeneity of microenvironments (e.g., soil nutrient availability) and the resultant partitioning of niches within root systems might have produced the observed aggregated patterns. Therefore, further studies are required to confirm whether the observed aggregated patterns reflect actual fungus–fungus ecological interactions in root systems. For example, experimental-inoculation methods will help to examine what kinds of mechanisms are responsible for the co-occurrence of those specific pairs of fungi.

Despite the proposition that fungus–fungus competitive or facilitative interactions could vary depending on host plant species, no significant difference in the levels of spatial segregation or aggregation was observed between Q. serrata and Q. glauca samples. The lack of host effects on the segregated or aggregated patterns may be partly attributed to the phylogenetic closeness of the two host species. Given that the root-associated fungal community composition was highly similar between the congeneric plant species (Fig. 2; cf. [74]), the environmental conditions experienced by root-associated fungi might be almost identical between the two host species. Hence, the use of congeneric host plant species in the present study may have precluded the detection of hosts' effects on fungus–fungus interactions. Therefore, the present results do not necessarily mean that host plants generally have no impact on the nature or strength of fungus–fungus interactions in roots: comparative analyses of C or T scores on phylogenetically-distant host plant species are awaited to further discuss the potential effects of host plants on the assembly processes of root-associated fungi.

As shown in this study, pyrosequencing technologies allow high-throughput profiling of root-associated fungal communities in hundreds (or more) of host individuals and will offer a breakthrough in the investigation of fungus–fungus interactions in roots. In addition, sequence-based taxonomic assignment (identification) potentially allows us to assess intraspecific genotypic diversity [75], providing further opportunities to infer fungus–fungus ecological/evolutionary interactions. Meanwhile, this pyrosequencing approach is based entirely on observational data on the co-occurrence of fungal OTUs in roots, and hence, it only provides insights into the potential consequences of competitive or facilitative interactions between root-associated fungal species/taxa. Therefore, complementary experimental studies are important to reveal the nature of fungus–fungus ecological interactions. For example, some experimental approaches have separated the mechanisms and consequences of interspecific interactions between fungal species [23], [24], [30] and others have quantitatively evaluated the effect of fungus–fungus interactions on host's performance [75], [76]. In combination with such experimental studies, the pyrosequencing-based high-throughput profiling of fungal communities will contribute to our understanding of the assembly processes of ecologically diverse root-associated fungi.

Supporting Information

Histograms of C and T scores that were obtained from the randomization of simulated data. In each combination of sample size (N) and sample counts of OTUs (Ri and Rj), the histogram of C or T scores of randomized data matrix were obtained. C and T scores tend to take the maximum (1) and minimum (0) values, respectively, when sample sizes (N) or sample counts of fungal OTUs (Ri and Rj) are small. Likewise, when Ri and Rj are much smaller than N, C and T scores tend to take the extreme values. Thus, all of sample size, sample counts of fungal OTUs, and the balance between them should be carefully inspected when screening pairs of fungal OTUs prior to randomization tests of C or T scores