Fungal communities in sediments of subtropical Chinese seas as estimated by DNA metabarcoding.

Ribosomal RNA internal transcribed spacer-1 (ITS1) metabarcoding was used to investigate the distribution patterns of fungal communities and the factors influencing these patterns in subtropical Chinese seas, including the southern and northern Yellow Sea and the Bohai Sea. These seas were found to harbor high levels of fungal diversity, with 816 operational taxonomic units (OTUs) that span 130 known genera, 36 orders, 14 classes and 5 phyla. Ascomycota was the most abundant phylum, containing 72.18% and 79.61% of all OTUs and sequences, respectively, followed by Basidiomycota (19.98%, 18.64%), Zygomycota (1.10%, 0.11%), Chytridiomycota (0.25%, 0.04%) and Rozellomycota (0.12%, 0.006%). The compositions of fungal communities across these three sea regions were found to be vary, which may be attributed to sediment source, geographical distance, latitude and some environmental factors such as the temperature and salinity of bottom water, water depth, total nitrogen, and the ratio of total organic carbon to nitrogen. Among these environmental factors, the temperature of bottom water is the most important driver that governs the distribution patterns of fungal communities across the sampled seas. Our data also suggest that the cold-water mass of the Yellow Sea likely balances competitive relationships between fungal taxa rather than increasing species richness levels.

As subtropical Chinese seas, both the Yellow Sea (YS) and Bohai Sea (BHS) are semi-enclosed continental seas of the northwestern Pacific Ocean in northern China1 (Fig. 1). The YS, composed of the northern Yellow Sea (NYS) and the southern Yellow Sea (n class="Chemical">SYS), is commonly considered one of the broadest marginal seas in the world. The physical and biotic environments of these sea regions have been greatly affected by Korean and Chinese rivers that carry abundant organic materials and terrestrial organisms23. In addition, as marginal seas of China, these seas have been profoundly influenced by the economic development of coastal areas, and their ecosystems have clearly been degraded. As a result, the biodiversity of the YS and BHS regions is declining for various reasons, particularly overfishing and environmental deterioration (pollution and coastal construction) resulting from anthropogenic activity in recent years3.

Figure 1

Sampled sites in the SYS (the southern Yellow Sea), NYS (northern Yellow Sea) and BHS (Bohai Sea) of China.

ES = the East Sea. The figure was generated using Golden Software Surfer® 13 (http://www.goldensoftware.com/products/surfer) and Adobe Photoshop CC (https://www.adobe.com/products/photoshop.html).

Previously reported data on the distribution of several plant and animal species indicate that species diversity levels decrease distinctly from the SYS to the n class="Chemical">NYS and BHS145. In the central region of the YS, a cold water mass (CWM) formed in water remains within a bottom layer under the thermocline at 15–30 m depth in summer, playing a critical role in the maintenance of high levels of local bottom fauna biodiversity1. Recently, through the use of molecular tools, researchers have found that the diversity and community compositions of planktonic nanociliates5 and bacteria6 in these sea regions have been profoundly influenced by biogeographic distance and environmental factors, particularly the temperature of bottom water (TBW). Fungi represent an important group of microbes living in marine ecosystems that play a critical role in the decomposition and mineralization of organic matter78. Moreover, temperature has been regarded as an important driver of fungal biogeography in marine environments89. Thus, new questions have been raised regarding how the CWM of the YS influences fungal communities inhabiting the northern Chinese seas and how the YS is associated with ecological importance of fungi in marine ecosystems.

Over the past decade, molecular tools have been widely used to estimate fungal diversity levels and numerous fungal taxa have been retrieved from marine various environments10, such as deep-sea sediments111213141516, costal n class="Chemical">water and sediments17, and marine invertebrates18 among others. However, no data of molecular survey data on fungal diversity levels in these subtropical China seas, i.e., the SYS, NYS and BHS, are currently available. This lack of information on fungal diversity and distribution patterns has hindered the conservation and utilization of fungal resources. Using ITS1 metabarcoding, the current study aims to address the following questions: 1) What are the fungal diversity and distribution patterns characteristic of these seas? 2) Are the effects of spatial and environmental factors on fungal communities similar to those of other organisms (for CWM in particular)?

