Insect herbivory dampens Subarctic birch forest C sink response to warming

Field site and treatment plots

The experiment was established in the mountain birch forest-tundra ecotone at the Turku University Kevo Subarctic Research Institute (69°45.4′N, 27°00.5′E; altitude 104 m a.s.l., mean annual temperature −1.3 °C, mean annual precipitation 354 mm in reference years 1981–2010; data from Finnish Meteorological Institute) in 2016. Adult mountain birches were trimmed in the experimental area, but field layer vegetation was left intact and twenty experimental plots (each 0.75 m × 1 m) with no mountain birch stumps were established. To include Subarctic woody species in the experiment and to have controlled and properly replicated plant material to accurately follow treatment responses in plants, plantlets of four birch species (local B. nana and B. pubescens subsp. czerepanovii and currently slightly more southern B. pendula and B. pubescens) were planted, 20 cm apart, amongst the field layer vegetation in each plot (12 plantlets in each plot, three genotypes per species).

Plant material and plot properties

Birch material was collected from natural Subarctic populations growing between 67°43′N and 69°01′N with mean annual temperatures varying from −1.2 to −1.9 °C (reference years 1981–2010; data from Finnish Meteorological Institute). The plantlets were cloned from the collected twig material using the micropropagation technique at the Haapastensyrjä Unit of the Natural Resources Institute Finland (Luke) at the end of 2015, potted into nursery peat (Kekkilä OPM 420 W) in May 2016 and planted to the field site in mid-July 2016. Mean plantlet height was 6.2 cm (SD 4.5, range 1–29 cm) at the time of planting and 15.6 cm (SD 11.0, range 1.5–66 cm) at the end of the growing season 2018. Birches can produce thousands of seeds m⁻² in Subarctic areas⁴⁸ and our planting density of 16 plantlets m⁻² is well within the natural variation of density in young birch seedling patches⁴⁹. Birches were chosen as test species due to their high abundance in Subarctic Fennoscandia.

Variation in soil and vegetation properties among the experimental plots was recorded after planting birches. Two to three soil cores (diameter 3 cm, depth 10 cm, or to the nearest rock) were taken from each plot, the samples were pooled, dried (70 °C, 48 h), and the organic matter (OM) content was measured as a loss on ignition (550 °C, 4 h). Lichens and plants were identified to species or genus level and their areal cover was visually estimated for the area inside the flux chamber collar (55 cm × 55 cm). Estimates of soil OM (range 9–38% of soil dry mass among the plots) and the cover of vascular plants (29–83% of area covered; Empetrum nigrum, Vaccinium uliginosum, and Vaccinium vitis-idaea as common species; cover of planted birches was on average 7% of the cover of all vascular plants), mosses (3–85%, Pleurozium schreberi and Hylocomium splendens) and lichens (0–23%, Nephroma arcticum and Cladonia arbuscula) were later used as covariates in the statistical models to control the effects of soil and plant plot-to-plot variation on response variables.

Experimental design and treatments

Treatments included two levels of warming (ambient, +3 °C) and two levels of insect herbivory (natural, reduced) in a fully factorial 2 × 2 design. In the field, the 20 treatment plots were first divided into five replicate blocks and then in each block, the four treatment combinations were randomly allocated to the plots. The distance between adjacent plots was >1.2 m, which was sufficient to avoid thermal effects between warmed and control plots (ensured using a thermal camera Flir E8; Flir Systems AB, USA).

