Vascular plant species response to warming and elevated carbon dioxide in a boreal peatland

Arctic and boreal peatlands form in regions with cool, wet climates where annual rates of primary productivity are greater than rates of ecosystem respiration. Over millennia, this has led to a significant store of organic carbon (Gorham 1991, Yu et al 2010, Treat et al 2019). It is estimated that one third of the global pool of soil organic carbon is stored in peatlands (Gorham 1991, Tarnocai et al 2009). As rapid warming changes the climate of the north, rates of carbon cycling are expected to change dramatically (Dise 2009). While studies have shown an increase in aboveground net primary productivity (NPP) with warming (Weltzin et al 2000, Boelman et al 2003, 2005), soil respiration rates are also likely to increase (Davidson and Janssens 2006, Fenner and Freeman 2011). The composition of vegetation communities plays an important role in determining the carbon balance of peatlands as either sinks or sources of carbon to the atmosphere. For example, shifts from dominance by graminoid species to woody plants, or changes in the relative abundance of vascular plants versus Sphagnum mosses may substantially alter rates of carbon sequestration in peatlands (Potvin et al 2015, Norby et al 2019). Characterizing how warming will affect peatland plant community is important to predicting the effects of climate change on terrestrial carbon cycling.

Experimental studies on the effects of climate change on peatlands have yielded a range of species responses depending on initial community composition and hydrology (i.e. nutrient-rich fens versus nutrient-poor bogs), and the type of experimental manipulation. Studies found that a direct effect of warming was to increase shrub presence in both bogs and fens, coincident with a decrease in non-vascular moss and graminoid cover (Weltzin et al 2000, 2003, Churchill et al 2015). Warming has also been shown to have a universally negative effect on forbs and non-vascular Sphagnum mosses (Weltzin et al 2000, 2003, Walker et al 2006, Dieleman et al 2015, Norby et al 2019). Fewer studies have tested the effects of elevated CO₂ (eCO₂) in peatlands, but some studies have demonstrated an increase in graminoid abundance in response to eCO₂ (Berendse et al 2001, Dieleman et al 2015). This result is consistent with other studies in graminoid-dominated ecosystems that have shown a positive effect of eCO₂ on the growth of grasses and sedges (Owensby et al 1999, Dieleman et al 2015, Mueller et al 2016).

Vascular species composition affects ecosystem productivity, and is therefore an important aspect of the peatland carbon cycle. Increases in both shrub and sedge cover have been found to correlate with greater estimates of aboveground NPP (Boelman et al 2003, Ward et al 2013, McPartland et al 2018). Vegetation manipulation studies that have explicitly tested the relationship between plant functional types (PFTs) and peatland carbon cycling have demonstrated that PFT cover is a strong predictor of the CO₂ sink strength of peatlands (Ward et al 2009, Kuiper et al 2014, Potvin et al 2015, Goud et al 2018). In some cases, the carbon sink strength of peatlands has been shown to increase with a greater cover of evergreen and woody shrub species (Churchill et al 2015, Ward et al 2013). However, the increased sink potential of woody plants may be outweighed by the detrimental effects of increased shrub cover on graminoid species (Goud et al 2018), and Sphagnum mosses (Potvin et al 2015, Norby et al 2019). Together, these studies indicate the sensitivity of peatland carbon sequestration to climate warming.

Northern ecosystems, and in particular nutrient-poor bogs have extremely low levels of species diversity compared with other terrestrial ecosystems (Gaston 2000). Consequently, the loss of an individual species from a bog may indicate the loss of an entire plant functional group. Loss of functional diversity in other ecosystems has been shown to negatively affect primary productivity and ecosystem function (Tilman et al 1996, 1997, Reich et al 2012). Given that warming has already been shown to affect peatland species diversity (Walker et al 2006, McPartland et al 2019), studying diversity is an important aspect of the overall picture of peatland ecological response to global change.

This paper reports findings from the first five years of the Spruce and Peatland Responses Under Changing Environments (SPRUCE) experiment on the vascular plant community change. SPRUCE is a whole-ecosystem warming and eCO₂ experiment located in northern Minnesota, USA. Since 2015, a suite of ecological, physiological, and biogeochemical processes have been tracked in response to a wide range of temperature and CO₂ treatments. We monitored species cover, leaf area, and standing biomass of the ground and shrub layers annually. Our research was guided by the following hypotheses: (1) aboveground net primary productivity (ANPP) would increase with warming, coincident with an increase in shrub cover and leaf area index (LAI), (2) abundance of graminoids would respond positively to warming and eCO₂, and (3) forb cover would decrease as a result of shading by shrubs. We also predicted a decrease in species and functional diversity with warming and eCO₂ as dominance by specific PFTs increased and some species became less common.