Results

Taxonomic assignment and community compositions

In total, 122,448 raw reads with a median length of 321 bases were obtained from 30 samples with pyrosequencing on one 1/4-plate. After the filtering process, 102,990 quality-filtered ITS1 reads assigned to 3,404 OTUs with 1,056 singletons were found from 30 samples. After the OTUs with fewer than five reads were removed, 816 OTUs (58.67%) of 79,847 reads (80.37%) were assigned to the fungal kingdom and 575 OTUs (41.33%) of 19,498 reads (19.63%) failed to be assigned to any known kingdom (Tables S2 and S3).

For the fungal kingdom, the phylum-level assignment of 52 OTUs remained elusive, and the other 764 OTUs spanned five fungal phyla, 14 known classes, 36 orders, 77 families, and 130 genera (Fig. 2, Tables S4 and S5). For the SYS, n class="Chemical">NYS and BHS, neither OTU richness (F = 0.44, P = 0.648) nor the Shannon index (F = 0.637, P = 0.537) differed according to a one-way ANOVA (Fig. 3). In addition, no significant correlation was observed between any one of the variables tested for OTU richness (P > 0.05).

Figure 2

Taxonomic distributions of (a) 816 OTUs and (b) 79,847 sequences at the fungal phylum and order levels across the entire region sampled.

Figure 3

(a) OTU richness per site and (b) Shannon index in three sea regions as demonstrated by boxplot with median and 95% confidence intervals displayed. Bars without shared letters indicate significant differences at the level of P value = 0.05.

In Ascomycota, Sordariomycetes (23.78% of all of the OTUs, 28.27% of all of the sequences) was found to be the most abundant class, followed by Dothideomycetes (25.25%, 26.49%) and n class="Chemical">Eurotiomycetes (19.24%, 23.75%). These abundant classes showed coherent distributions in different sampled sites with relatively higher abundance compared to other classes (Figure S1, Table S4). At the order level, Pleosporales was found to be the most abundant order, accounting for 10.91% and 17.60% of all of the OTUs and sequences, respectively, followed by Chaetothyriales (11.76%, 16.61%), Capnodiales (12.50%, 8.25%), Eurotiales (6.25%, 6.24%) and Xylariales (6.50%, 5.28%) (Fig. 2).

Within Basidiomycota, Agaricomycetes (6.25% of all of the OTUs, 8.48% of all of the sequences) was the most abundant class, followed by Tremellomycetes (7.35%, 6.21%) and Cystobasidiomycetes (3.19%, 2.29%). The least abundant classes were Agaricostilbomycetes (0.12%, 0.014%), n class="Chemical">Wallemiomycetes (0.25%, 0.041%) and Pucciniomycetes (0.25%, 0.015%) (Table S4). The most abundant order was found to be Polyporales (3.43%, 7.67%), followed by Tremellales (7.11%, 6.17%), Erythrobasidiales (2.94%, 2.11%) and Sporidiobolales (2.57%, 1.59%) (Fig. 2).

The following 16 genera were found to be the most common (>10 OTUs were present in all of the samples) in the sampled regions: Penicillium (31 OTUs), Trichomerium (30 OTUs), Nigrospora (28 OTUs), Shiraia (25 OTUs), Erythrobasidium (22 OTUs), Mycosphaerella (19 OTUs), Cryptococcus (18 OTUs), Rhodotorula (17 OTUs), Fusarium (15 OTUs), Eutypa (14 OTUs), Ganoderma (14 OTUs), Perenniporia (14 OTUs), Phialemonium (14 OTUs), Metschnikowia (11 OTUs), Sarocladium (11 OTUs) and Phialophora (11 OTUs) (Table S5). Other genera were frequently isolated from marine n class="Chemical">sediments, including Aspergillus (10 OTUs in 25 samples), Hypoxylon (10 OTUs in 30 samples), Libertella (10 OTUs in 30 samples) and Xylaria (10 OTUs in 30 samples).