Warming was realized using two 240 mm × 60 mm ceramic heaters (Elstein-Werk M. Steinmetz GmbH & Co., Germany; dummy heaters in ambient plots), installed 80 cm above the ground and controlled using a microprocessor-based feedback system. The feedback system maintained a fixed temperature differential between the ambient and heated plots using real-time temperature data from Pt100 sensors (Gräff GmbH, Germany) attached to a RMD680 multichannel transmitter (Nokeval Oy, Finland). The sensors were installed 30 cm above the ground (i.e., on the top of the vegetation layer) within green metal plates, mimicking plant leaf surfaces, in five ambient and five heated plots. The ceramic heaters warm surfaces more than the air, so using the temperature of the green metal plates for controlling heating we ensured that vegetation was not overheated. Warming was turned on in spring when the snow-cover melted below 20 cm, and turned off in late autumn when mean daily air temperatures in ambient plots remained permanently below −3 °C. While warming was on, air temperatures were continuously recorded using Pt100 sensors installed under white plastic plates 30 cm above the ground in three control and three warmed plots. Soil temperature and moisture were regularly, but not continuously, measured in 3–5 spots in each plot at the depth of 5 cm using SM150T soil moisture kit (Delta-T Devices Ltd., UK) and Testo 735 thermometer (Testo SE & Co., Germany) attached to a Pt100 sensor. On average, the plates mimicking plant leaves were 3.3 °C, the air 2.3 °C and the soil 1.2 °C warmer in heated than ambient plots across the warming periods in 2016–2018 (in 2016, warming did not cover the entire growing season but was turned on in early July).

The herbivore treatment was started in the beginning of the 2017 growing season and was accomplished by spraying the herbivore reduction plots with 0.1% solution of synthetic pyrethrin (Decis EC25, Bayer Crop-Science, Germany) and the control plots with tap water weekly, using two portable garden sprayers and a protective tent to eliminate wind drift. Decis EC25 has not found to have side-effects on plant growth or chemistry⁵⁰ and although deltamethrin, the active ingredient, contains N, the quantities of N that could enter soil during sprayings are negligible. In our experimental plots, the instantaneous top soil (0–10 cm) mineral N (sum of NH₄⁺–N and NO₃⁻–N) availability is, on average, 1.5 µg N per g dry soil (soil sampled in summer 2017). If all N in one spraying of a plot (69 µg N) entered the soil, the instantaneous mineral N availability in the soil would be increased by 0.1%. Similarly, the yearly maximum rate of N addition through sprayings is only 0.1% of the rate (1 g N m⁻² yr⁻¹) that was not found to affect B. glandulosa growth in a tundra experiment with comparable shrub vegetation³⁶.

Warming and herbivore treatments were also carried out through the 2019 growing season, which allowed us to provide supplementary leaf damage data for 2019 (Supplementary Tables 2 and 3).

Measuring and calculating CO₂ fluxes and VPD

Aluminum collars (outside dimensions 60 cm × 60 cm) were assembled in each plot to enable CO₂ exchange measurements using the closed chamber technique with transparent polycarbonate chambers (59 cm × 59 cm × 50 cm and 59 cm × 59 cm × 40 cm)⁵¹. As the collars were squares and the birch plantlets were planted 20 cm apart in three rows and four columns (following the shape of the plot), nine plantlets (3 × 3) were always included in the collars and three excluded (with the exception of one plot, where only six plantlets were included due to difficulties in the placement of the collar into the rocky soil surface). Each collar had grooves that were filled with moist quartz sand to provide an airtight seal between the collar and the chamber. The CO₂ concentration and air temperature inside the chamber were recorded continuously during measurements. A Vaisala CARBOCAP GMP343 (Vaisala Oyj, Finland) and a Picarro G2401 (Picarro Inc., CA, USA) online gas analyzer were used for taking CO₂ concentration measurements. The chamber closure time was 6 min for the GMP343 (used during 4 and 7 days in 2017 and 2018, respectively) and 2 min for the G2401 (1 day in 2017, 5 days in 2018). In all cases, the net CO₂ exchange was measured under the prevailing light conditions and with a dark hood. To determine the radiation response of photosynthesis, one or two additional shading levels were generated using meshes with 30 and 70% transparency when possible (160 of the 340 cases). The air inside the chamber was mixed with a battery-driven fan.

The solar radiation levels were measured with a photosynthetically active radiation (PAR) sensor (PQS1, Kipp & Zonen) on top of the chamber. In the GMP34₃-based system, there was a lid, made of the chamber material, above the PAR sensor to emulate the conditions inside the chamber, while the sensor was uncovered in the G2401-based system. The influence of the chamber wall on the measured PAR was tested afterwards with a similar transparent chamber, which indicated that the difference between the measurements taken inside and outside the chamber was <2%. In addition, it is important to note that the measurements used in the present study were each day conducted with the same system for all the chamber collars, so no systematic error was introduced into the comparison of treatments.

