The role of fungi in heterogeneous sediment microbial networks

While prokaryote community diversity and function have been extensively studied in soils and sediments, the functional role of fungi, despite their huge diversity, is widely unexplored. Several studies have, nonetheless, revealed the importance of fungi in provisioning services to prokaryote communities. Here, we hypothesise that the fungal community plays a key role in coordinating entire microbial communities by controlling the structure of functional networks in sediment. We selected a sediment environment with high niche diversity due to prevalent macrofaunal bioturbation, namely intertidal mangrove sediment, and explored the assembly of bacteria, archaea and fungi in different sediment niches, which we characterised by biogeochemical analysis, around the burrow of a herbivorous crab. We detected a high level of heterogeneity in sediment biogeochemical conditions, and diverse niches harboured distinct communities of bacteria, fungi and archaea. Saprotrophic fungi were a pivotal component of microbial networks throughout and we invariably found fungi to act as keystone species in all the examined niches and possibly acting synergistically with other environmental variables to determine the overall microbial community structure. In consideration of the importance of microbial-based nutrient cycling on overall sediment ecosystem functioning, we underline that the fungal microbiome and its role in the functional interactome cannot be overlooked.

. (a) The sesarmid crab Neosarmatium africanum and (b) the typical burrow structure with a hood at the burrow entrance. (c) Sampling design adopted to dissect the microbial structure at three depths and at five distances away from the burrow wall. (d-f ) qPCR for the total bacterial, archaeal and fungal communities. All 16S rRNA and ITS1 copies are normalized per g of sediment.
Identifying discriminately abundant taxa throughout the burrow wall. Linear discriminant analysis effect size (LEfSe) was used to identify taxa that were significantly more abundant at each depth and distance from the burrow wall ('Fraction').
Keystone nodes analysis (Fig. 3c-e) revealed that fungi in each 'Fraction' and 'Depth' showed a significantly higher level of radiality (χ2 deviance 16,10245 = 154.75, P < 0.01), eigenvector (χ2 deviance 16, 10153 = 258.06, P < 0.01) and number of directed edges (χ2 deviance 16,10245 = 793459, P < 0.01). This was confirmed by the significance of the relationship between betweenness centrality and degree of connection (Fig. 3b), where fungi consistently had the highest level of degree of connection with the highest level of betweenness centrality compared to bacteria and archaea (GAM; F 8,10316 = 132.2; P < 0.01). Fungal OTUs had the highest degree of connectivity and formed the majority of the keystone nodes in surface, subsurface and deep sediment (having an overall highest degree of connectivity, betweenness centrality and closeness centrality; Table 1). In the surface, in Fraction 1, the largest keystone node was the saprotrophic fungal family Trichocomaceae followed by the bacterial keystone node Pelobacteraceae. Fungi also formed the majority of keystone nodes in Fractions 2 to 4 (Sordariomycetes being a consistently highly connected node). The green non-sulphur bacteria Thermomicrobia and the anaerobic sulphate-reducing bacteria Desulfobulbaceae formed important nodes in Fractions 2 and 4, respectively (Table 1). In the subsurface (Table 1), fungi such as Pluteus (a wood-rotting saprobe) formed significant keystone nodes throughout the burrow wall. Large bacterial nodes included the purple sulphur bacteria Chromatiales (the second most important node in Fraction 2) and the ammonia oxidizing bacteria Candidatus Nitrososphaera. In the deep (Table 1), fungi formed the majority of the most important keystone nodes. In Fraction 1 these included OTUs belonging to the orders Sordariomycetes and Dothideomycetes. The bacteria Piscirickettsiaceae was also a significant node in Fraction 1 at the burrow wall. In Fraction 2 and bulk sediment, several keystone nodes were formed by members of the archaeal order Crenarchaeota. A bacterial OTU belonging to the family Desulfobacteraceae was a highly connected central node in Fraction 3. In Fraction 4, a member of the strictly aerobic bacteria family Hyphomonadaceae formed an important node.

