Increased biocrust cover and activity in the highlands of Iceland after five growing seasons of experimental warming

Background and Aims

One of the most important questions of our time is how ecosystems will be transformed by climate change. Here, we used a five-year field experiment to investigate the effects of climate warming on the cover and function of a sub-Arctic alpine ecosystem in the highlands of Iceland dominated by biocrust, mosses and vascular plants.

Methods

We used Open Top Chambers (OTCs) to simulate warming; standard surface and NDVI analyses to measure plant cover and function; gas analyzers to monitor biocrust respiration; and the Tea Bag Index approach to estimate mass loss, decomposition and soil carbon stabilization rates.

Results

Contrary to our initial hypothesis of warming accelerating an ecological succession of plants growing on biocrust, we observed a warming-induced decreased abundance of vascular plants and mosses —possibly caused by high temperature summer peaks that resemble heat waves— and an increase in the cover of biocrust. The functional responses of biocrust to warming, including increased litter mass loss and respiration rates and a lower soil carbon stabilization rates, may suggest climate-driven depletion of soil nutrients in the future.

Conclusion

It remains to be studied how the effects of warming on biocrusts from high northern regions could interact with other drivers of ecosystem change, such as grazing; and if in the long-term global change could favor the growth of vascular plants on biocrust in the highlands of Iceland and similar ecosystems. For the moment, our experiment points to a warming-induced increase in the cover and activity of biocrust.

Methods
Experimental design and abiotic measurements

In June 2018 we set up the Climate Research Unit at Subarctic Temperatures (CRUST) experiment (Fig. S1b; Salazar et al., 2022) in the highlands of Iceland (64°02' N, 19°13' W; 590 m.a.s.l.). The site is primarily covered with Anthelia juratzkana biocrust, mosses, and vascular plants (e.g. Salix herbacea; Fig. S1a), on soils classified as andosol Vitricryands (US Soil Taxonomy) in a cryic soil temperature regime with an abundance of relatively fresh volcanic glass of basaltic and andesitic composition (vitric materials). The area is subjected to periodic volcanic activity and intense (0.1-2 mm yr-1) dust deposition (Arnalds et al., 2015). The soils are coarse grained (sandy loam/loamy sand) with coarse sub-surface tephra layers (volcanic ash) that negatively affect soil water conductivity and storage capacity. Mean annual temperature and precipitation (1971-2000; Icelandic Met Office, 2024) at the site are ca 1.6 °C (Fig. S10) and 1600 mm, respectively. The snow-free season generally starts in late June and ends in October. Flooding is common in June, when the snow melts and the ice under the surface is still frozen. Anthelia juratzkana biocrust is often found in this type of late-snow bed ecosystems (Belland, 1983; Ottósson et al., 2016; Smáradóttir, 2020).

Experimental warming was simulated with OTCs (Hollister et al., 2023) that were built according to the protocols of the International Tundra Experiment (ITEX; Henry and Molau, 1997; Henry et al., 2022). We used a randomized block design. First, we visually selected 16 areas of at least 4x4 meters covered mainly by biocrust and separated by at least 10 meters. Then we randomly selected 8 areas (blocks, n=8). In each of these blocks, we set up two 1.5 m x 1.5 m experimental plots and randomly assigned treatments to them, an OTC and a control (i.e. ambient temperature) plot. The OTCs were left in place all year round.

We measured temperature, light intensity and moisture in OTCs and control plots at different times during five consecutive growing seasons from 2018 to 2022. We measured surface temperature and light intensity every 2 hours using HOBO pendant temperature/light data loggers UA-002 64 (Onset Computer, Bourne, Massachusetts, USA). We used one logger per treatment (OTC and control) in 4 of the 8 blocks (n=4). In the growing seasons of 2018, 2019, 2020 and 2022 we took paired measurements of soil moisture (2-3 cm depth) in all blocks (n=8) using an ECH2O EC-5 moisture sensor (Decagon Devices, Pullman, WA) attached to a HOBO micro station (Onset Computer, Bourne, Massachusetts, USA).

