Methods
Study Area. All data were collected from July 2019 to March 2020 on the Qilian Mountains (94°10′E- 103°04′E, 35°50′N-39°19′N), a mountain on the northeastern border of the Tibetan Plateau (China)40,41. The environmental characteristics of the Qilian Mountains vary significantly with the elevation gradients (maximum elevation range of 5806 m). The vertical precipitation distribution shows a quasi-linear manner with elevation. The main rainy season is normally from April to September, with average annual precipitation from 150 to 800 mm42,43. The eastern part is affected by the southeast monsoon and is relatively humid, with more precipitation than the western part. The mean annual temperature (MAT) decreases with elevation, starting at 8°C in the foothills and decreasing to 2 °C in the alpine meadows.
The Qilian Mountains are rich in minerals (i.e., asbestos, pyrite, chromite, copper, lead, and zinc)44. They mainly consisted of Cu-Ni, Cu-Ni-Pt deposits (mafic igneous rocks) related to volcanic magma, massive Zn-Pb-Cu and Zn-Cu sulfide deposits related to marine volcanic rocks and intrusive rocks (ophiolite), and several types of skarn-quartz vein tungsten (Cu, Mo) deposits45. Large-scale mining commenced in the 1980s and ceased at the beginning of the 21st century. Up to now, only several large mines are still in operation (over 1500 m. a.s.l). Owing to a long history of mining activities, natural habitats in the lowland and sub-montane zones of the Qilian Mountain were largely damaged. Some regions are protected as a national nature reserve of the Qilian Mountains (97°25′-103°46′ E, 36°43′-39°36′ N).
Experimental Design and Sample collection . We investigated eight abandoned mine areas of the Qilian Mountains spanned an elevational gradient of 1500 to 3,600 m a.s.l, on the basis of mineralization types, distribution of major minerals, periods of mining exploitation, habitat types, as well as elevations (Supplementary Materials Fig. 7). We excluded mine areas in climate zones like alpine and subalpine meadows (≥ 4000 m) because the mineral resources and the PTE-polluted habitats are mostly distributed in the lower elevation zones. The eight abandoned mine areas occurred in the four major natural habitats of the Qilian Mountains:
(1) Low-elevation deserts (arid and semi-arid desert, 1578-2183 m a.s.l), containing JYG, metal industry (salt desert, halophyte); SN, multi-metal industrial park (shrub desert, bushes);
(2) Low-elevation grasslands (semi-arid grasslands, 2365-3079 m a.s.l), containing HB, Chrome chemical factory (flat grasslands, vegetations); AR, an asbestos factory (mountain grasslands, shrubs and grass); XF, Chrome chemical factory (field, vegetations);
(3) High-elevation deserts (alpine desert steppes, 3043-3203 m a.s.l), only XTS, Solder iron mountain (dry/cold desert, shrubs and grass);
(4) High-elevation grasslands (semi-damp grasslands, 3427-3618 m a.s.l), containing HX, multimetal mine (mountain grasslands, shrubs and grass); SBG, multimetal mine (mountain grasslands, shrubs and grass).
Among the eight abandoned mine areas, we established sixty-four study sites (Supplementary Materials Fig. 7). A total number of 48 PTE-polluted habitats and 16 natural habitats (e.g., not contaminated or lightly contaminated) were selected. This arrangement was aimed at comparing the plant species diversity and the ecosystem functions of the habitats. All the study sites (n=64) were chosen from the core zones of the corresponding habitats to eliminate the disturbance of the transition zones. From 1 to 6 study plots (50 × 50 m2) per study sites (n=48, the PTE-polluted habitats) were arranged with a total number of 174 plots to reflect a within-habitat-type elevation gradient showing the fine-scale changes in biodiversity at changing environments. The soil and plant samples were collected from each plot. The soil samples were properly packed in plastic bags (500 g per piece), and plant samples (more than 3 individuals) were packed and stored in a low temperature (-4°C). Six topsoil samples (0-20 cm) were collected from the Qilian Mountains Natural Reserve where the ecological system had never been affected by mining activities. Then, the samples were taken back to the laboratory of the Physical and Chemical Center of the Institute of Geography, Chinese Academy of Sciences for further analysis.
PTEs, topography and climate data. PTEs. The PTEs was calculated to characterize the strength of potentially toxic elements by using four evaluation criteria (1) ecological risk index, (2) mining years, (3) closure years, (4) the relative proportion of mining area/plant/factory in the surrounding landscape46, 47,48. The four criteria were calculated as follows:
(1) Ecological risk index (RI ), introduced by Hakanson (1980)49, was calculated as the sum of the Ei values for each trace metal element according to the following equation:
\begin{equation} Ri=\sum E_{i}\nonumber \\ \end{equation}\begin{equation} E_{i}=T_{i}\frac{C_{i}}{B_{i}}\nonumber \\ \end{equation}
where Ei is the single risk index for the trace metal element i , dimensionless; Ti is the toxic-response factor for a given metal (1 for Zn, 5 for Cu, 5 for Ni, 5 for Pb, 2 for Cr, 30 for Cd); dimensionless; Ciis the concentration of metals in soils, mg kg-1;Bi is the background value for metals, mg kg-1.
(2) The mining years, from the beginning of mining activities to closure (open ~ shutdown). This indicator shows when the toxic elements started affecting the surrounding environment during the mining process and how long the process lasted.
(3) Closure years, from the mine shutdown to 2020 (shutdown ~ 2020). This indicator reflects the long-term effects of potentially toxic elements on ecosystems after pollution cessation.
(4) The relative proportion of mining area/plant/factory in the surrounding landscape, which is the relative proportion of natural habitats and mine areas. This indicator reflects the relationship between mine areas and the disturbed ecosystem areas around mine areas.
All four criteria of the PTEs were standardized before averaging them to final composite PTEs. According to the following equation:
y = (yiȳ )/(maximum(y ) − minimum(y ))
where y is the measured value of each component. The produced PTEs value is between 0 (lowest PTE intensity) and 1 (highest PTE intensity).
Climate data . The 1980-2018 climate data were collected from 8 weather stations close to the eight studied mine areas of the Qilian Mountains. In each mine area, the mean annual climate (MAC) data were calculated to characterize the climate change among different elevation gradients. Data included monthly surface evapotranspiration (ET0), average temperature, average precipitation, and average wind speed (win10) (http://www.cma.gov.cn/).
Topography data . Coordinates and elevation of the 64 study sites were recorded by hand-held GPS (China, Yili-X28). Data on the slope and aspect were extracted from the 90 m × 90 m Digital Elevation Model (DEM) of Gansu and Qinghai province (SRTM 90 m Digital Elevation Data; http://www.gscloud.cn/). ArcGIS 10.2 was used to edit and filter the 90 m resolution DEM data with the Qilian Mountains’ boundary and acquired the study area map.
Species diversity and indicators of ecosystem functions.Species richness, evenness, and richness. For each study site (n = 64), the total number of species (cumulative) was recorded and vascular plant species richness was calculated (Supplementary Materials Fig. 7). In each study site, 1 to 6 study plots (50 × 50 m2) were arranged depending on local terrain and plant distribution, with a total number of 174 plots. For each plot, 5 subplots with 25 sampling points (1 × 1 m2) in total were set up with the standard sampling method. Data on species coverage, richness, and evenness of dominant herbaceous species in each point were recorded and their mean values were calculated for further analysis50,51.
Soil physical and chemical properties (pH, N, P, K, EC and soil TOC ). Soil samples were taken from soil pits on a horizon basis, and in 10-cm intervals, using a standard soil auger. The samples were sieved (2-mm mesh), and roots and plant materials were removed. Field samples were dried at room temperature for 72 h, and ground for further analysis. For all study sites, results were standardized and averaged to 0-40 cm depth. Subsequently, chemical analyses were carried out in the Physical and Chemical Analysis Center of the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (Beijing, China). Soil pH was determined with an inoLab Level 1 pH meter in suspension with 1:1.25 (w/v) soil, deionized water, and 0.01 M CaCl2. Total organic carbon (TOC) was measured by wet oxidation with K2Cr2O7and the absorbance at 590 nm was measured with Spectrophotometer (Liqui TOC Ⅱ, Elementar, Germany). Thereafter, total nitrogen (N) was measured by a dry combustion automated C: N analyzer (Vario EL, Elementa). Available P extracted by 0.002 N H2SO4solutions was measured with the Truog method, while available K was determined with NH4OAc using gas chromatography-flame photometry.
Soil water content. We extracted the monthly global 0.5° × 0.5° soil moisture data from 1980 to 2004 for our studied sites in the Qilian Mountains from the 1948-2004 dataset originally calculated by Fan et al. (2004)52. Afterward, the average value of annual soil moisture was calculated for each of the 64 studied sites.
Plant biomass. We determined plant biomass per study site by summing up the biomass of shrubs and herbs. Shrub biomass was estimated in subplots of 5 × 20 m in the center of each study site and then extrapolated to the scale of study sites. Herbaceous biomass was determined by removing all plant biomass of the herbaceous layer in four representative subplots of 0.25 × 0.25 m per study site when the herbaceous biomass was maximal. Herbaceous samples were dried in a drying oven at 65 °C for 72 h and then weighed.
Root fine biomass. At each plot, 10 to 15 soil cores (inner diameter 3.5 cm) were taken randomly down to a depth of 40 cm. Soil samples were kept in plastic bags and stored below 4 °C. In the laboratory, samples were soaked and washed via a sieve of 200 μm mesh size. All root fragments were extracted by hand with tweezers (root lengths ≤ 0.2 cm was ignored). Root fractions were dried at 70 °C for 48 h, and weight and fine root biomass were expressed as mg ha-1.
Plant physiological properties (C/N, C/P, K, and TOC). All plant physical and chemical properties analyses were carried out in the Physical and Chemical Analysis Center of the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (Shandong, China). Total organic carbon (TOC) was measured by wet oxidation with K2Cr2O7and the absorbance at 590 nm was measured with Spectrophotometer (Liqui TOC Ⅱ, Elementar, Germany). Total nitrogen (N) was measured with the Kjeldahl method. Available P extracted by 0.002 N H2SO4 solutions was measured with the Truog method while available K was determined with NH4OAc using gas chromatography-flame photometry.
Plant bio-concentration factor (BCF) . Metal (Al, Fe, Mg, Cd, Cr, Cu, Ni, Pb, Sr, Zn, Mn, and Ba) concentration of both soil and plant samples was determined by using a mixture of acids (HCl-HNO3-HF-HClO4) to digest and the extraction solution was analyzed through ICP-OES (Optima 5300DV, PerkinElmer, USA). A quality assurance/quality check (QA/QC) program was established to ensure the accuracy of the results of metal concentrations in all samples. Furthermore, three samples from sixty-four studied sites were taken for inter-laboratory analysis to further ensure the validity of the data. Finally, the bio-concentration factor (BCF) of plants was calculated as follows53,54.
\begin{equation} \text{BCF}_{\text{metal\ j}}=\frac{1}{n}\sum_{i=1}^{n}{\frac{\text{CPT}}{\text{CE}}\text{\ \ \ \ \ \ \ \ }\left(i=1,\ 2,\ \ldots n\right)}\nonumber \\ \end{equation}
where CPT and CE are the concentration (mg kg-1) of metal j in plant tissue and soil, respectively, i is the i-th plant in the plot, andn is the total number of plants in the plot. Then, the factor analysis method was used to calculate the normalized BCF values.
Normalized Difference Vegetation Index (NDVI). NDVI dataset for the 64 study sites was extracted from Landsat 5 and 8 series (horizontal resolution of 30 × 30 m) and MYD13Q1 (horizontal resolution of 250 × 250 m). Cloud pollution was a prominent feature for the remote sensing data extraction, limiting the accuracy of environmental indices computation. To solve this problem, we first identified all pixels marked as 3 of MYD13Q1 quality and then deleted them along with eight surrounding pixels55. Using a ”Whitaker smoother” based on three iterations of 6,000λ, we calculated the overall average NDVI from 1980 to 2018 and extracted the pixel values corresponding to the study sites.
Leaf Area Index (LAI) . We measured the LAI at 25 points per study site and calculated the mean value from these measures. The leaf area was measured at ground level using a plant canopy analyzer in combination with a remote ‘above-canopy’ sensor (LAI 2200, LI-COR Bioscience). The above canopy reading was measured at ground level in an open area or forest gap as close to the site as possible, and LAI values were calculated using the program FV-2200 (LI-COR Bioscience).
Statistical Analysis. Generalized additive models (GAMs) was employed to analyze the correlation between species biodiversity (richness, evenness, and cover) or ecosystem function (response variables) and elevation (continuous predictor)56. All the ecosystem function indicators were transformed into z-axes before the analysis. Considering the consistency of data analysis, pollution types were added into the model as a binary factor variable (natural habitats versus PTE-polluted habitats). The binary factor variable was based on our main statistical analyses (i.e., the tests of seven major PTEs, climate, and topography models described below) on models in which the potential toxic element intensity was integrated as a continuous predictor variable (that is, the PTEs ). Poisson or Gaussian distributions were applied in the case of species richness and ecosystem functions and the basic dimension of the smoothing function was set to 4 to 6. Here, pollution type was added to the GAMs as an interaction factor to determine the correlations of response variables and elevation between PTE-polluted and natural habitats. If the interaction term was not significant (P > 0.05), the corresponding term was removed from the model but the additive effect of elevation and PTE dispersal was anyway tested (i.e., we tested for differences in the intercept between natural and PTE-polluted habitats). If the effect of PTE dispersal was not significant, the corresponding term was removed from the prediction variable and the response variable was modeled only by elevation, assuming that the mean distributions of response variables in the natural and PTE-polluted habitats were the same. Otherwise, we assumed the null model to be best supported by the data, and the response variable was modeled only by the intercept (i.e., the average values of the data).
To assess the effects on biodiversity and ecosystem functions determined by PTE dispersal and interacted by climate and topography, we evaluated the support for seven different hypotheses with linear models within a multimodel inference approach57. Multimodel inference takes both uncertainty in parameter estimation and uncertainty in model selection into account. Moreover, the multimodel inference is good at dealing with related predictor variables. In the model, strongly correlated predictors with the same explanatory power reduce support for the multimodel. The seven hypotheses or models we tested were:
(1) PTEs only (PTEs model); (2) PTEs + Climate (additive model); (3) PTEs + Topography (additive model); (4) PTEs + Climate + Topology (additive model); (5) PTEs × Climate (interaction model); (6) PTEs × Topography (interaction model); (7) PTEs × Climate × Topography (interaction model).
In the PTEs model (1), we assumed that the PTE impact alone could fully explain the variations of biodiversity and ecosystem functions. In the additive model (2) to (4), we assumed that PTEs, topography, and/or climate had additive effects. The best model was then selected and at least one PTEs, one climate, and/or one topography variable were included. In the interaction model (5) and (6), we tested the hypothesis that climate (or topography) could influence the impact of PTEs on biodiversity and ecosystem functions. The interaction terms of PTEs-climate (or PTEs-topography) were added in the GAMs. The best model (with the lowest value of Akaike information criterion with correction for small sample sizes (AICC )) was selected and at least one interaction term of PTEs-climate variables (or PTEs-topography variables) was included. In the interaction model (7), instead of testing the interaction terms of the three indicators, we run the model which involved both the interaction terms of PTEs-climate and those of PTEs-topography. The best model was selected (with the lowestAICC value) and at least one PTEs-topography and one PTEs-climate interaction terms were included.
For each response variable (plant richness, evenness, coverage, and ecosystem functions), we calculated the support for each of the above-mentioned models by calculating the AIC-based model weight (also known as ‘AIC weight’ or ‘model probability’)58. The lower the AIC value was, the better the fitting effect was. Meanwhile, we calculated the R2 value of the optimal model to represent the relative importance and significance of each prediction variable. For multimodel averaging analyses, the R package MuMIn (https://CRAN.R-project.org/package=MuMIn) was used. To analysis the robustness of the model, we analyzed the uncertainty of ecological risk and the influence on the stability of the model from other environmental variables (slope direction, atmospheric pressure, etc.) that could change the modeled results61.
To visualize and analyze the relative importance of PTEs, climate, and topography on species turnover (β-diversity), we calculated a dissimilarity matrix of plant diversity based on the Sørensen dissimilarity index. The nonmetric multidimensional scaling (NMDS) using Meta MDS function of the vegan R-package was used to ordinate the dissimilarity matrices in two-dimensional ordination space. For all biodiversity datasets, stress levels were low, which justified using this approach. The effect of PTEs, climate, and topography on the dissimilarity of plant communities was examined through the permutational multivariate analysis of variance (PERMANOVA) method (the function adonis in the vegan package). We started with the most complex model and deleted the least important predictor variables in sequence to obtain the optimal model for each of the seven assumptions.
To quantify the overall differences in ecosystem functions of multiple mine areas of the study area, we calculated two measures of ecosystem multifunctionality: (a) a multivariate index of multifunctionality, and (b) a measure of the average change of ecosystem functions in PTE-polluted, compared to the natural habitats.
The multivariate index of multifunctionality (a) was calculated by removing the default values in pairs and obtaining the Euclidean dissimilarity matrix based on the z-transformed ecosystem functions. Classical multidimensional scaling was used to visualize the dissimilarity matrix in the sorting space, while PERMANOVA was used to determine the significance of climate, PTEs, and topography on multidimensional scaling axis 1 and 2 scores.
We quantified the dissimilarity in ecosystem functions of PTE-polluted habitats to the average conditions in the natural habitats (b) by modeling ecosystem functions in natural ecosystems in a GAM with elevation as a predictor variable. Before analysis, all the ecosystem functions were transformed into z-axes. We then modeled the average dissimilarity in ecosystem functions as a function of PTE intensity and m-PTE intensity, and compared differences in PTE intensity among elevation zones by using simple ordinary linear models and AIC-based model inference. The values of m-PTEs were the multi-dimensional factors (axis1) obtained from multidimensional scaling of a total of 12 variables of climate, PTEs, and topography variables in the final model. This analysis was done for all 20 ecosystem functions and separately for the soil- and plant-mediated ecosystem functions.