Discussion
Acting as organic matter traps, burrows belonging to benthic herbivores of mangroves fuel microbial activity, promote niche diversification and habitat heterogeneity. As of yet we know little about mangrove sediment fungi other than their high diversity 23,24 , here we show that in these heterogeneous zones, fungi invariably play the major role in microbial network interactions from the internal burrow surfaces to the undisturbed sediments far from the burrow. Despite the varied sediment conditions from the burrow wall towards the bulk, and from the surface to the deep, the importance of fungi in the inter-kingdom networks was notable; fungal nodes were highly connected in every fraction of sediment (with the exception of deep bulk sediment), revealing a persistent central role of fungi in shaping the structure of the network (Table 1). In particular, the relationships among betweenness centrality and degree of connection (Fig. 3b) show that fungi persistently formed the keystone nodes essential for the interactome topology having, simultaneously, the highest levels of both of these network features. On the contrary, archaea and bacteria retained a high betweenness centrality but had a lower level of degree of connection revealing only a local effect on the overall topology of the network. These results are corroborated by other studies that observed fungi, as keystone species, to stabilize the network properties in other complex systems such as host-microbiome interactions 25 . Moreover, fungi showed a constant significant radiality compared to archaea and bacteria, which can be interpreted as a higher probability of an OTU to be functionally relevant for several other OTUs; eigenvector, that allows an immediate and informative evaluation of the interacting relevance of the OTU with the rest of the network, and the number of direct edges, indicating how the influence of each fungal node can extend to the rest of the network. Many studies have described fungi as superhighways, supporting their central role, that can enhance below ground communication, increasing the networking capabilities of the sediment microbiome and also networking with plants 26 . Certain members can establish long hyphae that boost the growth of microbial biofilms that are involved in the protection of the hyphae itself and in enhancing fungal communication 27 .
N. africanum is a highly efficient herbivore that removes up to 79% of mangrove ground litter 28,29 , storing more than half of this in their burrows 30 allowing fungi to proliferate thereby increasing the palatability and nutritional value of mangrove litter 31 (Supplementary File 5). Acting as hotspots of microbial activity, burrow walls were significantly enriched in POC and, correspondingly, higher numbers of bacteria, archaea and fungi.  Table 1. Keystones nodes obtained from the network analysis for surface (A), subsurface (B) and deep (C) sediment across each 'Fraction' considering the highest level of node degree distribution, closeness centrality and betweenness centrality. The lowest taxonomic resolution of the hub is shown.
While POC availability was found to be a significant driver of bacterial and archaeal community composition, instead fungal community composition was unaffected. This may be explained by the broad metabolic plasticity of fungi 32 and their dominance in the decomposition role in intertidal ecosystems 10,14 . At all depths, saprotrophic fungal nodes key to the microbial network were detected (e.g. the genus Pluteus) and fungi were the dominant network components throughout burrow sediment. By accumulating leaves in the walls of their burrows, N. africanum can facilitate and enrich the sediment fungal community with that of the leaves. Recently, the www.nature.com/scientificreports www.nature.com/scientificreports/ ability of fungi in the phyllosphere of the fallen leaves in a temperate beech forest to nurture the soil was demonstrated, by supplying a constant input of new fungal strains that establish in the forest litter and are involved in the litter decomposition process 33 . Even though we have not evaluated this aspect in the present study, we speculate that the burial of leaves by N. africanum can provide a constant input of fungal species to the burrow sediment, contributing to support the centrality of fungal taxa that shape the interactions and the assemblage of the entire sediment microbiome.
Although burrowing allows oxygen to penetrate into subsurface sediment, this oxygen is rapidly consumed by sediment microorganisms and organic-rich sediments characteristically have narrow zones of oxygen available as an electron acceptor. Indeed, anaerobic taxa were detected in the first layer of the burrow wall. Cellulose degradation in anaerobic sediment requires the interaction of microorganisms with diverse metabolism, with 5-10% degraded under anaerobic conditions by cellulose fermenting microbes to CH 4 , CO 2 and H 2 O 34 . Fungi may have an important role in decomposition of organic matter in anoxic sediment that provide substrates to fermentative bacteria, such as the anaerobic bacterium Pelobacteraceae, which feed the archaeal methanogenic community 35 . Moreover, saprotrophic fungi from the Basidiomycetes are known to produce methane under oxic conditions 36 , promoting co-occurrence with methanotrophic microbes. While this study did not aim to resolve the steep gradients of carbon flow and electron acceptors in sediment, we detected a larger scale pattern around the burrow.
Burrows of intertidal benthic fauna are also hotspots of organic nitrogen, enhancing nitrogen cycling processes such as nitrification and denitrification due to the increase in oxic sediment zones created by burrow walls 37,38 . Fungi are predominantly aerobic heterotrophs, however they are also metabolically capable of utilizing nitrate and/or nitrite, thus they are key players in anaerobic denitrification in intertidal sediments with low oxygen availability 39 . While denitrifying bacteria are obligate anaerobes and dominate under strongly reducing sediment conditions, denitrifying fungi can make use of sub-oxic conditions (300-900 uM O 2 39,40 ). Denitrification by fungi under anaerobic conditions has particular ecological impacts, since rather than producing N 2 as the final respiratory product, as do bacteria, the process ends with nitrous oxide, which is a significant greenhouse gas 41 . One of the functions to which bacteria were assigned in our study was nitrous oxide respiration, which may indicate a positive interaction between fungi and bacteria in N. africanum burrow sediment. Ingesting double the amount of litter they are able to assimilate, sesarmid crabs deposit highly nitrogen-rich faeces 42 . The availability of nitrate was a significant driver of the bacterial community, while PON was a significant driver of archaeal community composition. Nitrite influenced both archaeal and bacterial communities, while fungi were independent of both nitrite and nitrate. Due to the low availability of nitrate at the burrow wall (which increased away from the burrow) at every depth, we hypothesise that nitrate is rapidly consumed by the burrow wall microbiome and/or readily released into burrow water and flushed out at high tide 43 .
The high microbial community variability and turnover we observed in different sediment fractions, but also the patterns of sharing and abundance of certain taxa, may be explained by the large volume of sediment excavated from burrows by the host onto surface sediment around the burrow (sesarmid crabs are estimated to excavate in the range of 80-210 cm 3 m −2 d −1 30 ). Fungal diversity was higher in surface sediment unaffected by excavated material (in particular the orders Diaporthales, Sordariomycetes and Agaricomycetes), as was the abundance of photosynthetic bacteria such as Cyanobacteria and purple non-sulphur bacteria. In sediment closest to the burrow funnel and on the surrounding surface, we detected a number of sulphur compound reducing bacteria, indicating sulphur-compound rich and predominantly anaerobic sediment. This, along with the occurrence of similar taxa in Fraction 1 of the subsurface and surface of the burrow, can be explained by uni-directional excavation of burrow sediment, and thus microbes, and burrow maintenance. We detected a larger number of microbial OTUs to be shared between the surface, subsurface and deep sediment in the sediment adjacent to the burrow wall. Overall, however, deep sediment was found to have more unique OTUs which were not shared with other depth levels. Indeed, the deeper parts of N. africanum burrows are known to be highly stable over generations but the burrow openings are repeatedly repaired after high tides 44 , which would identify the absence of sediment mixing from the deep and the occurrence of a larger number of unique microbial OTUs at that depth. In a recent study of the smaller mangrove fiddler crab, which displays burrow plugging behaviour at high tide, sulphate reducing taxa were less abundant and aerobic taxa were more prevalent in sediment immediately around the burrow funnel 45 . This result may reflect the different behaviours of burrow hosts. Indeed, there was no indication of augmentation of oxygen through burrow sediment in this study.
Environmental variability guides the assembly of sediment bacteria and archaea by creating a halo of sequential physio-chemical characteristics around macrofaunal burrows, but here we suggest that key species detected by the study of the interactome also have an important effect on community assembly by driving co-occurrence patterns among microbial species. Understanding this complex relationship could be of pivotal importance in the comprehension of ecosystem functioning and restoration, especially for intertidal systems 46 . Combined, environmental factors and microbial species interactions can create new niches that boost microbial community diversity, stability and resilience through functional redundance 47,48 . For example, in riverine sediment contaminated with polycyclic aromatic hydrocarbons microbial co-occurrence patterns, interactions and key species varied with the extent of environmental contamination, thus demonstrating the combination of both the environmental factors and key network species in shaping sediment microbial communities 49 and in this case revealing important microbial species for effective restoration/rehabilitation action. In other systems, such as grassland, fungal networks were found to be more resistant to changing environmental conditions in the form of drought stress than bacterial networks 50 , indicating that fungi may be more resilient than bacteria under environmental variability which could be essential for ecosystem resilience. (2019) 9:7537 | https://doi.org/10.1038/s41598-019-43980-3 www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusions
Through their bioturbation activities and trapping mangrove leaves and organic matter inside their burrows, benthic fauna increase habitat heterogeneity and functional microbial diversity. At all the sediments depths and radial fractions that we investigated, the large proportion of the individuated keystone taxa were fungi ( Fig. 3; Table 1), and only few archaea and bacteria provided central roles, but without consistent patterns. Our results indicate the synergistic contribution of fungi and environmental variability in structuring the microbial communities and their interactome throughout intertidal sediment. Since microbial networking has consequences on rates of organic matter decomposition and nutrient and blue carbon cycling processes, we highlight the necessity of a wholly inclusive approach in which the picture of the fungal microbiome and its role in the functional interactome cannot be overlooked.

