Philippa Kaur

Our general strategy was to compare the performance of four approaches for inferring microbial associations from abundance data with overlying time-series signals. The approaches were (1) pairwise spearman correlation analysis (SCC) [1, 29], (2) Graphical lasso analysis (Glasso) [30, 31], (3) pairwise SCC analysis with a pre-processing step where seasonal and long-term splines were fit to and subtracted from each variable using a GAM (GAM-SCC), and (4) Glasso with the same GAM subtraction approach (GAM-Glasso). Our validation strategy for the GAM transformation consisted of generating mock datasets with underlying associations, masking those associations by adding seasonal and long-term signals to the abundance data, and comparing the predicted associations obtained from each network inference method to the true species-species associations.Data simulation: generating mock abundance data with time-series propertiesWe generated mock abundance datasets that had a predetermined, underlying network structure and contained long-term and seasonal species abundance patterns. First, a covariance matrix was generated to describe the relationships between species in a mock dataset (Fig. S1, Panel 1). The covariance matrices were constructed with underlying network structures that followed either a scale-free Barabási-Albert model, a random Erdős-Rényi model, or a model of network topology based on a real microbial dataset (American Gut dataset; Fig. S1) [32, 33]. The Erdős-Rényi and Barabási-Albert model datasets were generated so that each dataset contained 400 species and 200 samples, and the American Gut datasets were created so that each dataset contained 127 species and 200 samples. A random Bernoulli distribution was used to simulate the covariance matrix for the Erdős-Rényi networks. We set the probability of interactions occurring between species in a given Erdős-Rényi network to 1%. The Barabási-Albert networks were generated using the “sample_pa” function in the igraph package [34]. The “graph2prec” function in the SpiecEasi package was used to predict the covariance matrix of the American Gut dataset [33]. The covariance between species in a dataset was considered “high” or “low” when the true associations in the covariance matrix were set to 100 or 10 respectively (Fig. S1, Panel 1). These covariance matrices describe the “real”, underlying species interactions in our mock datasets.After generating a covariance matrix, the mean abundance for each species was generated from a normal distribution with a mean of 10 and a variance of 1. These mean abundance values and the covariance matrix were used to parameterize a multivariate normal distribution from which species abundance values for all 200 samples in a dataset were drawn (Fig. S1, Panel 2). The values generated from this multivariate normal distribution were the species abundance values without time-series features confounding the relationship between two associated species (Fig. S1, Panel 2).“Gradual” or “abrupt” seasonal trends were added to 0%, 25%, 50% or 100% of the species in each mock dataset. The gradual seasonal trend increased over 5 months, peaked at a specific month, and decreased over 5 months. Conversely, the abrupt seasonal signal increased over 2 months, peaked at a specific month, and decreased over 2 months (Fig. S1, Panel 3). These seasonal signals were generated by plugging a vector of consecutive integers of length 200 (Nt) into the gradual (Eq. (1)) or abrupt (Eq. (2)) seasonal equations (Fig. S1, Panel 3)…$$Gradual:S_t = left( {frac{{cos left( {N_t ast 2 ast frac{pi }{{12}}} right)}}{2}} right) + 0.5$$
(1)
$$Abrupt:,S_t = left( {left( {frac{{cos left( {N_t ast 2 ast frac{pi }{{12}}} right)}}{2}} right) + 0.5} right)^{10}$$
(2)
where N is the random vector of consecutive integers, S is the output seasonal vector, and t is the index of vectors N and S. The starting value of vector Nt was drawn at random for each species to allow the seasonal peaks to be centered at different months. Each element in the seasonal vector (St) was then multiplied by the corresponding element in the abundance vector (Xt) of a specific species to obtain mock species abundance values with a gradual or abrupt seasonal trend (Fig. S1, Panel 3).A long-term time-series trend was added to the abundance values of 0% or 50% of the species in the mock datasets (Fig. S1, Panel 4). When a long-term signal was applied to 50% of the species in a dataset, half of the species were randomly selected to have this long-term trend. Then, a vector of linear values was generated following Eq. (3) such that…$$Long – term,trend:,L_t = pm mleft( {L_{t – 1}} right) + 0.01$$
(3)
where Lt is the point in the line at the next time point and m is the slope of the line. The slope parameter (m) was generated from a random normal distribution with a mean of 0.01 and a variance of 0.01. The slope parameter (m) was also multiplied by −1 half of the time to ensure that half of the long-term trends increased over time and half decreased over time (Fig. S1, Panel 4). After generating the vector of linear values (Lt), each element of this vector was added to each element of the abundance vector (Xt) of a specific species to simulate long-term time-series trends (Fig. S1, Panel 4).Time-series predictor columns were added to each dataset after applying monthly and long-term abundance trends to a portion of the species in the mock datasets. The predictors that were used in the downstream GAM-based data transformation were the month of the year (i.e., 1–12) and the day of the time-series (i.e., 1–200). In total, we generated 100 mock datasets for every combination of conditions (84 combinations total; Table S1), resulting in 8400 mock time-series datasets that were used in the downstream count data transformation, GAM subtraction, and network analysis procedures.Data simulation: Simulating count data from abundance valuesThe 8400 time-series datasets that were generated using the methods described above were transformed to make the abundance values resemble high-throughput sequencing data because microbial time-series sampling efforts are often processed using such molecular methods (e.g., tag-sequencing, meta-omics). Analysis of high-throughput sequencing data is complicated by the compositional (i.e., relative) nature of the data and by the high number of zeros that may be prevalent in a dataset (i.e., zero-inflation; see Supplementary Information) [35, 36]. Relative abundances of different species in natural communities are also highly skewed, so that relatively few species constitute most of the organisms in a sample although many rare species are also present [37, 38]. Therefore, species abundances were first exponentiated to increase the prevalence of abundant species and to decrease the prevalence of rare species (Fig. S1, Panel 5). The exponentiated species abundance values were then converted to relative abundance values by dividing each species count by the sum of all species counts in a sample (Fig. S1, Panel 6). The resulting relative abundance values and time-series predictor variables were used in data normalization and GAM-transformation steps prior to carrying out the network analyses.Network inference: Count data normalization and GAM transformationSeveral steps were taken to back out the species-species relationships in the mock datasets. We advocate these steps to infer network structure from a real time-series dataset. A centered log-ratio (CLR) transformation was first applied to the species relative abundance values to normalize the mock species abundance data across samples using the “clr” function in the compositions package in R (Fig. 1) [39]. This transformation step is important to avoid spurious inferences induced by the inherent compositionality of relative abundance data [31, 33, 36]. In addition to the CLR transformation used in our main network iterations, we carried out additional network iterations using the modified CLR [40], cumulative sum scaling [41], and total sum scaling [42] transformations (see Supplementary Information). In all cases, the normalized dataset was copied, with one copy subjected to a subsequent GAM transformation, and the other one not GAM-transformed.Fig. 1: Steps used to carry out the GAM-based transformation of time-series species abundance data prior to carrying out pairwise spearman correlation (SCC) and graphical lasso (Glasso) ecological network analyses.The raw, species abundance data were first CLR-transformed (1). Generalized additive models (GAMs) were then fit to each species in the dataset (2) and the residuals of each GAM were checked for significant autocorrelation (3). The residuals of each GAM were extracted (4) and were used as input in the SCC and Glasso network analysis methods (5). Finally, the GAM-transformed network outputs were obtained (6; see text for additional details).Full size imageThe GAM transformation was carried out by fitting GAMs to each individual species in the dataset to remove monthly signals, long-term trends, and autocorrelation from the species abundance data. These GAMs were fit using the “gamm” function in the mgcv package in R [43, 44]. The GAMs that were used included the “month of year” parameter as a cyclical spline predictor and the “day of time-series” parameter as a penalized thin-plate spline predictor (“ts” in the mgcv package; Fig. 1), which given our one-dimensional data is analogous to a natural cubic spline. In addition, the first GAM included a continuous AR1 (“corCAR1” in the mgcv package) correlation structure term in the model. This corCAR1 model was revised for specific species when the GAM could not be resolved or when significant autocorrelation was detected in the GAM residuals (Fig. 1). The GAM revision step fit 4 new GAMs with different correlation structure terms (i.e., “AR1”, “CompSymm”, “Exp”, and “Gaus”) to the species that could not be fit using the corCAR1 model or that contained significant autocorrelation in the corCAR1 GAM residuals. Then, the correlation structure term that addressed these issues for the largest number of individuals was used as the GAM model for this group of species. After fitting a GAM to all of the species in the input dataset, the residuals of each GAM were extracted and were used as the new, GAM-transformed abundance values (Fig. 1). These GAM residuals represent species abundance values with a reduced influence of time (Fig. 2) and were used as input in the downstream GAM-SCC and GAM-Glasso network analyses.Fig. 2: A conceptual figure that demonstrates how the GAM transformation can remove seasonal signals while preserving ecologically relevant species co-occurrence patterns.In this example, the co-occurrence pattern between Species A and Species B persists even after the seasonal signals are removed by the GAM transformation.Full size imageNetwork inference: Network runs and statistical analysesThe pre-processed species abundance data with and without the GAM-removal of time-series signals were used in SCC and Glasso networks in order to compare the outputs of the SCC, GAM-SCC, Glasso, and GAM-Glasso network inference approaches (Fig. 1). Additional network iterations were also carried out using the CCLasso [45] and SPRING [40] network inference approaches (see Supplementary Information). For the SCC and Glasso network iterations, a nonparanormal transformation was applied to the species abundance datasets with and without the GAM transformation using the “huge.npn” function in the huge package in R [46]. Spearman correlation networks were then constructed by calculating the correlation between every pair of species in the mock abundance datasets. A Bonferroni-corrected p value of 0.01 was used as a cutoff to identify edges in these SCC networks. The Glasso networks were constructed by testing 30 regularization parameter values (i.e., lambdas) in each network using the “batch.pulsar” (criterion = “stars”; rep.num = 50) function in the pulsar package in R [47]. The lambda that resulted in the most stable network output was selected using the StARS method [48]. Finally, the graph that resulted from the StARS output was used to obtain a species adjacency matrix for the Glasso networks.The species-species associations predicted by the SCC, GAM-SCC, Glasso, and GAM-Glasso networks were compared to the true species-species associations and the F1 scores of the network predictions were calculated. The F1 score is a measure of classification performance (presence or absence of an edge) that accounts for uneven classes, which is essential when dealing with sparse networks. The F1 scores of the GAM-transformed networks were compared to the networks that did not undergo GAM transformation using paired Wilcoxon tests with Bonferroni correction. An adjusted p value of 0.01 was used as a cutoff to identify under what circumstances the GAM significantly improved the F1 score of a Glasso or SCC network.Network inference: Comparison of predicted network structuresAdditional networks were generated using the methods described above to compare the predicted network structures obtained from the GAM-Glasso, Glasso, GAM-SCC, and SCC approaches to the real network structures. These additional networks were constructed using smaller mock datasets to allow for better visualization of the network outputs and contained species with a gradual seasonal signal and high species-species covariance (see Supplementary Information). The average clustering coefficient and the degree distribution of these additional network outputs were calculated and used for the network structure comparisons. The average clustering coefficient of a network describes the likelihood that two species that are both associated with a third species are also associated with each other [49], and in a sense describes the “clumpiness” of a network. The network degree distributions describe the probability distribution of the number of interactions per node in a network [50]. More

Hughes, T. P. et al. Coral reefs in the Anthropocene. Nature 546, 82–90 (2017).ADS 
CAS 
PubMed 

Google Scholar 
Hughes, T. P. et al. Global warming and recurrent mass bleaching of corals. Nature 543, 373–377 (2017).ADS 
CAS 
PubMed 

Google Scholar 
Hughes, T. P. et al. Spatial and temporal patterns of mass bleaching of corals in the Anthropocene. Science (80- ). 359, 80–83 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Eakin, C. M., Sweatman, H. P. A. & Brainard, R. E. The 2014–2017 global-scale coral bleaching event: Insights and impacts. Coral Reefs 38, 539–545 (2019).ADS 

Google Scholar 
Glynn. Coral reef bleaching: Facts, hypotheses and implications. Glob. Chang. Biol. 2, 495–509 (1996).ADS 