The SPRUCE project is located at the USDA Forest Service's Marcell Experimental Forest (47.30°N, 93.29°W), within the Chippewa National Forest of Minnesota (figure 1). The Marcell Experimental Forest has been a research site of the USDA Forest Service since the early 1960s (Bay 1967, Verry and Timmons 1982, Harriss et al 1985). The mean annual air temperature from 1961 to 2005 was 3.3 °C and average annual precipitation was 768 mm. Mean annual temperatures have increased approximately 0.4 °C per decade since the 1980s (Kolka et al 2011). The SPRUCE experiment was designed following a regression‐based approach in which ecological changes could be measured in response to a broad range of increasing temperature treatments (Hanson et al 2017). The treatments are intended to explore the effects of increased above‐ and belowground temperatures and eCO₂ levels on ecological and biogeochemical processes at the southern edge of the boreal biome (Tfaily et al 2014). Belowground warming began in June 2014, aboveground warming started in August 2015, and eCO₂ treatments were initiated in spring of 2016. All subsequent years represent the whole-ecosystem warming and eCO₂ treatments.

Zoom In Zoom Out Reset image size

SPRUCE is located in an ombrotrophic bog dominated by Picea mariana (Mill.) B.S.P. (black spruce) and Larix laricina (Du Roi) K. Koch (tamarack), with an understory of ericacous shrubs, graminoids, and forbs. The experiment consists of ten large octagonal enclosures that are 12.8 m in diameter (114.8 m² in area), 8 m tall, with a belowground corral that penetrates through the peat to mineral material below, which averages 2.9 m deep across the experimental plots. The enclosures are maintained at a series of increasing temperatures by forced air warming and deep peat resistance heating (Hanson et al 2017). The treatments consist of five temperatures (+0, +2.5, +4.5, +6.75 and +9 °C) and two CO₂ levels (ambient and ~900 ppm). The temperature treatments are implemented continuously by maintaining a fixed differential from ambient temperatures in the two +0 °C reference plots. Mean air temperature measurements are made every half hour with thermistors (model HMP-155; Vaisala, Inc.) installed at the center of each plot at 0.5, 1, 2 and 4 m above the surface of the peat. In addition, we monitored three unchambered control plots, also 12.8 m in diameter, without the enclosed structure but with similar construction disturbances.

We measured the vegetation composition at SPRUCE from 2014 to 2019 between late June and mid-July of each year. Three 2 m² subplots were established in 2014 within each treatment plot by driving PVC stakes into the peat to mark the area for repeated measurement. Each subplot was placed to be representative of bog microtopography and included areas with both hummocks and hollows. A 1 × 2 m rectangular frame with 50 grid cells 20 × 20 cm in size was positioned on the PVC pipes above the shrub layer, and the presence and absence of all species in each grid cell was recorded. If a species was present in half of the 50 grid cells, it was considered to occupy 50% of the plot. We acknowledge that this method overestimates true cover because a small species if present in every cell would be measured at 100% cover despite not physically covering the entire plot. We classified each species into PFT for analysis of cover change (table S1 (available online at stacks.iop.org/ERL/12/124066/mmedia)).

Shrub layer growth assessments were obtained annually by destructive harvest of two small 0.25 m² plots within the chambers, or immediately outside of the unchambered control enclosures, separate from the survey plots. Sampling occurred annually from 2015 to 2019. Clipped plots were sampled at the peak of annual net primary production in the middle of August. In 2018 the number of plots per chamber was doubled to four. All vegetation above the peatland moss layer was clipped at the surface of the moss and placed in plastic bags for sorting. Sorting included the separation by species, followed by a separation of the tissues into current-year vs. older tissues. Current-year ANPP for each shrub species were evaluated via separation of all new stem and foliar growth from prior year tissues. Incremental changes in stem diameter were not accounted for. Graminoid and forb ANPP were equal to the total amount of species-specific tissues. The cumulative mass of current-year growth for leaves and stem tissues was interpreted as vascular species ANPP for the designated plot. Older tissues are included in an overall estimate of aboveground carbon stored in vascular plant biomass.

