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 China 1 (Fig. 1). The YS, composed of the northern Yellow Sea (NYS) and the southern Yellow Sea (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 organisms 2,3 . 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 years 3 .
Previously reported data on the distribution of several plant and animal species indicate that species diversity levels decrease distinctly from the SYS to the NYS and BHS 1,4,5 . 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 biodiversity 1 . Recently, through the use of molecular tools, researchers have found that the diversity and community compositions of planktonic nanociliates 5 and bacteria 6 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 matter 7,8 . Moreover, temperature has been regarded as an important driver of fungal biogeography in marine environments 8,9 . 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 environments 10 , such as deep-sea sediments [11][12][13][14][15][16] , costal water and sediments 17 , and marine invertebrates 18 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, 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).

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 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 (R 2 = 0.4791, P = 0.001) and NYS (R 2 = 0.2932, P = 0.015). However, no significant difference between the NYS and BHS (R 2 = 0.2291, P = 0.055) was found.
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). 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).

Discussion
Fungal diversity in Chinese seas sediments is poorly understood 19 , 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 ecosystems 20 .
Of the 130 fungal genera uncovered in our study, most of them are widespread in terrestrial, freshwater 9,19 and marine habitats 8,21 . For example, Aspergillus, Fusarium, Penicillium, Talaromyces, and Trichoderma were found to be dominant in wetland sediments along the Yangtze River 22 , which discharges large amounts of sediment into the SYS 2 . Bionectria, Cercophora, Didymella, Gliocladium, Hypocrea, Hypoxylon, Nectria, Nigrospora, Oidiodendron, Ramichloridium and Trichoderma, which have recently been found in Chinese freshwater 19 , and the typical EcM fungal genera Rhizopogon and Sebacina 23 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 Eurotiomycetes found in large proportions in subtropical Chinese seas have also been frequently found in deep-sea sediments in India 24 , the Pacific Ocean 15 and Arctic fjords 16 , 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 habitats 21 . Alternaria 16 , Aureobasidium 21 , Cladosporium 16,21 and Phoma 13,24 , 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 invertebrates 18 .
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 nutrients 12 . Members of Agaricomycetes are dominant in anoxic mangrove sediments of Saint Vincent Bay 25 . 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 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 sediments 12,13,21,24 and European coastal sediments 17 , 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 environments 21,26,27 and coastal regions 17,28 . 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 fungi 20,29,30 . 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 sequences 13,20,30 . Interestingly, chytrid-like sequences   Table 1. Correlations between Bray-Curtis distances of fungal community and spatial or environmental variables checked through a Mantel test. 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. recovered from freshwater environments are most closely related to fungi that are known to parasitize algae and protists 20 , 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 samples 31,32 , show considerable mismatching with Chytridiomycota, Glomeromycota, Microsporidia, Saccharomycetes, several taxa within Dothideomycetes, and Tremellomycetes 29 . 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 region 33,34 . Although several new primers designed for the amplification of ITS1 32 and ITS2 32,35 would allow for the selective investigation of fungal communities from various environmental samples, no universally accepted approach to covering all fungi has been developed 29 .
Sediment sources of the BHS, NYS and 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 × 10 8 and 10 × 10 8 tons of sediment per year, respectively 2 . 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 wetlands 22,36,37 , and of bacterial biogeography in marine sediments 6,38 . As in these previous studies, the present study found that the beta-diversity of fungal communities in 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 assemblages 39 . 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 algae 4 and animals 1,5 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 bacteria 6,40,41 , viruses 42 and nanociliates 5 . For fungi living in marine environments, pH levels were found to significantly affect the number of colony forming units in certain North Sea locations 43 . Both pH and C/N were shown to affect the mycobiomes of soils sampled in the southern maritime Antarctic 44 . 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 environments 8,9 . 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 summer 1 . 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 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 levels 45 , 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 exclusion 46,47 . Methods Sediment sample collection. The YS, which includes the SYS and NYS, has a surface area of 380 × 10 9 m 2 and an average depth of 44 m with a maximum depth of 100 m in the center. The BHS covers a 77 × 10 9 m 2 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 m 3 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, 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 Scientific RepoRts | 6:26528 | DOI: 10.1038/srep26528 between the sites (with kilometer units) (Table S6) were calculated using the distm function of the geosphere package in R 48 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 FastDNA ® 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 rRNA 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 51 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 settings (Tedersoo, personal communication). The Uchime_ref command in USEARCH 52 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 reference 53,54 . 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 artifactual 55 .
The longest sequence of each OTU was selected as the representative for local BLASTn 2.2.30+ searches 56 with a conservative setting (word_size = 7, penalty = − 3, reward = 1) against UNITE + INSDC database. Rules employed for the taxonomic placement of OTUs 37 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, NYS and BHS) and three temperature categories (< 10 °C, 10 °C -15 °C and > 15 °C) using Cytoscape 3.0.1 57 .
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 SYS, NYS and BHS).
NMDS was used to compare dissimilarities among the fungal communities across three sea regions based on Bray-Curtis dissimilarity 58 . Prior to our NMDS analyses, sequence counts were normalized using the DESeq2 59 package in R 48 and were then transformed using a lg (x + 1) algorithm 60 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 variables 60 . 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 ANOVA 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 test 22 .