INTRODUCTION

Anthropogenic climate warming is projected to increase global surface temperatures by 1.0-5.7 C by 2100 (IPCC, 2021). Climate warming studies in plant communities have primarily focused on the direct, individual effects of warming, including changes in the timing of phenological events and community structure and function (Parmesan, 2006;Parmesan & Yohe, 2003;Peñuelas & Filella, 2001;Renner & Zohner, 2018;Root et al., 2003;Walther et al., 2002). However, the interactive effects of warming with biotic factors could also be important in how these communities respond to climate change. For example, warming can affect insect herbivores' preferences and feeding patterns through changes in leaf chemistry and increased insect metabolic rates, and these warming-induced changes in feeding patterns affect the amount of herbivory experienced by plants (Hamann et al., 2021;Kharouba & Yang, 2021;Pinheiro et al., 2016;Welshofer, Zarnetske, Lany, & Read, 2018). Changes in herbivory levels can furthermore affect plants by altering their chemical composition, productivity, and phenology (Lemoine et al., 2017;Post & Pedersen, 2008;Ritchie et al., 1998;Welshofer, Zarnetske, Lany, & Read, 2018). In this case, the plants are indirectly affected by warming, mediated through changes in insect herbivory (Blois et al., 2013;Parmesan, 2006;Post, 2013;Zarnetske et al., 2012). Many such interactive effects, however, are not well understood (Blois et al., 2013;Parmesan, 2006;Post, 2013;Zarnetske et al., 2012). Some of the most commonly observed consequences of climate change are phenological shifts. Shifts in phenology may alter biotic interactions if they alter the relative activity periods of interacting partners. For example, Liu et al. (2011) found that under ambient conditions, gentian flowers typically bloom after the peak density of an insect herbivore. However, experimental warming advanced gentian flowering and delayed the emergence of the herbivore, leading to increased overlap between the two species. These phenological shifts, combined with increased herbivore densities and reduced densities of an alternate host plant, resulted in 100-fold greater damage to gentian flowers and fruits in warmed plots compared to ambient plots.

There can be substantial variation within and among species in the magnitude and direction of warming effects (Primack et al., 2009;Sherry et al., 2007;Youngflesh et al., 2021). Species unable to shift their phenology sufficiently in response to warming may experience negative fitness effects, potentially leading to decreased abundances or even local extinction (Willis et al., 2008). It remains unknown whether the association between the propensity for phenological shifts and extinction suggested by Willis et al. (2008) is driven by the direct effects of warming or indirect effects mediated through temporal mismatches with pollinators, herbivores, or competitors. Interacting organisms can also have their own, separate direct effects on plant phenology. For example, herbivory itself has been found to delay phenology, likely due to plants redirecting resources to repair tissue damage (Lemoine et al., 2017;Welshofer, Zarnetske, Lany, & Read, 2018).

Plant leaf traits, including morphological characteristics and chemical composition, are important cues that herbivores use to find quality food sources. Specific leaf area (SLA), which is the ratio of total leaf area to total leaf dry mass, is an important trait that can reflect whole plant growth (Liu et al., 2017). Warming treatments have been shown to have variable effects on SLA (Descombes et al., 2020;Hudson et al., 2011), but insect herbivores have been found to prefer plants whose leaves have smaller SLA (Dost alek et al., 2020;Pereira et al., 2020). Leaf palatability to herbivores has also been found to be positively correlated with N content and negatively correlated with C content (Schädler et al., 2003). Warming treatments have been found to decrease leaf contents of both C and N (Hudson et al., 2011;Yang et al., 2011). If there are warming-induced reductions in plant food quality, herbivores may need to consume greater quantities of plant material in order to meet nutritional demands (Hamann et al., 2021;Welshofer, Zarnetske, Lany, & Read, 2018). Paleontological records also document increased herbivory during periods of global warming in past geological times (Pinheiro et al., 2016). Herbivory itself can influence foliar C and N content, as insects may prefer to eat nutrient-rich leaf tissue, leading to declines in overall nutrient content for plants (Ritchie et al., 1998).

In terms of plant community composition, experimental warming treatments in alpine systems and temperate grasslands show reductions in species evenness and richness, with up to 25% of species lost, although the magnitude and the direction of diversity effects vary regionally (Elmendorf et al., 2012;Morecroft et al., 2009;White et al., 2014). Experimental warming may also lead to complex changes in community composition, as some species increase in abundance while other species decline (de Valpine & Harte, 2001;Li et al., 2011;Morecroft et al., 2009;Rudgers et al., 2014;Welshofer, Zarnetske, Lany, & Read, 2018). For example, Wangchuk et al. (2021) demonstrated that experimental warming treatments decrease overall plant diversity and richness through the promotion of grass species.

The interaction between climate warming and herbivory can mediate the impact of warming on plant community composition. In particular, herbivory can reduce the impacts of warming on plant diversity and richness (Kaarlejärvi et al., 2017;Post, 2013). Even in the absence of warming, herbivory can affect plant community composition by reducing dominant species and increasing light availability at ground level, therefore helping to maintain plant richness and diversity (Borer et al., 2014;Brown & Gange, 1989;Koerner et al., 2018;Mortensen et al., 2018;Post, 2013;Post & Pedersen, 2008;Price et al., 2022).

