The scale dependence of fungal community distribution in paddy soil driven by stochastic and deterministic processes

Fungal Ecology 42 (2019) 100856

Contents lists available at ScienceDirect

Fungal Ecology journal homepage: www.elsevier.com/locate/funeco

The scale dependence of fungal community distribution in paddy soil driven by stochastic and deterministic processes Jianshu Zhao a, Qun Gao a, Jizhong Zhou a, c, d, e, **, Mengmeng Wang a, Yuting Liang b, Bo Sun b, Haiyan Chu b, Yunfeng Yang a, * a

State Key Joint Laboratory of Environment Simulation and Pollution Control, School of Environment, Tsinghua University, Beijing, China State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing, China Institute for Environmental Genomics and Department of Microbiology and Plant Biology, University of Oklahoma, Norman, OK, USA d Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, USA e School of Civil Engineering and Environmental Sciences, University of Oklahoma, Norman, OK, USA b c

a r t i c l e i n f o

a b s t r a c t

Article history: Received 7 November 2018 Received in revised form 16 July 2019 Accepted 28 July 2019 Available online xxx

Fungal communities play important roles in terrestrial ecosystem functioning. Unraveling the relative importance of stochastic versus deterministic processes in shaping biogeographic patterns of fungal communities has long been a challenge in microbial ecology, owing to high biodiversity and difficulties in identifying fungal taxa. Using a unique anthropogenic system of geographically isolated paddy ‘islands’, we collected 198 soil samples with a spatially explicit design to examine how ecological processes shaped fungal biogeographic patterns. Fungal community structure showed scale-dependent distancedecay relationships. Stochastic processes (dispersal and drift) contributed more to community assembly than deterministic processes (selection) at the local scale, which was largely attributed to drift. In contrast, deterministic processes contributed more to community assembly than stochastic processes at the regional scale, with soil dissolved organic carbon being the most important measured factor. Collectively, scale dependence of fungal biogeographical patterns in paddy soils is influenced by differential contribution of deterministic and stochastic processes. © 2019 Elsevier Ltd and British Mycological Society. All rights reserved.

Corresponding Editor: Petr Baldrian Keywords: Soil fungal biogeography Scale dependence Community assembly Environmental selection Dispersal limitation Drift

1. Introduction Information of biodiversity patterns can provide valuable insights into the mechanisms influencing community assembly (Hanson et al., 2012). However, biodiversity patterns are especially difficult to be assessed in the microbial world owing to the overwhelming diversity (Widder et al., 2016), which can be 20,000e50,000 species per gram of soil (Roesch et al., 2007), and ~1012 species when taking all the soils on Earth into consideration (Locey and Lennon, 2016). Among various soil microbial groups, fungi play important roles in ecosystem functions such as soil organic carbon storage, decomposition and primary productivity of plants (Van Der Heijden et al., 2006; Cheng et al., 2012; Clemmensen et al., 2013). Therefore, it is essential to elucidate the

biogeographical patterns and underlying community assembly mechanisms of soil fungal communities. Global fungal diversity is estimated to be 2.2e3.8 million species (Hawksworth and Lücking, 2017). The diversity of most fungal groups, except for ectomycorrhizal fungi and several fungal classes, peaks in tropical ecosystems, which is partly attributed to particularly strong fungal endemicity in tropical regions (Tedersoo et al., 2014). Although there was strong spatial variability (also termed as spatial structure or distance decay) within fungal communities at the continental scale of North America (Talbot et al., 2014), the distance-decay relationship of fungal communities was weak at the local scale (Bahram et al., 2016). The distance-decay relationship of ectomycorrhizal fungi varied in different ecosystems and was scaledependent (Bahram and Tedersoo, 2013). However,

* Correponding author. ** Corresponding author. State Key Joint Laboratory of Environment Simulation and Pollution Control, School of Environment, Tsinghua University, Beijing, China. E-mail addresses: [email protected] (J. Zhou), [email protected] (Y. Yang). https://doi.org/10.1016/j.funeco.2019.07.010 1754-5048/© 2019 Elsevier Ltd and British Mycological Society. All rights reserved.

2

J. Zhao et al. / Fungal Ecology 42 (2019) 100856

