Environmental factors driving fungal distribution in freshwater lake sediments across the Headwater Region of the Yellow River, China

Dispersal limitation and environmental filtering are two primary processes involved in shaping microbial community structure. The pristine environmental and geographical relatively isolation of small lakes distributed in the Headwater Region of Yellow River (HRYR) offer a unique opportunity to test the relative roles of these two processes on fungal communities. Here, we investigated the fungal community in sediment samples from 10 lakes located in the HRYR using high-throughput sequencing. The results showed that the fungal community was dominated by Sordariomycetes, Leotiomycetes, Dothideomycetes, Pezizomycetes and Agaricomycetes. The results revealed that altitude, mean annual temperature, C/N ration, dissolve organic carbon and total nitrogen were the best predictors for shaping fungal community structure in these lakes. Significant spatial and environmental distance decay relationships in the fungal community were detected. The partial Mantel test indicated that the fungal community structure was significantly correlated with environmental distance but not with geographic distance. Overall, environmental filtering plays a more important role than dispersal limitation in fungal community structure at a local scale in such an pristine and isolated region.

Growing evidence indicates that fungal diversity and communities can alter a number of crucial ecosystem processes such as decomposition, nutrient cycling and productivity 1 . Understanding fungal biogeography and their underlying processes are essential for predicting the quality of ecosystem services under climate change 2,3 . Many studies have revealed that fungi display biogeographic patterns across spatial scales from centimeters to continents [4][5][6][7][8] . It has been reported that patterns of fungal diversity are controlled by dispersal limitation and environmental filtering 5,9,10 . However, the relative importance of these two factors and how they interact on different spatial scales with fungal communities remains unclear.
Several recent studies demonstrated that the processes controlling fungal communities vary by fungal groups, ecosystem types, and landscape heterogeneity or integrity. While dispersal limitation is a crucial in determining ectomycorrhizal fungal diversity at 10 m 2 to 10,000 m 2 scale 11 , environmental heterogeneity drives mycorrhizal fungal community variation at a finer scale 5 . On the other hand, the biogeography of arbuscular mycorrhizal fungi is driven by a combination of dispersal limitation and environment filtering [12][13][14] . Furthermore, the spatial influence on a fungal community is more pronounced in a contiguous site than in a disjointed site 4 . Therefore, to study integral fungal communities in varied ecosystems is essential to parse out the relative contributions of dispersal limitation and environmental filtering factors to fungal diversity. 1 State Key Laboratory of Mycology, Institute of Microbiology, Chinese Academy of Sciences, Beijing, 100101, China. 2

Results
Sediment geochemical characteristics. The major geographical and physicochemical characteristics of the sediments are listed in Table 1. Sediment pH varied from 7.71 to 8.84. Pairwise distances between sampled lakes were ranged from 4 to 148 km. There were significant unimodal trends of TN and dissolved organic carbon (DOC) with varying altitude (R 2 = 0.24 and R 2 = 0.28, respectively, P < 0.05). The C:N ratio were ranged from 16.3 to 23.0 across all study lakes.