Google Scholar 
Brown, B. E. Coral bleaching: Causes and consequences. Coral Reefs 16, 129–138 (1997).
Google Scholar 
Maynard, J. A. et al. Projections of climate conditions that increase coral disease susceptibility and pathogen abundance and virulence. Nat. Clim. Chang. 5, 688–694 (2015).ADS 

Google Scholar 
Hughes, T. P. et al. Global warming transforms coral reef assemblages. Nature 556, 492–496 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Anthony, K. R. N., Kline, D. I., Diaz-Pulido, G., Dove, S. & Hoegh-Guldberg, O. Ocean acidification causes bleaching and productivity loss in coral reef builders. Proc. Natl. Acad. Sci. U. S. A. 105, 17442–17446 (2008).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Huang, H. et al. Positive and negative responses of coral calcification to elevated pCO2: Case studies of two coral species and the implications of their responses. Mar. Ecol. Prog. Ser. 502, 145–156 (2014).ADS 
CAS 

Google Scholar 
Hoadley, K. D. et al. Physiological response to elevated temperature and pCO2 varies across four Pacific coral species: Understanding the unique host + symbiont response. Sci. Rep. 5, 1–15 (2015).
Google Scholar 
Schoepf, V. et al. Coral energy reserves and calcification in a high-CO2 world at two temperatures. PLoS One. 8, e75049 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
IPCC. In IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, (eds. Pörtner, H.-O. et al.) 1–36 (Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2019).Bahr, K. D., Jokiel, P. L. & Rodgers, K. S. Relative sensitivity of five Hawaiian coral species to high temperature under high-pCO2 conditions. Coral Reefs 35, 729–738 (2016).ADS 

Google Scholar 
Dove, S. G., Brown, K. T., Van Den Heuvel, A., Chai, A. & Hoegh-Guldberg, O. Ocean warming and acidification uncouple calcification from calcifier biomass which accelerates coral reef decline. Commun. Earth Environ. 1, 1–9 (2020).
Google Scholar 
Chow, M. H., Tsang, R. H. L., Lam, E. K. Y. & Ang, P. O. Quantifying the degree of coral bleaching using digital photographic technique. J. Exp. Mar. Bio. Ecol. 479, 60–68 (2016).
Google Scholar 
Amid, C. et al. Additive effects of the herbicide glyphosate and elevated temperature on the branched coral Acropora formosa in Nha Trang, Vietnam. Environ. Sci. Pollut. Res. 25, 13360–13372 (2018).CAS 

Google Scholar 
Anthony, K. R. N., Connolly, S. R. & Willis, B. L. Comparative analysis of energy allocation to tissue and skeletal growth in corals. Limnol. Oceanogr. 47, 1417–1429 (2002).ADS 

Google Scholar 
Edmunds, P. J. & Davies, P. S. An energy budget for Porites porites (Scleractinia). Mar. Biol. 92, 339–347 (1986).
Google Scholar 
Stimson, J. S. Location, quantity and rate of change in quantity of lipids in tissue of Hawaiian hermatypic corals. Bull. Mar. Sci. 41, 889–904 (1987).ADS 

Google Scholar 
Harland, A. D., Navarro, J. C., Spencer Davies, P. & Fixter, L. M. Lipids of some Caribbean and Red Sea corals: Total lipid, wax esters, triglycerides and fatty acids. Mar. Biol. 117, 113–117 (1993).CAS 

Google Scholar 
Grottoli, A. G., Tchernov, D. & Winters, G. Physiological and biogeochemical responses of super-corals to thermal stress from the northern gulf of Aqaba, Red Sea. Front. Mar. Sci. 4, 1–12 (2017).
Google Scholar 
Rodrigues, L. J. & Grottoli, A. G. Energy reserves and metabolism as indicators of coral recovery from bleaching. Limnol. Oceanogr. 52, 1874–1882 (2007).ADS 

Google Scholar 
Anthony, K. R. N., Hoogenboom, M. O., Maynard, J. A., Grottoli, A. G. & Middlebrook, R. Energetics approach to predicting mortality risk from environmental stress: A case study of coral bleaching. Funct. Ecol. 23, 539–550 (2009).
Google Scholar 
Baumann, J. H., Grottoli, A. G., Hughes, A. D. & Matsui, Y. Photoautotrophic and heterotrophic carbon in bleached and non-bleached coral lipid acquisition and storage. J. Exp. Mar. Bio. Ecol. 461, 469–478 (2014).CAS 

Google Scholar 
Hughes, A. D. & Grottoli, A. G. Heterotrophic compensation: A possible mechanism for resilience of coral reefs to global warming or a sign of prolonged stress?. PLoS ONE 8, 1–10 (2013).
Google Scholar 
Grottoli, A. G. et al. The cumulative impact of annual coral bleaching can turn some coral species winners into losers. Glob. Chang. Biol. 20, 3823–3833 (2014).ADS 
PubMed 

Google Scholar 
Grottoli, A. G., Rodrigues, L. J. & Palardy, J. E. Heterotrophic plasticity and resilience in bleached corals. Nature 440, 1186–1189 (2006).ADS 
CAS 
PubMed 

Google Scholar 
Levas, S. J. et al. Can heterotrophic uptake of dissolved organic carbon and zooplankton mitigate carbon budget deficits in annually bleached corals?. Coral Reefs 35, 495–506 (2016).ADS 

Google Scholar 
Jury, C. P., Delano, M. N. & Toonen, R. J. High heritability of coral calcification rates and evolutionary potential under ocean acidification. Sci. Rep. 9, 1–13 (2019).
Google Scholar 
Jury, C. P. & Toonen, R. J. Adaptive responses and local stressor mitigation drive coral resilience in warmer, more acidic oceans. Proc. R. Soc. B Biol. Sci. 286, 20190614 (2019).
Google Scholar 
Concepcion, G. T., Polato, N. R., Baums, I. B. & Toonen, R. J. Development of microsatellite markers from four Hawaiian corals: Acropora cytherea, Fungia scutaria, Montipora capitata and Porites lobata. Conserv. Genet. Resour. 2, 11–15 (2010).

Google Scholar 
Gorospe, K. D. & Karl, S. A. Genetic relatedness does not retain spatial pattern across multiple spatial scales: Dispersal and colonization in the coral, Pocillopora damicornis. Mol. Ecol. 22, 3721–3736 (2013).PubMed 

Google Scholar 
Wall, C. B., Ritson-Williams, R., Popp, B. N. & Gates, R. D. Spatial variation in the biochemical and isotopic composition of corals during bleaching and recovery. Limnol. Oceanogr. 64, 2011–2028 (2019).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Bahr, K. D., Tran, T., Jury, C. P. & Toonen, R. J. Abundance, size, and survival of recruits of the reef coral Pocillopora acuta under ocean warming and acidification. PLoS ONE 15, 1–13 (2020).
Google Scholar 
Rogelj, J. et al. Paris agreement climate proposals need a boost to keep warming well below 2 °C. Nature 534, 631–639 (2016).ADS 
CAS 
PubMed 

Google Scholar 
McLachlan, R. H., Price, J. T., Solomon, S. L. & Grottoli, A. G. Thirty years of coral heat-stress experiments: A review of methods. Coral Reefs 39, 885–902 (2020).
Google Scholar 
Grottoli, A. G. et al. Increasing comparability among coral bleaching experiments. Ecol. Appl. 31, e02262 (2021).CAS 
PubMed 
PubMed Central 

Google Scholar 
Grottoli, A. G. Variability of stable isotopes and maximum linear extension in reef-coral skeletons at Kaneohe Bay, Hawaii. Mar. Biol. 135, 437–449 (1999).
Google Scholar 
McLachlan, R. H., Dobson, K. L., Grottoli, A. G. Quantification of Total Biomass in Ground Coral Samples. Protocols.io (2020). https://doi.org/10.17504/protocols.io.bdyai7se.McLachlan, R. H., Muñoz-Garcia, A., Grottoli, A. G. Extraction of Total Soluble Lipid from Ground Coral Samples. Protocols.io (2020). https://doi.org/10.17504/protocols.io.bc4qiyvw.McLachlan, R. H., Price, J. T., Dobson, K. L., Weisleder, N. & Grottoli, A. G. Microplate Assay for Quantification of Soluble Protein in Ground Coral Samples. Protocols.io (2020). https://doi.org/10.17504/protocols.io.bdc8i2zw.McLachlan, R. H., Juracka, C. & Grottoli, A. G. Symbiodiniaceae Enumeration in Ground Coral Samples Using Countess™ II FL Automated Cell Counter. Protocols.io (2020). https://doi.org/10.17504/protocols.io.bdc5i2y6.McLachlan, R. H. & Grottoli, A. G. Geometric Method for Estimating Coral Surface Area Using Image Analysis. Protocols.io https://doi.org/10.17504/protocols.io.bdyai7se(2021).Muscatine, L., McCloskey, L. R. & Marian, R. E. Estimating the daily contribution of carbon from zooxanthellae to coral animal respiration. Limnol. Oceanogr. 26, 601–611 (1981).ADS 
CAS 

Google Scholar 
Levas, S. J. et al. Organic carbon fluxes mediated by corals at elevated pCO2 and temperature. Mar. Ecol. Prog. Ser. 519, 153–164 (2015).ADS 
CAS 

Google Scholar 
Perry, C. T. et al. Loss of coral reef growth capacity to track future increases in sea level. Nature 558, 396–400 (2018).ADS 
CAS 
PubMed 

Google Scholar 
Woodley, C. M., Burnett, A. & Downs, C. A. Epidemiological Assessment of Reproductive Condition of ESA Priority Coral (2013).Logan, C. A., Dunne, J. P., Eakin, C. M. & Donner, S. D. Incorporating adaptive responses into future projections of coral bleaching. Glob. Chang. Biol. 20, 125–139 (2014).ADS 
PubMed 

Google Scholar 
Rodrigues, L. J., Grottoli, A. G. & Lesser, M. P. Long-term changes in the chlorophyll fluorescence of bleached and recovering corals from Hawaii. J. Exp. Biol. 211, 2502–2509 (2008).PubMed 

Google Scholar 
Rowan, H. et al. Environmental gradients drive physiological variation in Hawaiian corals. Coral Reefs 40(5), 1505–1523. https://doi.org/10.1007/s00338-021-02140-8 (2021).Article 

Google Scholar 
Houlbrèque, F. & Ferrier-Pagès, C. Heterotrophy in tropical scleractinian corals. Biol. Rev. 84, 1–17 (2009).PubMed 

Google Scholar 
J. T. Price, thesis, The Ohio State University (2020). More

Literature search and study selectionA systematic literature search was conducted in the ISI Web of Science database for observational and experimental studies published from 1975 to 13 January 2020 using the following search terms: TOPIC: (grassland* OR prairie* OR steppe* OR shrubland* OR scrubland* OR bushland*) AND TOPIC: (drought* OR ‘dry period*’ OR ‘dry condition*’ OR ‘dry year*’ OR ‘dry spell*’) AND TOPIC: (product* OR biomass OR cover OR abundance* OR phytomass). The search was refined to include the subject categories Ecology, Environmental Sciences, Plant Sciences, Biodiversity Conservation, Multidisciplinary Sciences and Biology, and the document types Article, Review and Letter. This yielded a total of 2,187 peer-reviewed papers (Supplementary Fig. 1). At first, these papers were screened by title and abstract, which resulted in 197 potentially relevant full-text articles. We then examined the full text of these papers for eligibility and selected 87 studies (43 experimental, 43 observational and 1 that included both types) on the basis of the following criteria:

(1)

The research was conducted in the field, in natural or semi-natural grasslands or shrublands (for example, artificially constructed (seeded or planted) plant communities or studies using monolith transplants were excluded). We used this restriction because most reports on observational droughts are from intact ecosystems, and experiments in disturbed sites or using artificial communities would thus not be comparable to observational drought studies.

(2)

In the case of observational studies, the drought year or a multi-year drought was clearly specified by the authors (that is, we did not arbitrarily extract dry years from a long-term dataset). Please note that some observational data points are from control plots of experiments (of any kind), where the authors reported that a drought had occurred during the study period. We did not involve gradient studies that compare sites of different climates, which are sometimes referred to as ‘observational studies’.

(3)

