Plant community surveys
Plant community composition was measured in each plot annually from
2006-2020 during the peak of the growing season with pre-treatment data
collected in 2006. A point-intercept method was used to estimate species
presence at 100 points per plot in the field and these raw species
counts were used in subsequent modeling with a censoring term of the
total number of vegetative hits (i.e. excluding rock, litter,
non-vascular species) in a plot in a given year (mean=90). Thirty-three
unique plant species were present in control plots across all years,
however we only included species (n=20) with at least total 20
observations in control plots across all years.
For our modeling approach (see below), we summed the cover data of these
20 species into four species groups based on natural breaks in their
relative abundance in control plots over time. First, the ‘dominant’
species, Deschampsia cespitosa (grass) had an average of 42 ± 1.2
(SE) plot hits (range: 20-67) in ambient conditions (control plots)
forming a standalone group. Three ‘subdominant’ species were combined
into one group: Geum rossii (forb), Artemisia scopulorum(forb), and Carex scopulorum (sedge) which had an average of 10 ±
0.8 plot hits (range: 14-44) in ambient conditions. Four species were
combined into one ‘moderate’ group: Gentiana algida (forb),Trifolium parryi (legume), Bistorta bistortoides (forb),
and Caltha leptosepala (forb) which had an average of 3 ± 0.3
plot hits (range: 3-29) in ambient conditions. Finally, we placed the
remaining 12 species into one ‘rare’ group which had an average of 0.4 ±
0.2 plot hits (range: 0-10) in ambient conditions (Fig S1). Raw cover of
all species over time in all treatment and control plots are shown in
Fig S2.
We calculated changes in relative cover (plot hits) of each species
group (Dominant, Subdominant, Moderate, and Rare) with respect to the
pretreatment (2006) data for each year over the 15 year period within
each experimental treatment using the abundance_change function
in the package CODYN in R (Hallett et al. 2016). We then modeled
these values using a linear mixed model with a fixed 3-way interaction
and a global intercept (0+ time (years since 2006) x species group x
treatment) and a random intercept of (calendar) year to determine
whether each group increased, decreased, or did not change in relative
plot cover over the time period within a given treatment. Models were
run using the lmer function in package lme4 in R (Bates et
al. 2014; R Core Team 2020).