Distribution of fungal taxa and phylotypes.
A total of 163,471 quality sequence reads were obtained from 33 samples and assigned to 479 operational taxonomic units (OTUs). These were representative of twenty-four fungal classes belonging to five fungal phyla (Fig. 2). The class that presented the highest diversity was Agaricomycetes (81 OTUs), followed by Dothideomycetes (68 OTUs) and Sordariomycetes (61 OTUs). The abundance of the OTUs in the dataset varied from 0.0006% to 13.1%. Among the total 479 OTUs, twenty-one were shared between all sites and accounted for a relative abundance between 46% and 86%. A class Sordariomycetes predominated across all sampled lakes with abundance ranging from 8.1-37.6%, followed by Leotiomycetes Dynamics of fungal diversity and community. The fungal Shannon index was not related to local environmental variables, while fungal diversity slightly increased with warming temperature (linear regression, R 2 = 0.09, P = 0.051). pH was positively correlated with fungal OTUs richness (linear regression, r = 0.33, A plot of community similarity against geographic distance and environmental distance for each pairwise set of fungal samples revealed a significant distance decay relationship (Fig. 4). Furthermore, the slope of geographic distance decay curve varied significantly among the three spatial scales. The distance-decay slope was steeper (slope = −0.87, P < 0.001) in an intermediate range between 50 to 100 km than a range between 100 to 150 km (slope = −0.75, P < 0.001) and the overall slope (slope = −0.12, P < 0.001) (Fig. 4). In contrast, the distance-decay curve did not differ from zero at the smaller scale (P = 0.09). Multivariate Mantel correlograms showed that pairwise comparisons of fungal community composition were autocorrelated at a distance of less than approximately 60 km ( Figure S1). Mantel tests indicated that the community of fungi was significantly correlated with either environmental or geographical distance (Table 2), where environmental distance accounted for a higher proportion of variation in community structure. However, partial Mantel tests indicated that the fungal community was significantly correlated with environmental distance, but not with geographic distance (Table 2).
Variance partitioning analysis (VPA) was performed to quantify the relative contributions of geographic distance and sediment environmental parameters on the taxonomic structure of the fungal communities. Combination of a subset of environmental parameters (DOC, MAT, TN and C/N ratio) and geographic distance explained 27% of the observed variation, leaving 73% of the variation unexplained. The sediment properties explained 20.0% (P = 0.001), and geographic distance alone explained 7.0% of the variations (P = 0.03), and no interaction effect was detected (Fig. 5). Although the sediment properties together explained more of the variation, geographic distance by itself explained 12.2% of the observed variation, more than any of the other four of the individual sediment variables (Fig. 5). Thus both environmental factors and geographic distance are important factors in shaping bacterial community structures.

Discussion
Sediment fungal diversity is intuitively thought to be low in aquatic habitats due to uniform vegetation and the physiological constraints of submersion in water 34 . In our study, unexpectedly high fungal diversity was detected Site code Latitude Longitude Altitude (m) pH TC (g/kg) DOC (g/kg) TN(g/kg) MAT (°C) Lake area (km 2 ) Shannon index Chao1 OTU number  in the sediments of HRYR which is consistence with recent study on fungal communities in sediments 35 . We found the fungal classes Sordariomycetes, Agaricomycetes, Dothideomycetes, Leotiomycetes and Pezizomycetes in large proportions, consistent with previous studies [35][36][37][38] . Most of these dominant classes are terrestrial in   natural, but also capable of being amphibious 9,20,25 . The majority of Sordariomycetes members are terrestrial and frequently detected near shores of tropical and subtropical freshwater ecosystems 39,40 . Most of these fungi break down lignin and cellulose from plant debris in the land-water interface 41 . The fungal order Polyporales, belonging to the class Agaricomycetes, was highly abundant; these are efficient terrestrial wood-decay fungi that play an important role in the carbon cycle 42 . Polyporales can be introduced into lakes through high fungal propagule loads via inflowing streams, rainwater, and wind 35,43 . The expansion of desaturase-encoding genes in the Polyporales may help them easily adapt to the substrate and colonize succession series 42,44 . Dothideomycetes, and particularly the order Pleosporales are found in many low-temperature lakes 45 , due to their resistance or adaptation to low temperature 37 , which is essential for survival in the harsh environment of the HRYR. These findings have implications for microbial biogeography, suggesting that in such habitat, environmental selection may be a more important determinant of community composition than dispersal constraints 2 .
At the spatial extent of our study, although both environmental characteristics and geographic proximity was significantly related to fungal community similarity, the partial correlation coefficient indicated that environmental filtering had a greater influence on fungal community. In our study, fungal community similarity was controlled by environmental factors such as TN and MAT. Nitrogen availability affects the fungal growth 46-48 and may limit the fungal diversity in the soil 49,50 . The QTPs' low temperatures lead to lake fungal community structure in the HRYR similar to what is observed in polar regions by Cox et al. 2 and Newsham et al. 51 . Marginal result of increasing fungal diversity with warming temperature may result from enhanced metabolic activity and enable a switch from survival to growth and dispersal strategies 51 . Tian et al. 52 identified temperature as one of major factor regulating fungal community structure in Qinghai-Tibetan Plateau. Hence, low temperature as an environmental constraint may determine fungal community composition in these extreme habitats. Interestingly, sediments pH significantly affected the fungal alpha diversity but not the fungal community similarity, which is consistence with report of Wang et al. 53 . With an increasing pH, some endurable taxa started to grow, resulting in an increase of the community richness. Although, we measured the impact of various environmental factors known or assumed to be an important in shaping fungal communities, but we were unable to explain all the variation. The unexplained variation may result from other unmeasured factors that could have an effect on fungal assemblages.
Distance-decay data can be used to understand the factors driving community turnover patterns such as dispersal limitation and environmental heterogeneity 54 . In this study, the slope of geographic distance-decay relationship (slope = −0.12, P < 0.001) was similar or steeper compared to the other reported studies 5,55,56 , indicating that dispersal limitation was also important in shaping fungal community. However, there was no significant distance-decay relationship observed at less than 50 km scale, indicates that dispersal limitation was less important in shaping fungal community. One plausible explanation is that dispersal should be more efficient at this spatial scales since long-distance dispersal of fungi mainly occur by wind 57 , resulting in weaker spatial structure of the communities and local control by environmental filtering 29 . Furthermore, in the QTP, there is more area dominated by saline lakes than freshwater lakes. Since the saline lakes are relatively isolated from large river basins and thus are less hydrologically connected 58 , it would be expected that dispersal limitation has a stronger control on fungal community in the saline lakes than in freshwater ones.
In conclusion, a large proportion of fungal sequences obtained in high frequency matched with terrestrial taxa, indicating that they are transported into the aquatic system and subsequently adapt to the local conditions. Mantel tests and variance partitioning analysis suggested that fungal richness and community composition are strongly influenced by environmental factors such as temperature and TN. Nonetheless, a significant spatial distance decay relationship and widespread locally abundant fungi indicated that dispersal limitation also plays a role in shaping fungal community. Different patterns of diversification might be observed depending on the ecology of the studied taxa. It is therefore important to investigate organisms that are typical for each biome occurring in the region of the QTP.

Methods and Materials
Sampling site. We selected 10 freshwater lakes with areas ranging from 3 to 611 km 2 in the HRYR, including Lake Erling, the largest lake in this region (Table 1). Sampling map was generated by ArcGIS software under license. The minimal and maximal distance between any two lakes were 4 and 148 km, respectively ( Fig. 1 and  Table 1). Most of the lakes were not readily interconnected to each other via water routes. The dominant vegetation type in the catchments of the 10 lakes was alpine grassland.
Sampling collection. We collected surface sediments up to 15 cm deep from inshore sites with water depths of approximately 1 m from each of the 10 freshwater lakes (Fig. 1). Three samples, spaced between 3 and 5 m, were collected from each lake. Samples were kept at 4 °C for sediment geochemical characterization and stored at −80 °C for genomic DNA extraction. Global positioning system coordinates were recorded at each sampling point (Table 1), and imported into the NOAA website (http://www.nhc.noaa.gov/gccalc.shtml) to calculate the pairwise geographic distance. The air temperature was measured in the field. Local mean annual temperature (MAT) was obtained from local weather stations. The pH was determined with a soil to water ratio of 1:5 wt/ vol. Dissolved organic carbon (DOC) and total nitrogen (TN) were measured with a vario TOC cube analyzer (ELEMENTAR Analysensysteme, Germany). Before analyzing sediment DOC, samples were first acidified with 1 N HCl overnight to remove carbonates and subsequently washed to neutralize pH, dried in an oven and then ground in a mortar. Sequence analysis. The raw sequence data were demultiplexed using the QIIME toolkit 59 . Sequence reads with a quality score less than 20 were filtered out (split_libraries_fastq.py; QIIME), and the remaining reads were assigned to samples according to the corresponding barcode. The extracted sequences were binned into operational taxonomic units (OTUs) at 97% match using USEARCH 60 , denoting a species-level similarity using closed-reference OTU picking using the UNITE ITS database 61 as a reference (pick_otus.py; QIIME). One representative of each OTU was selected for downstream analyses, and a taxonomy was assigned to each (pick_rep_set. py, assign_taxonomy.py; QIIME).

Statistical analyses.
Singletons were removed prior to statistical analyses. Diversity indices Chao1, ACE and Shannon (H′) were estimated based on the OTU table rarefied at 3229 counts per sample. We used linear and quadratic regression to test whether fungal diversity was related to lake area and environmental variables (TN, pH, DOC and MAT).
Distance decay of similarity has been widely used to infer the relative importance of dispersal limitation and environmental filtering in structuring microbial community 62 . To determine the distance decay relationship, pairwise fungal community similarities between samples were calculated using the Sørensen index, which is widely used for calculating distance decay in both micro-and macroecology. The spatial distance decay relationship was then assessed as the slope of logarithmic space to enhance the linear fitting, according to the description of Nekola & White (1999) 63 . All logarithmic transformations were performed using the natural logarithm. Mantel tests were performed to examine community turnover along environmental and geographical gradients. Euclidean distances were calculated for all site pairs separately based on the environmental variables and lake coordinates. Environmental data was comprised of total N, DOC, pH, water content and temperature. A significant correlation was observed between environmental and geographic distance in our study (Mantel r = 0.41, P = 0.001) (Fig. S2). Therefore, we performed partial Mantel tests (with 1000 permutations) to disentangle the effects of environmental distance on fungal community similarity while holding geographical distance constant, and vice versa 30 .
To fully characterize the exact scale of spatial clustering, a multivariate Mantel correlogram analysis was performed using the 'mantel.correlog' command in the R software. A Mantel correlogram draws a graph in which Mantel correlation values r are plotted as a function of the geographic distance classes among the sampling sites.
BioEnv and distance-based redundancy analysis (db-RDA) were also used to identify the abiotic factors influencing fungal community composition, and subsequently used to construct the sediment property matrix for variation partitioning analysis. Variance was partitioned to examine the relative effect of geographic distance and environmental variables (best environmental subset) on microbial community composition. All of the above analyses were performed in the vegan package in R software.