Pacific Biosciences Single-molecule Real-time (SMRT) Sequencing Reveals High Diversity of Basal Fungal Lineages and Stochastic Processes Controlled Fungal Community Assembly in Mangrove Sediments

doi:10.21203/rs.3.rs-97364/v1

Download PDF

Research

Pacific Biosciences Single-molecule Real-time (SMRT) Sequencing Reveals High Diversity of Basal Fungal Lineages and Stochastic Processes Controlled Fungal Community Assembly in Mangrove Sediments

https://doi.org/10.21203/rs.3.rs-97364/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Mangrove wetlands are unique ecosystems with specific environmental characteristics, and are a hotspot of biodiversity. Although they probably harbor a variety of mangrove-specific fungi, the compositions of mangrove fungal community has been rarely investigated in detail, except for few published culture-based studies. In addition, the fungal community assembly and interaction patterns that impact the community composition in mangroves have not been explored to date.

Results: We used the Pacific Biosciences single-molecule real-time sequencing approach targeting the entire internal transcribed spacer region, to systematically investigate the composition, biogeographical patterns, assembly processes, co-occurrence patterns and shaping factors of the fungal communities in sediments of seven representative mangroves across the Southeast China. We recovered 15 phyla, including some early diverging fungal lineages not previously reported in mangroves. Phylogenetic analysis revealed an incredibly high proportion of Rozellomycota and Chytridiomycota, as accounting for up to one-third of all fungal abundance. Although the neutral community model described a moderate portion of community variation, the similarity of fungal communities exhibited strong distance-decay patterns. Furthermore, the mean values and most beta nearest-taxon index fell between -2 to 2, with Bray–Curtis-based Raup–Crick value generally greater than 0.95, suggesting that stochastic processes strongly shape the fungal community composition. Consistently, nonmetric multidimensional scaling and permutational multivariate analysis of variance confirmed the geographical location as a crucial factor driving the distribution of both, the dominant and rare taxa of mangrove fungi. The db-RDA analyses indicated a minor role of environmental selections in shaping the fungal community. Network analyses revealed that the deep sediments harbor more complex fungal networks with highly connected taxa than surface sediments, and that rare fungal taxa might play important roles in microbial interactions and ecological functions in mangrove sediments.

Conclusions: The investigation revealed high fungal diversity in mangrove sediments, with incredibly high numbers of basal fungal lineages, stochastic processes driving the assembly of fungal community, and geographic location strongly shaping fungal community composition in mangroves. These discoveries therefore spur further studies of the utilization and protection of fungal resources and communities in mangrove sediments.

General Microbiology

Fungal community

PacBio SMRT

Neutral community model

Stochastic process

Community assembly

Co-occurrence network

Mangrove sediment

Mangrove trees are unique forest biotypes widely distributed along the coast in tropical and subtropical regions, providing ecosystem services, such as protection of the coastal area, and provision of food and shelter for fish and shellfish [1, 2]. Despite covering just 0.5–0.7% of the tropical forest area, these forests account for 11% of the total terrestrial carbon and 10% of the terrestrial dissolved organic carbon exported to the ocean [1, 3]. Therefore, mangroves have been labeled as the “blue carbon sink” [1, 2, 4]. Specific environmental conditions, such as high levels of nutrients, high salinity and low oxygen, render the mangrove ecosystem a hotspot of biodiversity. These conditions support a variety of unique microorganisms with abundant and diverse metabolic pathways, driving complex nutrient and biogeochemical cycling, such as carbon cycling, ammonia oxidation, sulfate reduction and methane metabolism [2, 4–9].

Fungi are the key components of microbial community in mangrove sediments. They play a major role in the transformation of organic matter and greatly contribute to the supply of nutrients and energy to the mangrove ecosystem [10]. Early studies of mangrove fungi mostly involved traditional cultivation techniques. Using these techniques, over 625 fungal species have been obtained from mangroves, most of which are ascomycetes and basidiomycetes [11]. Many of these fungi produce useful enzymes or metabolites, such as extracellular hydrolases for the degradation of lignin and cellulose [4, 12, 13]; peptides and lactones with antimicrobial activities [14]; and organic acids and phosphatases for the solubilization of phosphate [15]. Fungal communities in the mangrove sediments are significantly affected by the environmental conditions, mainly biotic and abiotic factors. For example, the age and species of mangrove plant are the two most important biotic factors that shape fungal community compositions [16, 17]. The abiotic factors include, rainfall [18], sediment depth [16], and physicochemical properties of the sediment, such as salinity, pH, concentration of nutrients, etc. [19, 20].

Considering the limitations of traditional approaches, the mangrove mycobiota have not been extensively studied to date [21]. As an alternative, high-throughput sequencing (HTS) methods, represented by metabarcoding, metagenomics, and metatranscriptomics, provide ways to effectively explore the mangrove microbial community, with the potential for recovery of previously undescribed rare populations and biochemical pathways that are otherwise difficult to elucidate [22, 23]. Nonetheless, studies investigating fungal communities in mangrove sediments are rare. Using 454 pyrosequencing, Arfi et al. [24, 25] investigated the fungal diversity and composition in the mangrove sediment in New Caledonia, and revealed a predominance of Ascomycota and Basidiomycota. By analyzing metagenomics datasets, a similar fungal composition was determined in the gray mangrove sediment in Red Sea, in which fungal communities are more stable within the rhizosphere than in the bulk soil [26]. Recently, Luis et al. [27] used HTS to investigate the prokaryotic and fungal community during the wet and dry seasons in New Caledonian mangrove sediments. The data suggest that the prokaryotic communities are mainly shaped by sediment depth, with the fungal community almost evenly distributed over different depths, vegetation covers and seasons. Since most fungi from a particular environment are uncultivable [28], a significant portion of the mangrove fungal diversity is overlooked using culture-dependent methods [27, 29]. In addition, only single mangrove and small numbers of samples have been analyzed previously using culture-independent methods, and fail to provide an overview of the entire fungal community in mangrove sediments.

Second-generation HTS, represented by Illumina, the most widely used sequencing platform, rapidly generates hundreds of thousands of high-quality sequences per sample, greatly improving the study of microbiome. However, the generated sequence length is too short for accurate microbial identification at species level. For example, while the universal fungal DNA barcode, the internal transcribed spacer (ITS) region, typically spans 500–700 bases [30], the paired-end approach of MiSeq Illumina sequencing covers amplicons up to 550 bases long, and is only suitable for the ITS1 or ITS2 subregion (typically 250–400 bases) [30, 31]. The benefits of targeting the full ITS region include increased taxonomic resolution and suppressed amplification of material from dead organisms [30–32]. Third-generation HTS using Pacific Biosciences (PacBio) and Oxford Nanopore platforms produce long reads averaging 20–25 kb (up to 100 kb), and enable targeting the full ITS region, as well as parts of or even the entire flanking rRNA genes [30, 32]. Unfortunately, the high single-pass error rates (13–15%) limit their application in DNA metabarcoding [31, 33]. However, the circular consensus sequencing (CCS) version of PacBio single-molecule real-time (SMRT) sequencing greatly reduces the error rate by enabling sequencing of the same reads multiple times, which enables application of PacBio sequencing in metabarcoding [31, 34].

In the current study, we used PacBio SMRT sequencing to explore the fungal community composition, shaping factors, assembly and co-occurrence patterns in 110 sediment samples collected from seven representative mangroves in China (Fig. S1, Table S1). We hypothesized that (1) mangrove sediment harbors a high diversity of fungi, including a number of ancient fungal lineages; (2) stochastic processes control the assembly of fungal community; and (3) geographic location and sediment depth play significant roles in driving the composition of fungal communities in mangrove ecosystems.

Mangrove locations and sediment sampling

Mangroves cover approximately 17,800 ha in China, and are distributed in the coastal area from Hainan to Zhejiang. Seven representative mangrove nature reserves from Eastern to Southern China, i.e., Ximendao National Nature Reserve (XMD), Zhangjiangkou Mangrove National Nature Reserve (ZJK), Futian Mangrove National Nature Reserve (FT), Zhanjiang Mangrove National Nature Reserve (ZJ), Shankou Mangrove National Nature Reserve (SK), Dongzhaigang Mangrove National Nature Reserve (DZG), and Yalongwan Qingmei Estuary Mangrove Nature reserve (YLW), were selected for the current study (Fig. S1). The geographical locations of the seven nature reserves are significantly different. XMD is an artificial mangrove, located at the northern boundary at which mangroves survive. ZJK is the northernmost national nature mangrove reserve. FT is the only national nature reserve located in an urban hinterland in China. ZJ is the southernmost national mangrove reserve on the mainland of China. SK is located at the northwest of ZJ, separated from ZJ by the Leizhou Peninsula. DZG is the first mangrove nature reserve in China and located at northern Hainan Island. YLW, the southernmost mangrove in China, is a provincial nature reserve and located at southern Hainan Island.

For each mangrove, six (for all but FT) or seven (for FT) sampling sites were selected, with the same distance (> 500 m) between the adjacent sites (Fig. S1). Every sampling site was at least 1 m away from the nearest tree. At each sampling site, sediment samples were collected using a stainless steel sampler and separated into three layers, corresponding to the depths of 0–10 cm, 10–20 cm and 20–30 cm. To reduce sampling bias, three replicates at the vertices of an equilateral triangle were collected and mixed together for each site. Overall, 127 sediment samples were collected. All samples were placed in sterile zip-locked plastic bags, and transferred on ice to the laboratory, where they were stored at − 40 °C until analysis.

Environmental and spatial variables

