Climatic warming affects ecosystem-scale carbon fluxes directly through its impact on photosynthesis and respiration, and indirectly by altering the plant community.

This study reports on a 10-year factorial warming and dominant plant species removal experiment, conducted in a high- and a low-elevation montane meadow, to explore how dominant plants modify the effect of warming on the carbon cycle over time and in different locations.

At the low-elevation site, warming increased peak-season net carbon uptake in most years, primarily due to higher primary productivity. This effect was observed only in plots where the dominant species was present. Net ecosystem carbon uptake was generally positive, but it often shifted from net carbon uptake to net carbon release when the dominant species was removed, particularly in dry years. Surprisingly, the high-elevation site showed no response to the warming or plant removal treatments.

Synthesis: These findings highlight that dominant plant species can modify the impacts of warming on carbon fluxes. However, the effects of warming and plant species removal on the carbon cycle vary spatially and temporally. This study provides valuable insights into how both abiotic and biotic factors influence ecosystem carbon cycling and source-sink dynamics.

We conducted this experiment in the Gunnison National Forest in southern Colorado, USA at a low-elevation (Fig. 1a; 38.71, -106.82; 2740 m), and high-elevation (Fig. 1b;38.99, -107.06; 3460 m) site. The climate at both sites is classified as subarctic, defined by cool summers and year-round precipitation. Both sites are montane meadows and have a mean summertime temperature of 14.9 and 10.9°C, respectively, and a mean summertime precipitation of 143 and 151 mm (Prager et al., 2022). The low-elevation site is dominated by the perennial herb Wyethia × magna (Fig. 1a), which is a stable hybrid of W. amplexicaulis and W. arizonica (Weaver, 1915). Common grasses at the site are in the genera Poa and Fetusca, and common herbaceous species include Taraxacum officinale, Erigeron speciosus, and Achillea millefolium. The high-elevation site is dominated by Juncus drummondii (Fig. 1b), a perennial rush that grows in thick bunches. Other common species include Sibbaldia procumbens, Arctostaphylos uva-ursi, Senecio crassulus, and several Poa species. At the low-elevation site, we installed an electric fence to exclude cattle grazing. Early and late in the growing season, soil moisture tends to be highest at the low-elevation site, whereas soil moisture is highest at the high-elevation site during the peak growing season (Spinella et al., 2024).

In 2013, we constructed a factorial warming × dominant species removal experiment at the low- and at a high-elevation site (Prager et al., 2022). We established 32 2 × 2 m plots at each site arranged into 8 blocks, with a buffer of at least 2-m between adjacent plots. We assigned plots in each block to one of four experimental treatments: ambient temperature and dominant plant species present (ambient intact), ambient temperature and dominant removed (ambient removal), warmed temperature and dominant present (warmed intact), and warmed temperature and dominant removed (warmed removal). To warm plots, we used hexagonal open-topped chambers (OTCs), made from six pieces of transparent polycarbonate, arranged with sloping sides to give an inside diameter of 1.5 m (Hollister & Webber, 2000). We maintained the removal plots by clipping to ground level the dominant species at the beginning of each growing season and, as needed, throughout the growing season (Prager et al., 2022). In 2013, the dominant species at the low-elevation site, W. × magna, had a mean cover of 19.9% and the dominant species at the high-elevation site,  J. drummondii, had a mean cover of 19.1% (see Rewcastle et al., 2022). On average (± standard error), we removed 388±215 g of W. × magna biomass from each plot at the low-elevation site and 91.99±66.35 g of J. drummondii biomass from each plot at the high-elevation site in 2013. By comparison, only 27±5 g and 8±2 g of biomass were removed in 2023 from the low- and high-elevation sites, respectively.

To assess the impact of the treatments on C fluxes, we estimated net ecosystem C exchange (NEE) once in each plot during the peak of the growing season in 2015, 2017, 2018, 2019 (low elevation only), 2021 and 2022. In 2019 there was an avalanche that covered our high-elevation sites with snow through the growing season and we were unable to take carbon flux measurements in 2020 because of COVID-19 protocols. We collected flux measurements in late June-early July at the low-elevation and in mid to late July at the high-elevation site. NEE was measured using the protocol described by Prager et al. (2021), whereby a Licor-7500 infra-red gas analyzer (Li-Cor, Lincoln, Nebraska, USA) was used to measure CO2 concentration inside a sealed plastic chamber. Over the course of the experiment, we used several different chambers which varied in volume from 0.51m3 to 1.64m3. For details on the volume and ground surface area of each chamber, see the Supplementary Methods. For each measurement, we constructed a linear function between light and NEE, and used this to standardize NEE measurements to a common PPFD of 800 μmol m-2 s-1.