The CO₂ flux, i.e., the net ecosystem exchange (NEE), was calculated as

where p is atmospheric pressure, M is the molecular mass of CO₂ (44.01 g mol⁻¹), R is the universal gas constant (8.314 J mol⁻¹ K⁻¹), T is the mean air temperature during chamber closure, V is the chamber volume, A is the chamber base area, and (frac{{{mathrm{d}{c}}}}{{{mathrm{d}{t}}}}) is the mean CO₂ mixing ratio change in time calculated with linear regression (see Supplementary Figs. 1 and 2). To allow for the stabilization of the CO₂ flux after the chamber closure, the last 5 and 1.5 min of the data series recorded with the GMP343 and G2401, respectively, were included. A micrometeorological sign convention was used: a positive flux indicates a flux from the ecosystem to the atmosphere (emission), and a negative flux indicates a flux from the atmosphere into the ecosystem (uptake).

The measured NEE was partitioned into gross primary productivity (GPP) and ecosystem respiration (R_e):

R_e was obtained from the measurements with a darkened chamber. GPP was modeled by a rectangular hyperbola⁵²:

where a is the initial slope between GPP and PAR, and GP_max is the theoretical maximum gross photosynthetic rate. The parameters α and GP_max were first estimated for those light response measurements that had at least three light levels available and at least one data point with PAR > 800 μmol m⁻² s⁻¹ (n = 68 with 20, 15, 16, and 17 measurements from control, warming, herbivory reduction, and warming × herbivory reduction treatment plots, respectively). Using these parameter values, GP_max/α ratios were calculated and the median value of 203 μmol m⁻² s⁻¹ was chosen to be used as a common GP_max/α ratio. This fixed value was employed in all GPP calculations, while α (or, equally, GP_max) was available as a free parameter. The reason for introducing a common GP_max/α was to obtain consistent data also for those days when irradiance was limited and the saturation level of photosynthesis was not attained. Furthermore, there were no statistically significant differences in GP_max/α among the treatment combinations (P = 0.094 for the Mood’s median test). Using the common GP_max/α, the α parameter was fitted separately for each light response measurement (i.e., for each chamber collar in each day, n = 340), and NEE₈₀₀ and GPP₈₀₀, i.e., NEE and GPP at PAR = 800 μmol m⁻² s⁻¹, were calculated. NEE₈₀₀ and GPP₈₀₀ represent the CO₂ exchange in conditions typical of the daily maximum PAR and serve as standardized metrics that allow comparison of the ecosystem carbon sink potential among different days and treatments. The different steps of the fitting procedure described above are validated in Supplementary Figs. 3–6. We also standardized NEE at a lower (300 μmol m⁻² s⁻¹) and higher (1200 μmol m⁻² s⁻¹) PAR level, but there were no marked differences in the outcome of the statistical analysis (Supplementary Table 8).

The water vapor pressure deficit, VPD (kPa), was calculated as

where T is air temperature in °C and RH is relative humidity in %.

Soil N and microbial biomass C

Plant availability of mineral N (NH₄⁺ and NO₃⁻) was estimated in the soil organic layer (at 3 cm depth) using ion-exchange resin bags (UNIBEST Ag Manager™). Three resin bags were placed in each plot in July 2016 and in the following autumns one (2017) or two (2018) bags were transferred from each plot to the laboratory for N extraction in 50 ml of 2 M KCl. The KCl solution was filtered through a glass microfiber filter (Whatman, Germany) and NH₄⁺ and NO₃⁻ concentrations were analyzed using a Lachat QuikChem 8000 analyser (Zellweger Analytics, Lachat Instruments Division, USA).

