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.