MATERIALS AND METHODS
Site and Census design. We examine six years of census data
from a site on a southeast facing ridgeline with ~20°
slope, on Mount Baldy, in the Gunnison National Forest of western
Colorado (Blonder et al. 2018). The site contains a sparsely populated
alpine plant community with a total of eighteen species in its most
species-rich year. In the year with greatest abundance, there was an
average of only nine individuals per square meter. For a more complete
site description, a species list, and an analysis of the relationships
among species traits, individual survival, fecundity, recruitment, and
abiotic microenvironmental factors, see Blonder et al. (2018).
The relative sparseness of the vegetation, compared to other alpine
plant communities with short growing seasons, is largely due to the
substrate, which is composed of exposed and undeveloped rock strata,
comprising Mancos shale (Upper Cretaceous) downslope and quartz
monzonite porphyry (Upper Eocene) upslope. Soil across the site consists
primarily of ~5-10 cm of fine weathered black scree over
bedrock. Snowfall usually begins in October with snowmelt completed in
June or July.
The data set used here comprises six years from 2014 to 2019 and derives
from a complete annual vegetation census. All individuals including
seedlings were recorded and located to the nearest centimeter. The site
is partitioned by two-meter-wide corridors into a 5 by 10 grid with 50
plots, each with dimension two meters by two meters. Data are collected
annually during peak flowering season in late July or early August of
each year. Individuals are considered dead and removed from the census
in subsequent years if no above-ground growth was recorded for more than
two years or if they germinated and died within the same year. Data and
field methods are originally published for 2014-2017 in Blonder et al.
(2018) and updated for 2014-2019 in Ray et al. (2020).
Species-area relationship (SAR). The SAR describes the
dependence of species richness on censused area. A complete nested SAR
design is derived by dividing the entire community into disjoint
equal-area quadrats and averaging the number of species over all
quadrats of equal area. To construct species-area curves we ignored the
corridors between the 50 plots and combined the data into one continuous
plot (see Discussion for justification). We then re-sectioned the data
into a 32 by 64 grid, in order to have a sufficient number of either
square quadrats or \(2^{n}\)x\(2^{n-1}\) rectangular quadrats. Both
horizontal and vertical orientations were considered for rectangular
quadrats by averaging over equal areas (see Table S1 in Supporting
Information).
This design allowed us to study areas of relative size: 32, 64, 128,
256, 512, 1024, and 2048 (=32x64), where the smallest scale was set by
the requirement that we have at least six species on average at that
scale. The criterion of at least six species is based on a necessary
requirement for the validity of METE: exp(-S )
<< 1 (Harte 2011). This requirement is also why we
were not able to construct 50 separate SARs each year using each of the
50 plots as a study system; the species richness would have been too low
within the plots to be able to examine the SAR across a wide range of
spatial scales.
METE predicts that the slope of the species-area curve at a given scale
is a function of the ratio of the number of individuals, N, to
the number of species, S, at that scale. Hence, if data are
plotted as local slope, z, at a given scale versus the ratio, or
log of the ratio, of N to S at that scale, then all SARs
are predicted to collapse onto a universal curve (Harte et al. 2009).
Thus, we express the SAR by plotting z, the local slope of the
species-area curve expressed as log(S ) versus log(A ),
against log(N /S ), where N and S are measured
at the same spatial scale at which z is calculated.
To an excellent approximation, METE predicts that the area-dependent
slope, z (A ), as a function of N (A ) andS (A ) is given by
\(z(A)\cong\left[\log\left(2\right)\log\left(\frac{1}{1-e^{-\beta\left(N\left(A\right),S\left(A\right)\right)}}\right)\right]^{-1}\)(1)
where \(\beta\) is given by the solution to
\(\frac{N(A)}{S(A)}=\frac{\sum_{n=1}^{N(A)}e^{-\beta n}}{\sum_{n=1}^{N(A)}\frac{e^{-\beta n}}{\text{\ n}}}\)(2)
and natural logarithms are used. A more exact form for the METE
predicted slope z (A ) in Eq. 1 is given in Harte (2011).
Because our average values of N at each scale can be non-integer,\(\beta\) is calculated using linear interpolation between two adjacent
integer values.
To calculate the observed slopes of the species-area curve at each area,
we used the slope of the line connecting the coordinates of the next
lowest and next highest area in the species-area curve, which is
consistent with the way in which Eq. 1 above was derived. This meant we
could not calculate observed slopes at the endpoints of our species-area
curve, leaving a total of six slopes.
To quantify the change in prediction performance over the years we
compared the mean squared error between the observed and predicted
slopes for each year.
\(\text{error}=\ \frac{1}{K}\sum_{i=1}^{K}\left(\text{observed\ slop}e_{i}-\text{predicted\ slop}e_{i}\right)^{2}\)(3)
where K , the number of spatial scales where we could predict the
slope, is six.
Species-abundance distribution (SAD). The SAD,\(\Phi\left(n\middle|S,N\right),\) is the probability that a species
has n individuals given that the community has S species
and N individuals. The empirical SAD for any community is sparse
because a sample will not give data for every abundance ranging from 1
to N . We instead plot rank versus log(N ) distributions,
with abundances of each species ranked by decreasing abundance (see
Table S2).
METE predicts that the SAD is a logseries distribution given by:
\(\Phi\left(n\middle|S,N\right)=\frac{e^{-\beta n}}{n\log\left(\frac{1}{1-e^{-\beta}}\right)}\)(4)
where \(\beta\) is given by Eq. 2.
We construct a predicted rank abundance curve from Eq. 4 by computing
cumulative probabilities for each n -value from the predicted SAD.
To assess the fit of predictions, we compute the mean squared error
between the predicted and observed rank abundance curves for each of the
years. There is a small uncertainty in the way we calculate the
predicted rank abundance curve as it could be shifted up or down by ½ of
a rank. Thus, for each S and N value in the data set, we
drew 1,000 samples from the predicted SAD and rank ordered each sample
to determine how likely the observed rank abundance curve was.
Stress factors. We examined mortality and recruitment to
potentially quantify stress over time (see Table S3). Mortality is
defined as the number of individuals who died in each year. Mortality
data is not available in 2014 and 2019 because the study did not include
individuals who died in 2014 and data from 2020 is required to determine
which individuals died in 2019. Recruitment, defined as the number of
seedlings in each year, was measured in all six years. By fitting linear
models to the data, we derived correlations between the SAR and SAD
error metrics with mortality, recruitment, net loss (the difference
between mortality and recruitment), number of individuals, and species
richness, and thereby determine if any of these metrics quantify the
stress driving changes in the SAR and SAD.