Spatial distribution of community compositions and the effects of variables

Among the 816 OTUs examined, 241 OTUs were only found in one or two regions (Fig. 4a), indicating a notable variability of species compositions across these regions. The plot of non-metric multidimensional scaling (NMDS) shows that the spatial distribution of fungal community compositions varied considerably among the n class="Chemical">SYS, NYS and BHS (Fig. 5). Our ANOSIM analysis results also show that the fungal community in the SYS differs significantly from that in the BHS (R2 = 0.4791, P = 0.001) and NYS (R2 = 0.2932, P = 0.015). However, no significant difference between the NYS and BHS (R2 = 0.2291, P = 0.055) was found.

Figure 4

Network diagram of (a) 241 OTUs occurring only within one or two regions of the BHS (red color), NYS (green color) and SYS (blue color), and (b) 174 OTUs occurring only within one or two temperature categories of <10 °C (blue color), 10–15 °C (green color) and >15 °C (red color).

Figure 5

Non-metric multidimensional scaling (NMDS) plot of Bray-Curtis dissimilarity on abundance data transformed by a lg (x + 1) algorithm where x was designed as the number of normalized counts by DESeq2 package.

The BHS, NYS and SYS are represented with the black, red and green circles, respectively. Dotted ellipses indicate 95% confidence intervals around centroids of three sea regions.

Community compositions in the sampled regions were heavily shaped by most of the variables tested according to our db-RDA analysis, with the exception of total organic carbon (TOC) and pH (Fig. 6a). n class="Chemical">Nine variables explained 39.24% of variances (F = 1.4353, P = 0.001), indicating that these variables can partially account for fungal community dynamics across our sampled regions. After variable selection by AIC, TBW and latitude were found to be the main variables associated with dissimilarity changes in community compositions (Fig. 6b), explaining 16.49% of the variance (F = 2.6653, P = 0.001). The total variances were significantly explained by TBW (accounting for a 7.43% proportion (F = 2.4026, P = 0.001)) and latitude (accounting for a 9.06% proportion (F = 2.9279, P = 0.001)). Interestingly, fungal communities in different regions presented various responses to the same variable (Table 1). For example, TBW was found to be significantly correlated with fungal communities in the YS (P = 0.02) but not with the entire region or the BHS (P > 0.05). Water depth (WD) and salinity showed significant correlations with the entire region (P = 0.011, P = 0.003) but not with the YS and BHS (P > 0.05). The ratio of total organic carbon to nitrogen (C/N) only significantly affected community compositions in the BHS (P = 0.026).

Figure 6

Distance based redundancy analysis (db-RDA) of Bray-Curtis distance on normalized and transformed abundance data.

Three sampled regions are indicated with the black (the BHS), red (the NYS) and green circles (the SYS), respectively. Rectangles represent spatial or environmental factors. (a) The db-RDA ordination without variable selection tests. (b) The db-RDA ordination after variable selection tests. Gray lines indicate surface fitting of TBW.

Table 1

Correlations between Bray-Curtis distances of fungal community and spatial or environmental variables checked through a Mantel test.

		Spatial variables			Environmental variables
		GD	Longitude	Latitude	TBW	WD	TOC	TN	C/N	pH	Salinity
Entire region	r	0.2959	0.2431	0.3220	0.0976	0.1535	−0.0032	0.0232	0.0647	−0.0591	0.2642
	P	0.001*	0.007*	0.004*	0.133	0.011*	0.464	0.313	0.221	0.666	0.003*
YS	r	0.6630	0.2285	0.6956	0.3355	0.0094	0.1847	0.0891	0.1228	−0.0200	0.1656
	P	0.001*	0.039*	0.001*	0.002*	0.415	0.053	0.246	0.095	0.550	0.064