Aboveground cover and plant community analysis

We estimated percentages of area covered by biocrust, moss and vascular plants in the growing seasons of 2018, 2019, 2021 and 2022. For this, we placed a 50 x 50 cm quadrant, divided in 16 squares of 12.5 x 12.5 cm, in the center of each plot and visually estimated the cover of biocrust, moss and vascular plants (example in Fig. S2). In 2022, we complemented our cover analysis with: 1) Groundbase (SKYE instrument SpectroSense 2+) measurements of Normalised Difference Vegetation Index (NDVI) as proxy for primary productivity. 2) A detailed plant community analysis using the point intercept method with 25 evenly distributed points within the 50x50 cm quadrant where all species intercepted (hits) were recorded for each point along with visual estimates of total vascular plant, moss and litter cover as well as a complete list of vascular plant and moss species providing plot species richness.

Litter decomposition and stabilization factor

To compare the potential fate of organic matter and decomposition rates in the soil, we buried green (Lipton, EAN: 87 22700 05552 5) and rooibos (Lipton, EAN: 87 22700 18843 8) tea bags in OTC and control plots. Tea bags were weighted and buried in the field on July 27th, 2018 and June 25th 2021. Bags were buried at two depths: at 8 cm, to facilitate comparisons with other studies following the protocol proposed by Keuskamp et al. (2013); and under the biocrust (2-3 cm), which is more relevant to our particular study. We buried 2 subsamples per plot and per depth. The bags were collected on June 5th, 2019 (313 d of incubation) and August 10th, 2022 (411 d incubation), respectively. Notice that the tea bags were incubated for ca. 1 yr (similar to other tundra studies, e.g., Björnsdóttir et al., 2021; von Oppen et al. 2024) and not for 90 days (as in Keuskamp et al., 2013). From the measured mass loss, we calculated litter decomposition rates (k) and stabilization factor (S; an estimate of the fraction of labile compounds that stabilize and become recalcitrant during decomposition) as in Keuskamp et al., (2013). Our decomposition data was submitted to the TBI network: http://www.teatime4science.org/.

Biocrust respiration

We measured biocrust respiration (RB) in the growing seasons of 2019 (June and July), 2021 (June and July) and 2022 (June and August). For that, we installed PVC collars (10 cm depth and diameter) in Anthelia juratzkana biocrust patches within the OTCs and control plots, and measured CO2 fluxes periodically using standard procedures with gas analyzers attached to dark soil chambers. All measurements were done approximately in the middle of the day, between 11:00 and 13:00 UTC. We used different gas analyzers in different years. In 2019, we used a Li-6400XT portable system (Li-COR Inc., Lincoln NE, USA); in 2021 a Senseair K-33 ELG sensor (based on a design by Harmon et al., 2015), calibrated with an EGM-4 (PP Systems, Amesbury, Massachusetts, USA); and in 2022 we used an EGM-5 (PP Systems, Amesbury, Massachusetts, USA).

Statistical analysis

We used a combination of statistical approaches: 1) Linear regression (lm function in the R package) to analyze climate trends near our study site (hourly measurements of air temperature at the Vatnsfell station of the Icelandic Meteorological Office, IMO), and the potential relationships between surface light and OTC-warming. 2) Mixed effect models (lmer function), selected based on the Bayesian information criterion (BIC), to analyze differences in litter mass loss, decomposition rates, stabilization factors, and RB. 3) Analysis of variances (aov function) for cover, soil moisture, and NDVI analyses. Surface cover and moisture analyses were further complemented with Tukey Honest Significant Differences (HSD) analyses (TukeyHSD function), to explore potential differences between control and OTC plots in individual years. In the mixed effect models, we analyzed the fixed effects of the OTCs, time and/or tea bag type (depending on the case) and the random effects of the blocks (paired OTC and control plots). To select the mixed effects models with the lowest BICs (i.e. best fit), we used the glmulti package (Calcagno and de Mazancourt, 2010). All statistical work was conducted in R, version 4.1.2 (R Core Team, 2023).