Experimental warming treatments have also been found to increase community productivity via increased plant biomass (Wangchuk et al., 2021), but this increase can be dependent upon plant functional type (Lin et al., 2010). Other studies have noted that warming treatments may cause decreases in biomass, potentially due to strong competition for resources under conditions of high stress (De Boeck et al., 2008). Furthermore, herbivory on plants has been found to have contrasting effects on plant biomass due to herbivores having varying preferences for certain species over others (Post & Pedersen, 2008).

The relative influence of the individual and interactive effects of abiotic and biotic factors can be assessed by experimentally manipulating both climate and the presence or abundance of interacting species. However, of the 126 studies on in situ warming experiments with open-top chambers (OTCs) reviewed by Dobson and Zarnetske (2024), only 14 (11%) included a treatment involving species interactions, and only 57 (45%) spanned more than 3 years. In this study, we explored these individual and interactive effects by manipulating temperature and insect herbivory, separately and in combination, for 7 years in two early successional plant communities. We tracked the plant communities' responses to our experimental manipulations by measuring leaf herbivory, phenology (green-up, flowering, flowering duration, and seed set), plant composition, biomass, and leaf traits (SLA, C, N). We hypothesized that:

1. Warmed plots would have greater amounts of insect herbivory than ambient plots (Hamann et al., 2021;Pinheiro et al., 2016;Welshofer, Zarnetske, Lany, & Read, 2018).

2. Warmed plots would experience earlier green-up and flowering, delayed seed set, and longer flowering duration, especially under reduced herbivory (Lemoine et al., 2017;Peñuelas & Filella, 2001;Walther et al., 2002;Welshofer, Zarnetske, Lany, & Read, 2018;Zhou et al., 2022). 3. Warmed plots would have increased percent cover but lower plant species richness and diversity when compared to ambient plots (Elmendorf et al., 2012;Morecroft et al., 2009;Wangchuk et al., 2021;White et al., 2014). The effects of warming would be lessened with herbivory because herbivores can decrease percent cover while increasing plant species richness and diversity by reducing dominant species (Brown & Gange, 1992, 1989;Kaarlejärvi et al., 2017;Post, 2013;Post & Pedersen, 2008;Ritchie et al., 1998;Wangchuk et al., 2021). 4. Warmed plots would have higher plant biomass, higher SLA, and lower foliar C and N content (Descombes et al., 2020;Lin et al., 2010;Yang et al., 2011), especially with herbivory (Post & Pedersen, 2008;Ritchie et al., 1998).

METHODS

Site description

The study system consists of sites in two early successional plant communities in Michigan, USA, separated by 354 km and approximately three degrees latitude. The southern site, located at Kellogg Biological Station's Long-Term Ecological Research Site (KBS-LTER, 42.41 N, 85.37 W), was previously an agricultural field and is now dominated by Solidago canadensis, Poaceae spp., and Hieracium spp. (Appendix S1: Table S1). The mean annual temperature and precipitation in Kalamazoo County, where KBS is located, are 9.33 C and 975.4 mm (30 year means, PRISM Climate Group, n.d.; Appendix S1: Figure S1). The northern site, located at the University of Michigan Biological Station (UMBS, 45.56 N, 84.71 W), is within an old forest clearing that was clear-cut in 1994 and is now dominated by Centaurea stoebe, Pteridium aquilinum, and Carex pensylvanica (Appendix S1: Table S1). The mean annual temperature and precipitation in Emmet County, where UMBS is located, are 6.42 C and 770.7 mm (30 year means, PRISM Climate Group; Appendix S1: Figure S1).

Experimental design

At each site, 24 1 × 1 m plots are contained within a 25 × 36 × 3 m fence that prevents herbivory by deer.

Each 1 m 2 plot is contained within a 3 × 3 m buffer zone, and all plots are separated by at least 4 m to minimize potential shading from the OTCs. The experiments were established in the spring of 2015 and consist of a fully factorial design with warming and insect herbivore reduction treatments (ambient, warming, reduced herbivory, warming + reduced herbivory; n = 6 per treatment). There were no significant initial differences in plant composition between the plots within each site in 2015 before treatments were applied (Welshofer, Zarnetske, Lany, & Read, 2018).

Warming was achieved through the use of hexagonal OTCs designed for taller stature plant communities (Marion et al., 1997;Welshofer, Zarnetske, Lany, & Thompson, 2018; Appendix S1: Figure S1). OTCs simulate climate warming by passively increasing air temperatures in situ while also allowing for natural levels of precipitation, gas exchange, and solar radiation (Marion et al., 1997;Welshofer, Zarnetske, Lany, & Thompson, 2018). These chambers remained on the plots year-round and were constructed with clear, UV-transmitting 1/8 00 Lexan Polycarbonate (ePlastics, San Diego, CA). Insect herbivory was manipulated through insecticide applications throughout the growing season (for details, see Appendix S1: Section S1). Insecticide plots are termed "reduced herbivory," and noninsecticide plots are termed "herbivory" plots.