BHS	r	0.3665	0.0185	0.5371	0.1091	−0.1256	−0.0746	0.2464	0.2943	−0.0550	−0.0454
	P	0.028*	0.299	0.016*	0.232	0.672	0.627	0.121	0.026*	0.552	0.601

GD = geographical distance (km), TBW = temperature of bottom water (°C), WD = water depth (m), TOC = total organic carbon (mg/g), TN = total nitrogen (mg/g), C/N = the ratio of total organic carbon to nitrogen. The asterisk indicates significance at the level of P value = 0.05.

Discussion

Fungal diversity in Chinese seas sediments is poorly understood19, which impedes our understanding of the ecological importance of fungi in the marine ecosystem of China. In this study, 816 fungal taxa at the OTU level were revealed, greatly advancing our understanding of fungal diversity in Chinese seas. Moreover, several sequences (more than one fifth of all of the sequences) matched the database poorly and failed to be assigned to any known taxa, indicating the existence of considerable unknown biodiversity. These poorly recognized groups revealed by molecular approaches may prove to be important components of these ecosystems20.

Of the 130 fungal genera uncovered in our study, most of them are widespread in terrestrial, freshwater919 and marine habitats821. For example, Aspergillus, Fusarium, Penicillium, Talaromyces, and Trichoderma were found to be dominant in wetland n class="Chemical">sediments along the Yangtze River22, which discharges large amounts of sediment into the SYS2. Bionectria, Cercophora, Didymella, Gliocladium, Hypocrea, Hypoxylon, Nectria, Nigrospora, Oidiodendron, Ramichloridium and Trichoderma, which have recently been found in Chinese freshwater19, and the typical EcM fungal genera Rhizopogon and Sebacina23 were also recovered from our samples. Moreover, several of the genera detected are plant pathogens, and some are reasonably assumed to be dead or dormant in sediments of the sampled seas. Our data suggest that compositions of the fungal community are profoundly affected by the passive dispersal of spores and other propagates from terrestrial environments through terrestrial runoff and riverine inputs.

In Ascomycota, members of Sordariomycetes, Dothideomycetes and n class="Chemical">Eurotiomycetes found in large proportions in subtropical Chinese seas have also been frequently found in deep-sea sediments in India24, the Pacific Ocean15 and Arctic fjords16, indicating that these classes are ubiquitous in marine sediments. Members of Dothideomycetes, and particularly the order Pleosporales, are believed to have adapted or developed resistance to low temperatures and high levels of osmotic pressure, which may be essential to survival in marine habitats21. Alternaria16, Aureobasidium21, Cladosporium1621 and Phoma1324, which have often been found in deep-sea environments, were consistently observed in our study. In Eurotiomycetes, Aspergillus and Penicillium are known to be ubiquitous in various marine substrates as well as in marine sponge invertebrates18.

Fungal phylotypes of Basidiomycota are the most frequently recovered eukaryotes from deep sediment samples of the Peru Margin, and some species may play important roles in the utilization and recycling of nutrients12. Members of Agaricomycetes are dominant in anoxic mangrove sediments of Saint Vincent Bay25. In the seas sampled in our study, the most abundant class of Basidiomycota was Agaricomycetes, primarily represented by Polyporales sequences. This finding is unusual, as most species of this order are typical terrestrial wood decomposers and do not live in ocean n class="Chemical">sediments. Therefore, these Polyporales sequences targeted by sediments are assumed to be dead or possibly dormant spores that fell from the air or that were perhaps flushed into seas via rivers or terrestrial runoff. The yeast genera Bullera, Cryptococcus, Hannaella, Rhodotorula and Sporobolomyces, which are often found in deep-sea sediments12132124 and European coastal sediments17, accounted for 37.63% of all of the Basidiomycota sequences examined in this study. Two yeast genera, Dioszegia and Erythrobasidium, were uncovered from marine environments for the first time.