ectomycorrhizal fungi generally exhibited stronger distance-decay relationships of community similarity than saprotrophic fungi and corresponded to aboveground plants, which were consistent with observations elsewhere (Pickles et al., 2010; Glassman et al., 2017). Traditionally, determinants of microbial biogeographic patterns are analyzed under the paradigm of niche theory, i.e., “Everything is everywhere, but the environment selects”. Considering their prodigious dispersal capabilities, microbes have been thought to be ubiquitous (De and Bouvier, 2010). Therefore, environment selection is primarily responsible for site-to-site variation in species composition. Important environmental drivers such as salinity (Mohamed and Martiny, 2010) and pH (Dumbrell et al., 2010; Glassman et al., 2017) have been shown to strongly affect fungal community composition. However, fungal communities are also strongly impacted by processes such as dispersal, ecological drift and stochastic recruitment (Peay et al., 2010, 2012b; Bahram and Tedersoo, 2013). Consequently, both deterministic and stochastic processes can influence fungal communities, and relative contributions of those processes may vary between different ecosystem types and properties (Beck et al., 2015). Understanding the relative contribution of each process under different circumstances remains challenging, which is further complicated by the possibility that fungal community assembly processes may be spatial-scale
lyi et al., 2016). Also, there are techdependent (Mcgill, 2010; Va nical issues, e.g., commonly used rRNA biomarkers for identifying fungal species (e.g. 18S Small Subunit Ribosomal RNA) are coarse in resolution (Bruns and Taylor, 2016). These challenges can lead to difficulty in generalizing consistent ecological patterns and understanding underlying mechanisms (Verbruggen et al., 2017). Biogeographic patterns and the underlying mechanisms are more likely to emerge by focusing on similar habitats across different spatial scales (Dolan, 2006; Bahram et al., 2016). Therefore, we examined fungal communities in paddy soils using a spatially explicit experimental design, which was used to examine fungal community in forest soil ecosystems (Pickles et al., 2010), from 1 m to 668 km (Supplemental Fig. S1). On one hand, dispersal rates in strongly connected habitats may be high so as to interfere with local environmental selection, leading to homogenization through mass effects, i.e., taxa occurring in unsuitable habitats because of continuous immigration (Leibold and Norberg, 2004). Studying paddy fields can resolve the issue because they are geographically isolated, unique anthropogenic ecosystems with strong dispersal limitation. On the other hand, typical paddy fields are subjected to a long history of anthropogenic activities such as flooding and application of fertilizer, leading to strong environmental selection. Environmental heterogeneity caused by different rice planting regimes (e.g. the use of different rice species and agricultural managements) can also enrich habitat-specific fungal species, leading to different soil fungal communities (Lumini et al., 2011; Verbruggen et al., 2017). Green et al. (2004) showed that environmental selection was a more important driver of biogeographic patterns at small spatial scales (1 me10 km), but stochastic drivers of dispersal limitation and/or ecological drift were more important drivers for large-scale spatial patterns (>100 km). However, recent studies analyzing fungi communities have shown that environmental selection is more important at the large scale while stochastic processes are more important at the local scale (Peay et al., 2010; Bahram and Tedersoo, 2013; Bahram et al., 2016; Tisthammer et al., 2016). To clarify, we aim to address the following scientific questions using rice paddies distributed throughout China as a model system: (1) Do biogeographical patterns of fungal communities display scale-dependent distance-decay relationships? (2) How do community assembly processes interactively shape fungal biogeography at different spatial scales? and (3) What is the most important environmental factor measured in shaping fungal community structure?

2. Materials and methods 2.1. Site description and sample collection We collected a total of 198 soil samples from 18 paddy fields in Eastern China. Three paddy fields were sampled from each of the six counties, including Changxing, Zhejiang Province, Jurong, Jiangsu Province, Hefei, Anhui Province, Tongcheng, Anhui Province, Huangshi, Hubei Province and Xiantao, Hubei Province (Supplemental Fig. S1 (a)). The sample location was documented by GPS. A total of 11 samples were taken in each paddy field. Specifically, six soil samples were collected along a 76 m-long transect, and another five samples were collected along a vertical 76 m-long transect (Supplemental Fig. S1 (b)). Distances between two adjacent samples along both transects were 1 m, 5 m, 10 m, 20 m and 40 m, respectively. Samples taken from the same paddy field were considered to be at the local scale (1e113 m). Samples taken from different counties were considered to be at the regional scale (103e668 km). Samples taken in the same county but not in the same paddy field were considered to be at the meso scale (3.4e39 km, defined as an intermediate scale between the local and regional scale). Each sample was a mixture of 3 top soil cores at the depth of 0e20 cm and the diameter of 5 cm, which were usually 0.1 m apart. Soil samples were transferred into a plastic bag and stored in a portable 4
C fridge. Samples were immediately transported to the laboratory for processing. 2.2. Measurements of soil geochemical factors We measured soil physico-chemical factors including soil temperature, pH, dissolved organic carbon (DOC), dissolved total nitrogen (DTN), ammonium, nitrate, and volumetric water content (VWC) (Sparks et al., 1996; Ding et al., 2015). Soil temperature at the depth of 10 cm was measured with a mercurial thermometer in a location nearby soil cores during sampling. Soil pH was measured with a pH meter in a soil-water suspension (1:5, fresh soil: deionized water) after shaking for 30 min. Dissolved organic carbon (DOC) and dissolved total nitrogen (DTN; the sum of ammonium, nitrate and dissolved organic nitrogen) were measured with a Multi N/C 2100 analyzer (Analytik Jena, Thuringia, Germany). 2.3. DNA extraction and sequencing experiments Soil DNA was extracted from 1.5 g of soil by freeze-grinding, followed by SDS-based lysis (Liu et al., 2015). DNA was purified twice using 0.5% low melting point agarose gel before phenolechloroformebutanol extraction. DNA quality was assessed with a NanoDrop spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA). Two rounds of PCRs with a common primer pair (gITS7F (5ʹ-GTGARTCATCGARTCTTTG-3ʹ) and ITS4R (5ʹTCCTCCGCTTATTGATATGC-3ʹ) targeting the ITS2 region between 5.8S and 28S rRNA genes were performed (Tedersoo et al., 2018). The ITS2 region has been shown to be less variable than the ITS1 region and overall ITS (Blaalid et al., 2013). Then 250 bp paired-end sequencing of PCR amplicons targeting ITS2 region was performed on a MiSeq platform (Illumina, San Diego, CA, USA) at the Institute for Environmental Genomics, University of Oklahoma, Norman, OK, USA. 2.4. Raw sequence data processing and statistical analyses Raw sequences were separated to samples via sample-specific barcodes and with permission of up to one mismatched base pair. UCLUST was used to classify the operational taxonomic units (OTUs) based on the UNITE database at the 97% identity level after

J. Zhao et al. / Fungal Ecology 42 (2019) 100856

removal of chimera sequences (Abarenkov et al., 2010). We also classified OTUs based on the Warcup database, which generated similar results to UNITE. Singletons were removed before taxonomic annotations. The taxonomic assignment was conducted by RDP classifier with the minimal 50% confidence estimates (Wang et al., 2007). Then we resampled at the depth of 2500 sequences per sample to avoid bias associated with the sequencing depth. All raw sequence analyses were performed in an in-house Galaxy pipeline (http://zhoulab5.rccc.ou.edu:8080/root). The slopes of distance-decay relationships of microbial community were calculated at three spatial scales: the local scale (1 me113 m), the meso scale (3.4 kme39 km) and the regional scale (103 kme668 km). The slope coefficient at each spatial scale was calculated based on equation below: Ln (S) ¼ Ln (a) þ z  Ln (G) where S is the community similarity, G is the geographic distance, a is an intercept parameter and z is the slope coefficient of the distance-decay curve (Martiny et al., 2011). Parametric statistics including Analysis of Variance were performed in the package stats while non-parametric statistics including the Kruskal-Wallis test were performed using the package PMCMR (Pohlert, 2014). To analyze community dissimilarity, non-metric multidimensional scaling (NMDS) ordinations were constructed based on Bray-Curtis distance. To examine the relative importance of spatial-related processes (e.g. dispersal limitation) and environmental selection in shaping fungal community, the variation of fungal community composition data was partitioned between the measured spatial and environmental factors using canonical redundancy analysis (RDA)(Legendre et al., 2009). First, Principal coordinates of neighbor matrices (PCNM) eigenfunctions were computed across the sampling locations. PCNM eigenfunctions represent a ‘spectral decomposition of the spatial relationship among sampling locations’ and can be considered as purely spatial variables in an ordination-based analysis. Forward selection of PCNM variables based on permutation tests (9999 randomizations) was then performed, which identified 4 out of the 10 extracted PCNM variables that significantly explained the spatial structure of the fungal sequence abundance data. Lastly, variation of the community composition data was partitioned between the extracted PCNM spatial variables and soil environmental factors using RDA (Peresneto et al., 2006). To compare the relative contribution of environmental and spatial factors at three spatial scales, we performed scale-specific MRM models. At each spatial scale, we calculated R2 and the partial regression coefficient rho for only those pairwise comparisons within the scale. To test their significance, we performed permutations with the dataset and all pairwise comparisons were included by permutations of the distance matrix of fungal community. After 999 permutations, we re-ran the MRM models to calculate the MRM parameters (Martiny et al., 2011). We compared the original, observed parameters to the distribution of these permuted values in order to calculate statistical significance. Statistical analyses, including PCNM, variation partitioning analysis and MRM were carried out in the R statistical language version 3.5.1 (R-Core-Team, 2018). The vegan, spacemakeR and packfor libraries were used. 2.5. Null model-based community assembly framework To disentangle the relative importance of stochastic and deterministic processes in fungal community assembly, we utilized a null model modified from the Raup-Crick (RC) index considering species abundances (Chase et al., 2011; Stegen et al., 2013).

3

Specifically, the null model was constructed by performing a probability-based randomization process, in which randomly generated fungal OTU composition and abundance were assembled for each sample by randomly sampling from the total fungal OTU pool following four rules as described by Daleo et al. (2018). For all possible pairs of samples, fungal OTU composition of each sample was probabilistically generated for 9999 times. For each iteration, Bray-Curtis dissimilarity index between samples was calculated, and the resulting metric was the proportion of iterations in which the index was smaller than or equal to the actually observed BrayCurtis dissimilarity index between those pair of samples (Chase et al., 2011; Stegen et al., 2013). We standardized the metric to range from 1 to 1 by subtracting 0.5 and multiplying by 2. The resulting metric, referred as RCbray, can be used to calculate the proportion of observed dissimilarities across samples that are higher or lower than those estimated from the null model (9999 randomized matrices), with values near 1 (0.95 to 1) suggesting homogenizing dispersal (that is, mass effect), values near 1 (0.95e1) suggesting environmental selection and other values (0.95 to 0.95) suggesting drift (Bahram et al., 2016). RCbray values near 1 (0.95 to 1) and 1 (0.95e1) are considered as deterministic processes while other values (0.95 to 0.95) are considered as stochastic processes (Alberti et al., 2017). To disentangle relative contributions of these processes at different spatial scales, pairwise pairs of RCbray at each scale were extracted from the overall matrix of RCbray according to distance between samples before calculating the relative contributions. 3. Results 3.1. Environmental factors We measured a number of climatic and soil geochemical factors to reveal the extent of environmental heterogeneity across sampling sites (Supplemental Fig. S2). Soil temperature, soil water content, soil pH, ammonia, nitrate, dissolved organic carbon and dissolved total N were significantly (p < 0.050) different at the regional scale. Samples from Jurong had the highest concentration of dissolved total N and nitrate but the lowest soil water content and soil temperature. Samples from Xiantao had the highest soil pH while those from Hefei had the lowest pH. Samples from Xiantao had the highest dissolved organic carbon concentration while those from Huangshi had the lowest concentration. Samples from Huangshi had the lowest ammonium concentration but both Tongcheng and Xiantao samples had the highest ammonium concentration. 3.2. Overall fungal community composition A total of 11,342 fungal operational taxonomic units (OTUs) were identified from the 198 samples after removing singletons. Among them, only 13.9% of OTUs (1,585) were not assigned to known fungal taxa, indicating that soil fungi in paddy soils were generally well annotated. Although rarefaction curves showed that the total fungal diversity was not fully recovered in some samples (Supplemental Fig. S3), Good's coverage was 97.9% ± 1.94%, suggesting that detected sequences were generally sufficient. The most abundant class is Sordariomycetes across all samples (22.3%, Supplemental Fig. S4) and the most abundant OTU belongs to the genus Mortierella (3.8%, Supplemental Table S2). 3.3. Spatial patterns of fungal communities Fungal communities in the same county were similar to each other based on the NMDS analysis of taxonomic beta-diversity

4

J. Zhao et al. / Fungal Ecology 42 (2019) 100856

metrics (PERMANOVA Fcounty ¼ 13.60, R2county ¼ 0.23, p < 0.001, Supplemental Fig. S5). In accordance, a plot of community similarity versus geographic distance for each pairwise set of samples revealed that fungal community displayed a significant, negative distance-decay relationship (Slope ¼ 0.161, p < 0.001, Fig. 1). The slope of this curve varied significantly among three spatial scales, with the slopes significantly deviating from 0 at local and meso scales (p < 0.010). The distance-decay slope at the local scale was significantly shallower than the overall slope (slope ¼ 0.075; p ¼ 0.0039, Fig. 1) and was steeper at the meso scale than the overall scale (slope ¼ 0.246, p < 0.010, Fig. 1). In contrast, we did not observe a negative distance-decay relationship at the regional scale (slope ¼ 0.073, Fig. 1). Mantel correlogram analysis revealed that fungal communities showed significant spatial structure, exhibiting strong positive autocorrelation at scales <50 km while weak negative autocorrelation at scales >250 km (Supplemental Fig. S6), further confirming the scale-dependent distance-decay relationships. 3.4. Drivers of spatial patterns in soil fungal community To quantify the relative contributions of processes such as environmental selection and ecological drift on fungal community assembly at three scales, we applied null model-based microbial community assembly framework. The proportion of 0.95 0.95 and RCbray < 0.95 accounted for 73.4% when all sample pairs were included, suggesting that deterministic processes contributed more to overall fungal community structure than stochastic processes. When pair values of RCbray at three different spatial scales were compared, the proportion of RCbray > 0.95 and RCbray < 0.95 at the regional scale was higher than that at meso and local scales (82.3% versus 49.2% and 21.6%)(Fig. 2). RCbray values at the regional scale were significantly larger than those at meso and local scales (the Kruskal-Wallis test, chi-squared ¼ 2465.3, df ¼ 2, p < 0.001, Supplemental Fig. S7),

Fig. 1. Distance-decay relationships of fungal communities at three different spatial scales. Lines represent fitted linear regression of the overall distance-decay relationship (slope ¼ 0.161, N ¼ 16344, R2 ¼ 0.17, p < 0.001), the distance-decay relationship at the local scale (slope ¼ 0.075, N ¼ 828, R2 ¼ 0.024, p < 0.001), the meso scale (slope ¼ 0.246, N ¼ 1823, R2 ¼ 0.067, p < 0.001), and the regional scale (slope ¼ 0.073, N ¼ 13693, R2 ¼ 0.003, p < 0.001). Slopes of distance-decay relationships at three scales were tested against zero based on the permutation test. The significant test of slopes between different scales and with overall slope was performed based on the permutation test. The shadowed area represents 95% confidential interval of the fitted linear regression.

suggesting that deterministic processes (selection) were more important in shaping communities at the regional scale than at meso and local scales. Specifically, most communities were assembled by selection (82.3% of all overall pairwise comparisons) at the regional scale, whereas only 17.6% were assembled by ecological drift (Fig. 2). At the meso scale, 50.2% of fungal communities were assembled by drift, while 0.6% of fungal communities were assembled by homogenizing dispersal. At the local scale, 64.0% of fungal communities were assembled by drift, while 14.5% of fungal communities were assembled by homogenizing dispersal. To further verify those results, we performed redundancy analysis-based variation partitioning, which showed that environmental variables alone (9.0%) significantly explained more community variation than that of spatial variables alone (5.4%) at the regional scale (Fig. 3 (a)). When variation partitioning was run at the meso scale, spatial variables had higher explanatory power than environmental factors in all cases (Fig. 3 (b)e(g)), suggesting that spatial related stochastic processes (dispersal and drift) were more important at a smaller scale while environmental factors were more important at the regional scale. We then ran MRM models at three different spatial scales. Only geographic distance showed significant partial coefficients at the local scale (rho ¼ 0.008, p ¼ 0.03, Table S1). In contrast, dissolved total nitrogen and dissolved organic carbon showed significant partial coefficients at meso and regional scales (rho ¼ 0.01, p ¼ 0.027; rho ¼ 0.006, p ¼ 0.034, Supplemental Table S1), suggesting that they were the most important environmental factors in influencing fungal community composition.

4. Discussion Here, we show that there is a scale-dependent distance-decay relationship in fungal communities, similar to observations in bacterial and plant communities (Martiny et al., 2011; Nekola and Mcgill, 2014). Negative distance-decay relationships of soil fungal communities were observed at local and meso scales but not at the regional scale (Fig. 1). The steepest slope was observed at the meso

Fig. 2. Proportions of community pairs assembled by drift, homogenizing dispersal and selection at different spatial scales.

J. Zhao et al. / Fungal Ecology 42 (2019) 100856

5

Fig. 3. The Venn diagrams of variation partitioning analysis showing the effects of spatial and environmental variables (PCNM2~PCNM5) on community variation at regional (a) and meso scales (beg). Values indicate percentages of variation explained by each fraction, including alone and shared explained percentages and unexplained percentages. Forward selection was used to identify the best environmental and spatial variables explaining community variation. The ANOVA permutation test based on redundancy analysis was calculated to reveal variation explained by each set without the effect of the other. CX: Changxing; HF: Hefei; HS: Huangshi; JR: Jurong; TC: Tongcheng; and XT: Xiantao. Significant levels are shown as *** for p < 0.001, ** for p < 0.01 and * for p < 0.05.

scale, differing from what was observed in the landmark study by (Martiny et al., 2011), where the steepest distance-decay relationship was at the regional scale for marsh ammonium-oxidizing bacterial community (Martiny et al., 2011). There are several reasons that may account for this difference. First, frequent cultivation

can result in relatively homogeneous environments within a paddy field, so environmental selection may contribute less to the observed distance-decay pattern. Patchy species aggregation in bulk or rhizosphere soil of rice paddies might lead to divergent composition of fungal communities (Faust and Raes, 2012); such an

6

J. Zhao et al. / Fungal Ecology 42 (2019) 100856

observation remains undetermined here. Second, the priority effect can affect fungal community composition in laboratory microcosms (Kennedy et al., 2009; Peay et al., 2012a). It is possible that a priority effect exists in the natural environment, i.e., locally adapted fungal individuals may have been able to exclude later colonists via competitive exclusion (Kennedy and Bruns, 2010). Consequently, locally adapted fungal individuals increase in abundance and would produce a distanceedecay relationship. However, such local priority effects may diminish over the regional scale (Hanson et al., 2012). Third, the fungal community usually has high taxonomic diversity and thus is prone to under-sampling (Sloan et al., 2006), which might result in under-estimation of distance-decay slopes at the regional scale. The most abundant class Sordariomycetes identified in this study belongs to phylum Ascomycota and is highly abundant in low C:N soil (Yuan et al., 2018), which might result from removal of aboveground biomass during harvesting and application of nitrogen fertilizer to paddy fields. Also, the most abundant OTU belongs to the genus Mortierella, which is also the most abundant in grasslands (Pellissier et al., 2015) and peatland soils (Thormann, 2006). This genus is recognized as containing important saprotrophs for decaying leaves and other organic materials in terrestrial ecosystems, which is in accordance with a crucial role in soil carbon cycling. However, we did not detect high sequence abundance of arbuscular mycorrhizal fungi, which rice is obligately associated with. A possible reason is that our choice of ITS2 primers could potentially underestimate abundance of arbuscular mycorrhizal fungi since ITS2 primers have PCR biases against amplifying Glomeromycota and trovský et al., 2016), which are predominantly Chytridiomycota (Ve comprised of arbuscular mycorrhizal fungi. Within the Glomeromycota phylum, the same microorganism often carries divergent copies of ITS2 sequences (Diniandreote et al., 2016), which do not pair well with our choice of ITS2 primers. Fungi cannot fix carbon and nitrogen. Therefore, fungi rely on external provision of organic matter and nitrogenous compounds for growth. Many fungi prefer warm, sugary, acidic, aerobic, and moist conditions (Walker and White, 2017), which explains significant correlations between fungal communities and soil nutrients (DOC and DTN) (Table S1). Flooding and application of fertilizer in paddy soil impose strong environmental selection. As a consequence, deterministic processes are more important than stochastic processes in shaping fungal community assembly in paddy soils at the regional scale (82.4% versus 17.6%, Fig. 2). Similarly, environmental selection was predominant in shaping fungal community at large scales (Kivlin et al., 2014; Tedersoo et al., 2014; Linde et al., 2018) while stochastic processes such as dispersal and drift play a dominant role in structuring fungal communities at the local scale (Bahram et al., 2016). By disentangling deterministic and stochastic processes at different spatial scales, we showed that community assembly mechanisms for fungi were scale-dependent. At large spatial scales, species habitat associations may be strong, which determines community composition when environmental conditions rapidly change. Meanwhile, probabilistic events such as birth, death and dispersal contribute to existence of sporadically distributed individuals in the unfavorable habitat. When sampling scales decline to include fewer individuals and less habitat heterogeneity, relative contributions of those stochastic events on community composition increase accordingly (Chase et al., 2011). Our results in Fig. 3 (a) were in agreement with the previous study showing that environmental conditions were more important than geographic distance at the regional scale (Martiny et al., 2011). However, these results contradicted with the recent finding that geographic distance, instead of environmental selection, was the major limiting factor in shaping fungal biogeographical pattern in paddy soils (Yuan et al., 2018). This

contradiction probably results from differences in sequencing depths, since a low sequence depth can lead to over-estimation of contribution from dispersal limitation (Zhou et al., 2013). In addition, PCR amplification introduces bias to species abundance, which in turns affects community assembly process analyses. The use of RCbray, which is absent in Yuan et al.’s study, also introduces bias since the Bray-Curtis index is not metric. There are two competing hypotheses for microbial biogeographic studies (Wu et al., 2017). The size-dispersal hypothesis argues that larger organisms are more likely to be affected by dispersal limitation, because dispersal is more difficult for larger organisms. The size-plasticity hypothesis argues that larger organisms are more affected by environment selection because they have less plasticity in metabolic abilities than smaller organisms, and therefore, exist in a narrower habitat range. Similar to previous observations in fungal communities (Kivlin et al., 2014; Tedersoo et al., 2014), our results suggest that dispersal limitation does not strongly affect the distribution of fungal taxa even at relatively large spatial scales. Interestingly, there is evidence that AM fungi at large scale are not dispersal limited (Davison et al., 2015) whereas other studies found strong evidences for dispersal limitation by investigating fungal spores (Peay et al., 2012b; Peay and Bruns, 2014). The lack of observed dispersal limitation might be a result of the overreliance on OTU-based studies, or the selection of unsuitable primers targeting ITS region. Nevertheless, our model analysis provides evidence that environmental selection increases with geographic distance, as the proportion of sample pairs assembled by selection is smaller at the local scale (21.6.0%) than those at the meso scale (49.2%) and at the regional scale (82.3%). In contrast, sample pairs assembled by ecological drift are larger at the local scale (63.9%) compared with those at the meso scale (50.2%) and at the regional scale (17.6%), which offers strong support for higher drift when small-scale sampling includes less environmental heterogeneity (Chase, 2014). However, the null model-based framework may be subjected to statistical artifacts. Communities are considered to be assembled by both environmental selection and dispersal limitation if RCbray lies between 0.95 and 1, thus the effects of dispersal limitation is not considered in our framework according to Bahram et al. (2016). Also, it is likely that biotic selection (e.g. fungi-fungi interaction or fungi-plant interaction) falls into those criteria and thus the effect of biotic selection on community spatial patterns may be misidentified as a result of dispersal limitation and environmental selection (Chase et al., 2011). Therefore, introducing RCbray in community assembly processes has substantial influence on quantitation of contributions of selection and dispersal limitation. Lastly, communities assembled by ecological drift could also incorrectly include selection exerted by unmeasured environmental factors, evolution and historical contingencies (Zhou and Ning, 2017). In summary, both deterministic and stochastic processes contribute to fungal community structure and assembly in paddy fields. Deterministic processes of environmental selection play a more important role than dispersal-related processes and stochastic processes in structuring fungal community at the regional scale. By quantifying the effects of an individual environmental factor, we found that soil organic carbon was the most important factor associated with environmental selection. However, spatial related processes and stochastic events such as drift are more important for local communities, suggesting that spatial scales have to be taken into account when community assembly processes are examined (Mcgill, 2010; Chase, 2014). Taken together, our study provides valuable insights into how community assembly processes shape biogeographic patterns of fungal communities in paddy fields.

J. Zhao et al. / Fungal Ecology 42 (2019) 100856

Author's contribution YY and JZhou designed the experiments. JZhao, MW and QG performed bioinformatic analysis and further analysis. QG took soil samples and extracted DNA from soil samples. JZhao, YY and JZhou wrote the manuscript. YL, BS and HC helped write and improve the manuscript. All authors have contributed to the final manuscript. Data and code availability All sequencing datasets have been deposited in the NCBI Sequence Read Archive (SRA) under accession number SRP135966, BioProject ID PRJNA438873 with biosamples accession SAMN10259995 e SAMN10260188. Code for calculating RCbray is in https://github.com/jianshu93/Fungi_BioDiv_CommAssembly. Acknowledgements We want to thank two anonymous reviewers for their valuable comments on the manuscript. We also thank Kristen N. Wyckoff from Department of Civil and Environmental Engineering, University of Tennessee, Knoxville, TN, USA for help with the language. This research was supported by grants to Jizhong Zhou from National Science Foundation of China (41430856), and Yunfeng Yang from the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB15010102), National Science Foundation of China (41471202). Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.funeco.2019.07.010. References Abarenkov, K., Henrik Nilsson, R., Larsson, K.H., Alexander, I.J., Eberhardt, U., Erland, S., Høiland, K., Kjøller, R., Larsson, E., Pennanen, T., 2010. The UNITE database for molecular identification of fungierecent updates and future perspectives. New Phytol. 186, 281e285. Alberti, J., Bakker, E.S., Van Klink, R., Olff, H., Smit, C., 2017. Herbivore exclusion promotes a more stochastic plant community assembly in a natural grassland. Ecology 98, 961e970. Bahram, M., Kohout, P., Anslan, S., Harend, H., Abarenkov, K., Tedersoo, L., 2016. Stochastic distribution of small soil eukaryotes resulting from high dispersal and drift in a local environment. ISME J. 10, 885e896. Bahram, M., Tedersoo, L., 2013. The distance decay of similarity in communities of ectomycorrhizal fungi in different ecosystems and scales. J. Ecol. 101, 1335e1344. Beck, S., Powell, J.R., Drigo, B., Cairney, J.W., Anderson, I.C., 2015. The role of stochasticity differs in the assembly of soil-and root-associated fungal communities. Soil Biol. Biochem. 80, 18e25. Blaalid, R., Kumar, S., Nilsson, R.H., Abarenkov, K., Kirk, P., Kauserud, H., 2013. ITS1 versus ITS2 as DNA metabarcodes for fungi. Mol. Ecol. Resour. 13, 218e224. Bruns, T.D., Taylor, J.W., 2016. Comment on “Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism”. Science 351, 826826. Chase, J.M., 2014. Spatial scale resolves the niche versus neutral theory debate. J. Veg. Sci. 25, 319e322. Chase, J.M., Kraft, N.J.B., Smith, K.G., Vellend, M., Inouye, B.D., 2011. Using null models to disentangle variation in community dissimilarity from variation in adiversity. Ecosphere 2, 1e11. Cheng, L., Booker, F.L., Tu, C., Burkey, K.O., Zhou, L., Shew, H.D., Rufty, T.W., Hu, S., 2012. Arbuscular mycorrhizal fungi increase organic carbon decomposition under elevated CO2. Science 337, 1084e1087. Clemmensen, K.E., Bahr, A., Ovaskainen, O., Dahlberg, A., Ekblad, A., Wallander, H., Stenlid, J., Finlay, R.D., Wardle, D.A., Lindahl, B.D., 2013. Roots and associated fungi drive long-term carbon sequestration in boreal forest. Science 339, 1615e1618. Daleo, P., Alberti, J., Jumpponen, A., Veach, A., Ialonardi, F., Iribarne, O., Silliman, B., 2018. Nitrogen enrichment suppresses other environmental drivers and homogenizes salt marsh leaf microbiome. Ecology 99, 1411e1418. € Davison, J., Moora, M., Opik, M., Adholeya, A., Ainsaar, L., B^ a, A., Burla, S., Diedhiou, A.G., Hiiesalu, I., Jairus, T., Johnson, N.C., Kane, A., Koorem, K., €rtel, M., Reier, Ü., Saks, Ü., Singh, R., Vasar, M., Kochar, M., Ndiaye, C., Pa

7

Zobel, M., 2015. Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism. Science 349, 970e973. De, W.R., Bouvier, T., 2010. 'Everything is everywhere, but, the environment selects'; what did Baas Becking and Beijerinck really say? Environ. Microbiol. 8, 755e758. Ding, J., Zhang, Y., Wang, M., Sun, X., Cong, J., Deng, Y., Lu, H., Yuan, T., Van Nostrand, J.D., LI, D., Zhou, J., Yang, Y., 2015. Soil organic matter quantity and quality shape microbial community compositions of subtropical broadleaved forests. Mol. Ecol. 24, 5175e5185. Diniandreote, F., Pylro, V.S., Baldrian, P., Elsas, J.D.V., Salles, J.F., 2016. Ecological succession reveals potential signatures of marine-terrestrial transition in salt marsh fungal communities. ISME J. 10, 1984e1997. Dolan, J.R., 2006. Microbial biogeography? J. Biogeogr. 33, 199e200. Dumbrell, A.J., Nelson, M., Helgason, T., Dytham, C., Fitter, A.H., 2010. Relative roles of niche and neutral processes in structuring a soil microbial community. ISME J. 4, 337e345. Faust, K., Raes, J., 2012. Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538e550. Glassman, S.I., Wang, I.J., Bruns, T.D., 2017. Environmental filtering by pH and soil nutrients drives community assembly in fungi at fine spatial scales. Mol. Ecol. 26, 6960e6973. Green, J.L., Holmes, A.J., Westoby, M., Oliver, I., Briscoe, D., Dangerfield, M., Gillings, M., Beattie, A.J., 2004. Spatial scaling of microbial eukaryote diversity. Nature 432, 747e750. Hanson, C.A., Fuhrman, J.A., Horner-Devine, M.C., Martiny, J.B.H., 2012. Beyond biogeographic patterns: processes shaping the microbial landscape. Nat. Rev. Microbiol. 10, 497e506. Hawksworth, D.L., Lücking, R., 2017. Fungal diversity revisited: 2.2 to 3.8 million species. Microbiol. Spectr. 5. Kennedy, P.G., Bruns, T.D., 2010. Priority effects determine the outcome of ectomycorrhizal competition between two Rhizopogon species colonizing Pinus muricata seedlings. New Phytol. 166, 631e638. Kennedy, P.G., Peay, K.G., Bruns, T.D., 2009. Root tip competition among ectomycorrhizal fungi: are priority effects a rule or an exception? Ecology 90, 2098e2107. Kivlin, S.N., Winston, G.C., Goulden, M.L., Treseder, K.K., 2014. Environmental filtering affects soil fungal community composition more than dispersal limitation at regional scales. Fungal Ecol. 12, 14e25. Legendre, P., Mi, X., Ren, H., Ma, K., Yu, M., Sun, I.F., He, F., 2009. Partitioning beta diversity in a subtropical broad-leaved forest of China. Ecology 90, 663e674. Leibold, M.A., Norberg, J., 2004. Biodiversity in metacommunities: plankton as complex adaptive systems? Limnol. Oceanogr. 49, 1278e1289. Linde, S.V.D., Suz, L.M., Orme, C.D.L., Cox, F., Andreae, H., Asi, E., Atkinson, B., Benham, S., Carroll, C., Cools, N., 2018. Environment and host as large-scale controls of ectomycorrhizal fungi. Nature 558, 243e248. Liu, S., Wang, F., Xue, K., Sun, B., Zhang, Y., He, Z., Van Nostrand, J.D., Zhou, J., Yang, Y., 2015. The interactive effects of soil transplant into colder regions and cropping on soil microbiology and biogeochemistry. Environ. Microbiol. 17, 566e576. Locey, K.J., Lennon, J.T., 2016. Scaling laws predict global microbial diversity. Proc. Natl. Acad. Sci. U. S. A 113, 5970e5975. Lumini, E., Vallino, M., Alguacil, M.M., Romani, M., Bianciotto, V., 2011. Different farming and water regimes in Italian rice fields affect arbuscular mycorrhizal fungal soil communities. Ecol. Appl. 21, 1696e1707. Martiny, J.B., Eisen, J.A., Penn, K., Allison, S.D., Horner-Devine, M.C., 2011. Drivers of bacterial beta-diversity depend on spatial scale. Proc. Natl. Acad. Sci. U. S. A 108, 7850-7584. Mcgill, B.J., 2010. Matters of scale. Science 328, 575e576. Mohamed, D.J., Martiny, J.B., 2010. Patterns of fungal diversity and composition along a salinity gradient. ISME J. 5, 379e388. Nekola, J.C., Mcgill, B.J., 2014. Scale dependency in the functional form of the distance decay relationship. Ecography 37, 309e320. Peay, K.G., Belisle, M., Fukami, T., 2012a. Phylogenetic relatedness predicts priority effects in nectar yeast communities. Proc. Biol. Sci. 279, 5066-5066. Peay, K.G., Bruns, T.D., 2014. Spore dispersal of basidiomycete fungi at the landscape scale is driven by stochastic and deterministic processes and generates variability in plantefungal interactions. New Phytol. 204, 180e191. Peay, K.G., Matteo, G., Bruns, T.D., 2010. Evidence of dispersal limitation in soil microorganisms: isolation reduces species richness on mycorrhizal tree islands. Ecology 91, 3631e3640. Peay, K.G., Schubert, M.G., Nguyen, N.H., Bruns, T.D., 2012b. Measuring ectomycorrhizal fungal dispersal: macroecological patterns driven by microscopic propagules. Mol. Ecol. 21, 4122e4136. Pellissier, L., Niculitahirzel, H., Dubuis, A., Pagni, M., Guex, N., Ndiribe, C., Salamin, N., Xenarios, I., Goudet, J., Sanders, I.R., 2015. Soil fungal communities of grasslands are environmentally structured at a regional scale in the Alps. Mol. Ecol. 23, 4274e4290. Peresneto, P.R., Legendre, P., Dray, S., Borcard, D., 2006. Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology 87, 2614e2625. Pickles, B.J., Genney, D.R., Potts, J.M., Lennon, J.J., Anderson, I.C., Alexander, I.J., 2010. Spatial and temporal ecology of Scots pine ectomycorrhizas. New Phytol. 186, 755e768. Pohlert, T., 2014. The pairwise multiple comparison of mean ranks package (PMCMR). R package. http://CRAN.R-project.org/package¼PMCMR. R-Core-Team, 2018. R: A Language and Environment for Statistical Computing. R

8

J. Zhao et al. / Fungal Ecology 42 (2019) 100856

Foundation for Statistical Computing, Vienna, Austria. https://www.R-project. org/. Roesch, L.F., RR, F., A, R., G, C., AK, H., AD, K., SH, D., FA, C., WG, F., EW, T., 2007. Pyrosequencing enumerates and contrasts soil microbial diversity. ISME J. 1, 283e290. Sloan, W., Lunn, M., Woodcock, S., Head, I., Nee, S., Curtis, T., 2006. Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environ. Microbiol. 8, 732e740. Sparks, D.L., Fendorf, S.E., Toner IV, C.V., Carski, T., 1996. Kinetic methods and measurements. Part Methods Soil Anal. 3, 1275e1307. Stegen, J.C., Freestone, A.L., Crist, T.O., Anderson, M.J., Chase, J.M., Comita, L.S., Cornell, H.V., Davies, K.F., Harrison, S.P., Hurlbert, A.H., 2013. Stochastic and deterministic drivers of spatial and temporal turnover in breeding bird communities. Glob. Ecol. Biogeogr. 22, 202e212. Talbot, J.M., Bruns, T.D., Taylor, J.W., Smith, D.P., Branco, S., Glassman, S.I., Erlandson, S., Vilgalys, R., Liao, H.L., Smith, M.E., 2014. Endemism and functional convergence across the North American soil mycobiome. Proc. Natl. Acad. Sci. U. S. A 111, 6341e6346. ~ lme, S., Ko ~ljalg, U., Yorou, N.S., Wijesundera, R., Tedersoo, L., Bahram, M., Po Ruiz, L.V., Vascopalacios, A.M., Thu, P.Q., Suija, A., 2014. Global diversity and geography of soil fungi. Science 346, 1256688-1-1256688-10. Tedersoo, L., Tooming-Klunderud, A., Anslan, S., 2018. PacBio metabarcoding of Fungi and other eukaryotes: errors, biases and perspectives. New Phytol. 217, 1370e1385. Thormann, M.N., 2006. Diversity and function of fungi in peatlands: a carbon cycling perspective. Can. J. Soil Sci. 86, 281e293. Tisthammer, K.H., Cobian, G.M., Amend, A.S., 2016. Global biogeography of marine fungi is shaped by the environment. Fungal Ecol. 19, 39e46. V
alyi, K., Mardhiah, U., Rillig, M.C., Hempel, S., 2016. Community assembly and coexistence in communities of arbuscular mycorrhizal fungi. ISME J. 10,

2341e2351. Van Der Heijden, M.G.A., Streitwolf-Engel, R., Riedl, R., Siegrist, S., Neudecker, A., ineichen, K., Boller, T., Wiemken, A., Sanders, I.R., 2006. The mycorrhizal contribution to plant productivity, plant nutrition and soil structure in experimental grassland. New Phytol. 172, 739e752. Verbruggen, E., Sheldrake, M., Bainard, L.D., Chen, B., Ceulemans, T., De Gruyter, J., Van Geel, M., 2017. Mycorrhizal fungi show regular community compositions in natural ecosystems. ISME J. 12, 380e385.  L., Zelenka, T., Baldrian, P., 2016. The rpb2 gene reptrovský, T., Kolarík, M., Z, Ve resents a viable alternative molecular marker for the analysis of environmental fungal communities. Mol. Ecol. Resour. 16, 388e401. Walker, G.M., White, N.A., 2017. Introduction to fungal physiology. Fungi 1e35. Wang, Q., Garrity, G.M., Tiedje, J.M., Cole, J.R., 2007. Naïve bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261e5267. Widder, S., Allen, R.J., Pfeiffer, T., Curtis, T.P., Wiuf, C., Sloan, W.T., Cordero, O.X., Brown, S.P., Momeni, B., Shou, W., 2016. Challenges in microbial ecology: building predictive understanding of community function and dynamics. ISME J. 10, 2557e2568. Wu, W., Lu, H.-P., Sastri, A., Yeh, Y.-C., Gong, G.-C., Chou, W.-C., Hsieh, C.-H., 2017. Contrasting the relative importance of species sorting and dispersal limitation in shaping marine bacterial versus protist communities. ISME J. 12, 485e494. Yuan, C., Zhang, L., HU, H., Wang, J., Shen, J., He, J., 2018. The biogeography of fungal communities in paddy soils is mainly driven by geographic distance. J. Soils Sediments 18, 1795e1805. Zhou, J., Jiang, Y.-H., Deng, Y., Shi, Z., Zhou, B.Y., Xue, K., Wu, L., He, Z., Yang, Y., 2013. Random sampling process leads to overestimation of b-diversity of microbial communities. mBio 4 e00324-13. Zhou, J., Ning, D., 2017. Stochastic community assembly: does it matter in microbial ecology? Microbiol. Mol. Biol. Rev. 81, e00002e17.