LAI was measured at SPRUCE from 2017 to 2019 using an LI‐COR LAI 2200 (Lincoln, NE) instrument that estimated leaf area per unit ground area as the amount of incoming light available below a canopy (Gower and Norman 1991). Measurements were done near mid‐day on a cloud‐free day. The reference sensor was positioned on a tripod in an open meadow with no nearby tree cover. Measurements were taken at two central sampling locations in each of the three vegetation survey plots in each chamber away from disturbances. At each sampling location, two scans were taken each at 1 m above the shrub canopy and beneath the shrub canopy at the surface of the peat.

We used the vegetation survey data to analyze changes in community composition using non-metric multidimensional scaling (NMDS) ordination using a Bray-Curtis dissimilarity index (Minchin 1987). NMDS examines changes in community composition based solely on observations of species composition without reference to the experimental units. Ordination was performed using subplot-level vegetation survey data on all species from 2014 to 2019, and chamber-level mean values were extracted from the ordination results for plotting. This was done to capture within-chamber variation. The relationship of the experimental units to species compositional change was analyzed on subplot-level data using post-hoc permutation testing (Legendre and Anderson 1999). Permutation tests resample the existing data against a random distribution to determine the statistical significance of the predictor variables. For the permutation testing, we used the experimental temperature differential values (i.e. +0, +2.25, +4.5, +6.75 and +9 C), rather than actual measurements of temperature because we lacked pre-treatment continuous temperature data for 2014. An alpha-level of 0.05 was used to determine whether treatment was a statistically significant predictor of variation within the NMDS.

We calculated subplot-level species diversity using the Shannon diversity index, which incorporates both species richness and evenness into a single measure of alpha diversity:

$$\begin{equation}H{^{^{\prime}}} = \, - \mathop \sum \limits_{i = 1}^s {p_i}\log \left( {{p_i}} \right)\end{equation} \tag{ 1 }$$