Chytridiomycota and Zygomycota detected at low levels in this study were also recovered from deep-sea environments212627 and coastal regions1728. Members of the two fungal phyla may be considerably more diverse and numerous in marine habitats than we initially imagined, as the PCR processes and primers used to generate most marine clone libraries are biased toward Dikarya fungi202930. Moreover, Chytridiomycota and groups accounting for mostly unicellular species are difficult to infer for their phylogeny based on ITS regions due to high levels of genetic divergence and the current paucity of reference sequences132030. Interestingly, chytrid-like sequences recovered from freshn class="Chemical">water environments are most closely related to fungi that are known to parasitize algae and protists20, indicating that the taxa of Chytridiomycota found in this study are likely pathogens of other benthic organisms.

It is well known that various commonly used ITS primers may introduce biases during the amplification of different parts of ITS regions. As fungus-specific primers, ITS1F and ITS2, which facilitate the selective amplification of fungal ITS1 fragments from mycorrhizal and other environmental samples3132, show considerable mismatching with Chytridiomycota, Glomeromycota, Microsporidia, Saccharomycetes, several taxa within Dothideomycetes, and Tremellomycetes29. n class="Chemical">Nevertheless, a few studies have showed that the ITS1 region targeted by ITS1F/ITS2 primers can present similar patterns in the fungal community structure as those exhibited by the ITS2 region3334. Although several new primers designed for the amplification of ITS132 and ITS23235 would allow for the selective investigation of fungal communities from various environmental samples, no universally accepted approach to covering all fungi has been developed29.

Sediment sources of the BHS, NYS and n class="Chemical">SYS are mainly from Korean and Chinese rivers. Among these rivers, the Yangtze River and Yellow River are the main sources, contributing 4.7–5 × 108 and 10 × 108 tons of sediment per year, respectively2. Sediments of the southern part of the YS mainly derive from the Yangtze River, whereas sediments of the southern area of the BHS and the northern part of the YS are mainly derived from the Yellow River. Therefore, it is not surprising that we found clear difference among fungal communities in the BHS and NYS and SYS (Figs 5 and 6). Moreover, the fact that sites located near estuaries of the Yellow (sites 27, 28, 29 and 30) and Yangtze Rivers (sites 1, 2, 4, 5 and 6) were grouped separately into two distinct clusters (Figure S2) clearly highlights the effects of river discharge.

Geographical distance is an important driver of effective fungal community structuring in soils, rhizospheres and wetlands223637, and of bacterial biogeography in marine sediments638. As in these previous studies, the present study found that the beta-diversity of fungal communities in n class="Chemical">sediments of subtropical Chinese seas was significantly correlated with geographic distance (Table 1), supporting the hypothesis that dispersal limitations of taxa within habitats constitute another critical factor that shapes space variations in species assemblages39. However, the SYS, NYS and BHS were not found to differ significantly in terms of OTU richness (P > 0.05) and the Shannon index (P > 0.05), which was found to differ between marine algae4 and animals15 in the sampled seas (biodiversity levels increased from the BHS to the YS).

In addition to spatial factors, the physiochemical compositions of marine habitats were found to heavily affect microbial communities, such as bacteria64041, viruses42 and nanociliates5. For fungi living in marine environments, pH levels were found to significantly affect the number of colony forming units in certain North Sea locations43. Both pH and C/n class="Chemical">N were shown to affect the mycobiomes of soils sampled in the southern maritime Antarctic44. In contrast to these results, in the present study, pH was found to have no significant effect on either fungal richness levels or distribution patterns in our sampled seas, whereas other environmental variables such as WD, TBW, salinity, total nitrogen (TN) and C/N significantly influenced the spatial distribution patterns of fungal communities (Fig. 6a, Table 1). In addition, the results of the db-RDA variable selection tests (Fig. 6b) are consistent with the conclusion that temperature plays a critical role in the determination of fungal biogeography in marine environments89.

