Low abundance of Archaeorhizomycetes among fungi in soil metatranscriptomes

The Archaeorhizomycetes are recently discovered fungi with poorly resolved ecology. Even their abundance in soil fungal communities is currently disputed. Here we applied a PCR-independent, RNA-based metatranscriptomic approach to determine their abundance among fungi in eleven different soils across Europe. Using small subunit (SSU) ribosomal RNA transcripts as marker, we detected Archaeorhizomycetes in 17 out of 28 soil metatranscriptomes. They had average relative SSU rRNA abundance of 2.0% with a maximum of 9.4% among fungal SSU rRNAs. Network analysis revealed that they co-occur with arbuscular mycorrhizal Glomerales, which is in line with their previously suggested association with plant roots. Moreover, Archaeorhizomycetes ranked among the potential keystone taxa. This metatranscriptomic survey exemplifies the usage of non-targeted molecular approaches for the study of soil fungi. It provides PCR- and DNA-independent evidence for the low abundance of Archaeorhizomycetes in soil fungal communities, although they might be non-negligible players despite their low abundance.

Archaeorhizomycetes SSU rRNAs were detected in 17 out of 28 metatranscriptomes and in 8 out of 11 soils ( Fig. 1 and Table 1). The three sites where no Archaeorhizomycetes were detected were predominantly hypoxic, i.e. peatland and mofette soils. These metatranscriptomes had on average the lowest number of fungal reads, thus the detection limit could have been too high for Archaeorhizomycetes.
The overall relative abundance of Archaeorhizomycetes transcripts across all samples was 2.0% ( Fig. 1). Exploring only metatranscriptomes containing Archaeorhizomycetes, this class represented on average 3.3% of SSU rRNAs. The maximum proportion recorded was 9.4%. Thus, these PCR-independent metatranscriptome results support the view that Archaeorhizomycetes are rather a low abundant fungal class in soils and that reports of their high abundance might be derived from a preferential amplification of their ITS or rRNA genes in PCR assays 13 . However, a particularly high abundance of Archaeorhizomycetes in PCR based studies was reported in alpine tundra 14 and coniferous forest soils 3 , which were not included in our study. Furthermore, metatranscriptomics is not unbiased and shares some biases with DNA methods, such as protocol-dependent preferential nucleic acid extraction for certain taxonomic groups.
The relative abundance of Archaeorhizomycetes varied strongly between sites, with no SSU rRNA transcripts being detected in the predominantly hypoxic arctic peatland (PsS, PsK) and mofette (MO) soils (Fig. 1). Common to these soils was the low abundance of vascular plant roots, supporting once more the association of Archaeorhizomycetes with plant roots, although an obligate aerobic lifestyle could also result in such a distribution pattern. Remarkably, the highest relative abundance of Archaeorhizomycetes (average 5.6%) was found in soils of a former lead and zinc mining site (Fig. 1); with increasing heavy metal concentrations (MiL < MiM < MiH) having no effect on their relative abundance. These data indicate that Archaeorhizomycetes are tolerant to high concentrations of these metals, adding one additional piece of knowledge to their autecology.
In the other surveyed sites, Archaeorhizomycetes abundance was minimal -up to ~1% for grassland (MR, RS, RB) and negligible < 0.01% for beech forest (FL, FS) soils (Fig. 1). Exploring seasonal patterns in the abundance of Archaeorhizomycetes SSU rRNAs, we observed a high relative abundance in samples collected in spring, which is in agreement with a study conducted in a Colorado alpine tundra that found a peak in abundance of Archaeorhizomycetes during spring 14 . Nevertheless, this conclusion is rather speculative as there was no seasonal sampling done for any of the soils surveyed in our study. To provide a well-founded description of the Archaeorhizomycetes seasonal variability more research is needed.
To identify possible interactions of Archaeorhizomycetes with other fungi we conducted network analysis (Fig. 2). We analysed the co-occurrence patterns of the 52 most abundant fungal orders in the metatranscriptomes using the CoNET algorithm 15 in Cytoscape 3.0.2 16 (for details see Methods section). Archaeorhizomycetes had a strong co-presence relationship with the arbuscular mycorrhizal Glomerales (Fig. 2; Table 2). However, we are not able to resolve whether both fungal groups interact directly (share of mycorrhizal niche or parasitism) or only share the affinity to root environment without explicit relationship. Also other positively interacting fungal orders have representatives with strong association to plants as either pathogens or endophytes ( Table 2). The results of network analysis therefore further supported the hypothesis that Archaeorhizomycetes prefer to live in close vicinity of plant roots.

Figure 1. Relative abundance of the 16 most abundant fungal classes (mean + SE) as average of all datasets (black columns) and in respective soils (open columns).
Co-occurrence analysis placed Archaeorhizomycetes among the top 10 potential keystone taxa (Table 3), thus Archaeorhizomycetes seem to play a non-negligible role in soil fungal communities despite their low relative abundance.
This study concludes Archaeorhizomycetes as fungal class with minor representation in transcriptionally active part of soil fungal communities in different terrestrial biomes. Yet, it highlights the need for more studies to elucidate more aspects of their ecology, as Archaeorhizomycetes were identified as putative key players in soil fungal communities.

Methods
Metatranscriptome sequences were filtered using PRINSEQ 17 discarding short and low quality sequences (< 150 bp, min average quality < 0.25; exception dataset RB: < 80 bp, min average quality < 0.25). Metatranscriptomes generated with Illumina HiSeq were subsampled to 1 million reads and screened for eukaryotic SSU rRNA transcripts using SortMeRNA 18 with default settings and the Silva SSUref 111 database 19 . Metatranscriptomes generated with 454 pyrosequencing were not subsampled and not subjected to SortMeRNA. For fungal community profiling, whole metatranscriptome datasets (454 pyrosequencing) or eukaryotic SSU Figure 2. Co-occurrence network of significantly interacting fungal orders. Interactions with Archaeorhizomycetales are highlighted: positively (co-occurrence) interacting fungal orders are connected with green lines, negatively (mutual exclusion) interacting with red lines. The thickness of lines is proportional to significance of the interaction (q-value). The size of circle is proportional to the average relative abundance of fungal order in all datasets.   15 . Only those fungal orders whose sum of sequences in all datasets was higher than 200 were included into the network analyses. The parameters and settings for network analyses in CoNET algorithm were: -parent_child_exclusion, -row_minocc 10, -correlations (Spearman, Pearson, mutual information, Bray-Curtis dissimilatory and Kullback-Leibler dissimilatory), threshold for edge selection was set to 1000 top and bottom. In randomization step, 100 iterations were calculated for edge scores. In following bootstrap step 100 iterations were calculated, unstable edges were filtered out, Brown method was chosen as P-value merging method and Benjaminihochberg method for multiple test correction. Resulting network for fungal community was visualized and analysed (i.e. degree of nodes, betweenness centrality, closeness centrality) in Cytoscape 3.0.2 and nodes tending to have high mean degree, low betweenness centrality, and high closeness centrality were identified as potential keystone orders 22 .  Table 3. Potential keystone orders, their assignment to phyla, degree, betweenness centrality and closeness centrality. Orders are sorted according their degree (i.e. number of direct interactions with other fungal orders).