Data collection

Abiotic measurements

Hourly abiotic conditions were recorded at the plot level at each site using HOBO products (Onset Computer Corporation, Bourne, MA). Three ambient and three warmed plots were equipped with four-channel external U12-008 data loggers that recorded air temperature at 10 cm above ground and soil temperature at 5 cm below ground. These plots also contained microstation H21-002 data loggers that recorded air temperature at 1 m above ground and soil moisture at 5 cm below ground. Plastic dish solar shields were installed above each air temperature sensor to mitigate the impact of solar radiation on air temperature readings.

Leaf herbivory

Leaf herbivory was measured once per season at peak biomass prior to senescence, typically July-August (methods similar to The Herbivory Variability Network, 2023). We haphazardly selected four random leaves vertically distributed across the stem of three individuals of each measured species in each plot (Appendix S1: Figures S2 and S3). We then visually estimated the percentage of the leaf eaten (0%-100%).

Plant phenology

The phenology of all plant species within all plots (Appendix S1: Figures S4-S7) was monitored every 3-4 days. The beginning of data collection at each site was determined by the last snow melt of the year in the spring. Phenology consisted of green-up, flowering, flowering duration, and seed set. Green-up was calculated as the date at which a plot reached half of its maximum percent cover (Appendix S1: Table S2) to account for early season differences in the depth of plant litter, which might affect the detection of plants when they first emerge. A species was recorded as flowering during the period between first flower bud break (anthers exposed) and final flower senescence. Date of first flower was calculated as the average minimum date that a plot recorded a species flowering. The duration of flowering was the number of days between the average date of first flower to the average date of last flower. Seed set was determined when an individual exhibited a mature seed that was ready to be dispersed (pappus/achene, florets dehiscent, etc.), and was calculated as the average minimum date of first occurring seed set per plot.

Plant community composition

Percent aerial cover was visually estimated within all 1 m 2 plots as the percentage of the total plot occupied by each species in each plot (0%-100%; Appendix S1: Figure S8). Because each species in each plot could be estimated up to 100% cover, it is possible that the total calculated percent cover of any given plot could exceed 100%. This measurement was taken every 3-4 days through green-up and once a month post green-up.

Leaf traits (C, N, SLA)

Prior to senescence, green leaves were harvested for measurements of foliar C and N content and SLA. The species selected for these measurements were commonly found across all plots at each respective site (Appendix S1: Table S3). We chose 3-5 plants of the same species in each plot and harvested 4-5 green, mature leaves with little to no obvious insect damage or disease. The selected leaves were haphazardly selected from the top to the bottom of the plant to be representative of the whole plant.

The youngest fully expanded leaf from each individual plant sample was chosen for SLA, while the remaining leaves were stored separately for C and N analysis. SLA leaves were scanned fresh with a LI-COR LI-3000A Portable Leaf Area Meter with conveyor belt LI-3050A at KBS and a LI-COR LI-3100c at UMBS. After SLA leaves were scanned, all leaf samples, including those harvested for C and N analysis, were placed in a drying oven at 60 C for 36-48 h and subsequently weighed. Combustion analysis was then performed for C and N (see Appendix S1: Section S2).

Plant community biomass

In 2021, all aboveground plant biomass was harvested in a 0.20 m 2 area (1 × 0.20 m) within all 1 m 2 plots at both sites. Plant material was sorted into species, placed in individual paper bags, dried at 60 C for 3-4 days, and weighed for a final dry biomass weight per species in each plot.

Statistical analysis