Typically, the CWM of the YS emerges in winter and remains within a bottom layer under the thermocline at a depth of 15–30 m in summer1. During our sampling period, the effects of CWM on the TBW in the local sea region were obvious and are supported by the fact that the sampled sites located in the central region of the YS showed lower n class="Chemical">TBW levels than the other sites (Table S1). As a result, these sites (site 8 to site 13, site 15, site 16 and site 18) with TBW levels of less than 10 °C are closely related within the fungal community (Figs 6 and S2). Generally, species richness will decline with a decrease in temperature levels45, however, no obvious decline in fungal OTU richness in sites with lower TBW levels was observed (P > 0.05). Our network analysis shows that 174 OTUs (21.32% of total OTUs) appeared in one or two temperature categories (Fig. 4b), revealing that a clear variation of species compositions occurred along temperature gradients. These 41 OTUs occurring only within a certain specific temperature category may be taxa with narrow temperature adaptation features. In this case, we hypothesize that CWM may play a role in balancing the competitive relationship between cold- and warm-adapted species, proving critical to the alternation of species composition rather than to the enhancement of species richness. For instance, microbiome experiments have suggested that this increase in temperature may reduce species abundance in some cases through competitive exclusion4647.

Methods

Sediment sample collection

The YS, which includes the SYS and n class="Chemical">NYS, has a surface area of 380 × 109 m2 and an average depth of 44 m with a maximum depth of 100 m in the center. The BHS covers a 77 × 109 m2 sea area and has an 18 m average water depth with a 70 m maximum depth in the northern section of the Bohai Strait. Thirty sediment samples (Fig. 1, Table S1) were collected from the SYS, NYS and BHS using a 0.05 m3 stainless steel grab sampler (Wildco®, Florida, USA) in June 2013 aboard the Dongfanghong 2# research vessel. The collected sediment was undisturbed and compact. Ten sub-cores of sediment were collected using 3.5 cm-inner-diameter pre-cleaned glass tubes sealed with sterile plastic film from each sampler. These glass tubes filled with sediment were immediately stored at −80 °C until total DNA extraction.

Environmental and spatial variables

The spatial and environmental variables included sediment longitude, latitude and physiochemical properties (e.g., TOC, Tn class="Chemical">N, C/N, pH, WD, TBW and salinity of bottom water) (Table S1). Salinity, TBW and WD were measured using a Seabird 911 conductivity-temperature-depth (CTD) instrument. TOC and TN in sediments were determined using an Elemental Analyzer (EA3000, Euro Vector SpA, Milan, Italy) in our laboratory. The pH was measured by adding 10 ml of distilled water to 4 g of sediment and recording pH levels using a pH electrode (STARTER 300, OHAUS, Beijing, China). Geographic distances between the sites (with kilometer units) (Table S6) were calculated using the distm function of the geosphere package in R48 based on the longitude and latitude of the sites.

DNA extraction and 454-pyrosequencing

In the laboratory, aliquots of sediment (2–8 cm in depth) from the central parts of the sub-cores of each sampler were removed using a sterile spatula and mixed thoroughly. Genomic DNA extraction was performed on 0.5 g of composite sediment using the FastDn class="Chemical">NA® Spin kit for soil (MP Biomedicals, Solon, OH, USA) according to the manufacturer’s instructions. The remaining parts of the pooled samples were stored for physicochemical analysis.

