Climate mediates continental scale patterns of stream microbial functional diversity
Microbiome volume 8, Article number: 92 (2020)
Understanding the large-scale patterns of microbial functional diversity is essential for anticipating climate change impacts on ecosystems worldwide. However, studies of functional biogeography remain scarce for microorganisms, especially in freshwater ecosystems. Here we study 15,289 functional genes of stream biofilm microbes along three elevational gradients in Norway, Spain and China.
We find that alpha diversity declines towards high elevations and assemblage composition shows increasing turnover with greater elevational distances. These elevational patterns are highly consistent across mountains, kingdoms and functional categories and exhibit the strongest trends in China due to its largest environmental gradients. Across mountains, functional gene assemblages differ in alpha diversity and composition between the mountains in Europe and Asia. Climate, such as mean temperature of the warmest quarter or mean precipitation of the coldest quarter, is the best predictor of alpha diversity and assemblage composition at both mountain and continental scales, with local non-climatic predictors gaining more importance at mountain scale. Under future climate, we project substantial variations in alpha diversity and assemblage composition across the Eurasian river network, primarily occurring in northern and central regions, respectively.
We conclude that climate controls microbial functional gene diversity in streams at large spatial scales; therefore, the underlying ecosystem processes are highly sensitive to climate variations, especially at high latitudes. This biogeographical framework for microbial functional diversity serves as a baseline to anticipate ecosystem responses and biogeochemical feedback to ongoing climate change.
Functional gene alpha diversity monotonically declines towards high elevations and assemblage composition shows increasing turnover with greater elevational distances, these patterns being consistent across mountains, kingdoms and functional gene families.
Climate primarily explains the alpha diversity and compositional changes at mountain and continental scales, with local non-climatic predictors gaining more importance at mountain scale.
Under future climate scenarios, alpha diversity and assemblage composition show substantial changes across Eurasia, especially at mid- and high-latitudes.
Across terrestrial and aquatic habitats, microorganisms mediate many ecosystem processes [1, 2] crucially linked to climate through complex interactions and feedbacks . Gaining insight into the large-scale patterns of microbial diversity is essential for predicting climate change impacts on ecosystems worldwide. The assessment of biodiversity responses to climate change by functional approaches can improve the quantitative and predictive power of ecological research . However, functional biogeography remains largely unexplored for microbes , especially in freshwater ecosystems. The increasing availability of low-cost molecular methods, such as high-throughput sequencing and gene arrays [6, 7], allows linking microbial assemblages with ecosystem processes and enables the application of novel functional perspectives to classic microbial biogeography questions. Among these questions, elucidating the drivers that underlie the biogeographical patterns constitutes a key task.
Climate is among the main drivers shaping large-scale gradients in microbial diversity. For instance, temperature has been shown to drive the latitudinal gradient in diversity of planktonic marine bacteria  or the continental scale gradients in diversity and assemblage composition of functional genes involved in denitrification and nitrogen fixation across forest soils [9, 10]. Moreover, microbial functional gene diversity and composition can respond rapidly under warming and cooling scenarios, as demonstrated by experiments on soils [11, 12]. These findings are framed in the metabolic theory of ecology , and thus higher mutation and speciation rates are expected towards warmer areas in accordance with the kinetics of biological processes. Furthermore, precipitation can also drive large-scale gradients in microbial diversity by constraining above ground vegetation and soil moisture across bioclimatic zones [14, 15]. Nevertheless, the role of climate in shaping microbial functional diversity patterns remains scarcely addressed at large scales for a wide range of functional genes. This is especially true for freshwaters, the most threatened ecosystems on Earth , where studies over small spatial scales have demonstrated temperature to drive microbial diversity and composition in lakes  and streams , as well as to be positively related to functional diversity of biofilm bacteria in streams . Meanwhile, the relationship between precipitation and freshwater microbial assemblages has been understudied although, for instance, rainfall seasonality is linked to functional diversity of stream biofilms . At large geographical scales, precipitation may impact lake microbial assemblages as climate change-related phenomena are altering the nitrogen and phosphorus contents of water rainfall .
Here, we provide evidence from three mountainsides across Eurasia (Additional file 1: Figure S1) for climate to mediate large-scale gradients in stream microbial functional gene diversity. Firstly, we reveal the biogeographical patterns of functional diversity within and across mountainsides. We secondly determine the role of climatic and local non-climatic predictors, i.e. in situ abiotic conditions, in shaping such patterns. Finally, based on the anticipated link between functional diversity and climate [9, 10, 22, 23], we perform predictive models  to project the changes in functional diversity under future climate scenarios  across the Eurasian river network. Towards these aims, we focus on microbial functional genes from kingdoms archaea, bacteria and fungi which are involved in the cycling of carbon (C), nitrogen (N), phosphorus (P) and sulphur (S) as well as in stress-related processes (St). The metabolic potential of these functional genes was assessed by the DNA-based gene array GeoChip 4.0 , which contains approximately 82,000 probes and covers 141,995 gene sequences from 410 functional gene families. In total, we detected 15,289 functional genes linked to 88 functional gene families (Additional file 1: Table S1) included in these five functional categories (C, N, P, S and St). Microbial functional genes are ideal candidates for biodiversity projections under future climates as they meet two essential requirements of spatial distribution modelling , such as a high equilibrium with the environment because of their rapid response to changing conditions [11, 12], and the relatively high dispersal capacity reported for microbes . Moreover, literature indicates the utility in using elevational gradients for anticipating biodiversity responses to climate change over large geographical extents . Our study constitutes a remarkable contribution to the emerging field of functional biogeography , as most previous studies examined the large-scale patterns of microbial diversity from taxonomic or phylogenetic perspectives [8, 29].
Results and discussion
Within mountains, functional gene alpha diversity presented a significant (linear model, LM; P < 0.05) monotonic decline towards high elevations (Fig. 1a, Additional file 1: Table S2) and assemblage composition showed a significant (P < 0.05) increase in turnover with greater elevational distances (Fig. 1b, Additional file 1: Table S2). These patterns were highly consistent across mountains and kingdoms and, interestingly, concur with recent reports showing that tree functional diversity declines with elevation in forests worldwide . The uniformity in our alpha diversity gradients contrasts with the variety of elevational trends reported for microorganisms on taxonomic or phylogenetic approaches in both soil  and freshwater  environments. For instance, species richness of stream bacteria showed hump-shaped, monotonically increasing and U-shaped patterns on the same mountainsides in Norway, Spain and China, respectively . Such disparities between microbial taxonomic and functional diversities are an interesting phenomenon likely caused by functional redundancy which seems an inevitable outcome in open microbial systems where diversity is not limited by low immigration rates . Meanwhile, the consistency in our compositional turnover gradients agrees with the general spatial distance-decay patterns documented for microbial taxonomic composition  and with the elevational distance-decay patterns observed for microbial functional gene assemblages in soils . When scaling down to functional gene families, most of them also showed a significant (LM; P < 0.05) decline in alpha diversity with elevation (Fig. 1c, Additional file 1: Figure S2a), and this pattern was congruent across the three mountains for 58% of all families (51 out of 88). Taken together, these findings highlight a general response of stream microbial functional diversity to elevation across mountains, kingdoms and functional gene families.
Despite the remarkable congruency in our elevational patterns, the mountain in China showed the strongest trends for both functional gene alpha diversity and assemblage composition. Such a finding was consistent across kingdoms and functional gene families. For instance, regarding gene families, China showed the highest proportion of significant (P < 0.05) and the strongest (Bonferroni-pairwise t-test; P < 0.05) elevational gradients in alpha diversity and composition (Fig. 1c, d and Additional file 1: Figures S2a and S2b). These cross-mountain variations in the elevational trends might be related to distinct functional gene regional pools and latitudes  as the tropospheric temperature lapse rate decreases from the equator to the poles . Beyond the differences in slope steepness of our diversity gradients, the high congruency in the elevational trends both across kingdoms and gene families constitutes an interesting finding. For instance, the cross-family consistency in the elevational trends suggests that the microbial-mediated stream processes related to biogeochemical cycling and stress are constrained by similar driving forces.
Across mountains, functional gene assemblages exhibited continental scale patterns and regional differences on the overall gene pool were detected in alpha diversity (analysis of variance, ANOVA; P = 0.001; Additional file 1: Table S3) and assemblage composition (permutational analysis of variance; R2 = 0.387, P = 0.001; Table S3). The streams in China showed the highest mean alpha diversity (Bonferroni-corrected pairwise t-test; P < 0.05; Additional file 1: Table S3) as well as the most different composition (Additional file 1: Table S3 and Fig. S3a) for the overall functional gene pool, such a result being consistent across the three kingdoms. These findings expand the large-scale gradients in functional diversity documented for both macro-  and microorganisms [6, 9, 10] from the terrestrial to the freshwater realm and, namely for microbial functional gene diversity, from genes linked to denitrification  and nitrogen fixation  to genes involved in other biogeochemical cycles (e.g. C, P and S) and St-related processes.
To elucidate the drivers shaping our functional diversity gradients, we examined a set of 19 climatic predictors comprising the average conditions for years 1960–1990 and 11 local non-climatic predictors reported as important drivers of microbial diversity [37, 38] (see “Methods” for details). Regarding alpha diversity, we found that mean temperature of the warmest quarter (TWQ) was the most important predictor at both continental and mountain scales, revealed by hierarchical partitioning  (HP; Fig. 2a), LMs (Additional file 1: Table S4) and random forest  (RF; Additional file 1: Figure S4a). Such a finding is in line with the biological importance of the growing season conditions in mountains . Importantly, the alpha diversity response to the TWQ at continental scale was consistent across mountains; that is, such a response was not affected by site location (TWQ, P < 0.001; latitude, P = 0.109; TWQ × latitude, P = 0.573). Such a positive relationship between temperature and stream microbial diversity have been also reported using taxonomic approaches, e.g. in a recent study on streams from the Rocky Mountains . In accordance with the metabolic theory of ecology , this large-scale temperature dependence of functional diversity could be associated with increasing mutation and speciation rates towards warmer areas [8,9,10]. This result may be also compatible with the energy hypothesis, which proposes that higher diversity coexists with higher available energy . The consistent relationships between elevation and alpha diversity for most of the functional gene families further suggest a general temperature dependence for multiple stream processes. This finding was especially noticeable for families associated with stress processes (Additional file 1: Figure S5a) and highlights the need for future research on temperature dependence to consider a wider range of ecosystem processes in addition to mere C-cycling .
Regarding assemblage composition, temperature-related variables constituted the best predictors at mountain scale (Fig. 2c; Additional file 1: Table S4). This result is in line with previous studies on terrestrial habitats, which show changes in functional gene composition across temperature gradients [9, 10]. Across mountains, however, mean precipitation of the coldest quarter (PCQ) was consistently retained as the best continental scale predictor in HP (Fig. 2c), LMs (Additional file 1: Table S4), generalized dissimilarity models  (GDMs; Additional file 1: Figure S6a) and RF (Additional file 1: Figure S4b). The PCQ, i.e. winter precipitation, mediates the hydrology in high mountains through the accumulation of snow, consequently affecting local factors shaping the spatial patterns of stream microbial assemblages [31, 45, 46]. However, the high predictive power of the PCQ only at continental scale suggests a link between this factor and functional gene composition via hydrological mechanisms and agrees with the importance of precipitation when the spatial scale spans multiple biomes . Such a finding was further supported by the effect of site locations on assemblage composition (LM; PCQ, P < 0.001; latitude, P < 0.001, PCQ × latitude, P = 0.183). These results agree with previous observational and experimental studies in soils, which show precipitation to be a large-scale driver of microbial assemblage composition by constraining vegetation types across biomes [14, 15].
Compared to climate, local non-climatic predictors explained a lower amount of geographical variation in alpha diversity and assemblage composition at both mountain and continental scales, with local non-climatic predictors gaining more importance at mountain scale as revealed by various statistical methods, such as variation partitioning (Fig. 2b, d). These results contrast with the importance of local non-climatic factors, such as nutrients or pH, in shaping freshwater microbial taxonomic diversity . It suggests different driving forces for microbial taxonomic and functional gene diversities, which is in line with the decoupling observed between taxonomy and function in microbial systems . For example, no local predictor was consistently meaningful for alpha diversity at neither continental nor mountain scales under the HP and LM analyses (Fig. 2a; Additional file 1: Table S4). Regarding composition, only total nitrogen was identified as a significant driver in China. Interestingly, nitrogen availability has been shown to elicit changes in microbial assemblage composition, but not in alpha diversity, under elevated nutrient inputs in a global experiment on grasslands . Other potentially important predictors of freshwater microbial diversity, such as those informing about carbon availability , may also play a role. For instance, dissolved organic matter poorly predicts taxonomic elevational diversity and the biogeography of freshwater bacteria [31, 50], but is associated with bacterial trait structure over large spatial scales as a response to the terrestrial influence gradient in a continuum of small streams to large lakes . We therefore encourage further studies to assess the role of carbon availability in explaining functional gene diversity of stream microbes.
Collinearity between climatic and local non-climatic predictors is frequent for elevational gradients. By including three different mountainsides, we could test the role of all local non-climatic predictors across different models, either based on individual mountains or the overall dataset, and found climatic predictors to primarily mediate both alpha diversity and assemblage composition (see RF results on all predictors included in the full models at Additional file 1: Figure S7). For instance, among all local non-climatic predictors, only conductivity was revealed as a candidate to importantly drive alpha diversity when the entire dataset was analysed (see RF results prior to the removal of highly correlated predictors at Additional file 1: Figure S8). However, when each mountain was analysed separately, alpha diversity showed a significant and cross-mountain consistent response only to TWQ (LM; P < 0.05 for all mountains; Additional file 1: Figure S9), but not to conductivity (LM; P = 0.732 for Norway; Additional file 1: Figure S9) or other relevant local non-climatic predictors (Additional file 1: Fig. S9). Furthermore, we assessed all possible models including either the best climatic predictors or conductivity and their respective interactions with latitude to control for cross-mountain consistency and found that best models based on Akaike weights  were those containing the climatic predictors and the interaction with latitude (weight = 0.725 for alpha diversity; weight = 1.000 for assemblage composition; Additional file 1: Table S5). Despite these evidences, we encourage further complementary manipulative experiments  to elucidate how functional gene diversity response to those local non-climatic predictors that covary or interact with climate, such as conductivity or nutrient loads . Collectively, our findings show climate to mediate overall functional gene diversity in freshwater ecosystems as large spatial scale, suggesting that ongoing climate change can affect stream processes worldwide.
Given that elevational gradients offer excellent possibilities to conduct natural experiments for assessing climate change impacts , we employed LMs and GDMs to project the trends of alpha diversity and assemblage composition, respectively, under future climate scenarios (period 2061–2080) across the Eurasian river network. We selected three scenarios from the Fifth Assessment Report of the Intergovernmental Panel on Climate Change . They are based on the representative concentration pathways (RCP) 2.6, 4.5 and 8.5 and symbolize the most “optimistic”, one intermediate and the most “pessimistic” expectations of greenhouse emission rates, respectively, by the end of this century . For modelling, we combined together the TWQ and PCQ as predictors, i.e. the variables that best explained the continental scale gradients in overall alpha diversity and assemblage composition, respectively. These predictors jointly explained a large amount of variance for both overall alpha diversity (LMs; R2 = 0.543; P < 0.001) and composition (GDMs; D2 = 65.6%; P < 0.001; Additional file 1: Figure S6b). Our projections show a general increase in alpha diversity (Fig. 3a, Additional file 1: Figure S10) and elevated compositional turnover rates (Fig. 3b, Additional file 1: Figure S11) as a response to future climates. Compared to the baseline, the median increase in alpha diversity and the median turnover rates were predicted to be 11.0% and 29.6%, respectively, under the moderate RCP 4.5 scenario (Fig. 3c, e). Across kingdoms, i.e. archaea, bacteria and fungi, the changes in alpha diversity and composition followed similar trends, but significant differences in their magnitude (Additional file 1: Figure S12). Fungal functional genes showed the highest projected increase in alpha diversity, while bacteria accounted for the greatest changes in composition (Additional file 1: Figure S12). Such a finding indicates cross-kingdom variations in the magnitude of functional diversity response to climate change, which constitutes a significant step beyond previous reports on taxonomic approaches which show distinct biogeographical patterns and responses to nutrient availability for archaea, bacteria and fungi [50, 54]. Interestingly, the projected changes showed highly congruent trends across functional gene families (Fig. 3d, f). These results are consistent with the rapid response documented for macro- and microbial genes under climate change experiments [12, 55] and suggest that temperature and precipitation changes may strongly influence functional gene assemblages mediating multiple ecosystem processes, likely producing biogeochemical feedbacks to the climate [56, 57].
Across Eurasia, the largest increase in alpha diversity will occur in regions North of 60° N, while assemblage composition will mainly shift between 45° and 75° N. These projections show for the first time the highest sensitivity of northern regions to functional diversity changes affecting multiple stream processes, such as biogeochemical cycling. Despite no previous studies have anticipated the future trends of microbial functional diversity, neither in terrestrial nor freshwater ecosystems, there are interesting reports about plant traits and certain ecosystem functions [58, 59]. For instance, temperature sensitivity of soil respiration is predicted to mainly increase with warming at high latitudes by affecting microbial assemblage compositions .
The predicted changes can be further visualized through the lens of key biogeochemical processes mediated by gene families that showed a consistent alpha diversity response to elevation. Such a response was especially evident for gene families associated with the decomposition of recalcitrant C, such as chitin, aromatics and lignin (Additional file 1: Figure S13). This finding suggests an increase in CO2 emissions through the stimulation of old recalcitrant C upstream and northward, likely producing a positive feedback to global warming . In addition, all gene families linked to denitrification exhibited a similar alpha diversity response to elevation (Additional file 1: Figure S13). This finding suggests accelerated rates of nitrification with global warming that could increase anthropogenic emissions of N2O, a potent greenhouse gas contributing to climate change and stratospheric ozone destruction . Other gene families associated with essential processes, such as nitrogen fixation, assimilatory nitrogen reduction and nitrification, also showed a consistent decline in alpha diversity upstream, therefore suggesting a high sensitivity of these processes towards climate change.
Although we here show a big picture on the large-scale functional implications of microbial diversity changes, two major caveats could be acknowledged regarding the spatial modelling framework. First, our projections comprise considerable extrapolated areas due to the relatively limited climate gradient covered by the studied mountains. Further studies are encouraged to broaden the climate envelope by including additional elevational gradients. Second, our modelling informs on the variance in functional gene alpha diversity and assemblage composition exclusively linked to climatic predictors while it does not consider local- and landscape-level predictors. The reasons behind this strategy are mainly connected to the primary role of climate in driving large-scale microbial functional diversity as shown by our study and previous reports [9, 10] and also to the availability of future scenarios for climate but not for local variables. In this sense, the employed climate data set is the best one currently available for large-spatial scales and remote areas where there is no locally measured climatic data accessible. Despite our approach shows clear trends of increasing diversity and compositional change under future climate scenarios, there are uncertainties related to the exact magnitude of such changes which should be adjusted when more sophisticated current and future climate data become available [61, 62]. This recommendation is especially applicable to compositional change because of the uncertainty associated with interpolated precipitation provided by the climate gridded dataset . We note, however, that local physicochemical predictors vary within and among streams due to geological origins, climate variations, hydrological factors or anthropogenic disturbances, and have been demonstrated, together with landscape position, of substantial importance for microbial taxonomic diversity [64,65,66]. Although the inclusion of key local predictors will surely increase the accuracy of biodiversity projections under future climate scenarios, it currently constitutes a great challenge for large-scale studies. By assessing climate change effects on biodiversity from elevational gradients to large geographical extents and bearing in mind that the conclusions of our approach are subjected to the mentioned caveats, we here present a novel spatial exercise based on “elevation-for-latitude” substitution to provide the first general overview of climate change impacts on microbial functional diversity at the Eurasian scale.
In summary, our results demonstrate how climate constrains microbial functional diversity in lotic fresh waters and contribute to a better understanding of ecosystem responses to climate change. Warming and altered precipitation regimes under future climate will affect alpha diversity and lead to compositional turnover, these phenomena being especially evident in mid- and high-latitude regions across Eurasia. The predicted changes indicate a functional response of microorganisms linked to crucial biogeochemical processes, such as nutrient cycling, organic matter decomposition and greenhouse gas emissions, and anticipate biogeochemical feedback to ongoing climate change. Taken together, our findings serve as a baseline for estimating climate change impacts on stream processes through the study of microbial functional gene assemblages. We finally recommend further tests based on metatranscriptomics, metaproteomics and metabolomics approaches to support the validity of such findings and encourage future research to include field observations covering wider climate gradients and manipulative experiments across ecosystems and biogeochemical cycles, especially focusing on fresh waters and other nutrients in addition to carbon.
To span a wide climate gradient across Eurasia, we selected 52 stream sites located along three mountainsides from the subarctic, Mediterranean and subtropical regions (Additional file 1: Figure S1). For the subarctic region, we sampled the Bálggesvárri mountain in the Lyngen Alps Landscape Protected Area (Norway). This area presents a collection of glaciers and rocky peaks reaching a maximum elevation of 1833 m.a.s.l. For the Mediterranean region, samples were collected in Aigüestortes and Estany de Sant Maurici National Park. This area is known by its extensive set of glacial lakes located in the Pyrenees (Spain), a massive mountain range with its highest altitude reaching 3404 m.a.s.l. For the subtropical region, we selected the Laojun Mountain National Park, which is located at the junction of the Tibetan and Yungui Plateaus (China). This area presents a maximum elevation of 4513 m.a.s.l. and harbours typical mountain streams which, below 1800 m.a.s.l., turn into a deep, slowly discharging large river (Jinsha River, Upper Yangtze River). All these areas represent a landscape dominated by mixed evergreen and deciduous forests in lowlands, shrubs and grasses above the tree-line and bare rock peaks. Because landscape position has been reported of substantial importance for microbial taxonomic diversity [66,67,68], we selected for sampling only first, second at most, order streams to minimize the influence of such an additional confounding factor. Detailed information on the fieldwork procedures are described in Wang et al. . In short, the sampling sites were evenly spaced according to elevation and ranged from 18 to 771 m.a.s.l. in Bálggesvárri (n = 17), from 850 to 2500 m.a.s.l. in the Pyrenees (n = 17), and from 1828 to 4045 m.a.s.l. in Laojun (n = 18). Fieldwork was carried out in Autumn 2009 in China and Summer and Autumn 2012 in Norway and Spain, respectively. We started the sampling from the highest elevation and ended in the lowest parts of the streams, where the elevation stopped varying substantially. Based on stream width, five to ten cross sections were established at each location. Along the transects, we randomly selected 10 stones from riffle/run habitats and collected biofilm subsamples from the surface of each stone by scraping off a 9-cm2 area using a sterilized sponge. Subsamples were then pooled into a composite sample for each site and immediately frozen at − 18 °C.
We collected data on latitude, longitude and elevation as well as 11 spatial and physicochemical variables in situ or in the laboratory for each site (Additional file 1: Table S6) following the methods described in Wang et al. . Briefly, the latitude, longitude and altitude were taken using a GPS unit. Shading (% of canopy cover) was estimated from 10 locations in evenly spaced perpendicular transects that covered the whole study stretch. The water depth, current velocity, width and substratum grain size were obtained from 10 random locations within the same transects. We also measured the water temperature, conductivity and pH at each site. In the laboratory, we analysed chlorophyll-a by the 90% acetone extraction method  as well as total nitrogen and total phosphorus by peroxodisulphate oxidation and the spectrophotometry method , respectively. For Bálggesvárri Mountain, total nitrogen data were not available because the concentration level did not reach the minimum detection threshold required for the method employed. Additionally, we extracted 19 climatic variables with a spatial resolution of 30 arc seconds (~ 1 km) from WorldClim v.1.4. (www.worldclim.org/current) for each site. This information is based on interpolated climate station data and represents the average conditions for years 1960–1990  and constitutes the best climatic dataset currently available for large-spatial scales and remote zones where locally measured climatic data are not accessible.
We measured functional gene diversity and composition by applying GeoChip 4.0, a highly comprehensive (~ 82,000 probes covering 141,995 coding sequences) functional gene array widely used in biogeochemical, ecological and environmental analyses  because of its high optimization. GeoChip method presents several advantages compared to open-format technologies (e.g. metatranscriptomics, metaproteomics and metabolomics), such as high throughput, low detection limits, high reproducibility and potential for quantification, enabling a rapid performance with good correlations between target DNA concentrations and hybridization signal intensities . Microbial DNA was extracted from the biofilm samples using the phenol chloroform method  and purified using the QIAquick Gel Extraction Kit (QIAGEN Sciences, Germantown, MD, USA). The purified DNA was then quantified with a PicoGreen Kit (Eugene, OR, USA) and used for the GeoChip 4.0 hybridization. The DNA from each sample (500 ng) was labelled with the fluorescent dye Cy-3 (GE Healthcare, California, USA) by random priming . The DNA was purified as explained above and dried in a SpeedVac (Thermo Savant, New York, USA). Then, each dried and labelled DNA sample was resuspended in 42 μL of hybridization solution. This solution consisted of 1× HI-RPM hybridization buffer, 1× comparative genome hybridization blocking agent, 0.05 μg μL−1 Cot-1 DNA, 10 pM universal standard and 10% formamide (final concentrations). The solution was later denatured by remaining at 95 °C for 3 min. To remove the bubbles created during the denaturation process, the conditions were maintained at 37 °C for 30 min. The hybridizations were carried out at 67 °C for 24 h. Finally, the scanned images of the GeoChip hybridizations were obtained and converted by means of the Agilent Feature Extraction 11.5 software (Agilent Technologies, California, USA). For GeoChip data analyses, we employed signal intensity-based matrices as variation in functional gene abundances has an effect on diversity metrics and can also inform itself about environmental change . The signal intensities were quantified and processed based on previous pipelines  as follows: (i) removing probes with a signal-to-noise ratio of less than 2.0; (ii) normalizing the signal intensity of each probe by dividing the total signal intensity of a sample and multiplying by a constant;( iii) removing singletons, i.e. genes detected only once on each mountain; and (iv) selecting those genes involved in the carbon (C), nitrogen (N), phosphorus (P) and sulphur (S) cycles and stress-related (St) processes.
Mountain and continental scale patterns
Microbial functional genes are hierarchically organized, i.e. they are included in 88 functional gene families which are in turn grouped into five functional categories. Thus, alpha diversity and assemblage composition dissimilarity were calculated based on these two hierarchical levels. Firstly, we considered functional gene assemblages for the three kingdoms at the sampling site level which included all functional genes irrespective of their functional gene family and their functional category. Secondly, we repeated the procedure by grouping the functional genes into functional gene families. Then, we calculated the diversity and compositional dissimilarities of these functional gene assemblages. Given that various diversity metrics, such as observed richness, Chao1 estimated richness, Shannon-Wiener index, Simpson index and Inverse Simpson index, showed high correlations with each other (Pearson, r > 0.97 in all cases), we thus finally used alpha diversity based on the Shannon-Wiener index (base = exp(1)) because it is perhaps the most commonly employed diversity index in ecological research  and combines species richness (that is, functional gene richness in our case) and abundances (that is, signal intensities). For a more intuitive interpretation of diversity, we used true diversity, i.e. the effective number of attributes derived from the exponent-transformed Shannon entropy . Cross-site and cross-mountain differences in functional gene compositions (referred as turnover  in our study) were evaluated by Bray-Curtis dissimilarity, which accounts for the variation in assemblage structure due to changes in identities and abundances. We considered this measure of beta diversity because we chose to use quantitative data (i.e. signal intensities) given the importance of abundance for ecosystem functioning  and to exclude joint absences as this is an appropriate strategy for analysing communities along environmental gradients . We log-transformed signal intensities in all cases to emphasize the role of less-abundant genes as rare items can be essential in terms of ecosystem functioning . Across mountains, we tested for overall and pairwise differences in mean alpha diversity by ANOVA and Bonferroni-corrected t-tests, respectively. Global and pairwise differences in assemblage composition were assessed by PERMANOVA (1000 permutations), with Bonferroni-adjusted P values in pairwise comparisons. Both analyses were performed on the overall functional gene pool as well as on functional genes grouped into kingdoms and functional categories (i.e. C, N, P, S and St). Within mountains, the elevational gradients in alpha diversity and composition were examined on functional gene assemblages grouped into kingdoms and functional gene families (grouped by functional categories). For better comparison of model coefficients, elevation was z-standardized. The relationships between alpha diversity and elevation were tested by LMs, and the significance was determined by F-statistics. For kingdoms, we tested models including linear and quadratic terms. When both models were significant (P < 0.05) and significantly different (ANOVA; P < 0.05), we selected as the best model the one that minimized the Akaike’s information criterion adjusted for small sample size (AICc). For gene families, we considered only linear relationships as the alpha diversity-elevation relationships were estimated through model slope values. The relationships between assemblage dissimilarity and elevational Euclidean distances were examined by LMs for both kingdoms and functional gene families, with their significance determined by the Mantel test (Pearson correlations, 1000 permutations). Pairwise comparisons of mean slopes from LMs assessing the alpha diversity and compositional turnover response to elevation across functional gene families were tested by Bonferroni-corrected t-tests. For this analysis, gene families were grouped into functional categories and pairwise differences checked firstly with and secondly without accounting for the mountain. Because spatial patterns are most likely weak for those gene families containing low richness of functional genes, they were classified into richness-based quartiles and those within the lower quartile were excluded from the analyses. Finally, we examined the elevational patterns for 88 gene families (their distribution across the 5 targeted functional categories can be seen at Additional file 1: Table S1).
Pre-selection of climatic and local non-climatic predictors
To elucidate the drivers shaping the continental and mountain scale gradients in diversity and assemblage composition, we firstly selected ecologically relevant predictors that minimized collinearity to be included in the analyses. Prior to the initial selection, we ranked the individual influence of all predictors on alpha diversity and assemblage composition by RF analyses . We then included in further analyses the best ranked predictors that were not strongly correlated with the others (Pearson; r < 0.7; Additional file 1: Figure S14) prioritizing those predictors with a clear influence on freshwater microbial diversity at high elevations , such as those related to the growing season conditions, i.e. linked to the kinetics of biological processes, and the winter precipitation, i.e. mediating the hydrology in mountains through snow accumulation. This procedure, as well as the analyses themselves (see below), were applied to each dataset, i.e. all sites and the sites from individual mountains, on the overall functional gene assemblages.
Drivers underlying the biogeographical patterns
In order to get a comprehensive picture on how climatic and local non-climatic predictors contribute to the alpha diversity gradients, we employed an approach based on multiple statistical methods, including the HP, LMs and RFs. The combination of the HP and LMs provides robust evidence for meaningful predictors when both methods retain the same variables, and minimizes the probability of spurious results due to multicollinearity .
First, we used a multi-model approach to select the candidate model by running LMs on the full models (those containing all potential combinations of pre-selected single predictors) with alpha diversity as response variable. The models were ranked according to their AICc values and those with minimum AICc were chosen as candidate models. These models were further checked for spatial autocorrelation by Moran’s tests  and validated by visually checking their residuals for normality and homoscedasticity .
Second, we perform the HP, LM and RF analyses to evaluate the relationships between the predictors included in the candidate model and the response variable. For the HP analyses, the significance of the independent effects on the response variable was tested through a 1000 randomization-based procedure . For LMs, we z-standardized the predictors in order to better compare the coefficients within models. For all RFs (see also above), we set 3000 as the number of trees, and the number of predictors used in each split was one third of the included predictors; those variables showing negative or zero importance values were considered not predictive . For better comparison of the relative contributions of predictors, we transformed the variable importance scores into percentages relative to the sum of the importance values. To assess how predictors contributed to assemblage composition, we followed the same procedure employed for alpha diversity but using as a response variable the first axis coordinates from principal coordinate analysis (PCoA; Additional file 1: Figure S3). When assessing the cross-mountain gradient in composition, this procedure was further supported by GDMs, which predict spatial variation in assemblage composition between site pairs by modelling the response of biological dissimilarity matrices to environmental predictors. We run GDMs, plotted the I-splines for each predictor and obtained the impact of all individual predictors on the response dissimilarities (estimated as the variance explained by each predictor when the rest were kept constant). The uncertainty in the fitted I-splines was checked by plotting I-splines with error bands from 100 interactions-based bootstrapping , each data subsample retaining 70% of the sites from the full site-pair table. The significances of the full model and the individual predictors were assessed by a 100 permutation-based procedure . For both the PCoA and GDMs, the biological distances were obtained by applying the Bray-Curtis dissimilarity on signal intensity-based matrices previously log-transformed.
Additionally, we divided the variance in alpha diversity and composition associated with the climatic and local non-climatic predictors by variation partitioning  based on linear regression for alpha diversity and redundancy analysis  for the compositional matrices. In the latter case, signal intensity-based matrices were previously Hellinger-transformed to meet the requirements of linear ordination methods . The selection of predictors was performed by the ordiR2step procedure , which bases the forward selection of variables on adjusted R2 and P values of the initial pre-selected variables. When no predictors were retained for the local non-climatic set, the variance in alpha diversity in relation to the climatic predictors was checked by LMs, whereas the variance in composition was checked by a Mantel test.
Projections under climate change scenarios
We considered future climate conditions based on the most “optimistic”, one intermediate and the most “pessimistic” expectations of greenhouse emission rate scenarios, that is, the representative concentration pathways (RCPs) 2.6, 4.5 and 8.5. They broadly equal to atmospheric CO2 concentrations of 490, 650 and 1,370 p.p.m., respectively, by the end of this century . These future climate conditions based on the mentioned three CO2 emission scenarios come from the Fifth Assessment Report of The Intergovernmental Panel on Climate Change  and are available at the WorldClim site (www.worldclim.org/CMIP5v1) . We downloaded the climatic variables from all individual downscaled and bias-corrected climate model outputs available (n = 15–19) in the database at 2.5 arc minute spatial resolutions. We only considered the data that represent future climatic conditions averaged over the period from 2061 to 2080. For subsequent analyses, we used the ensemble means of each climatic variable over the individual climate model outputs.
We projected the changes in alpha diversity and shifts in assemblage composition over the entire Eurasian domain at a spatial resolution of 2.5 arc minutes combining the TWQ and PCQ as predictor variables, i.e. the best predictors of alpha diversity and assemblage composition, respectively. The current and future alpha diversities were obtained using LMs based on linear relationships. The variation in functional gene composition was projected using GDMs [44, 89] with a “space-for-time” substitution  based on Bray-Curtis dissimilarity matrices that were previously log-transformed. GDMs were run on the whole functional gene pool and on functional gene subsets for each gene family. The projections were shown on fluvial systems by masking the outputs with the global river classification layer . We excluded rivers that were classified as “very small” and “small” in the mask and applied a 5-km buffer around the streams for an adequate visualization of resulting continental-scale patterns.
All statistical analyses and plots were run in R statistical software V3.5.3 by using the packages vegan V2.54 , ggplot2 V3.1.1 , MuMIn V1.43.6 , hier.part V1.0-4 , gdm V1.3.11 , randomForestSRC V2.9.0 , raster V2.8-19  and corrplot .
Availability of data and materials
The GeoChip 4.0 dataset generated during the current study is available in the Gene Expression Omnibus repository with the accession number GSE128826, [www.ncbi.nlm.nih.gov/geo/]. Any other relevant data for this study are available from the corresponding author on reasonable request.
Falkowski PG, Frenchel T, Delong EF. The microbial engines that drive Earth’s biogeochemical cycles. Science. 2008;320:1034–8.
Von Schiller D, Acuña V, Aristi I, Arroitia M, Basaguren A, Bellin A, et al. River ecosystem processes: a synthesis of approaches, criteria of use and sensitivity to environmental stressors. Sci Total Environ. 2017;596-7:465–80.
Singh BK, Bardgett RD, Smith P, Reay DS. Microorganisms and climate change: terrestrial feedbacks and mitigation options. Nat Rev Microbiol. 2010;8:779–90.
McGill BJ, Enquist BJ, Weiher E, Westoby M. Rebuilding community ecology from functional traits. Trends Ecol Evol. 2006;21:178–85.
Violle C, Reich PB, Pacala SW, Enquist BJ, Kattge J. The emergence and promise of functional biogeography. Proc Natl Acad Sci U S A. 2014;111:13690–6.
Fierer N, Leff JW, Adams BJ, Nielsen UN, Bates ST, Lauber CL, et al. Cross-biome metagenomic analyses of soil microbial communities and their functional attributes. Proc Natl Acad Sci U S A. 2012;109:21390–5.
Tu Q, Yu H, He Z, Deng Y, Wu L, Van Nostrand JD, et al. GeoChip 4: a functional gene-array-based high-throughput environmental technology for microbial community analysis. Mol Ecol Resour. 2014;14:914–28.
Fuhrman JA, Steele JA, Hewson I, Schwalbach MS, Brown MV, Green JL, et al. A latitudinal diversity gradient in plaktonic marine bacteria. Proc Natl Acad Sci U S A. 2008;105:7774–8.
Wu B, Liu F, Weiser MD, Ning D, Okie JG, Shen L, et al. Temperature determines the diversity and structure of N2O-reducing microbial assemblages. Funct Ecol. 2018;32:1867–78.
Zhou J, Deng Y, Shen L, Wen CW, Yan Q, Ning D, et al. Temperature mediates continental-scale diversity of microbes in forest soils. Nat Commun. 2016;7:12083.
Wu L, Yang Y, Wang S, Yue H, Lin Q, Hu Y, et al. Alpine soil carbon is vulnerable to rapid microbial decomposition under climate cooling. ISME J. 2017;11:2102–11.
Xue K, Yuan MM, Shi ZJ, Qin Y, Deng Y, Cheng L, et al. Tundra soil carbon is vulnerable to rapid microbial decomposition under climate warming. Nat Clim Chang. 2016;6:595–600.
Brown JH, Gillooly JF, Allen AP, Savage VM, West GB. Toward a metabolic theory of ecology. Ecology. 2004;85:1771–89.
Angel R, Soares MI, Ungar ED, Gillor O. Biogeography of soil archaea and bacteria along a steep precipitation gradient. ISME J. 2010;4:553–63.
Zhou Z, Wang C, Luo Y. Response of soil microbial communities to altered precipitation: a global synthesis. Glob Ecol Biogeogr. 2018;27:1121–36.
Wiens JJ. Climate-related local extinctions are already widespread among plant and animal species. PLoS Biol. 2016;14:e2001104.
Lindström ES, Agterveld MPK-V, Zwart G. Distribution of typical freshwater bacterial groups is associated with ph, temperature, and lake water retention time. Appl Environ Microbiol. 2005;71:8201–6.
Hotaling S, Foley ME, Zeglin LH, Finn DS, Tronstad LM, Giersch JJ, et al. Microbial assemblages reflect environmental heterogeneity in alpine streams. Glob Chang Biol. 2019;25:2576–90.
Ylla I, Canhoto C, Romaní AM. Effects of warming on stream biofilm organic matter use capabilities. Microb Ecol. 2014;68:132–45.
Timoner X, Acuña V, Frampton L, Pollard P, Sabater S, Bunn SE. Biofilm functional responses to the rehydration of a dry intermittent stream. Hydrobiologia. 2014;727:185–95.
Camarero L, Catalan J. Atmospheric phosphorus deposition may cause lakes to revert from phosphorus limitation back to nitrogen limitation. Nat Commun. 2012;3:1–5.
Wieczynski DJ, Boyle B, Buzzard V, Duran SM, Henderson AN, Hulshof CM, et al. Climate shapes and shifts functional biodiversity in forests worldwide. Proc Natl Acad Sci U S A. 2018;116:587–92.
Yang Y, Gao Y, Wang S, Xu D, Yu H, Wu L, Lin Q, Hu Y, Li X, He Z, Deng Y, Zhou J. The microbial gene diversity along an elevation gradient of the Tibetan grassland. ISME J. 2014;8:430–40.
McMahon SM, Harrison SP, Armbruster WS, Bartlein PJ, Beale CM, Edwards ME, et al. Improving assessment and modelling of climate change impacts on global terrestrial biodiversity. Trends Ecol Evol. 2011;26:249–59.
IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge, New York: Cambridge University Press; 2013.
Franklin J. Moving beyond static species distribution models in support of conservation biogeography. Divers Distrib. 2010;16:321–30.
Barberán A, Casamayor EO, Fierer N. The microbial contribution to macroecology. Front Microbiol. 2014;5:203.
Sundqvist MK, Sanders NJ, Wardle DA. Community and ecosystem responses to elevational gradients: processes, mechanisms, and Insights for global change. Annu Rev Ecol Evol Syst. 2013;44:261–80.
Martiny JB, Eisen JA, Penn K, Allison SD, Horner-Devine MC. Drivers of bacterial beta-diversity depend on spatial scale. Proc Natl Acad Sci U S A. 2011;108:7850–4.
Hendershot JN, Read QD, Henning JA, Sanders NJ, Classen AT. Consistently inconsistent drivers of microbial diversity and abundance at macroecological scales. Ecology. 2017;98:1757–63.
Wang J, Meier S, Soininen J, Casamayor EO, Pan F, Tang X, et al. Regional and global elevational patterns of microbial species richness and evenness. Ecography. 2017;40:393–402.
Louca S, Polz MF, Mazel F, Albright MBN, Huber JA, O'Connor MI, et al. Function and functional redundancy in microbial systems. Nat Ecol Evol. 2018;2:936–43.
Hanson CA, Fuhrman JA, Horner-Devine MC, Martiny JB. Beyond biogeographic patterns: processes shaping the microbial landscape. Nat Rev Microbiol. 2010;10:497–506.
Qi Q, Zhao M, Wang S, Ma X, Wang Y, Gao Y, et al. The biogeographic pattern of microbial functional genes along an altitudinal gradient of the Tibetan Pasture. Front Microbiol. 2017;8:976.
McCain CM. Elevational gradients in diversity of small mammals. Ecology. 2005;86:366–75.
Neumann J. Latitudinal variation of tropospheric temperature lapse rate. Archiv für Meteorologie Geophysik und Bioklimatologie Serie A. 1955;8:351–3.
Fierer N, Jackson RB. The diversity and biogeography of soil bacterial communities. Proc Natl Acad Sci U S A. 2006;103:626–31.
Tedersoo L, Bahram M, Põlme S, Kõljalg U, Yorou NS, Wijesundera R, et al. Global diversity and geography of soil fungi. Science. 2014;346:1052–3.
Mac Nally R. Multiple regression and inference in ecology and conservation biology: further comments on identifying important predictor variables. Biol Conserv. 2002;11:1397–401.
Ishwaran H, Kogalur UB, Blackstone EH, Lauer MS. Random survival forests. Ann Appl Stat. 2008;2:841–60.
Körner C. The use of ‘altitude’ in ecological research. Trends Ecol Evol. 2007;22:569–74.
Hawkins BA, Field R, Cornell HV, Currie DJ, Guégan J-F, Kaufman DM, et al. Energy, water, and broad-scale geographic patterns of species richness. Ecology. 2003;84:3105–17.
Wang Q, Zhao X, Chen L, Yang Q, Chen S, Zhang W, et al. Global synthesis of temperature sensitivity of soil organic carbon decomposition: latitudinal patterns and mechanisms. Funct Ecol. 2018;33:514–23.
Ferrier S, Manion G, Elith J, Richardson K. Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Divers Distrb. 2007;13:252–64.
Savio D, Sinclair L, Ijaz UZ, Parajka J, Reischer GH, Stadler P, et al. Bacterial diversity along a 2600 km river continuum. Environ Microbiol. 2015;17:4994–5007.
Zeglin LH. Stream microbial diversity in response to environmental changes: review and synthesis of existing research. Front Microbiol. 2015;6:454.
Woodward FI, Lomas MR, Kelly CK. Global climate and the distribution of plant biomes. Philos Trans R Soc B. 2004;359:1465–76.
Besemer K. Biodiversity, community structure and function of biofilms in stream ecosystems. Res Microbiol. 2015;166:774–81.
Leff JW, Jones SE, Prober SM, Barberán A, Borer ET, Firn JL, et al. Consistent responses of soil microbial communities to elevated nutrient inputs in grasslands across the globe. Proc Natl Acad Sci U S A. 2015;112:10967–72.
Ruiz-González C, Niño-García JP, Lapierre J-F, del Giorgio PA. The quality of organic matter shapes the functional biogeography of bacterioplankton across boreal freshwater ecosystems. Glob Ecol Biogeogr. 2015;24:1487–98.
Burnham KP, Anderson DR. Model selection and multimodel inference: a practical information-theoretic approach. 2nd ed. New York: Springer-Verlag; 2002.
Verbeek L, Gall A, Hillebrand H, Striebel M. Warming and oligotrophication cause shifts in freshwater phytoplankton communities. Glob Chang Biol. 2018;24:4532–43.
Moss RH, Edmonds JA, Hibbard KA, Manning MR, Rose SK, van Vuuren DP, et al. The next generation of scenarios for climate change research and assessment. Nature. 2010;463:747–56.
Ma B, Dai Z, Wang H, Dsouza M, Liu X, He Y, et al. Distinct biogeographic patterns for archaea, bacteria, and fungi along the vegetation gradient at the continental scale in Eastern China. mSystems. 2017;2.
Ravenscroft CH, Whitlock R, Fridley JD. Rapid genetic divergence in response to 15 years of simulated climate change. Glob Chang Biol. 2015;21:4165–76.
Bardgett RD, Freeman C, Ostle NJ. Microbial contributions to climate change through carbon cycle feedbacks. ISME J. 2008;2:805–14.
Bradford MA, McCulley RL, Crowther TW, Oldfield EE, Wood SA, Fierer N, et al. Cross-biome patterns in soil microbial respiration predictable from evolutionary theory on thermal adaptation. Nat Ecol Evol. 2019;3:223–31.
Bjorkman AD, Myers-Smith IH, Elmendorf SC, Normand S, Rüger N, Beck PSA, et al. Plant functional trait change across a warming tundra biome. Nature. 2018;562:57–62.
Johnston ASA, Sibly RM. The influence of soil communities on the temperature sensitivity of soil respiration. Nat Ecol Evol. 2018;2:1597–602.
Beaulieu JJ, Tank JL, Hamilton SK, Wollheim WM, Hall RO Jr, Mulholland PJ, et al. Nitrous oxide emission from denitrification in stream and river networks. Proc Natl Acad Sci U S A. 2011;108:214–9.
Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37:4302–15.
Zelinka MD, Myers TA, McCoy DT, Po-Chedley S, Caldwell PM, Ceppi P, Klein SA, Taylor KE. Causes of higher climate sensitivity in CMIP6 models. Geophys Res Lett. 2020;47:e2019GL085782.
Marchi M, Sinjur I, Bozzano M, Westergren M. Evaluating WorldClim Version 1 (1961–1990) as the baseline for sustainable use of forest and environmental resources in a changing climate. Sustainability. 2019;11:3043.
Esposito A, Engel M, Ciccazzo S, Daprà L, Penna D, Comiti F, et al. Spatial and temporal variability of bacterial communities in high alpine water spring sediments. Res Microbiol. 2016;167:325–33.
Fegel TS, Baron JS, Fountain AG, Johnson GF, Hall EK. The differing biogeochemical and microbial signatures of glaciers and rock glaciers. J Geophys Res-Biogeo. 2016;121:919–32.
Ruiz-González C, Niño-García JP, Del Giorgio PA. Terrestrial origin of bacterial communities in complex boreal freshwater networks. Ecol Lett. 2015;18:1198–206.
Besemer K, Singer G, Quince C, Bertuzzo E, Sloan W, Battin TJ. Headwaters are critical reservoirs of microbial diversity for fluvial networks. Proc R Soc Lond B. 2013;280:20131760.
Crump BC, Amaral-Zettler LA, Kling GW. Microbial diversity in arctic freshwaters is structured by inoculation of microbes from soils. ISME J. 2012;6:1629–39.
Dalsgaard T, Nielsen LP, Brotas V, Viaroli P, Underwood G, Nedwell D, et al. Protocol handbook for NICE - Nitrogen Cycling in Estuaries: a project under the EU research programme: Marine Science and Technology (MAST III). Silkeborg: National Environmental Research Institute; 2000.
Jin XC, Tu QY. The standard methods for observations and analysis in lake eutrophication. 2nd ed. Beijing: China Environmental Science Press; 1990.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.
Zhou J, He Z, Yang Y, Deng Y, Tringe SG, Alvarez-Cohen L. High-throughput metagenomic technologies for complex microbial community analysis: open and closed formats. mBio. 2015;6.
Zhou J. DNA recovery from soils of diverse composition. Appl Environ Microbiol. 1996;62:316–22.
He Z, Deng Y, Van Nostrand JD, Tu Q, Xu M, Hemme CL, et al. GeoChip 3.0 as a high-throughput tool for analyzing microbial community composition, structure and functional activity. ISME J. 2010;4:1167–79.
Spellerberg I, Fedor PJ. A tribute to Claude Shannon (1916-2001) and a plea for more rigorous use of species richness, species diversity and the ‘Shannon-Wiener’ Index. Glob Ecol Biogeogr. 2003;12:177–9.
Jost L. Entropy and diversity. Oikos. 2006;113:363–75.
Anderson MJ, Crist TO, Chase JM, Vellend M, Inouye BD, Freestone AL, et al. Navigating the multiple meanings of beta diversity: a roadmap for the practicing ecologist. Ecol Lett. 2011;14:19–28.
Winfree R, Fox JW, Williams NM, Reilly JR, Cariveau DP. Abundance of common species, not species richness, drives delivery of a real-world ecosystem service. Ecol Lett. 2015;18:626–35.
Jain M, Flynn DFB, Prager CM, Hart GM, DeVan CM, Ahrestani FS, et al. The importance of rare species: a trait-based assessment of rare species contributions to functional diversity and possible ecosystem function in tall-grass prairies. Ecol Evol. 2014;4:104–12.
Ishwaran H. Variable importance in binary regression trees and forests. Electron J Stat. 2007;1:519–37.
Bivand RS, Wong DWS. Comparing implementations of global and local indicators of spatial association. Test. 2018;27:716–48.
Zuur A, Ieno EN, Smith GM. Analyzing Ecological Data. New York: Springer; 2007.
Shryock DF, Havrilla CA, DeFalco LA, Esque TC, Custer NA, Wood TE. Landscape genomics of Sphaeralcea ambigua in the Mojave Desert: a multivariate, spatially-explicit approach to guide ecological restoration. Conserv Genet. 2015;16:1303–17.
Manion G, Lisk M, Ferrier S, Nieto-Lugilde D, Mokany K, Fitzpatrick MC. gdm: Generalized Dissimilarity Modeling. R package version 1.3.11; 2018.
Anderson MJ, Gribble NA. Partitioning the variation among spatial, temporal and environmental components in a multivariate data set. Aust J Ecol. 1998;23:158–67.
Rao CR. The use and interpretation of principal component analysis in applied research. Sankhyā. 1964;26:329–58.
Legendre P, Gallagher ED. Ecologically meaningful transformations for ordination of species data. Oecologia. 2001;129:271–80.
Blanchet FG, Legendre P, Borcard D. Forward selection of explanatory variables. Ecology. 2008;89:2623–32.
Lehman A, Overton JM, Austin MP. Regression models for spatial prediction: their role for biodiversity and conservation. Biodivers Conserv. 2002;11:2085–92.
Blois JL, Williams JW, Fitzpatrick MC, Jackson ST, Ferrier S. Space can substitute for time in predicting climate-change effects on biodiversity. Proc Natl Acad Sci U S A. 2013;110:9374–9.
Dallaire CO, Lehner B, Sayre R, Thieme M. A multidisciplinary framework to derive global river reach classifications at high spatial resolution. Environ Res Lett. 2019;14:024003.
Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.5-2. 2018. https://CRAN.R-project.org/package=vegan.
Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag; 2016.
Bartoń, K. MuMIn: Multi-Model Inference. R package version 1.40.44. 2018. https://CRAN.R-project.org/package=MuMIn.
Walsh, C. & Mac Nally, R. hier.part: Hierarchical Partitioning. R package version 1.0-4. 2013. https://CRAN.R-project.org/package=hier.part.
Ishwaran H, Kogalur UB. Random Forests for Survival, Regression, and Classification (RF-SRC). R package version 2.6.1; 2018.
Hijmans, R. J. raster: Geographic Data Analysis and Modeling. R package version 2.7-15. 2018. https://CRAN.R-project.org/package=raster.
Wei T, Simko V. corrplot: Visualization of a correlation matrix. R package version 0.77. 2017. https://github.com/taiyun/corrplot.
We thank T. Yuan, Y. Zhang, K. Yang, X. Triadó-Margarit and J. Eskelinen for their support on field sampling or lab analyses. We also thank P. Abellán and C. Gutiérrez-Cánovas for comments on data analyses.
The study was supported by National Natural Science Foundation of China (91851117), the Program of Global Change and Mitigation (2017YFA0605200), CAS Key Research Program of Frontier Sciences (QYZDB-SSW-DQC043), CAS Strategic Pilot Science and Technology (XDA20050101), and National Natural Science Foundation of China (41571058, 41871048). AV was supported by the Chinese Academy of Sciences President’s International Fellowship Initiative (2018PS0007). EOC was supported by grant BRIDGES CGL2015-69043-P (Spanish Office for Science-MINECO-ERDF).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Picazo, F., Vilmi, A., Aalto, J. et al. Climate mediates continental scale patterns of stream microbial functional diversity. Microbiome 8, 92 (2020). https://doi.org/10.1186/s40168-020-00873-2