In this equation, Shannon diversity (H') is calculated at the subplot-level in which p_i represents the proportion of the population represented by species i (Gotelli and Ellison 2013). NMDS, permutation testing, and analysis of Shannon diversity were all performed using the R Vegan package for community ecology (Oksanen et al 2018).

We employed a linear mixed effects model comparison framework to examine the responses of plant abundance (for both species and PFT percent cover and biomass), LAI, and Shannon diversity to warming and eCO₂. We used the 2015–2019 survey and ANPP data for this analysis. We analyzed cover change for shrubs, graminoids, and forbs. We used growing-season temperature data collected by the sensor located in each chamber at 0.5 m above the peat surface. We calculated average chamber temperatures for the entire growing season, beginning in mid-April and extending through to the end of October (DOY 110–305). To determine the appropriate growing season window to use in our analysis we visually examined phenology photo series captured daily in each chamber by Richardson et al (2018), to establish when spring greening and autumn senescence occurred. Temperature was added to our model as a continuous variable, and CO₂ level was included as a factor. We structured the model to examine the interaction between the predictors in addition to the effect of each individually. To isolate the effects of experimental manipulation and control for change over time, we included year as a random effect (equation (2)):

$$\begin{equation}y = {\beta _{Temp\,}} + \,{\beta _{CO2}} + \,{\beta _{Temp\,}}*\,{\beta _{CO2}} + \,{u_i} + \,{\varepsilon _i}\end{equation} \tag{ 2 }$$

Where y represents change in the abundance of a particular species or PFT, measured either using our vegetation survey or destructive harvest method, ${\beta _{Temp\,}}\,$represents a continuous temperature variable, ${\beta _{CO2}}$ represents a binary factor of CO₂ level (ambient, or ~900 ppm), $\,{u_i}$ represents a random effect associated with sample year, and ${\varepsilon _i}$ is an error term. We compared models by sequentially removing the fixed effects for either temperature or CO₂ to determine the relative explanatory power of each predictor variable. We report only the results for the best model, determined by the Akaike Information Criterion (AIC). We used marginal and conditional R² values to report the amount of variance explained by the fixed effects and full model, respectively. We used change in AIC (ΔAIC) to assess the increases in model parsimony associated with omitting fixed effects in order to determine the best model. We also performed basic significance testing using the Satterthwaite's method for calculating p-values for fixed effects in linear mixed models (Kuznetsova et al 2017). All modeling was done at the chamber level using the average of the replicate subplots. Linear mixed effects modeling was done in R using the lme4 package (Bates et al 2015). In presenting our results, we opted to fit regression lines to the models that had a marginal R² value of greater than 0.2, indicating that the best model captured a reasonable amount of variance in the data.

The results of our NMDS analysis indicate a shift in species composition over the course of the experiment (figure 2). The ordination had a stress of 0.22 for two dimensions, indicating that the variance within the dataset was captured reasonably well. Results from the permutation testing indicate that temperature was strongly predictive of variation in species composition, but that there was a large amount of variation not accounted for by temperature (R² = 0.18, p < 0.001). Carbon dioxide level was not significantly correlated with compositional differences (R² = 0.01, p = 0.08). Sample year also emerged as a significant predictor (R² = 0.11, p < 0.001). Differences in species composition appeared slight until 2017 and increased substantially between 2017 and 2019 (figure 2). Analysis of changes in plot-level Shannon diversity indices indicated that species diversity is strongly negatively correlated with warming, which was cumulative across years (figure 3). The best model included only the effect of temperature (table 1, table S2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We analyzed the response of PFTs to warming and elevated CO₂ for both the vegetation survey and harvest datasets (table S1) (figure 4, table 1). For the sake of comparison, we left V. oxycoccos out of our survey dataset because it was not represented in the harvest dataset, but we did analyze it separately (See section 3.3). The harvest dataset also included only the forb Maianthemum trifolium. The survey dataset included observations of several forbs, but M. trifolium represented over 97% of forb observations, so the forb survey data are essentially also a characterization of cover change of that species.

Zoom In Zoom Out Reset image size

Shrub cover did not respond to treatment, but shrub ANPP increased with temperature (figure 4, top panel). The two datasets have contradictory results for graminoid abundance—the survey data show a weak decline, where the ANPP data shows a slight increase in graminoid biomass (figure 4, middle panel). The two datasets agree that both forb cover and ANPP have declined with warming (figure 4, bottom panel). Model fit was strongest for forbs in both datasets, and relatively weak for both shrubs and graminiods (table 1). In most cases, the most parsimonious model included only the fixed effect of temperature treatment (table 1, table S2).

We examined the individual responses of four common species—Rhododendron groenlandicum, Chamaedaphne calyculata, Vaccinium angustifolium, and Vaccinium oxycoccos (figure 5, table 1). V. oxycoccos was only represented in the survey data because it is small creeping species that did not contribute a significant amount of ANPP during the sample period. The survey data for R. groenlandicum and C. calyculata indicate no change over the study period, whereas the ANPP data show a weak increase and weak decrease for each of those species, respectively. The two datasets both show a marked increase in both cover and biomass increment of V. angustifolium, and the survey data shows a marked decrease in the abundance of V. oxycoccos (figure 5). In almost all cases, the most parsimonious model included only the fixed effect of temperature treatment, but CO₂ level emerged in some species such as R. groenlandicum that only exhibited weak trends, but those models were not significant (table 1, table S2). Those models captured very little variance in the data, and further visual examination of those cases did not reveal strong patterns associated with CO₂ level.

Zoom In Zoom Out Reset image size

In order to assess entire vascular plant community response to treatment we measured LAI and total aboveground carbon stocks, including all biomass and foliage for the ground layer (figure 6). LAI was measured in 2017–2019, and carbon mass data were available from 2015 to 2019. Both LAI and overall aboveground carbon stores showed an increase with warming, and no effect of CO₂, as determined via model selection (table 1, table S2).

Zoom In Zoom Out Reset image size

After five years of experimental manipulation of temperature and CO₂ levels at SPRUCE, we observed a series of changes to the composition, diversity, and productivity of the vascular plant community. Our results indicate that these changes are driven primarily by temperature, with only a minor effect of eCO₂ apparent at present. The two most striking findings from our study were the increase in aboveground vascular plant cover and the sharp decrease in species diversity. The increase in ANPP was driven by a few shrub species, while the decrease in diversity is the result of the elimination of uncommon species, and species that grow on the Sphagnum surface, in particular forbs. The decrease in species of low stature, may be due to a combination of the direct (warming) and indirect (shading) effects (Strauss 1991). Unlike in other types of ecosystems where diversity loss is associated with decreases in productivity (Tilman et al 1996, 1997), at SPRUCE diversity loss occurred despite an increase in NPP. The loss of diversity in peatlands is of concern given the already low levels of species richness present in these ecosystems, and the possibility that the loss of a single species may also represent the loss of an entire PFT (Reich et al 2012).

Shrubs are the most abundant PFT at SPRUCE, accounting for 65% of all survey observations. Shrubs showed the greatest overall increase in biomass, although this as not the case for all shrub species. Shrubs play a significant role in the carbon cycling of arctic and boreal ecosystems, including in peatlands. Numerous observational and experimental studies have documented increases in shrub presence across the far north (Elmendorf et al 2012, Bjorkman et al 2018, Myers-Smith et al 2020). In boreal and arctic ecosystems, shrubs are becoming more dominant as a direct response to warming (Walker et al 2006, Myers-Smith et al 2011, 2015). The increase in vascular species has been shown to alter ecosystem function by speeding up rates of carbon and nitrogen cycling (Hobbie 1996, Myers-Smith et al 2011, Mekonnen et al 2018). The shift to shrub dominance could lead to an increase in aboveground NPP, but also increase the overall rate of carbon cycling and contribute to a decrease in the long-term carbon storage of these ecosystems as rates of decomposition rise (Chapin et al 1996). These studies, which range from the high arctic to our study site at the southern edge of the boreal biome suggest widespread shift towards dominance by vascular woody plants and away from non-vascular species of mosses and bryophytes driven by climate warming.

Individual species dynamics were reflected somewhat differently by the survey and harvest datasets demonstrating the sensitivity of our assessment to the method of data collection. For example, the survey and ANPP datasets were not consistent regarding changes in graminoid abundance. Although neither response was particularly strong, the survey data indicated a weak decrease in abundance, whereas the harvest dataset indicated an increase in graminoid ANPP. These differences could be due to changes in leaf physiology due to warming and eCO₂ driving an increase in leaf mass per area which would be captured in biomass but not as percent cover. The two datasets also offer slightly different pictures of change in shrub abundance. In the survey dataset, shrub PFT cover is flat during the study period, whereas the biomass data show a marked increase in shrub ANPP. A likely explanation is that the survey data becomes saturated if a species is present in 100% of cells, regardless of its biomass. Therefore our method was unable to capture changes in cover in species that were very abundant at the beginning of the study.

Our results showed only minor effects of CO₂ treatment, and we are thus unable to draw conclusions regarding its effects on the vascular plant community at this time. Although our data were noisy, our model selection indicated that the best mixed effects models included only the term for temperature. In handful of cases, the best model included both terms, or included only the effect of CO₂. However, when we visually examined those cases, there was no separation evident between plots with ambient and eCO₂, and those models captured very low levels of variance overall. Further, there was only one instance where CO₂ emerged as a significant model predictor. Our findings are consistent with other recent studies from the SPRUCE project. Malhotra et al (2020) examined the belowground component of vascular productivity, and found no significant CO₂ effect despite large overall increases in fine-root growth associated with temperature. Hanson et al (2020) similarly found only a minor effect of CO₂ in a whole-ecosystem assessment of the SPRUCE carbon budget. Ward et al (2019) measured leaf-level physiology in the two most common shrub species (R. groenlandicum and C. calyculata) and did not observe significant gains in photosynthetic capacity. In the case of C. calyculata they observed a decrease in net photosynthesis in plants growing in eCO₂ treatments. These studies, along with our results, suggest that effects of eCO₂ on the structure and carbon budget of the ecosystem have been only marginally observed. Elevated CO₂ has been shown to increase the water-use efficiency in vascular plants (Chapin III et al 2011), particularly in xeric ecosystems such as grasslands (Owensby 1999, Morgan et al 2011). Despite evidence for drying at the soil surface (Norby et al 2019) SPRUCE remains a wetland site with no lack of water availability for deeply-rooted species such as the ones in our study. This could partially account for the lack of CO₂ effect at present.

The results of our study indicate an overall increase in vascular ANPP that should be considered in the context of changes to the whole ecosystem. Research by Hanson et al (2020) showed that warming is causing dramatic losses of carbon primarily as a result of increased soil respiration. These losses are many times greater than what is being offset by increased vascular plant uptake. In addition, results by Norby et al (2019) indicate that the crucial Sphagnum mosses that sequester the majority of carbon in peat soils are in decline as a result of warming, drying, and shading by shrubs, further reducing the carbon sink potential of the SPRUCE peatland. Although SPRUCE is located on the southern edge of the boreal biome, it is representative of peatlands around the circumpolar north (Tarnocai et al 2009). The total loss of carbon from this ecosystem even as warming stimulates plant growth corroborates many observational and experimental studies that have documented increases in plant growth at the same time as millennia-old carbon is being lost from northern ecosystems as the climate warms (Myneni et al 1997, Tape et al 2006, Schuur et al 2008, Myers-Smith et al 2020).