Over the course of the experiment, we used several different chambers which varied in volume and ground surface area covered. In 2015, 2017 and 2018, the chamber we used had a volume of 1.64m3 and an area of 1.43m2. The chamber we used in 2019, 2021 and 2022 had a volume of 0.512m3 and an area of 0.64m2, and the chamber we used in 2023 had a volume of 0.6m3 and an area of 1m2.

We measured CO2 concentration for 90 seconds under each of four light regimes, which were achieved by using a single- and double-layer plastic mesh coverings, each layer of which reduced the photosynthetic photon flux density (PPFD) by approximately 50%. PPFD was measured with an MQ-100 Apogee PAR meter (Apogee Instruments, Logan, UT, USA), which was mounted on top of or directly beside the CO2 sensor. The four regimes were full sun (PPFD>1000 μmol m-2 s-1), complete darkness (PPFD=0 μmol m-2 s-1), and two intermediate light levels. For each measurement, CO2 concentration was plotted against time, and any non-linear sections of data were removed. NEE (µmol·m-2·s-1) was calculated from the flux measurements using Eq. 1 (Prager et al., 2021), after converting CO2 concentrations to dry mole fractions. 

NEE=×V×(dC/dT)A

(1)

In Eq. 1, ρ is the air density in mol air/m3, calculated as P/RT, where P is the average pressure in Pa, R is the universal gas constant (8.314 J·air K-1·mol-1) and T is the mean air temperature in K. V is the chamber volume in m3, dC/dt is the slope of the concentration of CO2 in the chamber against time (µmol·CO2 mol-1·air s-1), and A is ground surface area inside the chamber in m2. In order to standardize measurements between plots measured under varying light conditions, we estimated NEE at PPFD=800 μmol m-2 s-1 (NEE800), by creating a light response curve. This was achieved by plotting NEE against measured PPFD and fitting a linear function to the four points. From this, we solved the equation for PPFD=800 μmol m-2 s-1. 

There were several flux measurements in 2015 and 2019 which did not have corresponding PPFD values recorded. For these missing values, we used averages of the PPFD measurements of other plots. In 2016, chamber temperature measurements were not recorded, so a corrected iButton temperature measurement was used instead, based on the regression between iButton and chamber temperatures from the previous year.

To determine whether treatment effects on NEE were driven more by their effects on photosynthesis or more by their effects on respiration, we estimated ER and GPP (at PPFD=800 μmol m-2 s-1) from each light response curve. Ecosystem respiration (ER) was defined as the x-intercept of the light response curve (i.e. where PPFD=0), and GPP was calculated by subtracting ER from NEE (for more details, see the Supplementary Methods).

To determine whether the treatments affected plant cover, as a proxy for biomass, we visually estimated percent cover of each species in the inner 1 m2 area of each plot each year at peak biomass, usually in early- to mid-July at the low-elevation site and late-July to early-August at the high-elevation site. For each plot, we calculated total cover and relative cover of the dominant plant species.

From 2015 – 2021, we measured volumetric soil moisture (0-15 cm) in each plot at the peak of the growing season. We randomly collected three measurements per plot using a TDR probe (Soilmoisture Equipment Corp, California, USA), and averaged them to calculate a plot-level estimate. From 2015-2021, the mean summer air temperature in each plot was recorded using iButtons (Maxim Integrated Corp, San Jose, California, USA) 5 cm above the soil surface (Spinella et al., 2024). We calculated a daytime average for each plot after removing outliers, and we defined mean summer air temperature as the mean of the daily measurements in the 10 days preceding the C flux measurements. In 2022, we installed TMS probes (TMS-4 standard; TOMST, Wilmington, USA) in each plot and logged soil moisture and air temperature 6 cm above the soil surface every 30 minutes. We averaged daily maximum air temperature across the week leading up to the flux measurements at each site in each year, as we determined, via linear models, that this was the best predictor of NEE.