The paper reported the amount or proportion of change in annual or growing-season precipitation (GSP) compared with control conditions. We consistently use the term ‘control’ for normal precipitation (non-drought) year or years in observational studies and for ambient precipitation (no treatment) in experimental studies hereafter. Similarly, we use the term ‘drought’ for both drought year or years in observational studies and drought treatment in experimental studies. In the case of multi-factor experiments, where precipitation reduction was combined with any other treatment (for example, warming), data from the plots receiving drought only and data from the control plots were used.

(4)

The paper contained raw data on plant production under both control and drought conditions, expressed in any of the following variables: ANPP, aboveground plant biomass (in grassland studies only) or percentage plant cover. In 79% of the studies that used ANPP as a production variable, ANPP was estimated by harvesting peak or end-of-season AGB. We therefore did not distinguish between ANPP and AGB, which are referred to as ‘biomass’ hereafter. We included the papers that reported the production of the whole plant community, or at least that of the dominant species or functional groups approximating the abundance of the whole community.

(5)

When multiple papers were published on the same experiment or natural drought event at the same study site, the most long-term study including the largest number of drought years was chosen.

In addition to the systematic literature search, we included 27 studies (9 experimental, 17 observational and 1 that included both types) meeting the above criteria from the cited references of the Web of Science records selected for our meta-analyses, and from previous meta-analyses and reviews on the topic. In total, this resulted in 114 studies (52 experimental, 60 observational and 2 that included both types; Supplementary Note 9, Supplementary Fig. 2 and ref. 25).Data compilationData were extracted from the text or tables, or were read from the figures using Web Plot Digitizer26. For each study, we collected the study site, latitude, longitude, mean annual temperature (MAT) and precipitation (MAP), study type (experimental or observational), and drought length (the number of consecutive drought years). When MAT or MAP was not documented in the paper, it was extracted from another published study conducted at the same study site (identified by site names and geographic coordinates) or from an online climate database cited in the respective paper. We also collected vegetation type—that is, grassland when it was dominated by grasses, or shrubland when the dominant species included one or more shrub species (involving communities co-dominated by grasses and shrubs). Data from the same study (that is, paper) but from different geographic locations or environmental conditions (for example, soil types, land uses or multiple levels of experimental drought) were collected as distinct data points (but see ‘Statistical analysis’ for how these points were handled). As a result, the 114 published papers provided 239 data points (112 experimental and 127 observational)25.For the observational studies, normal precipitation year or years specified by the authors was used as the control. If it was not specified in the paper, the year immediately preceding the drought year(s) was chosen as the control. When no data from the pre-drought year were available, the year immediately following the drought year(s) (14 data points) or a multi-year period given in the paper (22 data points) was used as the control. For the experimental studies, we also collected treatment size (that is, rainout shelter area or, if it was not reported in the paper, the experimental plot size).For the calculation of drought severity, we used yearly precipitation (YP), which was reported in a much higher number of studies than GSP. We extracted YP for both control (YPcontrol) and drought (YPdrought). For the observational studies, when a multi-year period was used as the control or the natural drought lasted for more than one year, precipitation values were averaged across the control or drought years, respectively. Consistently, in the case of multi-year drought experiments, YPcontrol and YPdrought were averaged across the treatment years. When only GSP was published in the paper (63 of 239 data points), we used this to obtain YP data as follows: we regarded MAP as YPcontrol, and YPdrought was calculated as YPdrought = MAP − (GSPcontrol − GSPdrought). From YPcontrol and YPdrought data, we calculated drought severity as follows: (YPdrought − YPcontrol)/YPcontrol × 100.For production, we compiled the mean, replication (N) and, if the study reported it, a variance estimate (s.d., s.e.m. or 95% CI) for both control and drought. In the case of multi-year droughts, data only from the last drought year were extracted, except in five studies (17 data points) where production data were given as an average for the drought years. When both biomass and cover data were presented in the paper, we chose biomass. For each study, we consistently considered replication as the number of the smallest independent study unit. When only the range of replications was reported in a study, we chose the smallest number.To quantify climatic aridity for each study site, we used an aridity index (AI), calculated as the ratio of MAP and mean annual PET (AI = MAP/PET). This is a frequently used index in recent climate change research27,28. AI values were extracted from the Global Aridity Index and Potential Evapotranspiration (ET0) Climate Database v.2 for the period of 1970–2000 (aggregated on annual basis)29.Because we wanted to prevent our analysis from being distorted by a strongly unequal distribution of studies between the two study types regarding some potentially important explanatory variables, we left out studies from our focal meta-analysis in three steps. First, we left out studies that were conducted at wet sites—that is, where site AI exceeded 1. The value of 1 was chosen for two reasons: above this value, the distribution of studies between the two study types was extremely uneven (22 experimental versus 2 observational data points with AI  > 1)25, and the AI value of 1 is a bioclimatically meaningful threshold, where MAP equals PET. Second, we left out shrublands, because we had only 14 shrubland studies (out of 105 studies with AI  More

IPCC. IPCC Fifth Assessment Report (AR5). 10–12 (IPCC, 2014).Cailleret, M. et al. A synthesis of radial growth patterns preceding tree mortality. Glob. Chang. Biol. 23, 1675–1690 (2017).PubMed 

Google Scholar 
Forzieri, G. et al. Emergent vulnerability to climate-driven disturbances in European forests. Nat. Commun. 12, 1–12 (2021).
Google Scholar 
Bonan, G. B. Forests and climate change: forcings, feedbacks, and the climate benefits of forests. Science https://doi.org/10.1126/science.1155121 (2008).Article 
PubMed 

Google Scholar 
Buras, A. & Menzel, A. Projecting tree species composition changes of European forests for 2061–2090 under RCP 4.5 and RCP 8.5 scenarios. Front. Plant Sci. 9, 1–13 (2019).
Google Scholar 
van der Maaten, E. et al. Species distribution models predict temporal but not spatial variation in forest growth. Ecol. Evol. 7, 2585–2594 (2017).PubMed 
PubMed Central 

Google Scholar 
Lebaube, S., Le Goff, N. L., Ottorini, J. M. & Granier, A. Carbon balance and tree growth in a Fagus sylvatica stand. Ann. Sci. 57, 49–61 (2000).
Google Scholar 
Dobbertin, M. Tree growth as indicator of tree vitality and of tree reaction to environmental stress: a review. Eur. J. For. Res. 124, 319–333 (2005).
Google Scholar 
Büntgen, U. Re-thinking the boundaries of dendrochronology. Dendrochronologia 53, 1–4 (2019).
Google Scholar 
Klesse, S. et al. Continental-scale tree-ring-based projection of Douglas-fir growth: Testing the limits of space-for-time substitution. Glob. Chang. Biol. 26, 5146–5163 (2020).PubMed 

Google Scholar 
Zhao, S. et al. The International Tree-Ring Data Bank (ITRDB) revisited: data availability and global ecological representativity. J. Biogeogr. 46, 355–368 (2019).
Google Scholar 
Babst, F. et al. When tree rings go global: challenges and opportunities for retro- and prospective insight. Quat. Sci. Rev. 197, 1–20 (2018).
Google Scholar 
Klesse, S. et al. Sampling bias overestimates climate change impacts on forest growth in the southwestern United States. Nat. Commun. 9, 1–9 (2018).
Google Scholar 
Yousefpour, R. et al. Realizing mitigation efficiency of European commercial forests by climate smart forestry. Sci. Rep. 8, 1–11 (2018).CAS 

Google Scholar 
Giesecke, T., Hickler, T., Kunkel, T., Sykes, M. T. & Bradshaw, R. H. W. Towards an understanding of the Holocene distribution of Fagus sylvatica L. J. Biogeogr. 34, 118–131 (2007).
Google Scholar 
Fang, J. & Lechowicz, M. J. Climatic limits for the present distribution of beech (Fagus L.) species in the world. J. Biogeogr. 33, 1804–1819 (2006).
Google Scholar 
Luterbacher, J., Dietrich, D., Xoplaki, E., Grosjean, M. & Wanner, H. European seasonal and annual temperature variability, trends, and extremes since 1500. Science 303, 1499–1503 (2004).CAS 
PubMed 

Google Scholar 
Luterbacher, J. et al. European summer temperatures since Roman times. Environ. Res. Lett. 11, 24001 (2016).Nabuurs, G. J. et al. By 2050 the mitigation effects of EU forests could nearly double through climate smart forestry. Forests 8, 1–14 (2017).
Google Scholar 
Walentowski, H. et al. Assessing future suitability of tree species under climate change by multiple methods: a case study in southern Germany. Ann. Res. 60, 101–126 (2017).
Google Scholar 
Mäkelä, A. et al. Process-based models for forest ecosystem management: current state of the art and challenges for practical implementation. Tree Physiol. 20, 289–298 (2000).PubMed 

Google Scholar 
Leech, S. M., Almuedo, P. L. & Neill, G. O. Assisted migration: adapting forest management to a changing climate. BC J. Ecosyst. Manag. 12, 18–34 (2011).
Google Scholar 
Sass-Klaassen, U. G. W. et al. A tree-centered approach to assess impacts of extreme climatic events on forests. Front. Plant Sci. 7, 1069 (2016).PubMed 
PubMed Central 

Google Scholar 
Bowman, D. M. J. S., Brienen, R. J. W., Gloor, E., Phillips, O. L. & Prior, L. D. Detecting trends in tree growth: not so simple. Trends Plant Sci. 18, 11–17 (2013).CAS 
PubMed 

Google Scholar 
Hacket-Pain, A. J. et al. Climatically controlled reproduction drives interannual growth variability in a temperate tree species. Ecol. Lett. 21, 1833–1844 (2018).PubMed 
PubMed Central 

Google Scholar 
Dorji, Y., Annighöfer, P., Ammer, C. & Seidel, D. Response of beech (Fagus sylvatica L.) trees to competition-new insights from using fractal analysis. Remote Sens. 11, 2656 (2019).Petit-Cailleux, C. et al. Combining statistical and mechanistic models to unravel the drivers of mortality within a rear-edge beech population. bioRxiv https://doi.org/10.1101/645747 (2019).Weigel, R., Gilles, J., Klisz, M., Manthey, M. & Kreyling, J. Forest understory vegetation is more related to soil than to climate towards the cold distribution margin of European beech. J. Veg. Sci. 30, 746–755 (2019).
Google Scholar 
Etzold, S. et al. Nitrogen deposition is the most important environmental driver of growth of pure, even-aged and managed European forests. Forest Ecol. Manag. 458, 117762 (2020).
Google Scholar 
Martínez-Sancho, E. et al. The GenTree dendroecological collection, tree-ring and wood density data from seven tree species across Europe. Sci. Data 7, 1–7 (2020).PubMed 
PubMed Central 

Google Scholar 
Hartl-Meier, C., Dittmar, C., Zang, C. & Rothe, A. Mountain forest growth response to climate change in the Northern Limestone Alps. Trees 28, 819–829 (2014).
Google Scholar 
Way, D. A. & Montgomery, R. A. Photoperiod constraints on tree phenology, performance and migration in a warming world. Plant Cell Environ. 38, 1725–1736 (2015).PubMed 

Google Scholar 
Martínez del Castillo, E. et al. Spatial patterns of climate – growth relationships across species distribution as a forest management tool in Moncayo Natural Park (Spain). Eur. J. Res. 138, 299 (2019).
Google Scholar 
Hacket-Pain, A. J., Cavin, L., Friend, A. D. & Jump, A. S. Consistent limitation of growth by high temperature and low precipitation from range core to southern edge of European beech indicates widespread vulnerability to changing climate. Eur. J. Res. 135, 897–909 (2016).
Google Scholar 
van der Maaten, E. Climate sensitivity of radial growth in European beech (Fagus sylvatica L.) at different aspects in southwestern Germany. Trees 26, 777–788 (2012).
Google Scholar 
Decuyper, M. et al. Spatio-temporal assessment of beech growth in relation to climate extremes in Slovenia – an integrated approach using remote sensing and tree-ring data. Agric. Meteorol. 287, 107925 (2020).
Google Scholar 
Kraus, C., Zang, C. & Menzel, A. Elevational response in leaf and xylem phenology reveals different prolongation of growing period of common beech and Norway spruce under warming conditions in the Bavarian Alps. Eur. J. Res. 135, 1011–1023 (2016).
Google Scholar 
Martínez del Castillo, E. et al. Living on the edge: contrasted wood-formation dynamics in Fagus sylvatica and Pinus sylvestris under mediterranean conditions. Front. Plant Sci. 7, 370 (2016).PubMed 
PubMed Central 

Google Scholar 
Čufar, K. et al. Temporal shifts in leaf phenology of beech (Fagus sylvatica) depend on elevation. Trees 26, 1091–1100 (2012).
Google Scholar 
Bontemps, J. D., Hervé, J. C. & Dhôte, J. F. Dominant radial and height growth reveal comparable historical variations for common beech in north-eastern France. Forest Ecol. Manag. 259, 1455–1463 (2010).
Google Scholar 
Latte, N., Lebourgeois, F. & Claessens, H. Increased tree-growth synchronization of beech (Fagus sylvatica L.) in response to climate change in northwestern Europe. Dendrochronologia 33, 69–77 (2015).
Google Scholar 
Zimmermann, J., Hauck, M., Dulamsuren, C. & Leuschner, C. Climate warming-related growth decline affects Fagus sylvatica, but not other broad-leaved tree species in central european mixed forests. Ecosystems 18, 560–572 (2015).CAS 

Google Scholar 
Tegel, W. et al. A recent growth increase of European beech (Fagus sylvatica L.) at its Mediterranean distribution limit contradicts drought stress. Eur. J. Res. 133, 61–71 (2014).
Google Scholar 
Hacket-Pain, A. J. & Friend, A. D. Increased growth and reduced summer drought limitation at the southern limit of Fagus sylvatica L., despite regionally warmer and drier conditions. Dendrochronologia 44, 22–30 (2017).
Google Scholar 
Dulamsuren, C., Hauck, M., Kopp, G., Ruff, M. & Leuschner, C. European beech responds to climate change with growth decline at lower, and growth increase at higher elevations in the center of its distribution range (SW Germany). Trees 31, 673–686 (2017).
Google Scholar 
Spiecker, H., Mielikäinen, K., Köhl, M. & Skovsgaard, J. P. Growth trends in European forests: studies from 12 countries. European Forest Institute Research Report (1996).Cavin, L. & Jump, A. S. Highest drought sensitivity and lowest resistance to growth suppression are found in the range core of the tree Fagus sylvatica L. not the equatorial range edge. Glob. Chang. Biol. 23, 1–18 (2016).
Google Scholar 
Mette, T. et al. Climatic turning point for beech and oak under climate change in Central Europe. Ecosphere 4, 1–19 (2013).
Google Scholar 
Michelot, A., Simard, S., Rathgeber, C. B. K., Dufrêne, E. & Damesin, C. Comparing the intra-annual wood formation of three European species (Fagus sylvatica, Quercus petraea and Pinus sylvestris) as related to leaf phenology and non-structural carbohydrate dynamics. Tree Physiol. 32, 1033–1045 (2012).PubMed 

Google Scholar 
Meier, I. C. & Leuschner, C. Belowground drought response of European beech: Fine root biomass and carbon partitioning in 14 mature stands across a precipitation gradient. Glob. Chang. Biol. 14, 2081–2095 (2008).
Google Scholar 
Leuschner, C. & Ellenberg, H. Ecology of Central European Forests. Vegetation Ecology of Central Europe. Vol. I (Springer, 2017).Allen, C. D., Breshears, D. D. & McDowell, N. G. On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere. 6, 1–55 (2015).
Google Scholar 
Pechanec, V., Purkyt, J., Benc, A., Nwaogu, C. & Lenka, Š. Ecological Informatics Modelling of the carbon sequestration and its prediction under climate change. https://doi.org/10.1016/j.ecoinf.2017.08.006 (2017).Speer, J. H. Fundamentals of Tree-Ring Research (University of Arizona Press, 2010).Biondi, F. & Qeadan, F. A theory-driven approach to tree-ring standardization: defining the biological trend from expected basal area increment. Tree-Ring Res. 64, 81–96 (2008).
Google Scholar 
Biondi, F. & Qeadan, F. Removing the tree-ring width biological trend using expected basal area increment. in USDA Forest Service RMRS-P-55 124–131 (2008).Karger, D. N. et al. Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4, 1–20 (2017).
Google Scholar 
De Martonne, E. Une nouvelle fonction climatologique: L’indice d’aridité. La Meteorol. 2, 449–458 (1926).Martínez del Castillo, E., Longares, L. A., Serrano-Notivoli, R. & de Luis, M. Modeling tree-growth: assessing climate suitability of temperate forests growing in Moncayo Natural Park (Spain). Ecol. Manag. 435, 128–137 (2019).
Google Scholar 
Bolker, B. M. et al. Generalized linear mixed models: a practical guide for ecology and evolution. Trends Ecol. Evol. 24, 127–135 (2009).PubMed 

Google Scholar 
Calcagno, V. & Mazancourt, C. De. glmulti: an R package for easy automated model selection with (generalized) linear models. J. Stat. Softw. 34, 1–29 (2010).
Google Scholar 
Detry, M. A. & Ma, Y. Analyzing repeated measurements using mixed models. JAMA J. Am. Med. Assoc. 315, 407 (2016).CAS 

Google Scholar 
Harrison, X. A. et al. A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ 2018, 1–32 (2018).
Google Scholar 
Johnson, J. B. & Omland, K. S. Model selection in ecology and evolution. Trends Ecol. Evol. 19, 101–108 (2004).PubMed 

Google Scholar 
Caudullo, G., Welk, E. & San-Miguel-Ayanz, J. Chorological maps for the main European woody species. Data Brief 12, 662–666 (2017).PubMed 
PubMed Central 

Google Scholar 
Meinshausen, M. et al. The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500. Geosci. Model Dev. 13, 3571–3605 (2020).CAS 

Google Scholar 
Karger, D. N. & Zimmermann, N. E. CHELSAcruts – High Resolution Temperature And Precipitation Timeseries For The 20th Century And Beyond. https://doi.org/10.16904/envidat.159 (2018).Norinder, U., Rybacka, A. & Andersson, P. L. Conformal prediction to define applicability domain – a case study on predicting ER and AR binding. SAR QSAR Environ. Res. 27, 303–316 (2016).CAS 
PubMed 

Google Scholar 
Metzger, M. J., Bunce, R. G. H., Jongman, R. H. G., Mücher, C. A. & Watkins, J. W. A climatic stratification of the environment of Europe. Glob. Ecol. Biogeogr. 14, 549–563 (2005).
Google Scholar  More

IPCC. in Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (eds Masson-Delmotte, V. et al) (Cambridge University Press. In Press, 2021).Otto, F. E. L. et al. Toward an inventory of the impacts of human-induced climate change. Bull. Am. Meteorol. Soc. 101, E1972–E1979 (2020).
Google Scholar 
Stanners, D. et al. in Sustainability Indicators. A Scientific Assessment (eds Moldan, B., Hak, T. & Dahl, A. L.) 127–144 (Island Press, 2007).Cohen-Shacham, E. et al. Core principles for successfully implementing and upscaling Nature-based Solutions. Environ. Sci. Policy 98, 20–29 (2019).
Google Scholar 
Seddon, N. et al. Getting the message right on nature-based solutions to climate change. Glob. Chang. Biol. https://doi.org/10.1111/gcb.15513 (2021).Keesstra, S. et al. The superior effect of nature based solutions in land management for enhancing ecosystem services. Sci. Total Environ. 610-611, 997–1009 (2018).CAS 

Google Scholar 
Seddon, N. et al. Understanding the value and limits of nature-based solutions to climate change and other global challenges. Philos. Trans. R. Soc. Lond. B Biol. Sci. 375, 20190120 (2020).
Google Scholar 
Gómez Martín, E., Máñez Costa, M. & Schwerdtner Máñez, K. An operationalized classification of Nature Based Solutions for water-related hazards: from theory to practice. Ecol. Econ. 167 https://doi.org/10.1016/j.ecolecon.2019.106460 (2020).Doswald, N. et al. Effectiveness of ecosystem-based approaches for adaptation: review of the evidence-base. Clim. Dev. 6, 185–201 (2014).
Google Scholar 
Chausson, A. et al. Mapping the effectiveness of nature-based solutions for climate change adaptation. Glob. Chang. Biol. https://doi.org/10.1111/gcb.15310 (2020).Rebelo, A. J., Holden, P. B., Esler, K. & New, M. G. Benefits of water-related ecological infrastructure investments to support sustainable land-use: a review of evidence from critically water-stressed catchments in South Africa. R. Soc. Open Sci. 8, 201402 (2021).
Google Scholar 
Berrang-Ford, L. et al. A systematic global stocktake of evidence on human adaptation to climate change. Nat. Clim. Change 11, 989–1000 (2021).
Google Scholar 
Griscom, B. W. et al. Natural climate solutions. Proc. Natl Acad. Sci. USA 114, 11645–11650 (2017).CAS 

Google Scholar 
Bastin, J.-F. et al. The global tree restoration potential. Science 365, 76–79 (2019).CAS 

Google Scholar 
Koch, A., Brierley, C. & Lewis, S. L. Effects of Earth system feedbacks on the potential mitigation of large-scale tropical forest restoration. Biogeosciences 18, 2627–2647 (2021).CAS 

Google Scholar 
Girardin, C. A. J. et al. Nature-based solutions can help cool the planet – if we act now. Nature 593, 191–194 (2021).CAS 

Google Scholar 
Sudmeier-Rieux, K. et al. Scientific evidence for ecosystem-based disaster risk reduction. Nat. Sustain. 4, 803–810 (2021).
Google Scholar 
Otto, F. E. L. Attribution of weather and climate events. Annu. Rev. Environ. Resources 42, 627–646 (2017).
Google Scholar 
Philip, S. et al. A protocol for probabilistic extreme event attribution analyses. Adv. Stat. Climatol. Meteorol. Oceanogr. 6, 177–203 (2020).
Google Scholar 
Herring, S. C., Christidis, N., Hoell, A., Hoerling, M. P. & Stott, P. A. Explaining extreme events of 2019 from a climate perspective. Bull. Amer. Meteorol. Soc. 102, S1–S112 (2021).Otto, F. E. L. et al. Challenges to understanding extreme weather changes in lower income countries. Bull. Am. Meteorol. Soc. https://doi.org/10.1175/bams-d-19-0317.1 (2020).Pall, P. et al. Anthropogenic greenhouse gas contribution to flood risk in England and Wales in autumn 2000. Nature 470, 382–385 (2011).CAS 

Google Scholar 
Kay, A. L., Crooks, S. M., Pall, P. & Stone, D. A. Attribution of Autumn/Winter 2000 flood risk in England to anthropogenic climate change: a catchment-based study. J. Hydrol. 406, 97–112 (2011).
Google Scholar 
Schaller, N. et al. Human influence on climate in the 2014 southern England winter floods and their impacts. Nat. Clim. Change 6, 627–634 (2016).
Google Scholar 
Wolski, P., Stone, D., Tadross, M., Wehner, M. & Hewitson, B. Attribution of floods in the Okavango basin, Southern Africa. J. Hydrol. 511, 350–358 (2014).
Google Scholar 
Ross, A. C. et al. Anthropogenic influences on extreme annual streamflow into Chesapeake Bay from the Susquehanna River. Bull. Am. Meteorol. Soc. 102, S25–S32 (2021).Mitchell, D. et al. Attributing human mortality during extreme heat waves to anthropogenic climate change. Environ. Res. Lett. 11, 074006 (2016).
Google Scholar 
Botai, C., Botai, J., de Wit, J., Ncongwane, K. & Adeola, A. Drought Characteristics over the Western Cape Province, South Africa. Water https://doi.org/10.3390/w9110876 (2017).Wolski, P. How severe is Cape Town’s “Day Zero” drought? Significance 15, 24–27 (2018).
Google Scholar 
Stafford, L., Shemie, D., Kroeger, T., Baker, T. & Apse, C. The Greater Cape Town Water Fund. Assessing the return on investment for Ecological Infrastructure restoration. Business case. (The Nature Conservancy, 2018).Otto, F. E. L. et al. Anthropogenic influence on the drivers of the Western Cape drought 2015–2017. Environ. Res. Lett. https://doi.org/10.1088/1748-9326/aae9f9 (2018).Pascale, S., Kapnick, S. B., Delworth, T. L. & Cooke, W. F. Increasing risk of another Cape Town “Day Zero” drought in the 21st century. Proc. Natl. Acad. Sci. USA https://doi.org/10.1073/pnas.2009144117 (2020).Van Wilgen, B. W., Measey, J., Richardson, D. M., Wilson, J. R. & Zengeya, T. A. Biological Invasions in South Africa (Springer Nature, 2020).Le Maitre, D. et al. Impacts of plant invasions on terrestrial water flows in South Africa in Biological Invasions in South Africa (eds van Wilgen, B. W., Measey. J., Richardson, D. M., Wilson, J. R. & Zengeya, T. A.) 431–457 (Springer, 2020).Brown, A. E., Zhang, L., McMahon, T. A., Western, A. W. & Vertessy, R. A. A review of paired catchment studies for determining changes in water yield resulting from alterations in vegetation. J. Hydrol. 310, 28–61 (2005).
Google Scholar 
Dennedy-Frank, P. J. & Gorelick, S. M. Insights from watershed simulations around the world: watershed service-based restoration does not significantly enhance streamflow. Glob. Environ. Change https://doi.org/10.1016/j.gloenvcha.2019.101938 (2019).Calder, I. D. & Dye, P. Hydrological impacts of invasive alien plants. Land Use Water Resour. Res. 7, 1–12 (2001).
Google Scholar 
Trabucco, A., Zomer, R. J., Bossio, D. A., van Straaten, O. & Verchot, L. V. Climate change mitigation through afforestation/reforestation: a global analysis of hydrologic impacts with four case studies. Agric. Ecosyst. Environ. 126, 81–97 (2008).
Google Scholar 
Farley, K. A., Jobbagy, E. G. & Jackson, R. B. Effects of afforestation on water yield: a global synthesis with implications for policy. Glob. Change Biol. 11, 1565–1576 (2005).
Google Scholar 
Jackson, R. B. Trading water for carbon with biological carbon sequestration. Science 310, 1944–1947 (2005).Filoso, S., Bezerra, M. O., Weiss, K. C. B. & Palmer, M. A. Impacts of forest restoration on water yield: a systematic review. PLoS ONE 12, e0183210 (2017).
Google Scholar 
Sitzia, T., Campagnaro, T., Kowarik, I. & Trentanovi, G. Using forest management to control invasive alien species: helping implement the new European regulation on invasive alien species. Biol. Invasions 18, 1–7 (2015).
Google Scholar 
Richardson, D. M. & Rejmánek, M. Trees and shrubs as invasive alien species – a global review. Divers. Distrib. 17, 788–809 (2011).
Google Scholar 
Everard, M. et al. Can control of invasive vegetation improve water and rural livelihood security in Nepal? Ecosyst. Serv. 32, 125–133 (2018).
Google Scholar 
Everard, M. Can management of ‘thirsty’ alien trees improve water security in semi-arid India? Sci. Total Environ. 704, 135451 (2020).CAS 

Google Scholar 
Archer, S. R. et al. Woody plant encroachment: causes and consequences in Rangeland Systems Springer Series on Environmental Management (ed. Briske, D. D.) Chapter 2, 25–84 (2017).Wood, M. Bootstrapped confidence intervals as an approach to statistical inference. Organ. Res. Methods 8, 454–470 (2016).
Google Scholar 
Tan, S. H. The correct interpretation of confidence intervals. Proc. Singapore Healthc. 19 (2010).Coetsee, C., Gray, E. F., Wakeling, J., Wigley, B. J. & Bond, W. J. Low gains in ecosystem carbon with woody plant encroachment in a South African savanna. J. Trop. Ecol. 29, 49–60 (2012).
Google Scholar 
Stevens, N., Erasmus, B. F., Archibald, S. & Bond, W. J. Woody encroachment over 70 years in South African savannahs: overgrazing, global change or extinction aftershock? Philos. Trans. R. Soc. Lond. B Biol. Sci. https://doi.org/10.1098/rstb.2015.0437 (2016).Venter, Z. S., Cramer, M. D. & Hawkins, H. J. Drivers of woody plant encroachment over Africa. Nat. Commun. 9, 2272 (2018).CAS 

Google Scholar 
Forsyth, G. G., Le Maitre, D. C., Smith, J. & Lotter, D. Upper Berg River Catchment (G10A) Management Unit Control Plan. (Natural Resources Management (NRM) Department of Environmental Affairs, 2016).Dirmeyer, P. A., Balsamo, G., Blyth, E. M., Morrison, R. & Cooper, H. M. Land‐atmosphere interactions exacerbated the drought and heatwave over northern Europe during summer 2018. AGU Adv. 2, e2020AV000283 (2021).
Google Scholar 
Rejmánek, M., Richardson, D. M. & Pysek, P. Trees and shrubs as invasive alien species – 2013 update of the global database. Divers. Distrib. 19, 1093–1094 (2013).
Google Scholar 
Terrer, C. et al. A trade-off between plant and soil carbon storage under elevated CO2. Nature 591, 599–603 (2021).CAS 

Google Scholar 
Ziervogel, G. et al. Climate change impacts and adaptation in South Africa. Wiley Interdiscip. Rev. Clim. Change 5, 605–620 (2014).
Google Scholar 
Thomas, A. et al. Global evidence of constraints and limits to human adaptation. Reg. Environ. Change https://doi.org/10.1007/s10113-021-01808-9 (2021).Dow, K., Berkhout, F. & Preston, B. L. Limits to adaptation to climate change: a risk approach. Curr. Opin. Environ. Sustain. 5, 384–391 (2013).
Google Scholar 
Manning, J. & Goldblatt, P. Plants of the greater Cape Floristic Region 1: the Core Cape Flora., (South African National Biodiversity Institute, 2012).Nel, J. L. et al. Strategic water source areas for urban water security: Making the connection between protecting ecosystems and benefiting from their services. Ecosyst. Serv. 28, 251–259 (2017).
Google Scholar 
Wolski, P. What Cape Town learned from its drought. Bull. At. Sci. https://thebulletin.org/2018/04/what-cape-town-learned-from-its-drought/ (2018).D. W. S. Cape Town River Systems State of Dams on 2021-08-16. Department of Water and Sanitation. Republic of South Africa. https://www.dws.gov.za/Hydrology/Weekly/RiverSystems.aspx?river=CT (2021).Rebelo, A. J. et al. The hydrological benefits of restoration: a modelling study of alien tree clearing in four mountain catchments in South Africa. Preprint at J. Hydrol. https://doi.org/10.21203/rs.3.rs-1316834/v1.DWAF. The Assessment of Water Availability in the Berg Catchment (WMA 19) by Means of Water Resource Related Models: Report 9 (Groundwater Model): Volume 9 (Breede River Alluvium Aquifer Model). (Department of Water Affairs and Forestry, 2008).DWAF. The Assessment of Water Availability in the Berg Catchment (WMA 19) by Means of Water Resource Related Models: Report 9 (Groundwater Model): Volume 3 (Regional Conceptual Model). (Department of Water Affairs and Forestry, 2008).Blake, D., Mlisa, A. & Hartnady, C. Large scale quantification of aquifer storage and volumes from the Peninsula and Skurweberg Formations in the southwestern Cape. Water SA 36, 177–184 (2010).
Google Scholar 
Holden, P. B., Rebelo, A. J. & New, M. G. Mapping invasive alien trees in water towers: a combined approach using satellite data fusion, drone technology and expert engagement. Remote Sens. Appl.: Soc. Environ. https://doi.org/10.1016/j.rsase.2020.100448 (2021).Midgley, J. & Scott, D. The use of stable isotopes of water in hydrological studies in the Jonkershoek Valley. Water SA 20, 151–154 (1994).
Google Scholar 
Van Genuchten, M. T. A closed‐form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898 (1980).
Google Scholar 
Moriasi, D. N., Gitau, M. W., Pai, N. & Daggupati, P. Hydrologic and water quality models: performance measures and evaluation criteria. Trans. ASABE 58, 1763–1785 (2015).
Google Scholar 
Stone, D. A. et al. A basis set for exploration of sensitivity to prescribed ocean conditions for estimating human contributions to extreme weather in CAM5.1-1degree. Weather Clim. Extremes 19, 10–19 (2018).
Google Scholar 
Risser, M. D., Stone, D. A., Paciorek, C. J., Wehner, M. F. & Angélil, O. Quantifying the effect of interannual ocean variability on the attribution of extreme climate events to human influence. Clim. Dyn. 49, 3051–3073 (2017).
Google Scholar 
Jones, G. S., Stott, P. A. & Christidis, N. Attribution of observed historical near-surface temperature variations to anthropogenic and natural causes using CMIP5 simulations. J. Geophys. Res. Atmos. 118, 4001–4024 (2013).
Google Scholar 
Sun, L. et al. Drivers of 2016 record Arctic warmth assessed using climate simulations subjected to factual and counterfactual forcing. Weather Clim. Extremes 19, 1–9 (2018).CAS 

Google Scholar 
Guillod, B. P. et al. weather@home 2: validation of an improved global–regional climate modelling system. Geosci. Model Dev. 10, 1849–1872 (2017).
Google Scholar 
Massey, N. et al. weather@home—development and validation of a very large ensemble modelling system for probabilistic event attribution. Q. J. R. Meteorol. Soc. 141, 1528–1545 (2014).
Google Scholar 
Taylor, K. E., Stouffer, R. J. & Meehl, G. A. An overview of CMIP5 and the experiment design. Bull. Am. Meteorol. Soc. 93, 485–498 (2012).
Google Scholar 
Flato, G. et al. in Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change 741–866 (Cambridge University Press, 2014).Hargreaves, G. H. & Samani, Z. A. Reference crop evapotranspiration from temperature. Appl. Eng. Agriculture 1, 96–99 (1985).
Google Scholar 
Cayan, D. R., Maurer, E. P., Dettinger, M. D., Tyree, M. & Hayhoe, K. Climate change scenarios for the California region. Clim. Change 87, 21–42 (2008).
Google Scholar 
Cannon, A. J., Sobie, S. R. & Murdock, T. Q. Bias correction of GCM precipitation by quantile mapping: how well do methods preserve changes in quantiles and extremes? J. Clim. 28, 6938–6959 (2015).
Google Scholar 
R Core Team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing. https://www.R-project.org/, 2020).Paciorek, C. J., Stone, D. A. & Wehner, M. F. Quantifying statistical uncertainty in the attribution of human influence on severe weather. Weather Clim. Extremes 20, 69–80 (2018).
Google Scholar 
Tadono, T. et al. Generation of the 30 M-Mesh Global Digital Surface Model by Alos Prism. ISPRS – International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences XLI-B4, 157–162, https://developers.google.com/earth-engine/datasets/catalog/JAXA_ALOS_AW3D30_V3_2#description (2016).
Google Scholar 
Takaku, J., Tadono, T., Tsutsui, K. & Ichikawa, M. Validation of “Aw3d” Global Dsm Generated from Alos Prism. ISPRS Ann. Photogramm. III-4, 25–31 (2016).
Google Scholar 
Viviroli, D. Increasing dependence of lowland population on mountain water resources. Nat. Sustain. 3, 917–928 (2020).
Google Scholar 
Meybeck, M. A New typology for mountains and other relief classes: an application to global continental water resources and population distribution. Mt. Res. Dev. 21, 34–45 (2001).DWS. Surface water home. Department of Water and Sanitation. Republic of South Africa. https://www.dws.gov.za/Hydrology/Unverified/UnverifiedDataFlowInfo.aspx (2021). More

Characteristics of heavy metal concentrationsOn the basis of the soil sample collection and chemical analysis, the concentration data for heavy metals in the urban soils of Wuxi were obtained. Through the statistical analysis of the soil heavy metal concentration data (Table 1), on the whole, the concentration of each heavy metal is as follows: Mn  > Zn  > Cr  > Ni  > Pb  > Cu  > Co  > Be  > Cd. Among these, the concentration range of Cr was 64.5–99 mg kg-1, and the average concentration was 72.9 mg kg−1. The concentration range of Ni was 31.4–67.5 mg kg−1, and the average concentration was 38.2 mg kg−1. The concentration range of Cu was 19.8–37.2 mg kg−1, and the average concentration was 25.5 mg kg−1. The concentration range of Zn was 72.4–1146 mg kg−1, and the average concentration was 90.2 mg kg−1. The concentration range of Cd was 0.34–1.06 mg kg−1, and the average concentration was 0.51 mg kg−1. The concentration range of Pb was 25.6–66.4 mg kg−1, and the average concentration was 37.6 mg kg−1. The variation coefficients of urban soil heavy metal concentration in Wuxi is between 0.09 and 0.33, which is less than 1. The spatial fluctuation of urban soil heavy metal concentration in Wuxi is small, indicating that the sources may be the same or similar.Table 1 Statistics of the heavy metal concentrations and Pb isotope ratios in the urban soils of Wuxi city (unit of heavy metal: mg kg−1; CV: coefficient of variation).Full size tableBy analysing the spatial distributions of the urban soil heavy metal concentrations in Wuxi, several obvious spatial distribution characteristics are found (Fig. 2). First, the heavy metals have high values in the central area of Wuxi, due to where has a high population density and various industries. The central aggregation of Pb is more obvious. Due to the dense roads in the city centre, vehicle traffic, bus stop signs and gas stations are mostly concentrated here, which will lead to Pb contents in this area that are significantly higher than those in other areas. In addition to the heavy metal concentrations, such as those for Cu, Zn and Cr in the downtown area, there are also areas with high values in western Wuxi and low values in eastern Wuxi. This phenomenon may be related to the land use types in Wuxi. In the western area of Wuxi, most land use types are urban and construction land, and the soils in this area are greatly disturbed by human activities. In the eastern region of Wuxi, woodland and grassland account for a large proportion of the land use types, which are less disturbed by human activities.Figure 2Spatial distribution characteristics of heavy metals in the urban soils of Wuxi city (unit: mg kg−1) [the figure was generated by Yan Li using the ArcGIS 10.2 (http:// https://developers.arcgis.com/)].Full size imageSource analysis of heavy metalsExploring for heavy metal pollution from emission sources is an important prerequisite for the study of urban soil pollution. By analysing the sources of heavy metals in soil environments, we can accurately determine which industries are major sources28,29,30 and whether there is homologous pollution. This is not only a theoretical basis for the study of lake sediment pollution and to clarify the risks brought by different pollution sources to the urban soil environment but also provides important guides for local government control of specific polluting industries and pollutant emissions. Based on this, the correlations and significance of heavy metals in the urban soils of Wuxi were analysed (Table 2). Generally, a heavy metal pollution source will emit multiple heavy metals at the same time. If the pollution source has a large emission, the concentration of these heavy metals in the environment will show a high level; on the contrary, if the emission of this pollution source is small, the concentration of these heavy metals in the environment will show a low level10. The correlations between the heavy metals Zn, Cr, Ni, Pb, Cu and Cd are between 0.655–0.907 and show strong correlations and significance at a level of 0.01. The strong significant correlations between different heavy metals indicate that these heavy metals have similar emission sources and transmission routes, which also means that they have consistent sources.Table 2 Correlations of Heavy Metals in the Urban Soils of Wuxi City.Full size tableTo further determine which industries are the sources of the heavy metals found in the urban soil of Wuxi, we analysed the Pb isotope data. The variation range of 208Pb/206Pb in soil is 2.09–2.12, and the average value is 2.10. The variation range of 206Pb/207Pb in soil is 1.17–1.18, and the average value is 1.177 (Table 1). After consulting relevant literature and materials, the main pollution sources of heavy metals in cities in eastern China include coal combustion, oil combustion, factory emissions, municipal wastes and so on3. Therefore, we collected the corresponding Pb isotope data in the emissions of heavy metal pollution sources. By collecting and comparatively analysing the Pb isotope data of known pollution sources (Fig. 3), it was determined that the Pb isotopes of the urban soil heavy metals in the soils of Wuxi have distinct characteristics. First, the Pb isotope distributions in the soils of Wuxi are relatively concentrated, and the ranges of variation are relatively small, which indicate that these heavy metals may have the same source or similar sets of sources. Second, the Pb isotopes in the urban soils of Wuxi city have few similarities with those of the uncontaminated soils and granites in eastern China; in contrast, the Pb isotopes in the urban soils of Wuxi are distributed in areas that are associated with coal combustion, automobile exhaust and urban waste (supplementary materials). The urban soil heavy metals in Wuxi generally have similar pollution sources and are greatly affected by human activities such as coal combustion and automobile exhaust emissions. Wuxi has a developed industrial economy and large numbers of factories. In the production and processing activities, the combustion of energy and fuel and the incomplete utilization of raw materials will lead to the enrichment of pollutants in the surrounding environment. By comparing other studies30,31, the Pb isotope analysis results in this study well indicate the source of soil heavy metals in Wuxi and make up for the Pb isotope data in this area. In the process of urban development, we should develop and apply clean energy, reduce the utilization of petroleum fossil fuels, and control the enrichment of heavy metals and other pollutants in the soil from the source.Figure 3Comparison of the Pb isotope compositions in the urban soils of Wuxi city with known sources.Full size imageEcological risk analysisBy calculating the potential ecological risk index for the heavy metals in the urban soils of Wuxi, the risks of heavy metals in the Wuxi soils were evaluated (Table 3). According to previous studies21, an Ei value lower than 40 indicates that a heavy metal is in a low-risk state at this location, and Ei values greater than or equal to 40 indicate that a heavy metal represents a high-risk state at this location. The average value of the potential ecological risk index of soil heavy metal Cd in Wuxi is 80.3, which represents a high-risk state. The average distributions of the potential ecological risk indexes of the heavy metals Cr, Cu, Zn, Pb and Ni are 1.8, 4.3, 1.1, 5.5 and 4.8, respectively, which all indicate a low-risk state. The risk statuses of different heavy metals may show certain correlations in space, which may be mutually complementary or antagonistic. Examining the spatial interactions of different heavy metal compound pollutants in urban soils plays an important role in the prevention and control of urban heavy metal pollution. Based on this, we used the Lisa analysis method to explore the spatial correlations of the different heavy metal risks in the urban soils of Wuxi (Fig. 4). The Moran scatter diagram can be divided into four quadrants that correspond to four different spatial patterns. High means that the variable value is higher than the average value, and Low means that the variable value is lower than the average value. In the upper right quadrant (High–High), a high-value area is surrounded by high-value neighbours; in the upper left quadrant (Low–High), a low-value area is surrounded by high-value neighbours; in the lower left quadrant (LL), a low-value area is surrounded by low-value neighbours; and in the lower right quadrant (High–Low), a high-value area is surrounded by low-value neighbours. High-High and Low-Low indicate that the differences between the region and its surrounding areas are small; that is, the regions with higher or lower values are concentrated, while the Low–High and High–Low quadrants indicate that the variable values between a region and its surrounding areas are different to a certain extent.Table 3 Ecological risk and health risk analysis of heavy metals in the urban soils of Wuxi (Cr-E represents the ecological risk of metal element Cr; Ni-E represents the ecological risk of metal element Ni; Cu-E represents the ecological risk of metal element Cu; Zn-E represents the ecological risk of metal element Zn; Cd-E represents the ecological risk of metal element Cd; Pb-E represents the ecological risk of metal element Pb; ADDderm-C is the average exposure to skin contact pathways for child; ADDderm-A is the average exposure to skin contact pathways for adult; ADDing-C is the average daily exposure to intake pathway for child; ADDing-A is the average daily exposure to intake pathway for adult; HI-C is the total health risk caused by accumulation of heavy metals in multiple ways in the same environmental medium for child; HI-A is the total health risk caused by accumulation of heavy metals in multiple ways in the same environmental medium for adult).Full size tableFigure 4LISA analysis of the ecological risks from different heavy metals [the figure was generated by Yan Li using the ArcGIS 10.2 (http:// https://developers.arcgis.com/)].Full size imageIn this study, two main results were obtained from spatial correlation Lisa analysis between different heavy metals. One is a High-High area, which is mainly distributed in the central and western regions of Wuxi city, which is consistent with the spatial distribution of the urban soil heavy metal concentrations in Wuxi city and is strongly disturbed by human activities. The other is the insignificant area, in which there are also large numbers of factories and enterprises and in which the forestland and grassland are distributed at intervals, which leads to an insignificant spatial correlation of soil heavy metal contents. Based on the above analysis, the high-risk areas for heavy metals in the urban soils of Wuxi are mainly concentrated in the central and western regions, and the relevant management activities need to be given great attention. In the eastern region, sporadic high-risk areas are also present, which should also receive due attention. Moran’s I is a method to measure the interdependence and degree of objects or phenomena by constructing statistics on certain characteristics or attributes for a certain spatial unit in the study area and the surrounding spatial units. It can be used to describe the spatial characteristics of spatial units such as aggregation or outliers in the distribution of certain attributes and is a very important technology in spatial data analysis33,34. However, few studies have applied it to the spatial relationship analysis of different heavy metals in urban soil.Health risk analysisBy using the health risk assessment model that is recommended by the U.S. EPA, this study calculated the health risks of soil heavy metals to adults and children through skin contact and ingestion. For both adults and children, the risk of soil heavy metals through ingestion was much higher than that caused by skin exposure (Table 3). For children, the total health risk that was caused by soil heavy metals is 0.078, which is four times that of adults. This may be related to children’s habits. Most children like to play with sand and climb around on the ground. These behaviours greatly increase the frequency of children contacting the soil, which thus increases the health risk caused by heavy metals in the soil. To further explore the spatial characteristics of the health risks of heavy metals in the soils of Wuxi, this study provides spatial predictions of the health risk values of soil heavy metals (Fig. 5). The total health risk values of soil heavy metals for children and adults have similar spatial distribution characteristics. High health risk values appear in the central area of Wuxi and decrease in a ring-shaped pattern. This is similar to the development degree of the city. The downtown area of Wuxi is densely populated, the pedestrian flow is very large, and the health risk of soil heavy metals in this area is very high, which poses a very serious potential threat. The health risk values for the western region of Wuxi are high, and there is also a potential threat. When compared with western Wuxi, eastern Wuxi has a lower risk.Figure 5Health risk analysis of heavy metals in the urban soils of Wuxi [the figure was generated by Yan Li using the ArcGIS 10.2 (http:// https://developers.arcgis.com/)].Full size image More

General characteristics of study soilsTable 2 presents the descriptive statistics regarding the soil characteristics. Significant changes were observed in the distribution of sand (110–850 g kg−1), silt (50–530 g kg−1), clay (100–610 g kg−1), and soil textural class (7 texture classes) showing the diversity of natural and human processes involved in the formation and development of these soils28. Almost all soil samples were alkaline (with reaction at a range of 7.4–8.1) and calcareous (with CCE at a range of 5.5–35%). The EC of some soils was  > 4 dS/m (about 7% of the soil samples), indicating the partial salinity of the study soils. The organic carbon and total N contents of the soils were, on average, 2% (0.8–3.1%) and 0.28% (0.05–0.51%), respectively, placing them within the range of the moderate class. Likewise, the mean CEC of the soil, which is an effective indicator of soil fertility and quality, was in the moderate class of 12–25 cmol kg−129. The CEC was found to be highly correlated with clay (r = 0.76 P  Pb (58 mg kg−1)  > Ni (55.4 mg kg−1)  > Cu (38.8 mg kg−1)  > Cd (0.88 mg kg−1). In most soil samples, these ranges are comparable with data reported for other urban soils around the world—e.g. Ref.30 in Poland, Ref.31 in China, and Ref.32 in Greece. The values of Cd, Cu, and Zn were below their acceptable ranges as per the international standards4 in all soil samples. Nonetheless, the Pb and Ni contents were higher than their acceptable ranges in 13.1% and 17.4% of the samples, respectively. Furthermore, the concentrations of the five elements were higher than their background values in all urban soil samples. This difference was considerable for Cd, Pb, and Ni. The heavy metals had CV in the order of Cd (53%)  > Pb (51%)  > Ni (46%)  > Zn (21%)  > Cu (18%). This CV variation implies great variations in Cd, Pb, and Ni, which is linked to anthropogenic activities33. The background values of the metals, estimated by the median absolute deviation method10,14, were 52.3, 18.7, 0.45, 29.1, and 30.8 mg kg−1 for Zn, Cu, Cd, Pb, and Ni, respectively.We compared the concentrations of the heavy metals between urban and non-urban soils and found significant increases in the concentration of the metals in most soil types (Fig. 2). The urban soils had 17–36%, 14–21%, 41–70%, 43–69%, and 13–24% higher Zn, Cu, Cd, Pb, and Ni contents than the non-urban soils. The effluent and waste entry from multiple food processing and storage units, dying plants, metal plating facilities, and plastic production in close proximity of the study area is believed to be the reason for the high concentration of these trace elements. Research in various parts of the world, e.g., Ref.34 in India, Ref.35 in Brazil, and Ref.36 in China, has documented that the facilities have introduced significant quantities of heavy metals to soils. However, traffic and agrochemicals also play a key role in the accumulation of heavy metals in this region10.Figure 2The comparison of the mean values of Zn (a), Cu (b), Cd (c), Pb (d), and Ni (e) between urban and non-urban soils in different soil types. Different letters indicate significant differences in metal content within each soil type at P  Ni  > Cu. These findings are comparable to the results reported by37 and12. The highest EF for all five elements was observed in the Fluvisols soil type, reflecting that this soil type had been exposed to element pollution induced by urban activities to a greater extent than the other soil types. In a study on the pollution potential of four soil types in Central Greece, Ref.38 reported different ranges of element pollution across different soil types.Figure 3The comparison of the mean enrichment factor of Zn (a), Cu (b), Cd (c), Pb (d), and Ni (e) between urban and non-urban soils in different soil types. Different letters indicate significant differences in enrichment factor within each soil type at P  Pb (1.89)  > Ni (1.86)  > Cu (1.73)  > Zn (1.51). Mean PI for non-urban soils followed the order Cd (1.5)  > Zn (1.4)  > Cu (1.33)  > Pb (1.31)  > Ni (1.29). Nearly 7% and 16% of the urban soils showed moderate pollution (MP, PI = 2–3) and high pollution classes (HP, PI  > 3) of PI for Cd and 39% and 4% showed the MP and HP class of PI for Pb, respectively. However, the PI class was low pollution (PI = 1–2) for all soil samples and soil types in the non-urban soils. The results on the pollution index indicate a widespread intensification of soil pollution in urban soils across all studied heavy metals.Table 3 The level and terminology of PI and Ei of the analyzed heavy metals in urban and non-urban soils.Full size tableEcological risk, Ei was similarly found to be significantly higher in the urban soils than in the non-urban soils, even though the concentration of all elements except Cd fell within the low-risk class (Ei ≤ 40) in both urban and non-urban soils (Table 3). The mean Ei for Cd was 58.7 (moderate-risk class) and 39.2 (low-risk class) in the urban and non-urban soils, respectively. This means that urban activities have enhanced the ecological risk class of Cd by one grade. Overall, Cd had the highest EF, PI, and Ei among all heavy metals and in all soil samples, indicating a greater risk potential by Cd than Zn, Cu, Pb, and Ni across the water-soil–plant-human domain. Elevated Cd pollution by anthropogenic activities has been widely reported in the literature10,12,39. Cadmium as a Group 1 carcinogen element40 can accumulate in plant tissue without exhibiting visual symptoms. Therefore, Cd generally transfers from soil to the food chain covertly. Cadmium pollution can also influence soil quality and reduce crop yields and grain quality3.Similar to EF, PI, and Ei, the mean ER was significantly elevated in all urban soil types than the non-urban soils (Fig. 4). Among different soil types, the ER magnitude was in the order of Fluvisols (66.6%)  > Regosols (66.1%)  > Cambisols (59.8%)  > Calcisols (47%). These results indicate that Fluvisols carry a higher ecological risk potential for heavy metal accumulations than other soil types. In the study region, Fluvisols due to higher fertility and productivity are subject to more intense and extensive agronomic operations than other soil types13. Heavy application of agrochemicals (e.g., pesticides, herbicides, insecticides, and chemical fertilizers), accelerate the heavy metal input to the Fluvisols. Widespread application of nitrogen fertilizers and subsequent reduction in average soil pH markedly increases the solubility of certain heavy metals (e.g., Zn, Cu, Cd) which can be another factor increasing the ecological risk of heavy metal contamination in Fluvisols41. In addition, these Fluvisols are located on the margin of open urban wastewater channels, which are sometimes used for irrigation. A combination of mentioned processes can be implicated for higher ER of Fluvisols than that of other soil types as for BF, PI, and Ei.Figure 4The comparison of the mean ecological risk of selected heavy metals between urban and non-urban soils in different soil types. Different letters indicate significant differences in ecological risk within each soil type at P  Cu  > Ni  > Cd  > Pb in the roots, partially differing from that of the grain—Zn  > Cu  > Pb  > Ni  > Cd. Heavy metals concentrations observed in the corn roots and grains are almost comparable with those reported by42 in China and43 in Peru.Table 4 Summary statistical attributes of the concentration of heavy metals in corn root (R) and grain (G) along with their BCF and TF.Full size tableThe accumulation of heavy metals in the edible parts of corn is of higher importance. In the present study, the concentrations of these metals were lower than the acceptable level in the corn grains based on international references44. So, the consumption of corns grown in the regions should not threaten human and animal health in the short term, but caution should be exercised in their long-term consumption because some of these elements, especially Cd and Pb, which have long decomposition half-lives, gradually accumulate in body organs, especially in kidneys and livers45. Besides, the ratio of Zn, Cu, Cd, Pb, and Ni of the corn grain to their acceptable standard concentration, known as the pollution index of crop heavy metals, Ref.12 was lower than 0.7 for most corn samples, indicating the unpolluted risk class.The mean concentrations of Cd, Pb, and Ni were 5, 3.1, and 9.2 times as great in the corn roots as in their grains. This observation exhibits a notable phytoremediatory function of corn roots through restriction of radial translocation of heavy metals to the xylems and eventually into the grains. A similar trend of heavy metal accumulation in different plant organs has been reported in previous observations46,47. Based on Kabata-Pendias4 and Adriano22, plant cells can use the defensive tools of the roots to cope with heavy metals, especially Cd and Pb—highly toxic metals to plant cytosols. Accordingly, plant cells can fix these elements in the root system by such approaches as precipitating on cell walls, storing in vacuoles, and/or chelating by phytochelatins, thereby alleviating their toxic effects and inhibiting their translocation to plant shoots. For Zn, Cu, and Cd metals, a significant correlation was observed between their concentration in corn roots and grains. But, a less significant correlation (P  Cu (0.17)  > Zn (0.12)  > Ni (0.02)  > Pb (0.01). This implies that Cd, and to a smaller extent Cu is taken up by corn roots from the soil more readily, but Pb and Ni are less absorbable. These results are consistent with the reports of48 and46. The greater value of BCF-Cd may be related to a combination of the specific factors e,g., Cd concentration and chemistry, as well as soil characteristics (e.g., soil texture, pH, and calcium carbonate content)4. As was already discussed, the examined soils were characterized by high alkaline (pH = 7.4–8.1) and calcareous properties (CCE = 5.5–35%) with a high concentration of Soluble salts (EC = 0.7–6.6 dS m−1). These characteristics can result in the formation of complex Cd ions, especially CdOH+, CdCl20, CdCl+, CdSO40, and CdHCO3+4,22. These ions are plant-available, resulting in a further increase in Cd BCF. Regarding Ni and Pb, the alkaline and calcareous properties of the soils may have motivated insoluble compounds such as NiHCO3+ and NiCO30 (for Ni) and Pb(OH)2, PbCO3, PbSO4, and PbO (for Pb)4,22. These compounds cannot be uptake by plant roots, which may have resulted in a significant decrease in the BCF of these metals versus the other analyzed elements.Like BCF, the heavy metals had TF of  Pb (0.21)  > Cd (0.2)  > Ni (0.15). This implies that Zn and Cu are translocated from roots to grains readily, about four times as great as the other metals, while Ni, Cd, and Pb are translocated in smaller concentrations.The comparison of BCF and TF of Cd showed that less than 30% of Cd, on average, accumulated in the corn roots were translocated to the grains. This states that Cd is immobilized by various mechanisms before it can find its way into the grains. Some of the important mechanisms include (i) the antagonistic effects of Cd with other equivalent elements, especially Zn, Fe, and Ca, in the vascular system of corn, which reduces its mobility in the corn root-stem-grain system22, (ii) Cd sequestration in active exchange sites on the cell wall in the corn root-stem pathway10, and (iii) the binding of Cd with some specific compounds, e.g., phytochelatins of root vacuoles, which immobilizes it before its translocation to grains4,22. Lin and Aarts52 remarked that Cd mostly tends to be trapped in root vacuoles, which reduces its translocation to the upper parts of the plants. In general, it was found that corn plants have a high potential to absorb and accumulate Cd in their roots and Zn in their grains, which is consistent with previous studies41. For the majority of heavy metals, the values of BCF and TF in different soil types were in the order of Fluvisols  > Regosols  > Cambisols  > Calcisols, indicating that the great variety of soil types for the uptake and translocation of heavy metals in the soil-root-grain of the corn (Fig. 5).Figure 5Effect of soil type on the mean bioconcentration factor (a) and translocation factor (b) of selected heavy metals in urban soils. Different letters indicate significant differences in bioconcentration and translocation factors among soil types for each metal at P  Zn  > Cu  > Pb  > Ni for children, differing from that for adults (Cu  > Cd  > Pb  > Zn  > Ni). The values of HQ was  1 in over 87% of the samples, implying the low non-carcinogenic risk of this metal for corn-consuming children in the study region53. Rapidly developing children’s nervous system are highly sensitive to environmental factors, including heavy metals, so even a relatively low concentration of Cd in children’s blood may irreversibly affect their mental growth and functioning54.The highest HI was observed in children (min = 1.16, max = 2.31, mean = 1.63) followed by women and men which was similar to the found pattern of HQ (Table 7). These data show a moderate non-carcinogenic health risk (1 ≤ HI  More

Foelix, R. F. Biology of Spiders (Oxford University Press, 2011).
Google Scholar 
World Spider Catalog. World Spider Catalog, Version 23.0. Natural History Museum Bern, online at http://wsc.nmbe.ch (2022).Nyffeler, M. & Sunderland, K. D. Composition, abundance and pest control potential of spider communities in agroecosystems: A comparison of European and US studies. Agric. Ecosyst. Environ. 95, 579–612 (2003).
Google Scholar 
Oldrati, V. et al. Peptidomic and transcriptomic profiling of four distinct spider venoms. PLoS ONE 12, e0172966 (2017).PubMed 
PubMed Central 

Google Scholar 
Herzig, V. et al. Animal toxins—Nature’s evolutionary-refined toolkit for basic research and drug discovery. Biochem. Pharmacol. 181, 114096 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Vollrath, F. & Knight, D. P. Liquid crystalline spinning of spider silk. Nature 410, 541–548 (2001).ADS 
CAS 
PubMed 

Google Scholar 
Moradmand, M. & Jäger, P. Taxonomic revision of the huntsman spider genus Eusparassus Simon, 1903 (Araneae: Sparassidae) in Eurasia. J. Nat. Hist. 46, 2439–2496 (2012).
Google Scholar 
Moradmand, M. The stone huntsman spider genus Eusparassus (Araneae: Sparassidae): Systematics and zoogeography with revision of the African and Arabian species. Zootaxa 3675, 1–108 (2013).PubMed 

Google Scholar 
Levy, G. The family of huntsman spiders in Israel with annotations on species of the Middle East (Araneae: Sparassidae). J. Zool. 217, 127–176 (1989).
Google Scholar 
Dunlop, J. A. et al. Computed tomography recovers data from historical amber: An example from huntsman spiders. Naturwissenschaften 98, 519–527 (2011).ADS 
CAS 
PubMed 

Google Scholar 
Moradmand, M., Schönhofer, A. L. & Jäger, P. Molecular phylogeny of the spider family Sparassidae with focus on the genus Eusparassus and notes on the RTA-clade and ‘Laterigradae’. Mol. Phylogenet. Evol. 74, 48–65 (2014).CAS 
PubMed 

Google Scholar 
Hutchinson, G. E. Cold spring harbor symposium on quantitative biology. Concl. Remarks 22, 415–427 (1957).
Google Scholar 
Pearman, P. B., Guisan, A., Broennimann, O. & Randin, C. F. Niche dynamics in space and time. Trends Ecol. Evol. 23, 149–158 (2008).PubMed 

Google Scholar 
Wake, D. B., Hadly, E. A. & Ackerlya, D. D. Biogeography, changing climates, and niche evolution. Proc. Natl. Acad. Sci. U. S. A. 106, 19631–19636 (2009).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Smith, A. B., Godsoe, W., Rodríguez-Sánchez, F., Wang, H. H. & Warren, D. Niche estimation above and below the species level. Trends Ecol. Evol. 34, 260–273 (2019).PubMed 

Google Scholar 
Peñalver-Alcázar, M., Jiménez-Valverde, A. & Aragón, P. Niche differentiation between deeply divergent phylogenetic lineages of an endemic newt: implications for Species Distribution Models. Zoology 144, 125852 (2021).PubMed 

Google Scholar 
Di Pasquale, G. et al. Coastal Pine-Oak Glacial Refugia in the mediterranean basin: A biogeographic approach based on charcoal analysis and spatial modelling. Forests 11, 673 (2020).
Google Scholar 
Du, Z., He, Y., Wang, H., Wang, C. & Duan, Y. Potential geographical distribution and habitat shift of the genus Ammopiptanthus in China under current and future climate change based on the MaxEnt model. J. Arid Environ. 184, 104328 (2021).ADS 

Google Scholar 
Kafash, A. et al. The Gray Toad-headed Agama, Phrynocephalus scutellatus, on the Iranian Plateau: The degree of niche overlap depends on the phylogenetic distance. Zool. Middle East 64, 47–54 (2018).
Google Scholar 
Namyatova, A. A. Climatic niche comparison between closely related trans-Palearctic species of the genus Orthocephalus (Insecta: Heteroptera: Miridae: Orthotylinae). PeerJ 8, e10517 (2020).PubMed 
PubMed Central 

Google Scholar 
Zhang, Z. et al. Lineage-level distribution models lead to more realistic climate change predictions for a threatened crayfish. Divers. Distrib. 27, 684–695 (2021).
Google Scholar 
Mammola, S. & Leroy, B. Applying species distribution models to caves and other subterranean habitats. Ecography (Cop.) 41, 1194–1208 (2018).
Google Scholar 
Mammola, S. et al. Challenges and opportunities of species distribution modelling of terrestrial arthropod predators. Divers. Distrib. 00, 1–19 (2021).
Google Scholar 
Saupe, E. E., Papes, M., Selden, P. A. & Vetter, R. S. Tracking a medically important spider: Climate change, ecological niche modeling, and the brown recluse (Loxosceles reclusa). PLoS ONE 6, 2 (2011).
Google Scholar 
Planas, E., Saupe, E. E., Lima-Ribeiro, M. S., Peterson, A. T. & Ribera, C. Ecological niche and phylogeography elucidate complex biogeographic patterns in Loxosceles rufescens (Araneae, Sicariidae) in the Mediterranean Basin. BMC Evol. Biol. https://doi.org/10.1186/s12862-014-0195-y (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
Taucare-Ríos, A., Nentwig, W., Bizama, G. & Bustamante, R. O. Matching global and regional distribution models of the recluse spider Loxosceles rufescens: to what extent do these reflect niche conservatism?. Med. Vet. Entomol. 32, 490–496 (2018).PubMed 

Google Scholar 
Wang, Y., Casajus, N., Buddle, C., Berteaux, D. & Larrivée, M. Predicting the distribution of poorly-documented species, Northern black widow (Latrodectus variolus) and Black purse-web spider (Sphodros Niger), using museum specimens and citizen science data. PLoS ONE 13, e0201094 (2018).PubMed 
PubMed Central 

Google Scholar 
Jiménez-Valverde, A., Decae, A. E. & Arnedo, M. A. Environmental suitability of new reported localities of the funnelweb spider Macrothele calpeiana: An assessment using potential distribution modelling with presence-only techniques. J. Biogeogr. 38, 1213–1223 (2011).
Google Scholar 
Monsimet, J., Devineau, O., Pétillon, J. & Lafage, D. Explicit integration of dispersal-related metrics improves predictions of SDM in predatory arthropods. Sci. Rep. https://doi.org/10.1038/s41598-020-73262-2 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
Salgado-Roa, F. C., Gamez, A., Sanchez-Herrera, M., Pardo-Diaz, C. & Salazar, C. Divergence promoted by the northern Andes in the giant fishing spider Ancylometes bogotensis (Araneae: Ctenidae). Biol. J. Linn. Soc. 132, 495–508 (2021).
Google Scholar 
Mammola, S., Goodacre, S. L. & Isaia, M. Climate change may drive cave spiders to extinction. Ecography (Cop.) 41, 233–243 (2018).
Google Scholar 
Ferretti, N. E., Soresi, D. S., González, A. & Arnedo, M. An integrative approach unveils speciation within the threatened spider Calathotarsus simoni (Araneae: Mygalomorphae: Migidae). Syst. Biodivers. 17, 439–457 (2019).
Google Scholar 
Pavlek, M. & Mammola, S. Niche-based processes explaining the distributions of closely related subterranean spiders. J. Biogeogr. 48, 118–133 (2021).
Google Scholar 
Bosso, L. et al. Nature protection areas of Europe are insufficient to preserve the threatened beetle Rosalia alpina (Coleoptera: Cerambycidae): evidence from species distribution models and conservation gap analysis. Ecol. Entomol. 43, 192–203 (2018).
Google Scholar 
Kafash, A. et al. Climate change produces winners and losers: Differential responses of amphibians in mountain forests of the Near East. Glob. Ecol. Conserv. 16, e00471 (2018).
Google Scholar 
Vásquez-Aguilar, A. A., Ornelas, J. F., Rodríguez-Gómez, F. & Cristina MacSwiney, G. Modeling future potential distribution of buff-bellied hummingbird (Amazilia yucatanensis) under climate change: species vs subspecies. Trop. Conserv. Sci. 25, 2 (2021).
Google Scholar 
Rosauer, D. F., Catullo, R. A., VanDerWal, J., Moussalli, A. & Moritz, C. Lineage range estimation method reveals fine-scale endemism linked to pleistocene stability in Australian rainforest herpetofauna. PLoS ONE 10, e0126274 (2015).PubMed 
PubMed Central 

Google Scholar 
Eyres, A., Eronen, J. T., Hagen, O., Böhning-Gaese, K. & Fritz, S. A. Climatic effects on niche evolution in a passerine bird clade depend on paleoclimate reconstruction method. Evolution 75, 1046–1060 (2021).PubMed 

Google Scholar 
Loyola, R. D., Lemes, P., Brum, F. T., Provete, D. B. & Duarte, L. D. S. Clade-specific consequences of climate change to amphibians in Atlantic Forest protected areas. Ecography (Cop.) 37, 65–72 (2014).
Google Scholar 
Muñoz, M. M. & Bodensteiner, B. L. Janzen’s hypothesis meets the bogert effect: Connecting climate variation, thermoregulatory behavior, and rates of physiological evolution. Integr. Org. Biol. 1, 1–12 (2019).
Google Scholar 
Entling, W., Schmidt, M. H., Bacher, S., Brandl, R. & Nentwig, W. Niche properties of Central European spiders: Shading, moisture and the evolution of the habitat niche. Glob. Ecol. Biogeogr. 16, 440–448 (2007).
Google Scholar 
Lafage, D., Maugenest, S., Bouzillé, J. B. & Pétillon, J. Disentangling the influence of local and landscape factors on alpha and beta diversities: opposite response of plants and ground-dwelling arthropods in wet meadows. Ecol. Res. 30, 1025–1035 (2015).
Google Scholar 
Peterson, A. T., Soberón, J. & Sánchez-Cordero, V. Conservatism of ecological niches in evolutionary time. Science 285, 1265–1267 (1999).CAS 
PubMed 

Google Scholar 
Wellenreuther, M., Larson, K. W. & Svensson, E. I. Climatic niche divergence or conservatism? Environmental niches and range limits in ecologically similar damselflies. Ecology 93, 1353–1366 (2012).PubMed 

Google Scholar 
Nosil, P. & Sandoval, C. P. Ecological niche dimensionality and the evolutionary diversification of stick insects. PLoS ONE 3, e1907 (2008).ADS 
PubMed 
PubMed Central 

Google Scholar 
McCormack, J. E., Zellmer, A. J. & Knowles, L. L. Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation?: Insights from tests with niche models. Evolution 64, 1231–1244 (2010).PubMed 

Google Scholar 
Goudarzi, F., Hemami, M. R., Malekian, M. & Fakheran-Esfahani, S. Ecological Characterization of the breeding habitat of Luristan newt (Neurergus kaiseri) at local scale. J. Nat. Environ. 72, 113–127 (2019).
Google Scholar 
Chase, J. M. & Leibold, M. Ecological Niches: Linking Classical and Contemporary Approaches (University of Chicago Press, 2003).
Google Scholar 
Bonte, D., Vandenbroecke, N., Lens, L. & Maelfait, J. P. Low propensity for aerial dispersal in specialist spiders from fragmented landscapes. Proc. R. Soc. B Biol. Sci. 270, 1601–1607 (2003).
Google Scholar 
GBIF.org. GBIF Occurrence Download. https://doi.org/10.15468/dl.2tc2ja (2021) doi:https://doi.org/10.15468/dl.2tc2ja.Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).
Google Scholar 
Jarvis, A., Reuter, H. I., Nelson, A. & Guevara, E. Hole-Filled SRTM for the Globe Version 4. Available from the CGIAR-CSI SRTM 90m Database. (2008) doi:https ://srtm.csi.cgiar .org.Hijmans, R. J. raster: Geographic Data Analysis and Modeling. R package version 3, 3–7 (2020).
Google Scholar 
Guisan, A., Thuiller, W. & Zimmermann, N. E. Habitat suitability and distribution models: With applications in R. (2017). doi:10.1017/ 9781139028271.Quinn, G. P. & Keough, M. J. Experimental Design and Data Analysis for Biologists (Cambridge University Press, 2002).
Google Scholar 
Naimi, B. Uncertainty Analysis for Species Distribution Models. R package version (2015).Phillips, S. J., Dudík, M. & Schapire, R. E. Maxent software for modeling species niches and distributions (Version 3.4.1). Available from url: http://biodiversityinformatics.amnh.org/open_source/maxent/. Accessed on 2022–2–12.Nǎpǎruş, M. & Kuntner, M. A GIS model predicting potential distributions of a lineage: a test case on hermit spiders (Nephilidae: Nephilengys). PLoS ONE 7, e30047 (2012).ADS 
PubMed 
PubMed Central 

Google Scholar 
Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).
Google Scholar 
Merow, C., Smith, M. J. & Silander, J. A. A practical guide to MaxEnt for modeling species’ distributions: What it does, and why inputs and settings matter. Ecography (Cop.) 36, 1058–1069 (2013).
Google Scholar 
Swets, J. A. Measuring the accuracy of diagnostic systems. Science 240, 1285–1293 (1988).ADS 
MathSciNet 
CAS 
PubMed 
MATH 

Google Scholar 
Schoener, T. W. The anolis lizards of Bimini: Resource partitioning in a complex fauna. Ecology 49, 704–726 (1968).
Google Scholar 
Warren, D. L., Glor, R. E. & Turelli, M. Environmental niche equivalency versus conservatism: Quantitative approaches to niche evolution. Evolution 62, 2868–2883 (2008).PubMed 

Google Scholar 
Warren, D. L. et al. ENMTools 1.0: an R package for comparative ecological biogeography. Ecography 44, 504–511 (2021).
Google Scholar 
Liu, C., Berry, P. M., Dawson, T. P. & Pearson, R. G. Selecting thresholds of occurrence in the prediction of species distributions. Ecography (Cop.) 28, 385–393 (2005).
Google Scholar 
Vale, C. G., Tarroso, P. & Brito, J. C. Predicting species distribution at range margins: Testing the effects of study area extent, resolution and threshold selection in the Sahara-Sahel transition zone. Divers. Distrib. 20, 20–33 (2014).
Google Scholar  More