The DNA samples were amplified using fungal-specific primer pairs ITS1F (5′-(A) CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2 (5′-(B) GCTGCGTTCTTCATCGATGC-3′) to generate PCR ITS1 rRn class="Chemical">NA fragments of approximately 400 bp, where A and B represent the two fusion primers (GCCTCCCTCGCGCCATCAG and GCCTTGCCAGCCCGCTCAG)49. The conditions for PCR included an initial hot start incubation (5 min at 94 °C) followed by 34 cycles of denaturation at 94 °C for 30 s, annealing at 55 °C for 30 s and extension at 72 °C for 1 min, followed by a final extension at 72 °C for 15 min (GS Titanium LV emPCR Purification Kit, TaKaRa DR500A and D4030A). The PCR products were pyrosequenced using a Genome Sequencer FLX 454 System (454 Life Sciences/Roche Applied Biosystems, Nutley, NJ, USA). The 454 SFF files and corresponding barcode files were deposited in the National Center for Biotechnology Information Sequence Read Archives (SRA) as BioProject ID SRP067807.

Bioinformatics analysis

Raw sequences were assigned to samples based on the barcodes and were trimmed (qaverage = 30, qwindowsize = 50, minlength = 200, maxambig = 1, maxhomop = 8, bdiffs = 1) with the trim.seqs command in Mothur 1.31.250. Then, ITSx 1.0.1151 was applied to trim out 18S and 5.8S rRNA genes and remove the compromised and nontarget sequences with a minlength setting for ITS1 of 50 bases and other parameters based on default settinn class="Disease">gs (Tedersoo, personal communication). The Uchime_ref command in USEARCH52 was used to detect and remove chimeras using the unified system for the DNA-based fungal species (UNITE) and international nucleotide sequence databases (INSDC) fungal ITS databases (the version released on March 11, 2015) as reference5354. These filtered sequences were clustered into operational taxonomic units (OTUs) at a 97% sequence similarity threshold using CD-Hit (www.cd-hit.org). OTUs represented by fewer than 5 sequences were removed, as these OTUs tend to be artifactual55.

The longest sequence of each OTU was selected as the representative for local BLASTn 2.2.30+ searches56 with a conservative setting (word_size = 7, penalty = −3, reward = 1) against UNITE + In class="Chemical">NSDC database. Rules employed for the taxonomic placement of OTUs37 were adopted here. The BLASTn e-value < e−50 was used as cutoff to assign a sequence to the fungal kingdom. However, query sequences with an e-value between e−20 and e−50 were manually checked against the 100 best-matching sequences for accurate assignment. These above queries were assigned to the fungal kingdom and were labeled by genus, family, order, or class at 90, 85, 80, and 75% sequence identity, respectively. The e-value > e−20 was used as threshold to classify sequences as “unknown”.

Networks

To display putative major OTUs associated with sea regions and temperatures, two networks were produced to visualize shared OTUs among three sea regions (the SYS, n class="Chemical">NYS and BHS) and three temperature categories (<10 °C, 10 °C –15 °C and >15 °C) using Cytoscape 3.0.157.

Statistical analyses

OTU richness (S), Shannon’s diversity (H) and evenness (H/lnS) values with sequence counts were calculated. A one-way ANOVA followed by a pairwise t-test was used to explore variations in the richness of OTUs and Shannon’s diversity levels across the sea regions (the n class="Chemical">SYS, NYS and BHS).

NMDS was used to compare dissimilarities among the fungal communities across three sea regions based on Bray-Curtis dissimilarity58. Prior to our n class="Chemical">NMDS analyses, sequence counts were normalized using the DESeq259 package in R48 and were then transformed using a lg (x + 1) algorithm60 where x was designated as the number of normalized counts. The anosim function in the vegan package of R was applied to evaluate changes in beta-diversity levels based on Bray-Curtis dissimilarity and permutations = 999.

A distance-based redundancy analysis (db-RDA) with Bray-Curtis distances was performed to visualize the relationship between compositional variations and environmental or spatial variables60. To further identify important variables associated with community structures, the step function in R with Akaike’s information criterion (AIC) was used to build an alternative db-RDA model. The significance of the models, axes and terms were tested using An class="Chemical">NOVA function in the vegan package. All samples were clustered using the Euclidean distance of scores for the first two axes of the two above listed db-RDA models. In addition, correlations between Bray-Curtis dissimilarities of the fungal community and spatial or environmental variables in different sea regions were checked through a Mantel test22.

Additional Information

How to cite this article: Li, W. et al. Fungal communities in sediments of subtropical Chinese seas as estimated by Dn class="Chemical">NA metabarcoding. Sci. Rep. 6, 26528; doi: 10.1038/srep26528 (2016).