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.