Materials and Methods
study site, study species, and sampling design. Sampling was carried out during April 2016 in a mixed Avicennia marina and Rhizophora mucronata riverine mangrove stand at the mouth of the temperate Mngazana estuary, South Africa (31°42′S, 29°25E). Sampling was performed 1 h before low tide, when burrows were uncovered by water, during the period of Spring Tide. We randomly sampled eight similar sized active N. africanum burrows along a 400 m transect, distinguished by fresh excavations and observation of the animal inside the burrow (Fig. 1a,b). Sediment sampling comprised of five sediment fractions along a horizontal gradient at the burrow wall (0-1 cm), mid-way along the excavated burrow content (13-14 cm), at the distal point of the excavated burrow content (26-27 cm) and a point away from the burrow (39-40 cm), which was repeated at three depths (surface: 0-1 cm, subsurface: 1-2 cm and deep: 20-21 cm; Fig. 1c). Bulk sediment was sampled at the three depths, >3 m away from the burrow and any visible root or pneumatophore or any other sign of bioturbation. From each burrow, we took 15 samples (total of 120 samples) stored on ice in the field and frozen at −20 °C within 3 h of collection.
DNA sequencing, metabarcoding and geochemical analysis. DNA was extracted from a 0.4 g sub-sample of each sediment sample using the PowerSoil Total DNA Isolation Kit (MoBio Inc., CA, USA) following the manufacturer's instructions. To identify bacterial and archaeal community composition we targeted the V4-V5 hypervariable region of the 16S rRNA gene using the primer pairs 341F-785R and 519F-806R respectively 51 . To identify fungal community composition, we amplified the ITS2 region of the internal transcribed spacer region using the primer pair ITS3-ITS4 52,53 . All libraries were prepared using the 96 Nextera XT Index Kit (Illumina ® ) and sequenced using the Illumina ® MiSeq platform with pair-end sequencing at the Bioscience Core Lab, King Abdullah University of Science and Technology. Details on raw read processing are provided in the Supplementary File 1. QIIME was used to assign taxonomy, using the Greengenes database for bacteria and archaea and the UNITE database for fungi 54,55 . Copies of the 16S small-subunit rRNA gene for bacteria and archaea and the ITS region for fungi were quantified following the qPCR protocol described by Fierer and Jackson 56 using the primer pairs Eub338-Eub518, Arc931F-Arc1100R and ITS1F-5.8 s for bacteria, archaea and fungi, respectively (see Supplementary File 1 for detail).
We tested differences in beta-diversity of bacteria, archaea and fungi amongst 'Depth' and 'Fraction' using a multivariate generalized linear model (GLM) with a negative binomial error distribution, in the R package "mvabund" 59 . In this test, and those following, we treated the burrow as a random effect by using the function 'manyany' where it is possible to specify the random factor. We used a GLM to test differences in bacterial, archaeal and fungal gene copies, using a log transformation for normality, among 'Depth' and 'Fraction' . OTU assemblages across 'Depth' were explored using Canonical Analysis of Principal coordinates (CAP 60 ). The Shannon diversity index and OTU richness were calculated using the function 'diverse' in PRIMER (v. 6.1) and 2-way PERMANOVA was used to test differences in diversity and richness among the factors 'Depth' and 'Fraction' .
Ternary plots, using the mean relative abundances of bacteria, archaea and fungi OTUs in surface, subsurface and deep sediment, were created using the R package "ggtern" 61 to examine the number of OTUs shared amongst different depths. The Sloan model 62 was used to test the homogenizing effect of sesarmid activity on the bacterial, archaeal and fungal communities in each 'Fraction' of the burrow. We determined if the structure of a given community fits a neutral assembly model where the abundances of taxa are driven by dispersal from a source community. For example, to test whether there is significant vertical mixing downwards, we examined if the abundances of the taxa in the subsurface or deep communities are driven by dispersal from the surface or subsurface communities, respectively. To identify taxa that were discriminately more abundant in each 'Fraction' at each 'Depth' (Wilcoxon P value: 0.05, LDA > 2), we used Linear discriminant analysis effect size (LEfSe, www. huttenhower.sph.harvard.edu/galaxy/) following Segata et al. 63  www.nature.com/scientificreports www.nature.com/scientificreports/ To identify co-existing or mutually exclusive OTUs amongst bacteria, archaea and fungi, an inter-kingdom co-occurrence network was built following Agler et al. 64 using the routine CoNet in Cytoscape 3.2.1 65 . Cytoscape was used to calculate the network topological parameters 66 and Gephi 1.9 to compute network modularity and visualization 67,68 . Network topological coefficients were calculated using the cytoscape plug-in Centiscape 69 and Network analyser 66 . Further details on the network analysis are provided in Supplementary File 1. Network centrality measures were compared across 'Depth' and 'Fraction' adding the explanatory variable 'Kingdom' (3 levels: archaea, bacteria, fungi; fixed and orthogonal) and 3-way analysis of variance (ANOVA) was used to test the effect of 'Kingdom' , 'Depth' and 'Fraction' on degree of connection and average path length. A GLM with a quasibinomial error distribution was used to test the effect of 'Kingdom' , 'Depth' and 'Fraction' on closeness, betweenness centrality measures, radiality, eigenvector and number of directed edges using the R package "MASS" 58,70 .
Keystone species were analysed by testing the relationship between degree of connection and betweenness centrality following de Vries 50 using a Generalized Additive Model included in the r packge mgcv 71 . Keystone species were detected by combining the highest level of three centrality measures: degree of connection, closeness centrality and betweenness centrality following Berry and Widder 72 .
After testing for multi-collinearity (see Supplementary File 1 for further detail), 2-way PERMANOVA (9999 permutations, Euclidean distance) was used to test differences in biochemistry, metals and grain size for 'Depth' and 'Fraction' . The contribution of different geochemical variables, metals and grain size classes to dissimilarity between each 'Depth' and 'Fraction' was assessed by SIMPER analysis. Grain size frequencies were calculated using the R package "G2sd" 73 . Distance-based multivariate analysis for a linear model 57 was used to determine the geochemical variables, metals and grain size classes that significantly explained the variation community composition amongst 'Kingdom' (using the corrected Akaike information criterion (AICc) to test significance 74 ). Differences in sediment phi were tested for 'Fraction' and 'Depth' with a 2-way ANOVA.
The FAPROTAX database was used to assign bacterial and archaeal OTUs to known metabolic or ecological functions (http://www.zoology.ubc.ca/louca/FAPROTAX75). We used the annotation tool FUNGuild v1.0 9 to assign functional guild, in terms of trophic mode, to fungal OTUs. 2-way PERMANOVA was used to test differences in functional group assignment amongst 'Depth' and 'Fraction' for each 'Kingdom' . The contribution of functional group assignment to dissimilarity was assessed by SIMPER analysis.