Microbial biomass samples were collected with an auger (diameter 2.5 cm) from three (or more if the soil was shallow) spots within each plot in August 2018. Three layers (organic layer, 0–5 cm layer of mineral soil, >5 cm layer of mineral soil) were separated from each soil core. Within each plot, the samples collected from the same layer were pooled and sieved through a 2-mm (mineral soil) or 6-mm mesh (organic soil) before analyzing microbial biomass carbon (MBC) using the chloroform-fumigation extraction method⁵³. In 2018, we further tested whether the insecticide sprayings could have direct effects on soil microbes. Soil from 10 random spots outside the experimental area was collected in early spring, homogenized (largest roots and rocks removed), and placed into ten 1.5-L pots. The pots were buried in the ground and covered with a 1–2 cm layer of dead, oven-dried (70 °C for 20 h, 100 °C for 4 h) Sphagnum moss to mimick the insulating moss layer in the experimental plots (where the mean moss depth was 1.2 cm). Along with the sprayings in the experimental plots, five of the pots were then sprayed with water and five with the insecticide. Soil was collected from the pots along with the microbial biomass C sampling and analyzed accordingly. No difference (t₈ = 0.389 and P = 0.708) between the mean MBC in soils sprayed with water (546.3 ± 85.0 mg MBC kg⁻¹ soil; mean ± SE) and insecticide (494.3 ± 103.4) was found (source data are provided as a Source Data file).

Phenology and growth of Betula plantlets

For estimating treatment effects on plant spring phenology, all birch plantlets were surveyed daily in spring 2017 and 2018. The date of the first bud of a plantlet opening was considered as the start of the growing season for that plantlet. A bud was considered open once the protective bud scales were completely separated and the emerging leaf was visible⁵⁴. During the peak growing season, plant performance was estimated two (2017) to three (2018) times by measuring the chlorophyll content of five full-grown topmost leaves in several branches in each plantlet, using the CCM300 non-destructive optical chlorophyll content meter (Opti-Sciences, USA). The autumn phenology was in turn estimated through chlorophyll breakdown: the chlorophyll content of the five topmost leaves of each plantlet was measured approximately every 3 days using the CCM300 until most of the plantlets had shed their leaves.

To assess plantlet shoot growth, their height was measured at planting and thereafter each autumn following growth cessation. The relative growth was then calculated as a height increment relative to the initial height.

Leaf herbivore damage of Betula plantlets

To estimate treatment effects on herbivore load, all leaves in each birch plantlet were surveyed for herbivore damage in the middle of August 2017. Leaf damage was illustrated using a Schreiner-type method^20,21,55 in which a damage index, ranging from 0 to 100, is produced for each plant individual by multiplying two scores, A and B. The score A is the plantlet mean of values (ranging from 0 to 25) that illustrate the size of the damaged area in a single leave: 0 = no damage, 1 = small damaged area (1–4% of leaf area damaged), 5 = medium damaged area (5–20%) and 25 = large damaged area (>20%). The score B (ranging from 0 to 4) tells the percentage of leaves of a plantlet that are damaged: 0 = 0%, 1 = 1–25%, 2 = 26–50%, 3 = 51–75%, and 4 = 76–100%. As no data was obtained for 2018, we repeated the leaf damage survey in August 2019 to ensure that results from 2017 can be generalized across years.

Statistical analysis

Treatment effects on response variables were tested using mixed models and Type I ANOVA (IBM SPSS Statistics 21), where the variance is allocated to explanatory variables in order of their appearance (in tables of statistics, order of explanatory variables follows their order in respective models). Soil OM content and the areal cover of vascular plants, lichens, and mosses were treated as covariates and included in the models first to remove plot-to-plot variation that might otherwise contribute to treatment effects. Warming and herbivory treatments were treated as fixed effects and field block as a random effect. For birch variables, species was included as a fixed effect and genotype, nested within species, as a random effect. Date, year, and soil layer (MBC was measured for three adjacent soil layers) were treated as repeated measurements and compound symmetry was used as a repeated covariance structure. Those response variables that were measured multiple times during a growing season were tested separately for 2 years. ANOVA assumptions were checked from the residuals and transformations were used where necessary (reported in table and figure legends).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com