The spatial variables included sampling depth, longitude, latitude and the numbers of tree species within 5 m of the sampling sites were documented (Table S1). The environmental variables included mean annual temperature (MAT), mean annual precipitation (MAP), salinity, pH, gravel proportion, total carbon (TC), total organic carbon (TOC), total nitrogen (TN), ammonium nitrogen (N/NH₄⁺), nitrate nitrogen (N/NO₃⁻), total phosphorus (TP) and total sulfur (TS) (Table S1). MAT and MAP data were obtained from the China Meteorological Administration (http://www.cma.gov.cn).

The Salinity, pH, N/NH₄⁺, and N/NO₃⁻ were measured in fresh sediments, and other parameters were determined in air-dried sediments. For the latter, sediments were air-dried for several days until the weights stabilized, and then ground to fine powder using a mortar and a pestle. The powder was subsequently filtered through a mesh (100 holes per inch) and weighed to calculate the gravel proportion. Coarse gravel was discarded and the fine powder was left. Salinity was measured by using an automatic compensation salinity refractometer (ATAGO CO., Japan). Sediment pH was determined using a pH meter (FE20, Mettler Toledo, Switzerland). Further, N/NO₃⁻, N/NH₄⁺, and TP levels were measured using a spectrophotometric method according to the National Environmental Protection Standards of the People’s Republic of China (http://www.mee.gov.cn/ywgz/fgbz/bz/bzwb/jcffbz/). TC, TOC, and TN were determined following the method of Zhang et al. [9]. TS was measured using the method of Zhang et al. [35].

Molecular analysis

For each sample, DNA was extracted from 300 mg of sediments using a DNeasy PowerSoil kit (Qiagen, Germany) according to the manufacturer’s instructions. The quantity and quality of the extracted DNA were examined using NanoDrop ND-2000c UV-Vis spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The DNA samples were stored at -40 °C before amplification by polymerase chain reaction (PCR).

After several preliminary trials, a two-step PCR (i.e., nested PCR) approach was selected to amplify the entire ITS region, using the primer pairs ITS9Munngs/ITS4ngs [30, 36] and ITS1F/ITS4 [37, 38] for the first and second PCR steps, respectively. For the first PCR step, the reaction volume was 30 µL; each reaction tube contained 1.5 U Takara Ex Taq HS (TaKaRa Biotechnology, Japan), 3 µL 10 × PCR Ex Taq buffer, 150 µM dNTPs, 0.12 µM of each primer, 2–10 ng genomic DNA, and water. Triplicate reactions were performed per sample to minimize the stochastic PCR effect of individual reactions. The PCR parameters were as follows: 94 °C for 5 min; 35 cycles of 94 °C for 30 s, 57 °C for 30 s, and 72 °C for 1 min 20 s; and a final elongation step at 72 °C for 10 min. The PCR products from triplicate reactions were then pooled equally, and diluted 10–40 folds to approximately the same DNA concentrations. Next, 1.5 µL of the resulting mixture was used as the template for the second PCR step. The amplification program was the same as that for the first PCR step, except that the annealing stage was performed at 58 °C for 30 s. Meanwhile, unique official 16-nt barcodes were added to the 5’-end of the primer pair ITS1F/ITS4 to enable multiplexing of samples within a single sequencing run.

Amplicon library building and SMRT sequencing were performed by the Annoroad Gene Technology Co. Ltd. (Beijing, China). Briefly, equimolar amounts of barcoded PCR products were pooled. Then, the hairpin sequencing adapters (SMRTbell template, SMRTbell Express Template Prep Kit 2.0) were ligated with the pooled amplicon libraries following the blunt-end ligation protocol of PacBio. Sequencing libraries were subsequently purified using the Enzyme Clean Up kit (Pacific Biosciences). The libraries were finalized for sequencing by annealing the sequencing primers to the SMRTbells and binding of the DNA polymerase to the template complex. The sequences were generated using the PacBio Sequal system.

Bioinformatics and phylogenetic analysis

CCS were extracted from the raw Sequal data using pbccs v4.02 software (https://github.com/PacificBiosciences/ccs) with default parameters except for -min-passes = 5, and then converted to fastq format using BAM2fastx tool (https://github.com/pacificbiosciences/bam2fastx/). CCS reads were demultiplexed based on the barcode-primer sequences, allowing for 0.1 mismatch using flexbar [39]. Flexbar was set to trim the barcode sequences from quality reads upon demultiplexing. To avoid unwanted multi-primer artifacts, reads with a full-length sequencing primer detected within them were removed [32]. The following quality-filter parameters were implemented in Mothur v1.44.1 [40]: qwindowaverage = 20, qwindowsize = 50, minlength = 100, maxambig = 0, maxhomop = 12. Full ITS region was extracted from filtered reads by targeting the identifiers of all eukaryotes based on hidden Markov models (HMMs) using ITSx v1.1.2 [41]. Putative chimeric sequences were then detected using uchime [42] de novo and reference-based (with UCHIME reference dataset as a reference, https://unite.ut.ee/repository.php) and removed using vsearch v2.15 [43]. After chimeras checking, filtered reads ranging from 300 bp to 900 bp in length [44] were selected using vsearch v2.15. The remaining sequences were de-replicated, sorted and clustered into operational taxonomic units (OTUs) at a 97% similarity threshold using USEARCH [45, 46]. Representative sequences were chosen using the abundance method [32].

Taxonomy was assigned to the representative OTU sequences by searching the UNITE + INSDC databases of all eukaryotes (version released on 3 February 2020) [47] using BLASTn [48]. The rules for robust taxonomic annotation followed those of Tedersoo et al. [32]. Briefly, the species, genus, family, order, class, and phylum levels were roughly matched based on 97%, 90%, 85%, 80%, 75% and 70% sequence similarity, respectively. For each OTU, 10 best BLAST matches were examined; in case of no specific assignment, 100 best matches were manually checked. However, since some sequences in the UNITE + INSDC databases for all eukaryotes are denoted as “Eukaryota_kgd_Incertae_sedis”, many OTUs could not be accurately assigned at the kingdom level after manually checking. To minimize this phenomenon, preliminary phylogenetic identification was used. Briefly, a phylogenetic tree was constructed in IQtree [49] using all representative sequences, and unassigned OTUs that clustered with the assigned fungal OTUs in a single clade were annotated as “Fungi_candidate” at kingdom level, and “Unidentified” at lower taxonomic levels (File S1). Subsequently, a phylogenetic tree comprising all the representative fungal sequences was constructed for more accurate assignments of fungal OTUs at phylum level. Several distinct clades were generated in the fungal phylogenetic tree (Fig. 1), representing different phyla, and OTUs clustered in each clade were thus assigned to the corresponding phylum but “Unidentified” at lower taxonomic levels. Representative sequences of all eukaryotic and fungal OTUs were aligned using MUSCLE [50] and trimmed using trimAl v1.2 [51], respectively. The maximum likelihood (ML) analyses were implemented using IQtree with 1000 replicates under the TIM + F + R10 model for both all eukaryotic and fungal sequences phylogeny, as selected by ModelFinder in IQtree according to Bayesian information criterion (BIC). The final consensus tree was visualized using iTol [52].

In the current study, the relative abundance thresholds were defined as 0.01% for rare taxa and 1% for dominant taxa, and all OTUs were classified into three categories: dominant taxa (DT), always rare taxa (ART), and conditionally rare taxa (CRT), according to the most recent publications [53]. DT were defined as the OTUs with the abundance of ≥ 0.01% in all samples or ≥ 1% in some samples; ART were defined as taxa with the abundance of < 0.01% in all samples; and CRT were defined as taxa with the abundance of < 1% in all samples and < 0.01% in some samples. However, only one OTU was classed as ART. To avoid confusion, ART and CRT were artificially combined into the rare taxa (RT) for further analyses.

Statistical analyses

Difference and indicator taxa of fungal communities

Fungal OTUs in the samples of 0–10 cm and other deeper layers were compared using edgeR [54] and limma package [55]. Multipatt command in the indicspecies package [56] was used to explore the indicator species in different sample depths. IndVal, the indicator value, was calculated to assess the extent to which a species was an indicator of a given habitat.

Alpha- and beta-diversity analysis

Rarefaction curves of all samples and different groups were calculated using rarecurve command in the vegan package [57]. The relative abundance of an individual taxon (e.g., phylum) within a sample was determined by comparing the read number of each taxon to the total reads number in that sample. α-diversity indices, such as OTU richness, Shannon-Wiener, Chao1, and evenness indices were calculated using vegan package [57]. The dataset was rarefied to the lowest reads number (1833) in the samples before the analysis of α-diversity and β-diversity. One-way analysis of variance (ANOVA) followed by Tukey’s honest significant difference (HSD) test was used to explore the variations of α-diversity across different group strategies. Pearson’s correlation coefficients and p-value were calculated to explore the associations between α-diversity, main fungal phyla and classes and environmental features. Prior to these analyses, the OTU richness and environmental feature values were transformed by log (x + 1).

The distance matrix of the fungal community (Hellinger transformation of the OTU abundance data) was constructed by calculating dissimilarity using the Bray˗Curtis method [57]. Fungal community compositions of all taxa, DT and RT in all samples and different groups of samples were ordinated using nonmetric multidimensional scaling (NMDS) using the metaMDS command in the vegan package [57]. Significant differences (p-value < 0.05) between groups were evaluated using the analysis of similarities (ANOSIM).

Fungal community assembly patterns

To determine the role played by stochastic processes in microbial community assembly, neutral community model (NCM) was used to predict the relationship between OTU detection frequency and relative abundance in the mycobiome [58]. The NCM used herein diverged from the neutral theory, which derives its name from the defining assumption of equivalent per-capita growth, death, and dispersal rates of species, and thus assumes that the ecological fitness of species is “neutral” [59, 60]. Briefly, because of the dispersal of dominant taxa at different sampling sites by chance, the model predicts that the dominant taxa are widespread, and that the rare taxa are more likely to be lost at different sites because of the ecological drift [53]. In the current study, the datasets for all and individual mangroves were analyzed separately. The sncm.fit_function.r code written by Burns et al. [60] was used to explore the relationship between the OTU detection frequency and relative abundance. Meanwhile, the geographic distance (in km) of samples, i.e., a straight-line distance between the sampling points, was calculated based on the longitude and latitude coordinates of each sampling site using the package geosphere [61]. Pearson’s correlation coefficients were calculated to determine the relationships between Bray-Curtis similarity of fungal community and geographic distance of the samples.

To further quantify the relative abundance of stochastic and deterministic processes that drive the fungal community assembly, beta nearest-taxon index (βNTI) was calculated using phylogenetic distance and OTU abundance [62]. βNTI is the number of standard deviations of the beta mean nearest taxon distance (βMNTD, i.e., the observed abundance-weighted-mean phylogenetic distance between closest relatives in two different communities) from the mean of the null distribution [62]. βNTI values between − 2 and 2 indicates a dominance of stochastic processes, whereas βNTI values smaller than − 2 or larger than 2 indicate that deterministic processes play a more important role in community assembly than stochastic processes [9, 62, 63]. In addition, the Bray–Curtis-based Raup–Crick (RC_Bray) value was calculated to further partition the pairwise comparisons that were assigned to stochastic processes. The fraction of pairwise comparisons with |βNTI| < 2 and RC_Bray < − 0.95 indicates a relative dominant influence of homogenizing dispersal on the community assembly. The fraction of pairwise comparisons with |βNTI| < 2 and RC_Bray > 0.95 suggests a crucial role of dispersal limitation. Finally, the fraction of pairwise comparisons with |βNTI| < 2 and |RC_Bray| < 0.95 indicates that“undominated” assembly, including weak selection, weak dispersal, diversification, and/or drift, may be the controlling processes [63, 64]. Rscript bNTI_Local_Machine.r written by Stegen et al. [62] was used to calculate the βNTI and RC_Bray values.

Fungal community driving factors

To explore the potential factors driving the fungal community composition, spatial and environmental variables were evaluated separately. The envfit function in vegan package was used to generate a preliminary explanation of the effects of spatial variables on the fungal community based on 999 permutations. Permutational multivariate analysis of variance (PERMANOVA, the Adonis command in the vegan package) was used to test the influence of spatial variables on community compositions.

A distance-based redundancy analysis (db-RDA) with Bray-Curtis dissimilarities was performed to explore the variations explained by environmental parameters [65]. The following steps were used to determine which parameters should be selected to build a db-RDA model. First, the collinearity of the explanatory variables was resolved by using the vif.cca command to find the values of variance inflation factors (VIFs) in the full models, including all variables [57], and the db-rda command was executed using the parameters with VIF values below 10. To identify the important variables associated with the community structure, forward selection was then performed using an ordiR2step function from a null model containing no explanatory variables to a full model. The ordiR2step function models forward with a maximum adjusted R² at every step, stopping when the adjusted R² starts to decrease, the adjusted R² of the scope is exceeded, or the selected permutation p-value (0.05) is exceeded [57]. Finally, the variables from step selection with the lowest Akaike’s Information Criterion (AIC) were used to build the final model. The vif.cca command was used again to test the collinearity of variables in the final models to ensure that the VIF values were below 10. The significance of the final model, axes, and terms was tested using the anova command.

Network patterns and keystone taxa of fungal communities

Network analysis reveals the interactions in a fungal community and keystone taxa within co-occurrence networks that are crucial for the composition and assembly of the community [66]. The co-occurrence patterns in fungal communities were assessed by performing network analysis using the psych [67] and igraph packages [68]. To reduce network complexity, only OTUs with a relative abundance of > 0.01% in the whole dataset were used, resulting in 593, 474, 589 and 493 OTUs for all samples, samples in 0–10 cm, 10–20 cm, and 20–30 cm layer samples, respectively. Pairwise correlations based on OTU abundance were performed using corr.test function in the psych package, and the p-value was adjusted using false discovery rate (FDR) correction. Spearman’s correlation coefficient between two OTUs with ρ value above 0.7 and p value below 0.01 was considered statically robust. Co-occurrence networks were built using the igraph package; calculation of network topological properties and visualization of networks were done using the interactive platform Gephi v0.92 [69]. In the networks, a degree represents the number of edges connected to a node. Closeness centrality (CC) is based on the average shortest paths and thus reflects the central importance of a node in disseminating information [70]. Betweenness centrality (BC) reveals the role of a node as a bridge between network components. Keystone taxa were therefore OTUs with the highest degree and CC scores, and the lowest BC scores [66]. For the overall network, OTUs with degree over 6, CC over 0.6, and BC below 0.12 were selected as the keystone taxa. For depth-specific networks, OTUs with degree over 15, CC over 0.4, and BC below 0.18 were selected as the keystone taxa.

In the current study, we generated 1,284,478 CCS reads from 110 mangrove sediment samples, with an average Q30 of 98.9%. Pre-processing resulted in 937,168 clean reads were obtained. The number of sequences in each sample ranged from 2186 to 31,416, with an average number of 8520 reads and an average length of 503 bp (Table S2). Overall, 844,885 sequences (90.2% of all sequences) were clustered into 2389 OTUs based on a 97% similarity threshold, of which 85.4% (2041 OTUs, accounting for 83.4% of clustered sequences) were assigned to fungi; the other sequences represented OTUs belonging to Alveolata (41 OTUs), Metazoa (39 OTUs), Protista (39 OTUs), Viridiplantae (20 OTUs), Stramenopila (12 OTUs), and unknown Eukaryota (197 OTUs) (File S1).

Diversity and composition of fungal communities in mangrove sediments

Although the rarefaction curves related to the detected OTUs indicated that most samples did not reach saturation (Fig. S2a), the fungal OTU accumulation curves for different mangroves were nearly asymptotic (Fig. S2b), indicating that the data were largely representative of the fungal diversity in mangrove sediments.

Comparative analyses of α-diversity revealed that the fungal diversity in sediments from 0–10 cm and 10–20 cm depth layers was greater than that in the 20–30 cm layer in terms of OTU richness (p < 0.05, Fig. S3a). The lowest OTU richness (83 OTUs) was detected in the 20–30 cm sediment layer in ZJK (ZJK6-3), and the highest richness (321 OTUs) was detected in the 10–20 cm sediment layer in DZG (DZG4-2). However, no significant differences in α-diversity were observed among the different mangrove locations (Table S3).

For a more accurate identification of fungal OTUs, a phylogenetic tree was constructed using representative sequences of all fungal OTUs (Fig. 1). The overall topology of the constructed phylogenetic tree was generally congruent with previous phylogenetic studies of fungi [71, 72]. In the phylogenetic tree, OTUs were obviously clustered into four representative clades, namely, Ascomycota, Basidiomycota, Chytridiomycota and Rozellomycota, with additional phyla scattered throughout (Fig. 1). For instance, OTUs of Blastocladiomycota, Mortierellomycota and Mucoromycota clustered together and located between Basidiomycota and Chytridiomycota. Aphelidiomycota, together with Kickxellomycota, were closely allied with a clade of unknown fungi within the Rozellomycota clade (Fig. 1). Further, a number of unidentified OTUs clustered into two distinct clades with few known OTUs. The first of these clades was closely related to the phyla Kickxellomycota, Aphelidiomycota and Olpidiomycota, and possibly represents a new fungal phylum. The second of these clades was clustered with Rozellomycota, suggesting that these unidentified OTUs might belong to Rozellomycota. Besides these two unknown fungal clades, most of the unidentified OTUs were close to Chytridiomycota and Rozellomycota, and have been assigned accordingly, indicating a high diversity of phyla Chytridiomycota and Rozellomycota in the mangrove ecosystem.

Based on the results of BLASTn and phylogenetic analyses (Fig. 1), 267 OTUs (13.1% of the fungal OTUs) remained unclassified at phylum level; additional 1774 OTUs were affiliated to 15 phyla, 43 classes, 88 orders, 187 families and 318 genera (Table S4). Ascomycota, Basidiomycota, Rozellomycota and Chytridiomycota were highly diverse and predominant (i.e., relative abundance of > 10%) in the fungal communities in mangrove sediments, accounting for 95.2% of all fungal sequences therein. Other phyla, namely Aphelidiomycota, Basidiobolomycota, Blastocladiomycota, Calcarisporiellomycota, Glomeromycota, Kickxellomycota, Monoblepharomycota, Mortierellomycota, Mucoromycota, Olpidiomycota, and Zoopagomycota, were recovered in small proportions, accounting for 3.4% of OTUs and 1.4% of all sequences (Fig. 2a, Fig. S4ab, and Table S4). At different sediment depth, Ascomycota, Rozellomycota, Basidiomycota, and Chytridiomycota accounted for the largest proportion of OTUs (Fig. S4c). According to the criteria described in Methods, 503 OTUs accounting for 93.7% of all fungal sequences were classified as DT and the remaining 1538 OTUs (6.3% of sequences) were classified as RT (Fig. 2a). DT and RT OTUs belonged to nine and 14 phyla, respectively.

Ascomycota (48.8% of sequences, 34.9% of OTUs) were represented by the classes Sordariomycetes (26.3% of sequences, 10.9% of OTUs) and Dothideomycetes (7.9% of sequences, 12.0% of OTUs). The two other dominant phyla, Basidiomycota (13.1% of sequences, 13.0% of OTUs) and Chytridiomycota (11.2% of sequences, 9.8% of OTUs) were represented by the classes Agaricomycetes (8.7% of sequences, 6.4% of OTUs) and Rhizophydiomycetes (1.8% of sequences, 2.5% of OTUs), respectively. However, OTUs from the dominant phylum Rozellomycota (22.1% of sequences, 25.7% of OTUs) were completely uncategorized at class level (Fig. 2b, Table S4).

Comparison of the fungal OTUs in deep (10–30 cm) and surface (0–10 cm) sediments revealed a significant enrichment of 164 OTUs and a significant depletion of 66 OTUs in the deeper layers. Most of the enriched and depleted OTUs belonged to Ascomycota (25 enriched and 56 depleted), Rozellomycota (11 enriched and 42 depleted), Basidiomycota (16 enriched and 21 depleted) and Chytridiomycota (seven enriched and 24 depleted), and 21 OTUs were unidentified at phylum level (five enriched and 16 depleted, Fig. 3). Depleted OTU with the highest relative abundance in the deep layers was OTU8, an unknown species of Psathyrella (Agaricomycetes, Basidiomycota). Enriched OTU with the highest relative abundance in the deep sediments was OTU6, an unidentified species of phylum Rozellomycota.

Based on IndVal and p-values, 327 indicator OTUs were identified by indicator species analysis, including 161 OTUs in the 0–10 cm layer; 68 OTUs in the 10–20 cm layer; and 98 OTUs in the 20–30 cm layer. Most of these belonged to Ascomycota (120 OTUs), Rozellomycota (58 OTUs), Basidiomycota (46 OTUs) and Chytridiomycota (33 OTUs). The OTUs with the highest IndVal values, which represent the indicator possibility, were OTU3 (unidentified fungus), in the 0–10 cm layer; OTU6 (unidentified Rozellomycota), in the 10–20 cm layer; and OTU63 (Malassezia sp.) in the 20–30 cm layer (Table S5).

Fungal community assembly patterns in mangrove ecosystems

The NCM did not fit well fungal community assembly in all and individual mangroves (R² < 0.6) and the estimated immigration rate (m = 0.01 in all mangroves, m = 0.02–0.05 in individual mangroves) was very low (Fig. 4a, Fig. S5), suggesting the occurrence of only a few dispersal events and low ecological drift in mangroves. Since the NCM explained the low portion of community assembly in mangroves, βNTI was used to explore the relative roles of stochastic and deterministic processes in shaping the microbial community assembly (Fig. 4b). The average βNTI values in all (seven) mangroves and the majority of βNTI (93.8%) in all samples were between − 2 and 2, suggesting that stochastic processes are more important for community assembly than deterministic processes. Further, the mean values (0.958) and the majority (92.1%) of RC_bray in all samples were greater than 0.95, indicating a crucial role of dispersal limitation on the community assembly (Fig. S6). βNTI values in all samples were significantly (p ≤ 0.01) but weakly (R < 0.1) correlated with changes in MAP and TP (Fig. S7). Further, Bray-Curtis similarity of fungal communities (including all, DT and RT) of all samples and sample depths was significantly and negatively correlated with the geographical distance, indicating a significantly strong distance-decay relationship (p < 0.01, Fig. 4c, Fig. S8). Furthermore, Bray-Curtis similarity of DT was significantly higher than that of RT in all samples and sample depths (Fig. S9), suggesting a strong environmental preference or limitation of ecological dispersal.

Spatial and environmental selection shaping the fungal community composition

NMDS and envfit analyses indicated that the composition of fungal community in mangrove sediments is significantly related to the geographical location (R² = 0.7902, p ≤ 0.001), longitude (R² = 0.6434, p ≤ 0.001) and latitude (R² = 0.7325, p ≤ 0.001), which was further supported by the PERMANOVA analysis (Fig. 5a, Fig. S10a). In addition, geographical location, latitude, and longitude significantly affected the fungal communities in different sediment depths (Fig. S10a and S11). Complementary NMDS, envfit and PERMANOVA analyses confirmed a pronounced strong effect of sediment depth on the community composition in different mangroves (p < 0.05, Fig. S10b and S12). For both DT and RT, the location (p ≤ 0.01), longitude (p ≤ 0.01) and latitude (p ≤ 0.01) strongly affected the fungal communities in all samples and sample depths (Fig. 5a, Fig. S10a and S11). Further PERMANOVA analysis also confirmed a significant (p ≤ 0.01) effect of depth on the community compositions of DT and RT in all mangroves (Fig. S10b).

For all taxa, pH, N/NO₃⁻ and N/NH₄⁺ were positively correlated with OTU richness, and pH was positively correlated with the Shannon and Evenness indices. For DT, TC, TOC, and TP were positively correlated with OTU richness, and pH showed positive correlations with Shannon and Evenness indices. For RT, gravel proportion, TN, N/NO₃⁻, and N/NH₄⁺ positively co-varied with OTU richness and Shannon index, while TC, TOC, and TP were only positively correlated with OTU richness (Table 1, Table S6).

Table 1

Correlation coefficients (Spearman’s rho) between OTU richness, alpha diversity, and the relative abundance of main phyla and classes with environmental factors^a
Parameter/phylum	Correlation coefficient (ρ)
Parameter/phylum	Depth	MAT	MAP	Gravel	Salinity	pH	TC	TOC	TN	N/NH₄⁺	N/NO₃⁻	TP	TS
OTU richness	-0.1474	0.0955	0.0224	0.1247	0.0354	0.2536**	0.0102	0.0217	0.1048	0.2494**	0.1984*	0.0782	0.0704
Alpha diversity (Shannon index)	-0.0481	-0.0681	-0.077	-0.0464	0.0956	0.246**	-0.1657	-0.1443	-0.0628	0.0002	0.0289	-0.048	-0.0617
Aphelidiomycota	0.0587	0.0823	-0.0921	0.1393	0.0415	0.0902	-0.1406	-0.1098	-0.1455	-0.0785	0.0409	0.2304	0.0056
Ascomycota	0.109	0.257**	0.0498	0.0813	-0.1011	-0.2081	0.3649***	0.386***	0.2456**	0.0279	-0.0515	-0.2768	0.1953*
Basidiomycota	0.1116	0.1376	-0.0348	0.0116	0.1049	-0.0295	-0.0341	0.0187	-0.0156	-0.0155	-0.0709	-0.1001	0.0272
Chytridiomycota	-0.1139	-0.7471	-0.2323	-0.4097	0.4211***	0.0733	-0.4722	-0.5564	-0.472	-0.4457	0.002	0.0736	-0.5592
Rozellomycota	-0.1556	0.0958	0.1273	0.1532	-0.2593	0.132	0.0471	0.015	0.1029	0.3292***	0.1243	0.2501**	0.1573*
Agaricomycetes	-0.0985	0.1539	0.1177	-0.0154	-0.1226	-0.1812	0.0914	0.1304	0.1317	0.0529	-0.1229	-0.1407	0.069
Chytridiomycetes	-0.1264	-0.2188	0.0449	-0.1635	0.0365	0.0232	-0.0944	-0.0891	-0.0512	-0.0934	-0.0836	0.0268	-0.0851
Cystobasidiomycetes	0.2844	-0.3682	-0.1744	-0.1824	0.299***	0.0774	-0.244	-0.3068	-0.2091	-0.3549	0.0874	-0.1428	-0.1984
Dothideomycetes	0.1558	-0.0144	0.3157***	0.1081	-0.1377	0.1505	0.1698*	0.1842*	0.0562	0.1138	-0.1965	0.1752*	0.2264**
Eurotiomycetes	-0.145	-0.0006	0.0138	-0.0899	0.0778	-0.003	-0.0186	0.0175	-0.0124	0.0208	-0.0572	-0.0469	-0.0847
Leotiomycetes	0.1086	-0.0032	0.155	0.0603	-0.1493	0.0984	-0.0284	-0.0125	0.0265	0.087	0.0362	0.023	0.053
Lobulomycetes	-0.0965	-0.6081	-0.1945	-0.3081	0.3197	0.0379	-0.3157	-0.3979	-0.3276	-0.3497	0.0343	0.0618	-0.5058
Rhizophydiomycetes	-0.1424	-0.0232	-0.0836	-0.0274	0.1622*	0.0228	-0.1819	-0.1543	-0.1455	-0.1107	-0.0495	-0.0888	-0.1171
Saccharomycetes	-0.0627	0.1004	-0.1228	-0.0007	0.0222	0.083	0.0589	0.0383	0.0129	0.1592*	-0.091	0.0588	0.0343
Sordariomycetes	0.026	0.3111***	-0.0224	0.0499	-0.0768	-0.3189	0.2612**	0.2893**	0.2374**	-0.0251	0.0045	-0.2943	0.1615*
Tremellomycetes	0.119	0.0649	-0.1347	0.0518	0.1837*	0.1058	-0.0854	-0.0468	-0.1493	0.119	-0.001	0.1389	-0.0214
^a MAT, mean annual temperature; MAP, mean annual precipitation; Gravel, Gravel proportion; TC, total carbon, TOC, total organic carbon; TN, total nitrogen; TP, total phosphorus ; TS, total sulfur. , p < 0.05; , p < 0.01; **, p < 0.001.

Correlation analyses were further performed to explore the relationship between the environmental variables and the relative abundance of main phyla and classes (Table 1). The relative abundance of Ascomycota was significantly positively correlated with MAT, TC, TOC, TN, and TS. Similarly, the relative abundance of Sordariomycetes was positively correlated with MAT, TC, TOC, TN, and TS. On the other hand, Eurotiomycetes, another dominant class of Ascomycota, was significantly correlated with MAP, TC, TOC, TP, and TS. Furthermore, a significantly positive correlation was observed between the relative abundance of Chytridiomycota and salinity. The relative abundances of the classes Cystobasidiomycetes, Rhizophydiomycetes and Tremellomycetes were also positively correlated with salinity. Finally, the relative abundance of Rozellomycota was positively correlated with N/NH₄⁺, TP, and TS.

The environmental variables selected to build the db-RDA models and their effects on the community compositions of different taxa are visualized in Fig. 5b. All designed models were statistically validated based on ANOVA. For all taxa, a model encompassing 11 significant variables (the depth, gravel proportion, MAT, MAP, N/NH4⁺, tree number, salinity, TOC, N/NH₄⁺, N/NO₃⁻, and TP) was built; it explained 27.5% of the variations in fungal communities (Fig. 5b). The explanations of each model build for the seven mangroves ranged from 24.4% (XMD) to 41.0% (SK). Variables in the models for DT communities in all samples were the same as model for all taxa, but exhibited a slightly higher explanation of variation (28.8%). The models for DT communities in different mangroves explained 26.0% (SK) to 48.2% (ZJ) of community variation. A 10-variable model was constructed for RT in all samples, but only 19.2% of the variations. The models for RT communities in different mangroves explained 10.6% (SK) to 23.0% (FT) of the community variations (Fig. 5b).

Co-occurrence patterns of fungal communities in mangrove sediment

Based on the Spearman correlation analysis, a co-occurrence network consisting of 124 nodes (OTUs) and 176 edges was generated for all samples, as depicted in Fig. 6. This overall network with low complexity and high modularity had a diameter of 5, average degree (AD) of 1.419, modularity of 0.884, and average path length (APL) of 1.89. The network comprised 24 modules, of which the top eight modules accounted for only 65.32% of the nodes (Fig. 6). The nodes in the network were assigned to seven phyla, accounting for 75.8% of all nodes, with 24.2% remaining unidentified at phylum level. Most nodes were annotated as DT. Network analysis revealed a 100% positive edge, indicating that the co-occurring relationships accounted for the entire fungal network. In the overall network, 16 OTUs were recognized as keystones; most of them belonged to modules 1–4 (Table S7). Six keystone taxa belonged to Ascomycota, four belonged to Basidiomycota, five belonged to Rozellomycota and Chytridiomycota, and one was unidentified fungus. The calculated topological properties are summarized in Table S8.

To better illustrate the fungal co-occurrence patterns in each sediment depth, fungal networks were then constructed separately for the three sampling depths (Fig. S13), and their topological properties were calculated correspondingly (Table S8). The network generated for the 0–10 cm layer encompassed 276 nodes and 610 edges; that for the 10–20 cm layer encompassed 331 nodes and 1136 edges; and that for the 20–30 cm layer encompassed 315 nodes and 1047 edges. The 0–10 cm layer network contained 28 modules, with the top eight modules accounting for 81.9% of all nodes. Five OTUs in the network were recognized as keystone taxa, with two OTUs belonging to RT (Table S7). The 10–20 cm layer network encompassed 31 modules, and the top eight modules accounted for 77.7% of nodes. Among these nodes, 11 OTUs were identified as keystone taxa, with one belonging to RT (Table S7). The 20–30 cm layer network comprised of 20 modules, with the top eight in the proportion of 90.2% of all nodes. Further, 32 OTUs in this network were recognized as keystone taxa, five of which were RT (Table S7). The numbers of edges, nodes, and average degrees in the deep-sediment fungal network were higher than these in the surface-sediment network. Further, the APL was lower in the deep-sediment fungal network, indicating a possibly higher frequency of co-occurrence of fungi in deep sediments, and their relatively close correlation. By contrast, fungi in the surface sediment were less connected and their correlations were relatively loose.

Although mangroves are generally considered as a hotspot of biodiversity and harbor a high diversity of specific fungi, studies, especially these based on culture-independent methods, have been poorly performed. The present study constitutes the first attempt to systematically investigate the fungal communities in mangrove sediments in China using PacBio SMRT sequencing. The reported discovery of a great number of early diverging fungal lineages, assembly patterns, and co-occurrence relationships within fungal communities will spur further studies into the utilization and protection of fungal resources and communities in mangrove sediments.

Mangrove ecosystem encompasses a high number of basal fungal linages

Fungi represent one of the most diverse groups of life on Earth, with an estimated 2.2–5.1 million species [28, 71, 73]. According to the latest taxonomy, the fungi described to date are affiliated with 19 fungal phyla [74]. Remarkably, 15 fungal phyla were identified in mangrove sediment samples in the current study, revealing the presence of highly divergent fungal lineages in mangrove sediments.

In the current study, Rozellomycota and Chytridiomycota, two phyla that have been rarely reported in mangroves to date, were predominant and incredibly highly diverse. These two zooporic phyla are distributed worldwide and represent the basal fungal clades [71, 75]. Chytridiomycota is highly diverse and widely distributed in the extreme environments characterized by high salinity, such as ocean sediments [76, 77], hydrothermal vents [78], and methane seeps [79]. Rozellomycota are also remarkably phylogenetically diverse [71] and found in both, aquatic ecosystems and terrestrial habitats [75, 80]. Since these two phyla have been established based mainly on the molecular data, their members are mostly novel and uncultured [71, 81]. Consistently, OTUs belonging to Rozellomycota and Chytridiomycota identified in the current study remain almost completely unidentified at genus and species levels. In addition, several other basal fungal phyla with low OTU proportion and relative abundance were recovered in the current study, most of which had never been reported in mangroves.

There are several possible explanations for the high proportion of basal fungal lineages uncovered in the current study. First, mangrove sediments harbor highly divergent and diverse fungal lineages because of a relatively high speciation but low extinction rate in tropical habitats [71]. However, while mangroves are transitional ecosystems between the land and sea that have unique characteristics, which may support distinct fungal species, the mangrove fungi are relatively poorly covered by biodiversity and taxonomic studies. Second, similar to the phyla Ascomycota and Basidiomycota, Rozellomycota and Chytridiomycota are comprised by a small proportion of DT and a large proportion of RT, which are difficult to detect because of methodological bias and artifacts [26, 82, 83]. Fortunately, nested PCR is highly sensitive and can detect low amount of DNA [84, 85]. Third, PacBio SMRT can yield homogeneous results even for GC-rich sequences, and more accurate annotations are supplied in the new version of fungal and eukaryotic UNITE databases.

In the current study, a number of fungal OTUs were found to be specific to surface or deep sediments; most of these were depleted in deep sediment layers compared with surface sediments (Fig. 3). That might be because of the different niches or preference for specific environmental conditions of the identified fungal taxa [27, 86]. Most of such specific fungal taxa were Ascomycota and Rozellomycota, with some Basidiomycota and Chytridiomycota, which were also predominant in mangrove ecosystems, indicating that Ascomycota and Basidiomycota might be more sensitive to the vertical variations of environmental characteristics in mangrove sediments than other fungi. Psathyrella sp. (OTU8) was a depleted fungal taxon with the highest relative abundance; the genus contains reported as a saprotrophic and mycorrhizal fungi from the forest ecosystem [87]. Similarly, Ganoderma gibbosum, a species of Agaricomycetes, was the best indicator of the 0–10 cm sediment layer, together with two unknown fungi with a higher IndVal value than G. gibbosum. Interestingly, Paraconiothyrium cyclothyrioides, previously isolated from a contaminated mangrove in Brazil [88] and reportedly a cause of human disease [89], was the best indicator of the 10–20 cm sediment layer. Species of Paraconiothyrium are commonly found in the soil or woody plants, as endophytes or agents of plant disease, or in the clinic as human pathomycetes [89, 90]. In the 20–30 cm sediment layer, Malassezia sp. was the indicator species with highest IndVal value. Members of Malassezia are likely the most widespread fungi on Earth and have been identified in a startling diversity of habitats and locations, from the human skin, through the polar regions, to deep-sea sediments [91].

Collectively, the presented data indicate an incredibly high divergence and diversity of fungi in mangrove sediments, especially early diverging and unknown fungal lineages.

Fungal community assembly is mainly controlled by stochastic processes

According to Zhang et al. [9], deterministic processes in the prokaryotic community assembly in mangroves across southeast China to a greater extent than stochastic processes. By contrast, the data obtained in the current study on the fungal community assembly in mangrove sediments indicate a relatively greater importance of stochastic processes. NCM accounted for a moderate portion of community variation (R² = 0.56), with an extremely low estimated immigration rate (m), suggesting that although stochastic processes are relatively important for the fungal community assembly, the dispersal and ecological drift are very low in mangrove sediments. The community variation (R²) in individual mangroves determined by NCM was much lower than that in all samples considered together, as also supported by the ensuing db-RDA analyses, which indicated that the environmental conditions generally influence the fungal community composition in an individual mangrove to a greater extent than in overall mangroves (Fig. 5b). These observations indicate that the effect of stochastic processes on the community assembly in an individual mangrove is weaker than that of overall mangroves, and might be accompanied by other community assembly mechanisms, including environmental selection and species interactions, which are difficult to evaluate [53, 92].

To explore the relative effects of stochastic and deterministic processes on fungal community assembly, βNTI was then calculated based on the OTU abundance and their phylogenetic distance. Unlike NCM fitting, βNTI analysis supported a crucial role of stochastic processes in the community assembly (Fig. 4b). Based on the distribution of βNTI scores, 93.8% of all comparisons were consistent with a random phylogenetic turnover; the fungal community assembly was therefore mainly controlled by stochastic processes in the investigated mangroves [62]. In addition, the RC_bray values suggested that dispersal limitation is more crucial for the community assembly than homogenizing dispersal and “undominated” assembly in mangrove sediments (Fig. S6) [63]. Consistently, only changes in MAP and TP exerted a significant (p < 0.01) but weak (-0.1 < R < 0.1) effect on βNTI values (Fig. S7), indicating only a weak effect of environmental variables on the variation of fungal community composition [54]. The significantly strong distance-decay correlation between the similarity of fungal communities and geographic distance further confirmed the importance of stochastic processes (Fig. 4c, Fig. S8). That is because community similarity is predicted to decrease along the geographic distance as a result of the dispersal limitation and ecological drift [86, 93]. In summary, the above observations suggest that stochastic processes, mainly dispersal limitation, strongly shape the fungal community composition and that deterministic processes play a minor role in community distribution in mangrove sediment.

Spatial and environmental selections on fungal communities in mangrove sediments

Spatial and environmental variables are critical elements determining the microbial diversity and community composition in various environments [9, 94–96]. However, the influence of spatial and environmental factors on the fungal community in mangrove sediments has rarely been studied, except for the presence of plants and sediment depth [27]. Further, the knowledge of many ecosystems is mostly based on the dominant and entire communities, while the role of rare taxa is unaccounted for [53]. Based on previous studies, some rare taxa are in fact metabolically active in the environment and may constitute keystone species that regulate the functions of aquatic ecosystem; the “rare biosphere” is hence of great importance to the metabolic and ecological functions of aquatic habitats [53, 83, 97]. Consequently, in the current study, the fungal community was separated into DT and RT, to explore the biogeography and potential controlling factors of the “rare biosphere” in mangrove sediments.

In the current study, the fungal all taxa, DT and RT communities significantly grouped together according to the mangrove location (Fig. 5a), suggesting a higher similarity of the communities within the same mangrove than between mangroves. This was supported by the distance-decay relationship of community similarity and geographical distance (Fig. 4c, Fig. S8). These observations suggest a similar biogeography of DT and RT fungi in mangrove sediments, which is consistent with previous studies of bacterial communities in coastal Antarctic lakes [98] and microeukaryotic communities in Tingjiang River [53]. In agreement with these observations on the effect of geographic location on fungal community, Zhang et al. [9] reported a similar distribution of archaea and bacteria in mangrove sediments in China, indicating that microbial communities in mangrove sediments are strongly shaped by the geographical location, possibly because of the combined effects of climate, niche conservatism, and rates of dispersal, evolutionary radiation and extinction [94, 99]. Depth is a known and important factor that strongly affects microbial communities in mangrove sediments [8, 16, 27]. Similar to Luis et al. [27], a significant but weak effect of the depth on fungal communities was observed herein for all taxa, DT, and RT fungi, and it impacted DT more so than RT (Fig. 5b, Fig. S10).

Although several significant environmental variables were identified and used to build the db-RDA model of fungal community composition, the explained variations indicated that the environmental and spatial factors play only a minor role, as relatively low variations were revealed in all samples (Fig. 5b). Specifically for RT, only 19.2% of community variations were explained by five parameters. Consistently, db-RDA or variance partitioning analysis (VPA) also revealed a large proportion of unexplained microeukaryotic or prokaryotic community variations in different habitats in several previous studies [53, 100–102]. There are several potential explanations of this phenomenon. First, additional important influencing factors exist that have not been included in the current study [53]. Second, while co-occurrence relationships among microbes significantly affect the community compositions, they cannot be quantified by db-RDA [22, 103]. Further, as suggested by previous studies, taxa that are rare under one condition may become prevalent once the conditions become suitable. In other words, a species might be recognized as rare may because of a dearth of suitable microenvironments or microniches [26, 53, 83]. In the current study, the ensuing db-RDA analyses of DT and RT communities yielded different models and explanations as both these communities were significantly impacted by different environmental parameters. PERMANOVA confirmed the different influences of spatial variables on DT and RT communities (Fig. S10). That could be because different taxonomic and functional groups of microbial communities may occupy different ecological niches and be shaped by contrasting underlying factors [86, 104]. In summary, the above findings indicate a crucial role of geographic location and a minor role of environmental selections in driving the fungal community in mangrove sediments.

Co-occurrence network patterns and keystone taxa in fungal communities in mangrove sediments

Instead of thriving in isolation, microorganisms form complex interaction networks [70]. Interactions between microorganisms such as antagonism and cooperation, significantly affect the microbial community composition and dynamics [22, 66, 103]. According to the analyses performed in the current study, the fungal network in deep sediment is made up of highly connected taxa (nodes) and is more complex than that in the surface sediment. That may be because of frequent interactions between microbial communities in the surface sediments and other ecosystems, promoted by biotic factors, i.e., animal activities, or abiotic factors, i.e., daily rising tide. Such frequent interactions and the appearance of microbes belonging to other ecosystems may result in a co-occurrence network of low complexity [105, 106]. Further, the above biotic and abiotic factors significantly affect the microenvironmental conditions in the surface sediment rather than those in the deep sediment, which greatly impacts the structure and reduces complexity of fungal ecological network in the surface sediment [107]. In the current study, in different depth networks, 34–44% of the nodes were RT, supporting the notion that RT may play important roles in the microbial interactions and ecological functions of aquatic ecosystems [53, 83, 97].

Keystone taxa are highly connected taxa that are particularly important for maintaining the connectivity of a community network [53, 66, 108]. In the current study, OTU365 and OTU243, two keystone taxa in the overall network, were identified as yeast-like Hannaella sp. and Tausonia pullulans, respectively. These two fungi reportedly greatly contribute to carbon metabolism and cycling in the mangrove sediments [27, 109]. Remarkably, an unknown species of Zygophlyctis (OTU242), a genus reported as an algal parasite, was identified as a keystone taxon in deep sediments [110]. Node OTU242 in the network had a degree of 14, with all correlations being significantly positive, indicating that members of Zygophlyctis may have a fabulously number of parasitic or symbiotic objects and play important roles in the aquatic ecosystem. Furthermore, the majority of identified keystone taxa in different networks were unclassified at genus levels. Consequently, the fungi that play key roles in mangrove fungal communities await further explorations and identification.

In conclusion, in the current study, we used PacBio-based metabarcoding approach for a comprehensive assessment of the compositions, controlling factors, assembly, and network patterns of the fungal communities in mangrove sediments in China. The investigation revealed a high level of fungal diversity, including a number of early diverging fungal lineages such as Rozellomycota and Chytridiomycota. Based on NCM modeling, βNTI analysis, and the strong distance-decay relationship of community similarity, stochastic processes, mainly dispersal limitation, are proposed to be the determinants of fungal community assembly. Geographical locations and depth were the key determinants of fungal community structure in mangrove sediments. Further, the db-RDA analysis indicated only a minor role of environmental selection on the fungal community. The fungal network was more complex in the deep sediment than in the surface sediment, and the identified keystone taxa may play important roles in nutrient cycling and maintaining network connectivity. Overall, the findings of the current study inform an overview of fungal community in mangrove sediment, and promote the cultivation of uncultivable fungi, and the utilization and protection of fungal resources in mangroves.

HTS: high-throughput sequencing; ITS: internal transcribed spacer; PacBio: Pacific Biosciences; CCS: circular consensus sequencing; SMRT: single-molecule real-time; MAT: mean annual temperature; MAP: mean annual precipitation; TC: total carbon; TOC: total organic carbon; TN: total nitrogen; N/NH4⁺: ammonium nitrogen; N/NO3⁻: nitrate nitrogen; TP: total phosphorus; TS: total sulfur; PCR: polymerase chain reaction; HMMs: hidden Markov models; OTU: operational taxonomic units; DT: dominant taxa; ART: always rare taxa; CRT; conditionally rare taxa; RT: rare taxa; ANOVA: analysis of variance; HSD: Tukey’s honest significant difference; NMDS;: nonmetric multidimensional scaling; NCM: neutral community model; βNTI: beta nearest-taxon index; βMNTD: beta mean nearest taxon distance; RC_Bray: Bray–Curtis-based Raup–Crick; PERMANOVA: permutational multivariate analysis of variance; db-RDA: distance-based redundancy analysis; VIF: variance inflation factor; AIC: Akaike’s Information Criterion; FDR: false discovery rate; CC: closeness centrality; BC: betweenness centrality; APL: average path length; VPA: variance partitioning analysis.

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and materials

All raw sequences from this study have been submitted to the NCBI Sequence Read Archive (SRA) database under the BioProject number PRJNA 668858. The rawdata has also been deposited in NODE (http://www.biosino.org/node) with accession number OEP 001211 as backup.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was financially supported by the National Science and Technology Fundamental Resources Investigation Program of China (Grant No. 2019FY100700), National Natural Science Foundation of China (Grant No. 91851105, 31970105) and China Postdoctoral Science Foundation (no. 2020M672779).

Authors' contributions

ZFZ and ML conceived the study. ZFZ and YPP collected sediment samples in filed. ZFZ performed the laboratory experiments with the help of YPP and YL. ZFZ analyzed data and wrote the paper. ML revised the manuscript. All authors read and approved the final manuscript.

Acknowledgements

We thank Dr. Cui-Xing Ye in Hainan Academy of Ocean and Fisheries Sciences, Dr. Wen-Lu Lan in Guangxi Marine Environment Monitoring Centre station, Dr. Chun-Fang Zheng in Wenzhou University and Dr. Jian-Biao Chou in Zhejiang Mariculture Research Institute for help with sample collection. We acknowledge Dr. Joanna Mackie for providing professional writing services.

Giri C, Ochieng E, Tieszen LL, Zhu Z, Singh A, Loveland T et al. Status and distribution of mangrove forests of the world using earth observation satellite data. Global Ecol Biogeogr. 2011; 20:154-159.
Alongi DM. Carbon cycling and storage in mangrove forests. Ann Rev Mar Sci. 2014; 6:195-219.
Donato DC, Kauffman JB, Murdiyarso D, Kurnianto S, Stidham M, Kanninen M. Mangroves among the most carbon-rich forests in the tropics. Nat Geosci. 2011; 4:293-297.
Friesen SD, Dunn C, Freeman C. Decomposition as a regulator of carbon accretion in mangroves: a review. Ecol Eng. 2018; 114:173-178.
Liu Y, Zhou Z, Pan J, Baker BJ, Gu JD, Li M. Comparative genomic inference suggests mixotrophic lifestyle for Thorarchaeota. ISME J. 2018; 12:1021-1031.
Liu X, Li M, Castelle CJ, Probst AJ, Zhou Z, Pan J et al. Insights into the ecology, evolution, and metabolism of the widespread Woesearchaeotal lineages. Microbiome. 2018; 6:102.
Cai M, Liu Y, Yin X, Zhou Z, Friedrich MW, Richter-Heitmann T et al. Diverse Asgard archaea including the novel phylum Gerdarchaeota participate in organic matter degradation. Sci China Life Sci. 2020; 63:886-897.
Pan J, Chen Y, Wang Y, Zhou Z, Li M. Vertical Distribution of Bathyarchaeotal Communities in Mangrove Wetlands Suggests Distinct Niche Preference of Bathyarchaeota Subgroup 6. Microb Ecol. 2019; 77:417-428.
Zhang CJ, Pan J, Duan CH, Wang YM, Liu Y, Sun J et al. Prokaryotic Diversity in Mangrove Sediments across Southeastern China Fundamentally Differs from That in Other Biomes. mSystems. 2019; 4.
Holguin G, Vazquez P, Bashan Y. The role of sediment microorganisms in the productivity, conservation, and rehabilitation of mangrove ecosystems: an overview. Biol Fert Soils. 2001; 33:265-278.
Schmit JP, Shearer CA. A checklist of mangrove-associated fungi, their geographical distribution and known host plant. Mycotaxon. 2003; 85:423-477.
Sahoo K, Dhal NK. Potential microbial diversity in mangrove ecosystems：A review. Indian Journal of Marine Sciences. 2009; 38:249-256.
Ravindran C, Naveenan T, Varatharajan GR. Optimization of alkaline cellulase production by the marine-derived fungus Chaetomium sp. using agricultural and industrial wastes as substrates. Bot Mar. 2010; 53:275-282.
Schlingmann G, Milne L, Carter GT. Isolation and identification of antifungal polyesters from the marine fungus Hypoxylon oceanicum LL-15G256. Tetrahedron. 2002; 58:6828-6835.
Behera BC, Singdevsachan SK, Mishra RR, Dutta SK, Thatoi HN. Diversity, mechanism and biotechnology of phosphate solubilising microorganism in mangrove—A review. Biocatalysis and Agricultural Biotechnology. 2014; 3:97-110.
Sarma VV, Hyde KD. A review on frequently occurring fungi in mangroves. Fungal Divers. 2001; 8:1-34.
Kumaresan V, Suryanarayanan TS. Endophytes assemblages in young mature and senescent leaves of Rhizophora apiculata: evidence for the role of endophytes in mangrove litter degradation. Fungal Divers. 2002; 9:81-91.
Costa IP, Maia LC, Cavalcanti MA. Diversity of leaf endophytic fungi in mangrove plants of northeast brazil. Braz J Microbiol. 2012; 43:1165-1173.
Cortes-Lorenzo C, Rodriguez-Diaz M, Lopez-Lopez C, Sanchez-Peinado M, Rodelas B, Gonzalez-Lopez J. Effect of salinity on enzymatic activities in a submerged fixed bed biofilm reactor for municipal sewage treatment. Bioresour Technol. 2012; 121:312-319.
Colares GB, Melo VMM. Relating microbial community structure and environmental variables in mangrove sediments inside Rhizophora mangle L. habitats. Appl Soil Ecol. 2013; 64:171-177.
Lin X, Hetharua B, Lin L, Xu H, Zheng T, He Z et al. Mangrove Sediment Microbiome: Adaptive Microbial Assemblages and Their Routed Biogeochemical Processes in Yunxiao Mangrove National Nature Reserve, China. Microb Ecol. 2019; 78:57-69.
Lima-Mendez G, Faust K, Henry N, Decelle J, Colin S, Carcillo F et al. Determinants of community structure in the global plankton interactome. Science. 2015; 348:1262073.
Cabral L, Junior GVL, Pereira de Sousa ST, Dias ACF, Lira Cadete L, Andreote FD et al. Anthropogenic impact on mangrove sediments triggers differential responses in the heavy metals and antibiotic resistomes of microbial communities. Environ Pollut. 2016; 216:460-469.
Arfi Y, Buee M, Marchand C, Levasseur A, Record E. Multiple markers pyrosequencing reveals highly diverse and host-specific fungal communities on the mangrove trees Avicennia marina and Rhizophora stylosa. FEMS Microbiol Ecol. 2012; 79:433-444.
Arfi Y, Marchand C, Wartel M, Record E. Fungal diversity in anoxic-sulfidic sediments in a mangrove soil. Fungal Ecol. 2012; 5:282-285.
Shade A, Jones SE, Caporaso JG, Handelsman J, Knight R, Fierer N et al. Conditionally rare taxa disproportionately contribute to temporal changes in microbial diversity. mBio. 2014; 5:e01371-01314.
Luis P, Saint-Genis G, Vallon L, Bourgeois C, Bruto M, Marchand C et al. Contrasted ecological niches shape fungal and prokaryotic community structure in mangroves sediments. Environ Microbiol. 2019; 21:1407-1424.
Hawksworth DL, Lucking R. Fungal Diversity Revisited: 2.2 to 3.8 Million Species. Microbiol Spectr. 2017; 5.
Singh P, Raghukumar C, Meena RM, Verma P, Shouche Y. Fungal diversity in deep-sea sediments revealed by culture-dependent and culture-independent approaches. Fungal Ecol. 2012; 5:543-553.
Nilsson RH, Anslan S, Bahram M, Wurzbacher C, Baldrian P, Tedersoo L. Mycobiome diversity: high-throughput sequencing and identification of fungi. Nat Rev Microbiol. 2019; 17:95-109.
Heeger F, Bourne EC, Baschien C, Yurkov A, Bunk B, Sproer C et al. Long-read DNA metabarcoding of ribosomal RNA in the analysis of fungi from aquatic environments. Mol Ecol Resour. 2018; 18:1500-1514.
Tedersoo L, Tooming-Klunderud A, Anslan S. PacBio metabarcoding of Fungi and other eukaryotes: errors, biases and perspectives. New Phytol. 2018; 217:1370-1385.
Goodwin S, McPherson JD, McCombie WR. Coming of age: ten years of next-generation sequencing technologies. Nat Rev Genet. 2016; 17:333-351.
Travers KJ, Chin CS, Rank DR, Eid JS, Turner SW. A flexible and efficient template format for circular consensus sequencing and SNP detection. Nucleic Acids Res. 2010; 38:e159.
Zhang XX, Zhang YJ, He XF, Chen CP, Chen NC. Improvement of Extraction and Determination Method of Total Sulfur in Soil. Journal of Anhui Agricultural Sciences. 2017; 45:91-93.
Tedersoo L, Lindahl B. Fungal identification biases in microbiome projects. Environ Microbiol Rep. 2016; 8:774-779.
White TJ, Bruns TD, B. LS, Taylor JW. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In: Innis MA, Gelfand DH, Sninsky JJWT, editors. PCR protocols: a guide to methods and applications New York: Academic Press; 1990. p. 315–322.
Gardes M, Bruns TD. ITS primers with enhanced specificity for basidiomycetes - application to the identification of mycorrhizae and rusts. Mol Ecol. 1993; 2:113-118.
Dodt M, Roehr JT, Ahmed R, Dieterich C. FLEXBAR-Flexible Barcode and Adapter Processing for Next-Generation Sequencing Platforms. Biology (Basel). 2012; 1:895-905.
Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009; 75:7537-7541.
Bengtsson-Palme J, Ryberg M, Hartmann M, Branco S, Wang Z, Godhe A et al. Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for analysis of environmental sequencing data. Methods Ecol Evol. 2013; 4:914-919.
Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics. 2011; 27:2194-2200.
Rognes T, Flouri T, Nichols B, Quince C, Mahe F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016; 4:e2584.
Schoch CL, Robbertse B, Robert V, Vu D, Cardinali G, Irinyi L et al. Finding needles in haystacks: linking scientific names, reference specimens and molecular data for Fungi. Database (Oxford). 2014:1-21.
Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010; 26:2460-2461.
Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013; 10:996-998.
Abarenkov K, Zirk A, Piirmann T, Pöhönen R, Ivanov F, Nilsson RH et al. Full UNITE+INSD dataset for eukaryotes. In.: UNITE Community; 2020: Version 04.02.2020.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009; 10:421.
Nguyen L, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies. Mol Biol Evol. 2015; 32:268-274.
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004; 32:1792-1797.
Capella-Gutierrez S, Silla-Martinez JM, Gabaldon T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009; 25:1972-1973.
Letunic I, Bork P. Interactive Tree Of Life (iTOL) v4: recent updates and new developments. Nucleic Acids Res. 2019; 47:W256-W259.
Chen W, Ren K, Isabwe A, Chen H, Liu M, Yang J. Stochastic processes shape microeukaryotic community assembly in a subtropical river across wet and dry seasons. Microbiome. 2019; 7:138.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010; 26:139-140.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47.
De Cáceres M. How to use the indicspecies package. 2020; https://cran.r-project.org/web/packages/indicspecies/vignettes/indicspeciesTutorial.pdf.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, D. M et al. vegan: Community Ecology Package. R package version 2.5-6. 2019; https://CRAN.R-project.org/package=vegan.
Sloan WT, Lunn M, Woodcock S, Head IM, Nee S, Curtis TP. Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environ Microbiol. 2006; 8:732-740.
Hubbell SP. A unified neutral theory of biodiversity and biogeography. Princeton: Princeton University Press; 2001.
Burns AR, Stephens WZ, Stagaman K, Wong S, Rawls JF, Guillemin K et al. Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J. 2016; 10:655-664.
Hijmans RJ. geosphere: Spherical Trigonometry. R package version 1.5-10. 2019; https://CRAN.R-project.org/package=geosphere.
Stegen JC, Lin X, Konopka AE, Fredrickson JK. Stochastic and deterministic assembly processes in subsurface microbial communities. ISME J. 2012; 6:1653-1664.
Jiao S, Yang Y, Xu Y, Zhang J, Lu Y. Balance between community assembly processes mediates species coexistence in agricultural soil microbiomes across eastern China. ISME J. 2020; 14:202-216.
Tripathi BM, Stegen JC, Kim M, Dong K, Adams JM, Lee YK. Soil pH mediates the balance between stochastic and deterministic assembly of bacteria. ISME J. 2018; 12:1072-1083.
Legendre P, Gallagher ED. Ecologically meaningful transformations for ordination of species data. Oecologia. 2001; 129:271-280.
Berry D, Widder S. Deciphering microbial interactions and detecting keystone species with co-occurrence networks. Front Microbiol. 2014; 5:219.
Revelle W. psych: Procedures for psychological, psychometric, and personality research. Northwestern University, Evanston, Illinois, Software. 2017.
Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal, complex systems. 2006; 1695:1-9.
Bastian M, Heymann S, Jacomy M. Gephi: an open source software for exploring and manipulating networks. International AAAI Conference on Weblogs and Social Media. 2009.
Banerjee S, Walder F, Buchi L, Meyer M, Held AY, Gattinger A et al. Agricultural intensification reduces microbial network complexity and the abundance of keystone taxa in roots. ISME J. 2019; 13:1722-1736.
Tedersoo L, Bahram M, Puusepp R, Nilsson RH, James TY. Novel soil-inhabiting clades fill gaps in the fungal tree of life. Microbiome. 2017; 5:42.
Wijayawardene NN, Pawłowska J, Letcher PM, Kirk PM, Humber RA, Schüßler A et al. Notes for genera: basal clades of Fungi (including Aphelidiomycota, Basidiobolomycota, Blastocladiomycota, Calcarisporiellomycota, Caulochytriomycota, Chytridiomycota, Entomophthoromycota, Glomeromycota, Kickxellomycota, Monoblepharomycota, Mortierellomycota, Mucoromycota, Neocallimastigomycota, Olpidiomycota, Rozellomycota and Zoopagomycota). Fungal Divers. 2018; 92:43-129.
Blackwell M. The fungi: 1, 2, 3 ... 5.1 million species? Am J Bot. 2011; 98:426-438.
Wijayawardene NN, Hyde KD, Al-Ani LKT, Tedersoo L, Haelewaters D, Aptroot A et al. Outline of Fungi and fungus-like taxa. Mycosphere. 2020; 11:1060-1456.
Grossart H, Wurzbacher C, James TY, Kagami M. Discovery of dark matter fungi in aquatic ecosystems demands a reappraisal of the phylogeny and ecology of zoosporic fungi. Fungal Ecol. 2016; 19:28-38.
Wang Y, Guo X, Zheng P, Zou S, Li G, Gong J. Distinct seasonality of chytrid-dominated benthic fungal communities in the neritic oceans (Bohai Sea and North Yellow Sea). Fungal Ecol. 2017; 30:55-66.
Chrismas N, Cunliffe M. Depth-dependent mycoplankton glycoside hydrolase gene activity in the open ocean-evidence from the Tara Oceans eukaryote metatranscriptomes. ISME J. 2020; 14:2361-2365.
Calvez TL, Burgaud G, Mahe S, Barbier G, Vandenkoornhuyse P. Fungal Diversity in Deep-Sea Hydrothermal Ecosystems. Appl Environ Microb. 2009; 75:6415-6421.
Nagahama T, Takahashi E, Nagano Y, Abdel-Wahab MA, Miyazaki M. Molecular evidence that deep-branching fungi are major fungal components in deep-sea methane cold-seep sediments. Environ Microbiol. 2011; 13:2359-2370.
Jones MD, Forn I, Gadelha C, Egan MJ, Bass D, Massana R et al. Discovery of novel intermediate forms redefines the fungal tree of life. Nature. 2011; 474:200-203.
Ivarsson M, Bengtson S, Neubeck A. The igneous oceanic crust – Earth's largest fungal habitat? Fungal Ecol. 2016; 20:249-255.
Pedrós-Alió C. Marine microbial diversity: can it be determined? Trends Microbiol. 2006; 14:257-263.
Lynch MD, Neufeld JD. Ecology and exploration of the rare biosphere. Nat Rev Microbiol. 2015; 13:217-229.
Fan ZY, Li XR, Mao DP, Zhu GF, Wang SY, Quan ZX. Could nested PCR be applicable for the study of microbial diversity? World J Microb Biot. 2009; 25:1447-1452.
Yu G, Fadrosh D, Goedert JJ, Ravel J, Goldstein AM. Nested PCR Biases in Interpreting Microbial Community Structure in 16S rRNA Gene Sequence Datasets. PLoS One. 2015; 10:e0132253.
Bahram M, Kohout P, Anslan S, Harend H, Abarenkov K, Tedersoo L. Stochastic distribution of small soil eukaryotes resulting from high dispersal and drift in a local environment. ISME J. 2016; 10:885-896.
Hobbie EA, Weber NS, Trapp JM. Mycorrhizal vs saprotrophic status of fungi: the isotopic evidence. New Phytol. 2001; 150:601-610.
Ghizelini AM, Martins KG, Giesselmann UC, Santoro E, Pasqualette L, Mendonca-Hagler LCS et al. Fungal communities in oil contaminated mangrove sediments - Who is in the mud? Mar Pollut Bull. 2019; 139:181-188.
Colombier MA, Alanio A, Denis B, Melica G, Garcia-Hermoso D, Levy B et al. Dual Invasive Infection with Phaeoacremonium parasiticum and Paraconiothyrium cyclothyrioides in a Renal Transplant Recipient: Case Report and Comprehensive Review of the Literature of Phaeoacremonium Phaeohyphomycosis. J Clin Microbiol. 2015; 53:2084-2094.
Mostert L, Groenewald JZ, Summerbell RC, Robert V, Sutton DA, Padhye AA et al. Species of Phaeoacremonium associated with infections in humans and environmental reservoirs in infected woody plants. J Clin Microbiol. 2005; 43:1752-1767.
Amend A. From dandruff to deep-sea vents: Malassezia-like fungi are ecologically hyper-diverse. PLoS Pathog. 2014; 10:e1004277.
Zhou JZ, Ning D, L. Stochastic Community Assembly: Does It Matter in Microbial Ecology. Microbiol Mol Biol R. 2017; 81:e00002-00017.
Chase JM, Myers JA. Disentangling the importance of ecological niches from stochastic processes across scales. Philos Trans R Soc Lond B Biol Sci. 2011; 366:2351-2363.
Tedersoo L, Bahram M, Polme S, Koljalg U, Yorou NS, Wijesundera R et al. Fungal biogeography. Global diversity and geography of soil fungi. Science. 2014; 346:1256688.
LeBlanc N, Kinkel LL, Kistler HC. Soil fungal communities respond to grassland plant community richness and soil edaphics. Microb Ecol. 2015; 70:188-195.
Li W, Wang M, Pan H, Burgaud G, Liang S, Guo J et al. Highlighting patterns of fungal diversity and composition shaped by ocean currents using the East China Sea as a model. Mol Ecol. 2018; 27:564-576.
Pedrós-Alió C. The rare bacterial biosphere. Ann Rev Mar Sci. 2012; 4:449-466.
Logares R, Lindström ES, Langenheder S, Logue JB, Paterson H, Laybourn-Parry J et al. Biogeography of bacterial communities exposed to progressive long-term environmental change. ISME J. 2013; 7:937-948.
Zhang ZF, Cai L. Substrate and spatial variables are major determinants of fungal community in karst caves in Southwest China. J Biogeogr. 2019; 46:1504–1518.
Zhang W, Pan Y, Yang J, Chen H, Holohan B, Vaudrey J et al. The diversity and biogeography of abundant and rare intertidal marine microeukaryotes explained by environment and dispersal limitation. Environ Microbiol. 2018; 20:462-476.
Mo Y, Zhang W, Yang J, Lin Y, Yu Z, Lin S. Biogeographic patterns of abundant and rare bacterioplankton in three subtropical bays resulting from selective and neutral processes. ISME J. 2018; 12:2198-2210.
Zhang H, Hou F, Xie W, Wang K, Zhou X, Zhang D et al. Interaction and assembly processes of abundant and rare microbial communities during a diatom bloom process. Environ Microbiol. 2020; 22:1707-1719.
Wei G, Li M, Li F, Li H, Gao Z. Distinct distribution patterns of prokaryotes between sediment and water in the Yellow River estuary. Appl Microbiol Biotechnol. 2016; 100:9683-9697.
Liu L, Yang J, Yu Z, Wilkinson DM. The biogeography of abundant and rare bacterioplankton in the lakes and reservoirs of China. ISME J. 2015; 9:2068-2077.
Chen X, Wei W, Wang J, Li H, Sun J, Ma R et al. Tide driven microbial dynamics through virus-host interactions in the estuarine ecosystem. Water Res. 2019; 160:118-129.
Booth JM, Fusi M, Marasco R, Michoud G, Fodelianakis S, Merlino G et al. The role of fungi in heterogeneous sediment microbial networks. Sci Rep-UK. 2019; 9:7537.
Li X, Zhang Q, Ma J, Yang Y, Wang Y, Fu C. Flooding Irrigation Weakens the Molecular Ecological Network Complexity of Soil Microbes during the Process of Dryland-to-Paddy Conversion. Int J Env Res Pub He. 2020; 17:561.
Zhu HZ, Zhang ZF, Zhou N, Jiang CY, Wang BJ, Cai L et al. Diversity, Distribution and Co-occurrence Patterns of Bacterial Communities in a Karst Cave System. Front Microbiol. 2019; 10:1726.
Fabian J, Zlatanovic S, Mutz M, Premke K. Fungal-bacterial dynamics and their contribution to terrigenous carbon turnover in relation to organic matter quality. ISME J. 2017; 11:415-425.
Seto K, Van den Wyngaert S, Degawa Y, Kagami M. Taxonomic revision of the genus Zygorhizidium: Zygorhizidiales and Zygophlyctidales ord. nov. (Chytridiomycetes, Chytridiomycota). Fungal Syst Evol. 2020; 5:17-38.

SupplementaryFileS1IQtree.contree
File S1 Maximum likelihood (ML) phylogenetic tree based on the representative ITS sequences of all OTUs for taxonomic identification.
TableS8.xlsx
Table S8 Topological properties of co-occurring networks constructed for overall network and networks in different depth.
TableS7.xlsx
Table S7 Keystone taxa in the overall network and networks in different depth. Keystone taxa were selected on the basis of high degree, high closeness centrality (CC) and low betweenness centrality (BC) (degree > 6 and CC > 0.6 and BC < 0.12 for overall network, degree > 15 and CC > 0.4 and BC < 0.18 for network in different depth)
TableS6.xlsx
Table S6 Correlations between α-diversity and environmental variables.
TableS5.xlsx
Table S5 Indicator species of mangrove sediments samples in different depth.
TableS4.xlsx
Table S4 Annotation OTUs and distribution in different samples.
TableS3.xlsx
Table S3 Pairwise comparison of -diversity indices between samples in different mangroves. Numbers in the table are the p-values generated by One-way ANOVA followed by Tukey’s honest significant difference (HSD) method.
TableS2.xlsx
Table S2 Statistics of intermediate results of raw data processing.
TableS1.xlsx
Table S1 Information of samples and α-diversity.
Fig.S13.pdf
Fig. S13 Co-occurrence networks of fungal communities of samples in (a) 0–10 cm, (b) 10–20 cm, and (c) 20–30 cm depth layer. Nodes of networks in each depth are colored by modularity class, phylum and dominant/rare taxa respectively. Node size is proportional to the betweenness centrality (BC) of each node, and edge thickness is proportional to the weight of each correlation. Positive correlations are displayed in red lines and negative are in gray.
Fig.S12.pdf
Fig. S12 NMDS plots of community compositions of all taxa (a), dominant taxa (b), and rare taxa (c) in seven sampled mangroves, respectively. Dotted ellipses indicate 95 % confidence intervals around the centroids of three sediment depth and stress values are labeled at the upper left corner. Similarity values among the samples in different depths (Depth) were examined via the ANOSIM test, which are shown at the bottom left corner.
Fig.S11.pdf
Fig. S11 NMDS plots of community compositions of all taxa (a), dominant taxa (b), and rare taxa (c) in 0–10 cm, 10–20 cm, and 20–30 cm depth layer, respectively. Dotted ellipses indicate 95 % confidence intervals around the centroids of seven mangroves and stress values are shown at the upper left corner. Similarity values among the samples of different mangroves (“Location”) were examined via the ANOSIM test, which are shown at the bottom left corner
Fig.S10.pdf
Fig. S10 Bubble plots show the influences of spatial variables on the community compositions of all taxa, dominant taxa, and rare taxa in (a) all samples and sample depths, and (b) samples in each mangrove.
Fig.S9.pdf
Fig. S9 Comparisons of similarity of fungal community compositions of (a) all samples and (b) samples in different depth. p-values are labeled above comparisons.
Fig.S8.pdf
Fig. S8 Distance-decay patterns based on the geographic distance and bray-Curtis similarity of fungal community compositions of all taxa, dominant taxa and rare taxa of all samples, the shaded area around the lines covers 95% confidence interval of the correlations. Associated correlation coefficients and p-values are shown on each panel.
Fig.S7.jpg
Fig. S7 The relationships between NTI and changes in (a) mean annual temperature (MAT), (b) mean annual precipitation (MAP), (c) Gravel proportion, (d) salinity, (e) pH, (f) total carbon (TC), (g) total organic carbon (TOC), (h) total nitrogen (TN), (i) N/MH4+, (j) N/NO3-, (k) total phosphorus (TP), and (l) total sulfur (TS) are plotted. Associated correlation coefficients and p values are provided on each panel.
Fig.S6.pdf
Fig. S6 Boxplots of Bray–Curtis-based Raup–Crick (RCbray) in all and individual mangroves, horizontal dash lines indicate RCbray value of 0.95 and -0.95.
Fig.S5.pdf
Fig. S5 Fitting of the neutral community model (NCM) of community assembly in different mangroves. The solid blue lines indicate the best fit to the NCM, and the dashed blue lines represent 95% confidence intervals around the predication. OTUs that occur more or less frequently than predicted by the NCM are shown in different colors, m indicates the immigration rate, and R2 indicates the fitness to this model.
Fig.S4.pdf
Fig. S4 (a) OTUs proportions and (b) relative abundance of phylum observed in different samples, OTU proportions of (c) phylum, (d) class, (e) genus, and (f) species observed in samples in different depth.
Fig.S3.pdf
Fig. S3 Comparisons of (a) operational taxonomic unit (OTU) richness, (b) Shannon indexes and (c) chao1 indexes in different depth, as demonstrated by boxplots with the median and 95% confidence interval displayed. p-values are labeled above each comparison.
Fig.S2.pdf
Fig. S2 Rarefaction curves for (a) each sample and (b) different mangroves.
Fig.S1.pdf
Fig. S1 Locations of the seven sampled mangroves along the coast of Southeast China and the sampling sites in each mangrove.

Download PDF

Version 1

posted

You are reading this latest preprint version

Pacific Biosciences Single-molecule Real-time (SMRT) Sequencing Reveals High Diversity of Basal Fungal Lineages and Stochastic Processes Controlled Fungal Community Assembly in Mangrove Sediments

Status:

Version 1

Abstract

Figures

Background

Methods

Mangrove locations and sediment sampling

Environmental and spatial variables

Molecular analysis

Bioinformatics and phylogenetic analysis

Statistical analyses

Difference and indicator taxa of fungal communities

Alpha- and beta-diversity analysis

Fungal community assembly patterns

Fungal community driving factors

Network patterns and keystone taxa of fungal communities

Results

Diversity and composition of fungal communities in mangrove sediments

Fungal community assembly patterns in mangrove ecosystems

Spatial and environmental selection shaping the fungal community composition

Co-occurrence patterns of fungal communities in mangrove sediment

Discussion

Spatial and environmental selections on fungal communities in mangrove sediments

Co-occurrence network patterns and keystone taxa in fungal communities in mangrove sediments

Conclusions

List Of Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1