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).