All analyses were conducted using R (R Core Team 2020). All response variables were calculated at the plot level to test for the overall effects of the treatments on the plant community; species-specific effects can be found in the supplement (Appendix S1: Figures S2-S10 and S13-S16, Tables S4-S19). We tested for the individual and interactive effects of warming level (warmed vs. ambient), herbivory level (herbivory vs. reduced herbivory), and year using linear mixed-effects models in R with the package lmerTest (Kuznetsova et al., 2017;R Core Team, 2020). Fixed effects included warming level, herbivory level, and year, with interactions between all three factors (Response

To test whether species or groups themselves differed for each response variable, we included species, plant origin (native or exotic), or plant growth habit (forb or graminoid) as fixed effects in separate species models (Response variable i = β0 + β1warmed + β2insecticide + β3year_factor + β43(warmed × insecticide) + β5(warmed × year_factor) + β6(insecticide × year_factor) + β7(warmed

). Plot number was included as a random effect for all models to account for inherent variation between plots. For leaf herbivory, SLA, C, and N models, individual plant ID was nested within species within plot number and included as a random effect.

To test for evidence of interactive effects for each response variable, we looked for a significant interaction between warming level and herbivory level (Appendix S1: Tables S4-S19). If significant, we then tested the pairwise comparisons of all treatments using the emmeans package (Lenth, 2022; Appendix S1: Tables S20-S32). We also used pairwise comparisons to determine treatment effects for specific years. If there was no significant interaction, we tested for the individual effects of each treatment. We confirmed the data fit the assumption of normality prior to running our models and that there were no outliers with Bonferroni-adjusted outlier tests. SLA data were transformed using a cubed root transformation, while percent cover was transformed with natural log transformations. For UMBS, we applied a log transformation to species richness. Leaf herbivory data did not fit the assumptions of normality, as they contained an excess of zeros and were over dispersed; therefore, we ran a negative binomial hurdle model using the glmmTMB package in R (Brooks et al., 2017; Appendix S1: Tables S19 and S32, Section S3). With this model, we evaluated the probability of a leaf being eaten (a binomial response), and if eaten, the amount of the leaf eaten (a truncated negative binomial response).

Data from 2015 were removed from green-up, first flower, and flowering duration analyses because the OTCs were not in place at that point. For plant composition, we calculated the average percent cover during the month with the greatest recorded percent cover (KBS: August; UMBS: July). We also calculated the average percent cover of forb and graminoid species (functional type) and native and exotic species (origin). Shannon diversity index and species richness were calculated from the plant composition data using the R package vegan (Oksanen et al., 2020). UMBS 2021 data were removed from green-up, flowering, seed set, and plant composition analyses due to infrequent data collection in that year.

We were also interested in quantifying the effect of natural temperature variation on our response variables, without consideration of our warming treatment and its effects. Therefore, we quantified the mean temperature in the ambient plots for each site per year. These temperature data were calculated independently for each response variable to match the date ranges of each variable (Appendix S1: Table S2). The models included mean temperature as a fixed effect and plot as a random effect (Response

). We did not include natural temperature variation models for SLA or biomass because we did not have at least 5 years of data compared to the other response variables. We also did not include natural temperature variation models for plant origin and plant functional type percent cover, nor species diversity metrics because we believed the overall percent cover results encompassed these metrics.

We compared hourly site-level warmed and ambient temperatures at 1 m and 10 cm above ground and at 5 cm below ground, and soil moisture at 5 cm below ground. We removed large outliers from the hourly data that were likely due to sensor malfunctions (e.g., temperatures recorded as >49 C or <-30 C). We tested for the effects of the OTCs on hourly temperature and moisture data using Welch's two-sample t-tests. At KBS, 2021 data from one set of paired sensors were removed due to a sensor malfunction, and 2021 data were removed at UMBS for one set of paired sensors due to a wasp nest covering the sensor. For 2018, 10-cm air temperature and 5-cm soil temperature data were removed at KBS due to sensor malfunctions.

RESULTS

Abiotic measurements

From 2016 to 2021, the OTCs increased 1 m air temperatures by an average of 1.9 C at KBS (t 66223 = -27.2, p < 0.001) and 3.0 C at UMBS (t 68232 = -40.2, p < 0.001) during daytime hours in the growing season (07:00-19:00, April-August; Figure 1a). The amount warmed by the chambers varied within and among years, but OTCs were consistently warmer than ambient plots at 1 m (Appendix S1: Tables S33 and S34). Air temperatures at 10 cm in the OTCs were 0.6 C cooler than ambient at KBS (t 18939 = 5.1, p < 0.001; Figure 1a). In contrast, OTCs at UMBS were 1.8 C warmer than ambient plots at 10 cm (t 68289 = -21.6, p < 0.001). Slight winter warming was also achieved, with chambers warming by 0.6 C at KBS (t 49539 = -9.7, p < 0.001) and 0.6 C at UMBS (t 52242 = -10.2, p < 0.001) from November to February. From 2016 to 2021, OTC soil temperatures at 5 cm belowground were 0.8 C cooler (t 38051 = 14.8, p < 0.001) and 0.9 C warmer (t 68371 = -16.0, p < 0.001) than ambient temperatures at KBS and UMBS, respectively, from April to August during daytime hours (07:00-19:00; Figure 1c,d). We found only small effects of warming on soil moisture at 5 cm, as moisture levels decreased by 1% (t 66885 = 20.7, p < 0.001) and 0.4% (t 67023 = 9.9, p < 0.001) in the OTCs at KBS and UMBS, respectively (Appendix S1: Figure S17).

Leaf herbivory

At both KBS and UMBS, the herbivory reduction treatment effectively reduced both the probability of a plant being eaten and the amount of leaf area eaten (Figure 2). The treatment appeared to be especially effective at UMBS, as that site contained more years with significant differences between herbivory and reduced herbivory treatments (Appendix S1: Table S32). We also found that warmed plants at UMBS typically had a greater probability of being eaten, as well as a slight increase in the amount eaten by herbivores (Figure 2b,d; Appendix S1: Tables S19 and S32). For example, regardless of herbivory treatment, warming increased the probability of a plant being eaten by 0.10-0.20 in 2018 at UMBS (herbivory: z = 4.12, p < 0.001, reduced herbivory: z = 3.89, p < 0.001; Appendix S1: Table S32). However, at KBS, we did not find a clear trend of warming effects on the probability of being eaten or the amount eaten by herbivores (Figure 2a,c; Appendix S1: Table S19). Certain plant types (e.g., native species) were more likely to be eaten than their counterparts (e.g., exotic species; Appendix S1: Section S3).

When considering the effect of natural temperature variation on herbivory at KBS, we found a temperature increase from 15 to 16 C increased the probability of a plant being eaten by 0.17 (z 2321 = -2.92, p = 0.004) but had no effect on the amount eaten (z 2321 = 0.04, p = 0.97; Appendix S1: Figure S18, Table S19). At UMBS, a temperature increase from 15 to 16 C did not have a significant effect on the probability of a plant being eaten (z 2121 = 0.45, p = 0.65), but if eaten, the amount of leaf area eaten decreased by 7.6% (z 2121 = -11.0, p < 0.001; Appendix S1: Figure S18, Table S19).

Plant phenology

Green-up

At KBS, the effect of warming on green-up depended upon the presence of herbivores (warming × herbivory interaction: F 1,24 = 3.30, p = 0.08; Appendix S1: Table S4). Overall, in plots with herbivores present, warming advanced green-up by 6.7 days (t 29 = -2.75, p = 0.01; Figure 3a). However, in plots with reduced herbivory, warming did not have a significant effect on green-up (t 29 = -0.42, p = 0.68; Figure 3a; Appendix S1: Table S20). We also found that herbivory only advanced green-up when plots were warmed (t 29 = -2.81, p = 0.009), whereas in ambient plots, there was no significant effect of herbivory on green-up (t 29 = -0.45, p = 0.66; Appendix S1: Table S20). There was a stronger effect of warming on green-up in some years compared to others (F 5,119 = 3.36, p = 0.007; Appendix S1: Table S4). For example, warming led to the advancement of green-up by 7 days in 2017 (t 171 = 1.76, p = 0.08) and 15 days in 2021 (t 172 = 3.68, p < 0.001; Appendix Table S20).

At UMBS, the effect of warming and herbivory on green-up differed between years (warming × herbivory × year interaction: F 4,96 = 4.63, p = 0.002; Appendix S1: Table S4). For example, in 2016, herbivory advanced green-up by 18.5 days, but only in ambient plots (t 131 = -4.18, p < 0.001); there was no significant effect of herbivory on green-up in warmed plots that year (t 131 = 1.26, p = 0.21; Figure 3b; Appendix S1: Table S20). However, most years did not demonstrate a significant warming or herbivory effect on green-up (Appendix S1: Table S20). We did not find significant evidence that natural temperature variation affected green-up at KBS (F 1,71 = 0.97, p = 0.33; Appendix S1: Figure S19, Table S4). At UMBS, however, for each unit increase in ambient temperature, green-up advanced by 5.2 days (F 1,48 = 17.8, p < 0.001; Appendix S1: Figure S19, Table S4).

Flowering

At KBS, the effect of warming on the date of the first flower depended on the year (F 5,118 = 4.68, p = 0.001; Appendix S1: Table S5). Warming advanced flowering by 7-10 days in 2017, 2018, 2019, and 2021, but did not have a significant effect in 2016 or 2020 (Figure 3c; Appendix S1: Table S21). There was no significant effect of herbivory on the date of first flower (F 1,24 = 0.99, p = 0.33). At UMBS, there was no significant effect of warming (F 1,24 = 0.31, p = 0.58; Figure 3d) or reduced herbivory (F 1,24 = <0.001, p = 0.99; Figure 3d) on the date of the first flower. Natural temperature variation did not significantly affect the date of first flower at KBS (F 1,99 = 1.64, p = 0.20; Appendix S1: Figure S20). However, at UMBS, for each unit increase in ambient temperature, the date of first flower advanced by 2.2 days (F 1,80 = 8.01, p = 0.006; Appendix S1: Figure S20). There was no significant effect of any treatment on the duration of flowering (Appendix S1: Section S4).

Seed set

At KBS, the effect of warming on the date of first seed set depended on the year (F 6,138 = 5.27, p = 0.002; Appendix S1: Table S7). Warming delayed seed set by 21 days in 2015 (t 191 = -4.75, p < 0.0001) and 8 days in 2019 (t 186 = -2.08, p = 0.04), but did not have a significant effect in other years (Figure 3g; Appendix S1: Table S22). The effect of warming also depended on the presence of herbivores (warming × herbivory interaction: F 1,24 = 4.43, p = 0.05; Appendix S1: Table S7). Warming delayed seed set by 8.7 days in reduced herbivory plots (t 29 = -3.25, p = 0.003; Appendix S1: Table S22), but had no significant effect in plots with herbivory present (t 29 = -1.44, p = 0.60). Similarly, reduced herbivory delayed seed set by 5.9 days in warmed plots (t 29 = -2.21, p = 0.04; Appendix S1: Table S22), but had no significant effect in ambient plots (t 29 = 0.50, p = 0.62). At UMBS, there was an interactive effect of warming and herbivory on date of seed set for some years (warming × herbivory × year interaction: F 4,96 = 3.95, p = 0.01; Appendix S1: Table S7). For example, in 2019 and 2020, warming advanced seed set in reduced herbivory plots (Figure 3h; Appendix S1: Table S22). When considering the effect of natural temperature variation on the date of first seed set at KBS, we found that for each unit increase in temperature, the date of first seed set was delayed by 4.2 days (F 1,100 = 7.55, p = 0.007; Appendix S1: Figure S22). At UMBS, for each unit increase in ambient temperature, date of first seed set advanced by 2.2 days (F 1,88 = 10.7, p = 0.002; Appendix S1: Figure S22).

Plant community composition

Percent cover

At KBS, the effect of warming on percent cover depended on year (F 6,140 = 2.94, p = 0.01; Appendix S1: Table S8). In particular, there was a significant increase in percent cover in 2020 and 2021 in warmed plots (Appendix S1: Table S23). Warmed plots had 1.4 times the percent cover of ambient plots in 2020 (t 175 = -2.1, p = 0.03) and 1.6 times the percent cover of ambient plots in 2021 (t 175 = -2.18, p = 0.001; Figure 4a; Appendix S1: Table S23). Depending on the year, reduced herbivory positively affected percent cover (F 6,140 = 3.21, p = 0.01). For example, the reduced herbivory treatment had 1.4 times the percent cover of the herbivory treatment in 2021 (t 167 = -2.32, p = 0.02), but this was the only year where a significant effect of reduced herbivory was found (Appendix S1: Table S23). Depending on the year, warming increased all four plant "types" (forb, graminoid, native, exotic; Appendix S1: Tables S9-S12, Sections S5 and S6). At UMBS, neither warming (F 1,24 = 1.71, p = 0.20) nor reduced herbivory (F 1,24 = 1.46, p = 0.24; Figure 4b) had a significant effect on percent cover. Natural temperature variation did not significantly affect percent cover at KBS (F 1,57 = 0.21, p = 0.65; Appendix S1: Figure S23). At UMBS, for each unit increase in ambient temperature, percent cover increased by 0.06% (F 1,48 = 17.36, p < 0.001; Appendix S1: Figure S23).

Plant diversity

At KBS, warming decreased species richness by an average of 1.0 species over the study period (F 1,24 = 12.89, p < 0.001; Figure 4c). Reduced herbivory did not significantly affect species richness overall (F 1,24 = 0.01, p = 0.92); instead, the effect depended on the year (F 6,140 = 4.07, p = 0.001; Figure 4c; Appendix S1: Table S13). Only in 2016 did reduced herbivory positively affect plant species richness (t 121 = -2.52, p = 0.01; Appendix S1: Table S28). At UMBS, there was no significant effect of warming (F 1,24 = 0.52, p = 0.48) or reduced herbivory (F 1,24 = 0.21, p = 0.65; Figure 4d; Appendix S1: Table S13) on species richness. For both sites, neither warming nor reduced herbivory significantly affected Shannon diversity over the study period (Appendix S1: Figure S11, Tables S14 and S29, Section S7).

Leaf traits (C, N, SLA)

At KBS, herbivory decreased nitrogen content across the two representative species by 0.10% (F 1,38 = 3.84, p = 0.06; Appendix S1: Figure S12, Table S16). Only in 2021 did warming significantly decrease N content (t 87 = 2.10, p = 0.04; Appendix S1: Table S31). At UMBS, warming only significantly decreased N content in 2018 (t 29 = 1.99, p = 0.06; Appendix S1: Table S31). There was no significant effect of herbivory reduction on N content (F 1,43 = 0.53, p = 0.47). When considering the effect of natural temperature variation on nitrogen content, increasing ambient temperatures increased N content by 0.11% at KBS (F 1,183 = 2.96, p = 0.09) and by 0.43% at UMBS (F 1,202 = 163, p < 0.001; Appendix S1: Figure S25). There was no significant effect of any treatment on leaf C content (Appendix S1: Section S8).

Warming led to an average increase in SLA of 11.2 cm 2 /g across the six representative species at KBS (F 1,116 = 4.52, p = 0.04; Appendix S1: Figure S12), but there was no significant overall effect of herbivory reduction (F 1,116 = 0.45, p = 0.50; Appendix S1: Figure S12, Table S17). SLA also increased over time (F 2,56 = 244, p < 0.001), with the highest SLA found in 2021 in the warmed and reduced herbivory treatment (Appendix S1: Figure S12). At UMBS, neither warming (F 1,88 = 0.002, p = 0.97) nor herbivory reduction (F 1,88 = 0.026, p = 0.87) had a significant effect on SLA, but SLA did decrease over time (F 2,756 = 23.8, p < 0.001; Appendix S1: Figure S12).

Biomass

At KBS, both warming and herbivory affected total plant biomass. Warming increased biomass by 30 g/0.20 m 2 (F 1,19 = 7.38, p = 0.014; Figure 5a), whereas herbivory decreased biomass by 25 g/0.20 m 2 (F 1,19 = 5.88, p = 0.026; Figure 5a; Appendix S1: Table S18). However, we did not find that the effect of warming varied significantly based on the presence or absence of herbivores (F 1,19 = 0.26, p = 0.62). At UMBS, neither warming (F 1,20 = 0.002, p = 0.97) nor herbivory reduction (F 1,20 = 0.82, p = 0.38) had a significant effect on biomass (Figure 5b).

DISC USS I ON

In this multiyear experiment, we found that plant responses to warming are largely driven by the separate effects of warming and herbivory. We found little evidence for herbivore-mediated interactive effects on plant traits and community composition. Plants were more responsive to the warming treatment at our southern site, KBS. Warming ultimately led to increased plant productivity and shifts in phenology and composition; in warmed plots at KBS, we found greater percent cover, increased biomass and SLA, earlier green-up and flowering, delayed seed-set, and reduced species richness. Aside from the lack of support for strong interactive effects, our findings generally support our hypotheses and align with past warming studies (Karimi et al., 2021;Wangchuk et al., 2021).

Previous studies have noted a positive effect of warming on percent cover, specifically with warming more strongly promoting the growth of exotic species and graminoids (Wangchuk et al., 2021;Willis et al., 2010), and often a coincident reduction in species richness due to fewer species benefiting from warmer conditions (Elmendorf et al., 2012;Morecroft et al., 2009;White et al., 2014). Those findings demonstrate that future plant communities may be dominated by more exotic and graminoid species, and fewer native and forb species. In support of Hypothesis 3, we also found that warming reduced species richness (Figure 4c), however, depending on the year at KBS, both native and exotic species and both graminoid and forb species increased under warming (Appendix S1: Figures S9 and S10). Therefore, in this system, exotic and graminoid species may not outcompete native and forb species under a new climate regime. Nonetheless, more research is needed on the effects of warming on these competitive interactions, as some species may benefit more than others (De Boeck et al., 2008), and these responses are likely more nuanced than broad exotic versus native or forb versus graminoid responses to warming. We did not find evidence that insect herbivory significantly mediated these outcomes, although herbivory did reduce percent cover and species richness in some years (Appendix S1: Figures S9 and S10).

Plant biomass at KBS was affected by both warming and herbivory reduction, but not their interaction. While warming increased plant biomass, herbivory reduced it (Figure 5a); this finding supports our Hypothesis 4. This suggests that herbivory may ameliorate warming effects on plant growth (Post & Pedersen 2008), especially in systems where herbivory also increases under warmer conditions (Hamann et al., 2021). The increased plant productivity in the OTCs at KBS may have also shaded our temperature sensors, hence leading to cooler 10-cm air temperatures and 5-cm soil temperatures in the OTC plots (Figure 1a,b).

The expansion of the growing season, due to advanced green-up and delayed seed set, may explain why plant productivity increased in warmed plots at KBS. We found support for Hypothesis 2 that warming will lead to earlier green-up, flowering, and delayed seed set at KBS (Figure 3). Interestingly, warming did not lead to significantly longer flowering duration (Figure 3e,f). In contrast, we saw the opposite effect of warming on seed set in some years at UMBS, where warming led to an advancement of the date of first seed set (Figure 3h). Because warming did not advance green-up or delay seed set at UMBS, and therefore did not expand the growing season, this might explain why there was no significant warming effect on biomass at that site (Figure 5b).

For green-up and seed set, we did find evidence for interactive effects between warming and herbivory. In general, there was potential for herbivory to advance green-up, but the effects of herbivory depended upon the site and warming treatment. At KBS, there was a strong effect of herbivory on green-up, but only in warmed plots (Figure 3a). In contrast, at UMBS, herbivory only affected green-up in ambient plots (Figure 3b). We also found that herbivory has the potential to alleviate the effects of warming on seed set at KBS; warming only delayed seed set in reduced herbivory plots, demonstrating that herbivory may be buffering the effect of warming on seed set. These findings contribute to an existing body of evidence that warming alters phenology (Peñuelas & Filella, 2001;Walther et al., 2002;Zhou et al., 2022), however, here we show that these particular phenological responses to warming depend upon the level of biotic interactions with insect herbivores. The opposing effects of herbivory on green-up at KBS versus UMBS may be due to underlying differences in species types and environmental conditions between the sites, which are discussed in more depth below. The support for Hypothesis 2 suggests that the timing of plant life cycle events may change with climate change, and this change has the potential to alter community dynamics and plant interactions with other organisms (Liu et al., 2011).

In support of Hypotheses 4, warming increased SLA and herbivory decreased N content for plants at KBS (Appendix S1: Figure S12). Interestingly, although warming decreased N content, we did not find a concurrent warming effect on leaf herbivory at KBS (Figure 2a,c). However, at UMBS, warming increased the probability of a plant being eaten, which supports Hypothesis 1 (Figure 2b), but warming ultimately did not affect leaf traits. Typically, we would expect decreased N content to lead to increased herbivory because insects require more leaf material to meet nutritional needs (Hamann et al., 2021). Our ability to capture plant responses to herbivory may have been limited by our method of measuring herbivory. For example, we only measured chewing damage on leaves for herbivory, leaving out herbivory by sap-sucking and stem-boring insects. Other forms of herbivory might lead to differing results when compared to outcomes due to insect-specific herbivory (Lebbink et al., 2023). Moreover, certain plant types (i.e., native plants and forbs) were more likely to be eaten than their counterparts (i.e., exotic plants and graminoids) (Appendix S1: Section S3). Other factors such as plant size, plant defenses, and plant relatedness can also affect herbivory levels (The Herbivory Variability Network, 2023). The differential effects of warming on different plant types suggest that some plants may be more vulnerable to herbivory under climate change, which could lead to overall changes in community composition (Primack et al., 2009;Sherry et al., 2007;Youngflesh et al., 2021). While we did not measure insect abundance and presence in this experiment, future studies would benefit from monitoring the insect community alongside the plant community to better link changes in herbivore presence with their effects on plants.

The results of our models using natural temperature variation as a continuous explanatory variable sometimes differed from our results using our warming treatment as a categorical explanatory variable (warmed vs. ambient). For example, our warming treatment models show that experimental warming had no significant effect on green-up and flowering at UMBS (Figure 3), while greater natural temperatures advanced both green-up and flowering at UMBS (Appendix S1: Figures S19 and S20). These models may differ because naturally warmer years increase temperatures at the whole community level, whereas our warming treatments warm at the plot level and likely have a larger effect on sessile organisms.

Warming at the community-level warms not just the plant community, but also affects the organisms that interact with that community, including primary and secondary consumers and plant species outside of a 1-m 2 plot. Because natural temperature variation affects communities in different ways than warming treatments, the results of climate change experiments may differ from potential future climate warming.

Throughout our analyses, we noted the high amount of both yearly and species-specific variability in plant responses to both of our treatments. This inherent variation demonstrates the importance of conducting longer-term climate studies on plant communities. For example, it may take many years for the composition of plant communities to respond to a change (Bahlai et al., 2021;Cusser et al., 2021;Dickson & Gross 2013). Dickson and Gross (2013) found N addition caused an increase in aboveground productivity within a few years, but it required 14 years for plant species richness to decrease. In our study, we also found that some response variables required multiple years to demonstrate a response; for example, we only found a treatment effect on percent cover in the final 2 years of the experiment (Figure 4a). Although our treatments were in place for seven consecutive years, we recognize that an even longer duration experiment may yield stronger plant compositional responses to warming and reduced herbivory and may be necessary to uncover the interactive effects of warming and herbivory.

We also saw that plants at the southern site (KBS) were more responsive to warming than at the northern site (UMBS). We expected the more northern site to have greater sensitivity to climate changes (Prevéy et al., 2017), however, other site contexts like plant community and soil type also play a role. In particular, the plant community differences between the sites may explain why the northern site was not as responsive to warming compared to the southern site. The most common species at UMBS is an exotic forb (Centaurea stoebe) (Appendix S1: Table S1), which may have a greater tolerance to changing temperatures compared to native species (Hahn et al., 2012). The OTCs also warmed by a much greater amount at UMBS (Figure 1b,d), and these hotter temperatures could have led to plant stress and mitigated any increased growth response to warming, which may explain why plant productivity variables (e.g., biomass, percent cover, SLA) did not change as a result of either the warming or the herbivory reduction treatments. Furthermore, the soil at UMBS differs from the soil at KBS; it is sandier and drier on average (Appendix S1: Figure S17). Previous studies have noted the importance of soil traits in mediating plant responses to climate change (Bjorkman et al., 2018;Collins et al., 2021;Elmendorf et al., 2012;Wolkovich et al., 2012), therefore the soil type at UMBS may be ameliorating some warming and herbivory effects. Future experiments could identify and test how site contexts contribute to differential responses to warming and herbivory, which would lead to a greater understanding of which environmental factors are the most important for determining community responses to warming.

These results may also depend upon our experimental design of using OTCs for warming. While OTCs are a common method for manipulating temperature in plant communities, they are not without their limitations. For example, the structure of the OTCs can limit wind and precipitation, and increase humidity (Ettinger et al., 2019;Hollister et al., 2022). The chambers themselves may also limit dispersal between plots in the community and therefore could affect plant composition and herbivory in unintended ways. Despite these potential limitations, OTCs are a well-known and effective method for manipulating the abiotic environment in situ.

Conclusions

Our study demonstrates that warming and herbivory can have strong direct effects on plant communities, but that warming and insect herbivore-mediated interactive effects may be more subtle in these early successional systems. Under current and future climate scenarios, warming is likely to affect biotic interactions alongside plant communities themselves, leading to complex responses to warming. Furthermore, the strength and direction of these effects can vary by ecological context. Thus, it is still beneficial to include levels of biotic interactions, multiple traits, and community type when studying climate change effects on plants and their communities, especially over multiple years. By including these biotic interactions in climate change experiments, we can gain a more holistic understanding of how communities may respond to a changing climate.