BVOC Emissions From a Subarctic Ecosystem, as Controlled by Insect Herbivore Pressure and Temperature

The biogenic volatile organic compounds, BVOCs have a central role in ecosystem–atmosphere interactions. High-latitude ecosystems are facing increasing temperatures and insect herbivore pressure, which may affect their BVOC emission rates, but evidence and predictions of changes remain scattered. We studied the long-term effects of + 3 °C warming and reduced insect herbivory (achieved through insecticide sprayings) on mid- and late summer BVOC emissions from field layer vegetation, supplemented with birch saplings, and the underlying soil in Subarctic mountain birch forest in Finland in 2017–2018. Reduced insect herbivory decreased leaf damage by 58–67% and total ecosystem BVOC emissions by 44–72%. Of the BVOC groups, total sesquiterpenes had 70–80% lower emissions with reduced herbivory, and in 2017 the decrease was greater in warmed plots (89% decrease) than in ambient plots (34% decrease). While non-standardized total BVOC, monoterpene, sesquiterpene and GLV emissions showed instant positive responses to increasing chamber air temperature in midsummer samplings, the long-term warming treatment effects on standardized emissions mainly appeared as changes in the compound structure of BVOC blends and varied with compounds and sampling times. Our results suggest that the effects of climate warming on the total quantity of BVOC emissions will in Subarctic ecosystems be, over and above the instant temperature effects, mediated through changes in insect herbivore pressure rather than plant growth. If insect herbivore numbers will increase as predicted under climate warming, our results forecast herbivory-induced increases in the quantity of Subarctic BVOC emissions.

(achieved through insecticide sprayings) on midand late summer BVOC emissions from field layer vegetation, supplemented with birch saplings, and the underlying soil in Subarctic mountain birch forest in Finland in 2017-2018. Reduced insect herbivory decreased leaf damage by 58-67% and total ecosystem BVOC emissions by 44-72%. Of the BVOC groups, total sesquiterpenes had 70-80% lower emissions with reduced herbivory, and in 2017 the decrease was greater in warmed plots (89% decrease) than in ambient plots (34% decrease). While non-standardized total BVOC, monoterpene, sesquiterpene and GLV emissions showed instant positive responses to increasing chamber air temperature in midsummer samplings, the long-term warming treatment effects on standardized emissions mainly appeared as changes in the compound structure of BVOC blends and varied with compounds and sampling times. Our results suggest that the effects of climate warming on the total quantity of BVOC emissions will in Subarctic ecosystems be, over and above the instant temperature effects, mediated through changes in insect herbivore pressure rather than plant growth. If insect herbivore numbers will increase as predicted under climate warming, our results forecast herbivory-induced increases in the quantity of Subarctic BVOC emissions.

INTRODUCTION
Ecosystems in Subarctic and Arctic regions, while known to be sensitive to warming (Walsh 2014), are currently experiencing temperature increments twice the rate of the global average (Pachauri and others 2014). Even if the global temperature increase could be limited to the 2°C target by the end of the century (IPCC 2018), annual mean temperatures in the Arctic would likely increase by 4°C by the mid-century and by over 5°C by the year 2100 (Overland and others 2019). Such warming will have consequences on vegetation distribution, plant species composition and plant growth in Arctic areas (Box and others 2019;Forkel and others 2016;Keenan and Riley 2018). For instance, already now, in the Subarctic Fennoscandia warming has increased plant growth and led to a northward expansion of shrub cover and the mountain birch (Betula pubescens ssp. czerepanovii) tree line (Forkel and others 2016;Myers-Smith and others 2011;Truong and others 2007).
Besides affecting vegetation, climate warming also plays a significant role in the species composition and population dynamics of insect herbivores. Warming can lead to species range expansions towards higher latitudes (Jepsen and others 2011), cause periodic increases and enhance the outbreak potential in insect populations (Dale and others 2001;Jepsen and others 2011). Warming will also likely increase the non-outbreak background densities of insects. At high latitudes, background insect herbivory is typically low (Barrio and others 2017; Kozlov 2008; Li and others 2019), but observations of downy birch (Betula pubescens) leaf damage in Fennoscandia show that the leaf area eaten increases from 1-2 to 5-7% along a latitudinal gradient from 70°N to 60°N (Kozlov 2008). As this gradient climatically roughly represents the predicted warming in the next 100 years in these areas (Ruosteenoja and others 2016), insect herbivore damages on plant leaves in Herbivory, Temperature and Subarctic BVOC Emissions the Arctic can be assumed to more than double in the near future. These predictions are supported by recent findings from the tundra, where dwarf birch (Betula nana) growing in experimentally warmed plots had fourfold leaf damages in comparison with plants growing in control plots (Li and others 2019).
Biogenic volatile organic compounds (BVOCs) are a major group of secondary metabolites released by all organisms, but especially by the terrestrial plants through their aboveground tissues (Laothawornkitkul and others 2009). Soils can be both sources and sinks of biogenic VOCs with varying levels of contribution to the ecosystem-scale emissions (Tang and others 2019). Of the aboveground emissions, isoprene (C 5 H 8 ) and monoterpenes (C 10 H 16 ) are the dominant and most widely studied BVOC groups, which together with the sesquiterpenes (C 15 H 24 ) constitute terpenoids. BVOCs are reactive gases, which in the atmosphere undergo chemical reactions with O 3 and hydroxyl (OH) and nitrate (NO 3 ) radicals (Atkinson and Arey 2003), and therefore, have a lifetime from minutes to days only (Kesselmeier and Staudt 1999). The low-volatile oxidation products of BVOCs enhance the formation and growth of secondary organic aerosols (SOAs) (Hao and others 2009;Paasonen and others 2013), which in the atmosphere scatter and absorb solar radiation and act as cloud condensation nuclei. As this in turn controls the radiative balance and climate on Earth (Carslaw and others 2009), BVOCs have a significant role in the ecosystem-atmosphere interactions and feedbacks (Kulmala and others 2013; Š impraga and others 2019).
BVOC emission rates at both tree canopy and ecosystem scales, including those from the soils, depend on temperature (Aaltonen and others 2011;Kramshøj and others 2016;Laothawornkitkul and others 2009). For instance, climate warming is known to increase BVOC emissions from Northern Hemisphere trees (Holopainen and others 2018), and the effects are much greater at the temperature-limited northern ecosystems compared to other ecosystems (Kramshøj and others 2016). Similarly, in a Subarctic tundra heath, simulated 2°C warming led to 1.5-1.8-fold ecosystem-scale isoprene emissions (Tiiva and others 2008) and doubled the emission rates of mono-and sesquiterpenes (Faubert and others 2010). In another experiment, exposing Arctic tundra to six years of experimental warming increased the total ecosystem emission rate by 260% (Kramshøj and others 2016). Soil BVOC emissions mainly originating from leaf litter, soil organic carbon and plant roots (Svendsen and others 2018; Tang and others 2019), were found to be less influenced by warming in these high-latitude studies (Faubert and others 2012;Kramshøj and others 2016).
In addition to being affected by temperature, plant BVOC emissions are known to be influenced by insect herbivory, which typically triggers higher emission rates with unique BVOC blends (Trowbridge and Stoy 2013). For instance, even a low herbivore pressure by geometrid moths can increase BVOC emissions from Subarctic mountain birch forests (Rieksta and others 2020). Herbivorytriggered BVOC blends are dominated by terpenoids (mono-and sesquiterpenes and hemiterpenes) (Irmisch and others 2014;McCormick and others 2014;Rowen and Kaplan 2016) and green leaf volatiles, GLVs (Ameye and others 2018) and can offer information on host plant location and the herbivore species, hence providing cues for herbivore enemies (Bruce and Pickett 2011;McCormick and others 2014;Ode 2013). BVOCs thus play a significant ecological role in terrestrial communities in mediating within-and betweenplant communication and plant-insect interactions (Holopainen and others 2013;Trowbridge and Stoy 2013).
In a recent study in the Arctic tundra, Li and others (2019) found that long-term warming and herbivory-mimicking methyl jasmonate application increased dwarf birch terpenoid emissions and that the increase was much higher when the treatments were combined. Although warming is likely to increase insect pressure in high latitudes (Kozlov 2008), evidence of responses of BVOC emissions to warming and insect herbivory remains scarce for these regions. Here, we present results from an open-air warming experiment established at a Subarctic mountain birch forest in northern Finland, where we measured ecosystem-scale BVOC emissions (thus including emissions from both the ground layer vegetation and the underlying soil) in response to long-term experimental 3°C warming and insect herbivory reduction. Using the same experimental set-up, Silfver and others (2020) have recently showed that plant growth and the ecosystem CO 2 uptake potential were significantly, and almost equally, increased by both warming and reduced leaf herbivory. In our study, we predicted that total BVOC emission rates (1) will show instant responses to air temperature variation, but (2) will also follow the increased plant growth and CO 2 uptake and show long-term positive responses to the warming treatment and (3) will decrease with reduced insect herbivory due to decreasing amounts of herbivory-induced BVOC emissions.

STUDY SITE AND TREATMENT PLOTS
The experiment was established in mountain birch forest at the Turku University Kevo Subarctic Research Station in Utsjoki, Finland (69°45¢ N, 27°01¢ E, 104 m a.s.l.) in spring 2016. In the site, the mean annual temperature is -1.3°C and the mean annual precipitation 354 mm (reference years 1981-2010; data from Finnish Meteorological Institute). When the site was established, adult mountain birches were trimmed, but the field layer vegetation (see below) was left intact. The experimental area of 220 m 2 was divided into five replicate blocks. Four treatment plots (each 100 cm 9 75 cm), with no mountain birch stumps included and always > 1.2 m apart, were established in each block.
Subarctic woody species were included in the experiment as plantlets of four birch species-dwarf birch (B. nana), mountain birch (B. pubescens ssp. czerepanovii), downy birch (B. pubescens) and silver birch (B. pendula)-that were planted 20 cm apart in each treatment plot amid the ground layer vegetation in mid-July 2016 (twelve plantlets in each plot, three genotypes for each species). The plantlets were cloned at the Haapastensyrjä Unit of the Natural Resources Institute Finland (Luke) at the end of 2015 using standard micropropagation techniques (Ryynä nen 1996) applied to twig material collected from natural Subarctic populations growing between 67°43¢ N and 69°01¢ N. After in vitro replication the plants were first grown on vermiculite in a greenhouse for several weeks and then transferred into plastic trays filled with fertilized peat. Plantlets in plastic trays were acclimatized to outdoor conditions before being transferred to the experimental site. Besides representing woody species, the plantlets provided well replicated and precise responses of plant growth (Silfver and others 2020) and leaf damages to the treatments applied in the experiment. Such data could not have been obtained using natural ground layer vegetation that varied in species composition and abundance from plot to plot. At the time of planting, the mean height of birch plantlets was 6.2 cm (range 1-29 cm) and altogether the twelve plantlets covered on average 3.8% of treatment plot area (Table S1).
Aluminium collars (outside dimensions 60 9 60 cm) were assembled into each treatment plot to enable CO 2 (Silfver and others 2020) and VOC measurements using the closed chamber technique. Nine of the twelve birch plantlets were included in the collar, except in one plot, where only six plantlets were included due to difficulties in assembling the collar into the rocky ground.

Plot Vegetation and Soil Properties
Baseline vegetation and soil properties were examined in treatment plots in the beginning of growing season 2017. On the site, the field layer vegetation is a spatially heterogeneous mixture of evergreen and deciduous shrubs, mosses and lichens, with few graminoids and herbs (Table S1). In each plot, plants and lichens were identified to species or genus level and their aerial cover was estimated for the collar area (inside dimensions 55 9 55 cm). Vascular plants (including the Betula plantlets) covered 29-83% of the area (with Empetrum nigrum, Vaccinium uliginosum, V. vitis-idaea, V. myrtillus and Rhododendron tomentosum as the most abundant species), mosses 3-85% (Pleurozium schreberi, Hylocomium splendens and Barbilophozia lycopodioides) and lichens 0-23% (Nephroma arcticum and Cladonia arbuscula) (Table S1). These values were later used as covariates in the statistical models to control the effects of vegetation plot-toplot variation on BVOC emission rates. At plot level, vascular plant cover correlated negatively with lichen cover (Spearman's q = -0.53, P = 0.016, n = 20), while moss cover did not clearly correlate with either vascular plant cover (q = 0.20, P = 0.405) or lichen cover (q = -0.25, P = 0.281).
During the experiment, the trimmed mountain birches sprouted and formed bushes in the experimental area. The foliage volume of each bush that located closer than 2 m from the centre of any plot was estimated by measuring two crossed radii and the height of the bush. The obtained foliage volume was then divided by the distance of the bush from the centre of a plot to form an index that described the strength of the bush effect on the plot. If several bushes grew in the vicinity of a plot, their indices were summed together. The bush effect index values calculated for each plot were used as a covariate in the statistical models to remove the potential effect of varying number, size and proximity of sprouted mountain birch bushes and their roots on the BVOC emission rates quantified on treatment plots.
The soil in the experimental area is podzol with a clear organic layer, sandy mineral soil and randomly placed rocks. Soil C concentration is on average 41 and 1.7% of dry mass, N concentration 1.6 and 0.14% of dry mass and the C/N ratio 26 and 12 in the organic layer and in the top 0-5 cm mineral layer, respectively (Nele Meyer, pers. comm.). Among the treatment plots, the depth of the organic layer varies from 2 to 6.5 cm (Nele Meyer, pers. comm.) and the soil organic matter (OM) content, in the top 10 cm of soil, from 9 to 38% of soil dry mass (Silfver and others 2020). Soil OM content was used as a covariate in the statistical models to control the effects of soil plot-to-plot variation on BVOC emissions.

Experimental Design and Treatments
Treatments included two levels of temperature (ambient, + 3°C) and two levels of insect herbivory (natural, reduced) in a fully factorial 2 9 2 design, leading to four treatment combinations: C = control (ambient temperature, natural herbivory), R = reduced insect herbivory (ambient temperature, reduced herbivory), W = warming (+ 3°C temperature, natural herbivory), W + R = warming combined with reduced insect herbivory (+ 3°C temperature, reduced herbivory). In each replicate block, treatment combinations were randomly allocated to the four plots.
Warming was realized in plots using two 240 mm 9 60 mm ceramic heaters (Elstein-Werk M. Steinmetz GmbH & Co, Germany; dummy heaters in ambient plots), installed 80 cm above the ground. The heaters were controlled by a microprocessor-based feedback system, which maintained a 3°C difference between the ambient and heated plots using real-time temperature data. Data were collected by sensors installed 30 cm above the ground (that is, on the top of the vegetation layer) within green metal plates, mimicking plant leaf surfaces, in five ambient and five heated plots (see more details in Silfver and others 2020). In 2016, when the experiment was established, warming was turned on in early July, but in later years, always in spring when the snow-cover melted below 20 cm. Warming was turned off in late autumn when mean daily air temperatures in ambient plots remained permanently below -3°C. Air temperatures ( Figure 1) were recorded using sensors installed under white plastic plates 30 cm above the ground in three ambient and three warmed plots. On average (± SE, n = 3), the plates mimicking plant leaves were 3.3 ± 0.37°C, the air 2.3 ± 0.35°C and the soil (at the depth of 5 cm) 1.2 ± 0.15°C warmer in heated than ambient plots across the warming periods in 2016-2018 (Silfver and others 2020).
The herbivore treatment was started in the beginning of the 2017 growing season and continued through 2017 and subsequent growing seasons. The herbivore reduction plots were sprayed with 0.2 L of 0.1% synthetic pyrethrin solution (Decis Mega EW50, Bayer Crop-Science, Germany) and the control plots with 0.2 L tap water weekly, using two portable garden sprayers and a protective tent to eliminate wind drift. Synthetic pyrethrin has not been found to have sideeffects on birch growth or leaf chemistry (Silfver and others 2013), neither on soil microbial growth (Silfver and others 2020) nor birch leaf litter decomposition rate (Silfver and others 2013). The quantity of nitrogen that could enter the soil while deltamethrin, the active ingredient, degrades after spraying is also negligible (Silfver and others 2020). The main feeding guilds of insects were leaf feeders, miners and rollers. To quantify the leaf damage in plots with natural herbivory and the effects of the insecticide in the plots with reduced herbivory, all leaves in each Betula plantlet were surveyed in each plot in the middle of August in 2017 and 2019 (results for summer 2018 are not available). Leaf damage was illustrated using a modified Schreiner-type method (Fritz and others 1998;Prittinen and others 2003), 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%. The low number of sufficiently clear classes of damaged area makes the damage index feasible and repeatable in field conditions, where a high number of leaves has to be examined and no samples can be collected for laboratory inspection. Variation in leaf damage was assessed using the birch plantlets as they included controlled genetic material of several species and were evenly distributed among and within treatment plots.

BVOC Sampling and Analysis
BVOCs were sampled from the vegetation headspace (thus including all aboveground and belowground parts of the entire ecosystem) using a transparent Teflon chamber (609 60 9 50 cm, wall thickness 1.5 mm) in each plot in late June and mid-August in both 2017 and 2018. BVOC samples were collected block by block, randomly from different treatment combinations and over 3-4 d in each sampling campaign. In the treatment plot, the chamber was placed on aluminium collar grooves filled with moist quartz sand to provide an airtight seal. At the time of the campaign days (Figure 1), the means of daily mean and daily max air temperatures, respectively, were 18 and 21°C in June 2017, 11 and 13°C in August 2017, 14 and 17°C in June 2018 and 10 and 16°C in August 2018.
Using a battery-operated sampling pump (Rietschle Thomas 5002 12 V DC, Puchheim, Germany), clean charcoal-filtered and MnO 2 scrubbed air was pumped into the chamber through Teflon tubing and a vent made at the upper part of the chamber. In the sampling pump, the flow rate of incoming air was set to 0.6 L min -1 using a mini-Buck calibrator (Model M-5, A.P. Buck, Inc., Orlando, FL, U.S.A). After flushing air into the chamber for 30 min, the headspace sample was pulled at a rate of 0.2 L min -1 via Teflon tubing through another vent made at the other side of the chamber into a stainless steel tube (ATD sample tubes; PerkinElmer, Boston, MA, USA), filled with Tenax TA and Carbograph 5TD adsorbents (Tenax TA, 100 mg and Carbopack B, 100 mg, mesh 60/ 80; Supelco, Bellefonte, PA, USA). The sampling period lasted for 45 min, and the total volume of the collected air sample was 9 L.
Air was circulated inside the chamber using a small battery-operated fan fixed at the lower part of the chamber. The temperature and air humidity inside the chamber were monitored using wireless data loggers (HygrochronDS1923-f5iButton, Maxim Integrated Products, San Jose, CA, USA) and the photosynthetically active radiation (PAR) using a quantum sensor (Li-Cor, Lincoln, NE, USA). Air temperatures within the chamber varied from 12 to 44°C across the four campaigns ( Figure 2). This variation was due to varying environmental conditions and was used for testing the instant responses of non-standardized BVOC emission rates to air temperature. Over all samplings, the chamber air temperature was on average higher in warmer (25.4 ± 1.5°C, mean ± SE, n = 40) than in ambient (24.3 ± 1.2°C) plots, but the difference was not statistically significant in either 2017 (F 1,14 = 1.67, P = 0.217) or 2018 (F 1,14 = 0.41, P = 0.532). Mean PAR levels inside the chamber were 880 (June 2017), 560 (August 2017), 550 (June 2018) and 393 lmol m -2 s -1 (August 2018). Once the sampling was over, the BVOC-sampled tube was immediately sealed with a Teflon-coated brass cap and stored in a refrigerator. A blank air sample was collected from an empty closed Teflon chamber-the same chamber as used for BVOC samplings from treatment plots-to examine the purity of the background air and to detect compounds possibly originating from the used materials or analysis system. BVOC samples were analysed within two weeks using a gas chromatography-mass spectrometry (GC -MS) (Hewlett Packard GC type 6890, Waldbronn, Germany; MSD 5973, Beaconsfield, UK) with thermal desorption as described by Ghimire and others (2016). BVOCs were identified by using pure standards and by comparing the compounds' mass spectra with the spectra of standards and the Wiley library (ver. WILEY275.L). The compounds that had no standards were quantified using a-pinene (used for non-oxygenated monoterpenes), 1,8-cineole (oxygenated monoterpenes) and longifolene (sesquiterpenes) as reference compounds.
The detection limits set in the chromatograms were 1 ng for mono-and sesquiterpenes and 3 ng for green leaf volatiles. To remove the effect of temperature and radiation variation among samplings on BVOC emission rates, isoprene emissions were standardized to a temperature of 30°C and a PAR level of 1000 lmol m -2 s -1 using the model by Guenther (1997). Emission rates of all other compounds were standardized to a temperature of 30°C using the algorithm developed by Guenther and others (1993) with a temperature coefficient (b) of 0.09 for monoterpenes (Guenther and others 1993), 0.16 for sesquiterpenes (Helmig and others 2006) and 0.1 for other VOCs (Guenther and others 2012). Standardized emission rates were used for testing the long-term responses of BVOC emission rates to the warming treatment (in contrast to instant responses to air temperature variation). Additionally, to check that these responses were not suppressed by standardization, the responses were also tested using non-standardized emission rates.
The concentrations of compounds detected in blank samples were subtracted from concentrations found in samples collected from the treatment plots. BVOC emission rates were then calculated in ng m -2 h -1 using the equation: where E = VOC emission rate (ng m -2 h -1 ), F = flow rate of incoming air (L h -1 ), C2 = concentra- Figure 2. Relationship between chamber air temperature and non-standardized total BVOC, monoterpene, sesquiterpene and GLV emissions during BVOC sampling campaigns in June and August 2017 and 2018. P-values are for temperature effects from linear mixed models (with baseline covariates and herbivory treatment also included), and regression lines are shown for temperature effects with P £ 0.1. tion of compound in the outgoing air (ng L -1 ), C1 = concentration of compound in the incoming air (considered as 0 ng L -1 as the incoming air was filtered) and A = chamber surface area (m 2 ) (Ortega and Helmig 2008). Standardized mean emission rates for the four treatment combinations and for the four samplings are given in Tables 1 and S2, respectively, and the non-standardized rates for the four treatment combinations are shown in Table S3 and used in Figure 2 to illustrate the effect of chamber air temperature on BVOC emission rates.

Statistics
The main and interaction effects of warming and reduced herbivory on standardized and non-standardized total BVOC, monoterpene, sesquiterpene and GLV emission rates, proportions of group emissions of all emissions, and birch leaf damage were tested using linear mixed model and Type I ANOVA (IBM SPSS Statistics 21), where the variance is allocated to explanatory variables in order of their appearance in the model. The baseline aerial cover of vascular plants, lichens and mosses, soil OM content and bush effect index were treated as covariates and included in the model before treatments to remove plot-to-plot variation that might otherwise hide or contribute to treatment effects. Those covariates that had P > 0.1 in the full model were removed from the final model unless this led to less clear treatment effects. Warming and herbivory treatments were treated as fixed effects, field block and birch genotype as random effects and the samplings within years as repeated measures. The effect of chamber air temperature variation on non-standardized total BVOC and group emission rates was tested for each sampling individually using a model, where baseline covariates were chosen as explained above, chamber temperature was included as another covariate and herbivory treatment as a fixed factor. Finally, the effect of warming treatment on chamber air temperature was tested for each year separately using a model with only the warming treatment, month and their interaction included. ANOVA assumption of normal distribution was checked using histograms and Shapiro-Wilk tests of model residuals. If residuals were not normally distributed, the data were square-root-or log 10transformed, the choice between these depending on the skewness of the distribution. For mixed models, the significant effects are presented in Figures. 3, 4 and 5 and the full statistical models with F-and P-statistics in Tables S4-S7. Leaf damage results were already presented by Silfver and others (2020), but since leaf damages are essential for interpreting BVOC emissions, they are also presented in this study.
The main and interaction effects of warming and reduced herbivory on monoterpene, sesquiterpene and GLV blends (combining information of compound presence/absence and emission rates) were tested using PERMANOVA. Baseline covariates for these analyses were chosen following the principles described above. For PERMANOVA results, the significant treatment effects are presented in text and the full statistical models with F-and P-statistics in Table S8. Those compounds that were found to significantly contribute (> 10% contribution) to the observed treatment effects in PERMANOVA are shown in Figure 6 as examples of single compound responses (with no statistical inference of treatment effects on single compounds, however).

PLANT LEAF DAMAGE
Insect damage indices calculated for birch leaves were on average 67% and 58% lower in plots with reduced insect herbivory than in plots with normal herbivory in 2017 and 2019, respectively (Figure 3). Moreover, marginally significant herbivory 9 warming interaction effects appeared in both years (Figure 3). These interactions were because reduced herbivory decreased leaf damage more under ambient temperature (on average by 78%) than in heated plots (46%) (Figure 3). Warming did not increase leaf damage, except for B. pubescens under reduced herbivory in 2017 (significant warming 9 herbivore reduction 9 Betula species interaction effect in Table S4).

Identified BVOCs and Emission Rates Across the Samplings
Altogether 18 monoterpenes, 20 sesquiterpenes, five GLVs, isoprene (a hemiterpene) and methyl salicylate were identified in the emissions during the four samplings (Table 1, S2). Monoterpenes and sesquiterpenes dominated the BVOCs and together comprised about 95% of all emissions (Table 1), except for June 2017 when the proportion of mono-and sesquiterpenes together was 81% and the proportion of GLVs was 18% (Table S2). Isoprene covered 0.9-5.6% and methyl salicylate 0.1-0.4% of emissions across the treatments and samplings ( and GLVs increased with increasing chamber air temperature in the June, but not in the August samplings ( Figure 2).

Treatment Effects on Total BVOC and Compound Group Emission Rates
When standardized emission rates were analysed using univariate linear mixed models (see complete statistical results in Table S5), total BVOC emission rates were on average 44 and 72% lower in plots with reduced herbivory than in plots with natural herbivory in 2017 and 2018, respectively, whereas the warming treatment had no statistically significant effect (Figure 4). Responses of sesquiterpene emission rates to treatments followed these trends with on average 70 and 80% lower emissions with reduced herbivory than with natural herbivory in 2017 and 2018 (Figure 4). A significant herbivory 9 warming interaction effect on sesquiterpene emission rates also appeared in 2017. This was because herbivory reduction led to much greater decrease in emissions in warmed plots (89% decrease) than in plots with ambient temperature (34% decrease) (Figure 4). The treatment response pattern in monoterpene emission rates followed the pattern in total BVOC and sesquiterpene emission rates, but the herbivory reduction effect was not statistically significant ( Figure 4, Table S5). Neither was there a warming treatment effect on monoterpene emissions (Figure 4, Table S5). Also, no effects on total GLV  (Table S4) fitted to log 10 -transformed data (n = 57-60 per treatment in 2017, n = 44-55 per treatment in 2019). These results were originally presented in Figure 2d and Suppl. Table 1 in Silfver and others (2020). Quantified using reference compounds. Mean (SE; n = 5) emission rates (ng m -2 h -1 ) of total BVOCs, monoterpenes (MTs), sesquiterpenes (SQTs) and green leaf volatiles (GLVs), with percentages of group emissions of all emissions, and individual compounds for each treatment combination (C = control, R = reduced herbivory, W = warming, W + R = warming with reduced herbivory) across the four samplings. All emission rates were standardized to a temperature of 30°C and the isoprene emission rate further to a PAR level of 1000 lmol m -2 s -1 .
emission rates were found in 2018, whereas a significant warming effect and warming 9 month interaction effect appeared in 2017 (Figure 4). The interaction effect was because warming increased emission rates tenfold in August, but not in June (Figure 4). The marginally significant warming 9 herbivory interaction effect in 2017 was because herbivory reduction increased emission rates in warmed plots, but not in plots with ambient temperature (Figure 4). The proportion of sesquiterpenes of total BVOCs was lower in plots with reduced herbivory than in plots with natural herbivory in both ambient and heated plots in 2018, whereas in 2017 this effect was found in heated plots only as suggested by the significant warming 9 herbivory interaction effect ( Figure 5; see complete statistical results in Table S6). As a mirror image, the proportion of monoterpenes of total BVOC emissions increased with reduced herbivory in heated plots as suggested by the significant warming 9 herbivory interaction effect in both years ( Figure 5). The proportion of GLVs of all BVOCs was affected by a significant warming 9 month interaction effect in 2017 and a warming 9 herbivory interaction effect in 2018 ( Figure 5). In 2017, the interaction is explained by warming increasing GLV proportion in August, but not in June ( Figure 5). In 2018, the interaction can be interpreted either as warming increasing GLV proportion in plots with natural herbivory, or as reduced herbivory lowering GLV proportion in heated plots ( Figure 5).
Using non-standardized emission rates (Table S3) instead of standardized rates (Table 1) did not change the results in terms of warming effects on BVOC emissions (Table S7). No warming effects emerged for total BVOCs, monoterpenes or sesquiterpenes, and for GLVs, the warming main effect (P = 0.013) and the warming 9 month (P = 0.075) and herbivory 9 warming (P = 0.021) interaction effects appeared in the same way as with the standardized data in 2017. The only new warming effect was the marginally significant (P = 0.081) effect on GLVs in 2018, when warming increased non-standardized GLV emissions on average by 25% (Table S7).
In both years, the compounds that had > 10% contribution to the observed warming effect were a-copaene (18 and 14% contribution in 2017 and 2018, respectively), c-gurjunene (12 and 11%) and palustrol (11 and 22%). The warming treatment decreased a-copaene emissions on average by 34 and 15%, but increased c-gurjunene emissions 2.1fold and 4.0-fold in 2017 and 2018, respectively ( Figure 6). Palustrol was emitted exclusively by those plots that had R. tomentosum with one exception: one plot with no R. tomentosum also emitted palustrol in August 2018. However, while all 14 plots with R. tomentosum emitted palustrol in 2018, only 6 emitted palustrol in 2017. When palustrol emission rates were standardized to mean (8%) R. tomentosum plot cover, emission rates were 5.8-fold in warmed plots in comparison with ambient plots in 2017, whereas in 2018 emissions were reduced by warming on average by 20% ( Figure 6).

Green Leaf Volatiles
The three most abundant GLVs (all with mean emission rates > 190 ng m -2 h -1 ) were nonanal (39-50% of GLV emissions in the four samplings), Figure 4. Total BVOC and compound group emission rates (mean + SE; n = 5) from treatment plots (C = control, R = reduced herbivory, W = warming, W + R = warming with reduced herbivory) in June and August 2017 and 2018. Pvalues are shown for treatment effects when P £ 0.1 (M = month effect). All emission rates were standardized to a temperature of 30°C and the isoprene emission rate further to a PAR level of 1000 lmol m -2 s -1 . Means and errors are estimated marginal means and standard errors produced by statistical models (Table S5) fitted to square-root-transformed (total BVOCs, SQTs and MTs) or log 10 -transformed (total GLVs) data.
cis-3-hexenyl-acetate (1-40%) and cis-3-hexen-1ol (7-21%) ( Table S2). The multivariate test found a warming 9 month interaction effect on the GLV blend in 2017 (P = 0.017), but no main or interaction effects of warming in 2018 (Table S8). The 2017 interaction was because the effect of warming on the GLV blend was much stronger in August than June. In August 2017, nonanal was the only GLV that was found in the majority of field plots (others were found in 0-4 plots) and contributed 76% to the observed warming effect by having 2.6fold increase in emission rates in warmed plots in comparison with ambient plots (Figure 6). No effects of herbivory reduction on the GLV blend were found in either 2017 (P = 0.600) or 2018 (P = 0.403).  (Table S6) fitted to non-transformed data. Figure 6. Emission rates (mean + SE; n = 5) from treatment plots (C = control, R = reduced herbivory, W = warming, W + R = warming with reduced herbivory) in June and August 2017 and 2018 for those SQTs (a-copaene, c-gurjunene, and palustrol) and GLVs (nonanal) that were found to contribute by > 10% to the treatment effects in the SQT and GLV blends in the multivariate PERMANOVA analyses (Table S8). All emission rates were standardized to a temperature of 30°C. Means and errors for a-copaene, c-gurjunene and nonanal are estimated marginal means and standard errors produced by statistical models fitted to square-root-transformed data. The means and errors for palustrol were calculated from data standardized to the mean (8%) R. tomentosum plot cover and are shown as square-root-transformed.

DISCUSSION
We predicted that BVOC emissions from the mountain birch forest field layer vegetation, with the underlying soil included, would have an instant response to temperature variation during sampling, and in the long-term, increase with the warming treatment and decrease with the reduced insect herbivory. Our results support these predictions in terms of the herbivory effect: the results show that reducing leaf damage by 58-67% led to a significant, general reduction in BVOC and in particular sesquiterpene emission rates, with reductions as high as 80% observed. Non-standardized BVOC emission rates had an instant response to temperature variation in the chamber, but this only occurred in June, not in August. The warming treatment mainly modified herbivory effects on relative BVOC proportions, with clearer herbivory effects in heated than ambient plots, and affected the BVOC group blend structure; these effects were compound specific or restricted to certain samplings only and contrasting responses to warming appeared within blends. These results suggest that the instant positive effect of warmer climate on total BVOC emission rates in Subarctic ecosystems may depend on plant phenological state and that the long-term effects may be mediated through changes in insect herbivore pressure rather than changes in plant growth. It appears that if background insect herbivory increases as predicted due to the warmer climate (Kozlov 2008;Li and others 2019), BVOC emissions from these ecosystems may increase in the future.

Herbivore Control of BVOC Emissions
Silfver and others (2020) reported in their Supplementary Table 1 that 25% and 16% of birch leaves in our plots with natural herbivory were wounded by insects in 2017 and 2019, respectively, and of the wounded leaves, 85% and 64% had only 1-4% of their area damaged. In plots with reduced herbivory, 12 and 10% of leaves were wounded, respectively, and of the wounded leaves, 88 and 83% had 1-4% of leaf area damaged. These values support the current belief that non-outbreak background herbivory is low at high-latitude ecosystems (Barrio and others 2017;Kozlov 2008;Li and others 2019). Considering this, it is remarkable that total BVOC emissions decreased by 44-72% with the insecticide treatment. That these reductions closely match the magnitude of 58-67% reductions in leaf damage indicates that BVOC emission rates in these ecosystems are under a strict control of herbivory and leaf damage, despite the low herbivore pressure. Reductions in BVOC emissions also closely match the magnitude of changes in CO 2 uptake potential (26-52% increase) and plant growth (50% increase) earlier reported for the same experiment (Silfver and others 2020). Altogether these results suggest that the seemingly low background insect herbivory may significantly affect the ecosystem-atmosphere interactions in Subarctic ecosystems.
Our results are in line with earlier studies reporting increased BVOC emissions from plants with higher herbivore load and greater feeding damages (Niinemets and others 2013;Yli-Pirilä and others 2016). For instance, Yli-Pirilä and others (2016) reported increased sesquiterpene emissions from mountain birch when damaged by autumnal moth Epirrita autumnata larvae in a laboratory study. Similarly, Li and others (2019) showed that total monoterpene emissions from dwarf birch branches were increased fourfold by herbivorymimicking methyl jasmonate application in an Arctic field study. In a recent field study in Subarctic and pan-Arctic regions, Rieksta and others (2020) reported fivefold total BVOC emissions and threefold sesquiterpene emissions from mountain birch shoots with low intensity herbivory by the geometrid moths Operophtera brumata and E. autumnata. Our results add to these findings by providing evidence that natural insect herbivore communities, at their background densities, significantly increase BVOC emissions from multispecies Subarctic vegetation.
Synthetic pyrethrin application has earlier been shown to have no side-effects on B. pendula growth, leaf chemistry, leaf litter decomposition rate and soil microbial growth others 2013, 2020). We have no data of lack of pyrethrin sideeffects on the entire Subarctic community that we manipulated, since communities with no herbivores, which are needed for side-effect tests, are not achievable at the field scale. However, considering that no side-effects were earlier found on a wide array of leaf metabolites in B. pendula (Silfver and others 2013) and that the magnitude of pyrethrin application effects on BVOC production closely match the magnitude of effects on leaf damage, plant growth and photosynthesis rate (Silfver and others 2020), It is unlikely that the effects on BVOC emissions were pyrethrin side-effects. Another aspect of the herbivory treatment that warrants contemplation is whether the effects of herbivory reduction can be used for predicting the consequences of increasing herbivory in the future. In the case of BVOC emission rates, it is particularly important that herbivory reduction does not completely erase leaf damages as this would turn off herbivory-induced BVOC emissions. In such case, the observed decrease in emissions would likely lead to predictions of excessively increasing emissions with equivalent increase in leaf damage. In our study, 16-25% of birch leaves were wounded under natural herbivory, 10-12% under reduced herbivory. All of the BVOC compounds that were collected from plots with natural herbivory were also found in plots with reduced herbivory. For instance, the emission rate of trans-b-caryophyllene, a typical stress-induced compound that is a common BVOC also in healthy plants (Mä ntylä and others 2008), was on average 437 ng m -2 h -1 in plots with natural herbivory and 190 ng m -2 h -1 in plots with reduced herbivory, thus exhibiting a 57% reduction in emission rates due to herbivory reduction. Together these findings suggest that the leaf damage threshold for induced emissions was also achieved and that the reductions in BVOC emission rates we observed with decreasing leaf damage should provide reliable estimates of emission rates also for increasing leaf damage.

Warming Effects on BVOC Emissions
Previous studies in Subarctic have reported twofourfold increases in total monoterpene and sesquiterpene emissions with 1 -2.5°C warming (Faubert and others 2010;Tang and others 2018) and up to sixfold increase in total BVOC emissions with 4°C warming (Lindwall and others 2016). Equally, isoprene emissions have been found to increase 1.5-6-fold with 1 -2°C warming (Tiiva and others 2008;Tang and others 2018). Supporting these findings, we found that the non-standardized BVOC emissions from our Subarctic vegetation showed an instant increase with increasing chamber air temperature, although this response also depended on plant phenological stage: while the positive response was clear in the June sampling in both years, no response was found in the August sampling in either year. Besides being explained by plant phenology, this contrast in response could be due to plant water stress in August samplings, but since the water vapour pressure deficit (VPD) was much lower during August than June samplings (Silfver and others 2020), this is not a likely explanation. Earlier studies have shown that plant phenology has a role in BVOC emission rates. For example, in mountain birch, high early season BVOC emission rates were linked to high BVOC-storing glandular trichome density in young leaves and the emissions declined with deterioration of glandular trichomes in the late growing season (Mofikoya and others 2018).
Due to earlier reported warming-induced increases in plant CO 2 capture and growth in our field plots (Silfver and others 2020), we expected to also find a positive long-term effect of warming on BVOC emissions. However, we found no statistically significant increase in either standardized or non-standardized emission rates of total BVOCs, monoterpenes and sesquiterpenes. Leaf damage by herbivores neither increased with warming. These results indicate that under warmer conditions, in case the herbivory pressure does not increase, plants may invest the enhanced C capture to growth rather than to BVOC-based chemical defence against the herbivores. That warming effects on BVOC emissions could be controlled by such growth-defence trade-off is in accord with earlier findings: for instance, when Li and others (2019) found a twofold increase in monoterpene emissions from dwarf birch due to warming in their field study, they also found a fourfold increase in leaf damage due to warming. As such the fact that insect herbivory did not increase with warming in our study cannot be used as an evidence that herbivory would not respond to warming. Treatment plots in field experiments are vulnerable to both avoidance and congregation of freely moving animals and may thus tell little of responses of herbivores under large-scale climatic changes (Moise and Henry 2010).
Our results also reveal that compounds within BVOC blends can have opposite responses to the warming treatment. Among sesquiterpenes, a-copaene emissions decreased and c-gurjunene emissions increased with warming, which obviously can lead to a neutral response to warming at the BVOC group level. Many significant warming effects were also found at some samplings only: sesquiterpene palustrol emissions were increased by warming in 2017, but not in 2018, and the GLV nonanal emission was increased mainly in the August 2017 sampling. Finally, isoprene emissions, which in earlier studies have responded to warming (Tang and others 2018;Tiiva and others 2008), were low in our study, contributing less than 2% to total BVOC emissions and found in a half or fewer field plots in three of the four samplings. Altogether these findings suggest that while the level of insect herbivore pressure has a consistent and significant effect on total BVOC emissions (the only exception that we found for the general herbivory effect were the marginally increased GLV emissions with herbivory reduction in warmed plots in 2017), long-term warming effects are weaker, compound specific and temporally variable.
Although we cannot distinguish between aboveground and belowground emissions in our study, it is good to realize that some of the contrasting responses of single BVOCs to warming may stem from different responses of soil and foliage emissions. For instance, sesquiterpenes can also be emitted by soil microorganisms (Bourtsoukidis and others 2018) and plant litter decomposed by microbes (Svendsen and others 2018). The warminginduced c-gurjunene emissions in our study might thus partly stem from increased activity of c-gurjunene releasing soil microbes (Bourtsoukidis and others 2018) in warmed plots. Another important aspect is that besides acting as a source, soils can also act as a sink for BVOCs, and particularly for monoterpenes through microbial uptake (Tang and others 2019). We know that warming increased soil microbial activity in our experiment as nitrogen mineralization increased (Silfver and others 2020). The reason why we found no positive response of monoterpenes to warming might thus be partially explained by monoterpene emissions from vegetation being increasingly absorbed by soil microbes in the warmed plots. However, these are speculative hypotheses and need to be verified in future investigations.
Using 13 CO 2 -labelling of Subarctic dwarf birch, Ghirardo and others (2020) recently showed that the contribution of de novo monoterpene biosynthesis to monoterpene emissions increased from 60-68% to 87-95% under warming. This shift in monoterpene allocation led to negative effects of warming on C sequestration as the net exchange of CO 2 decreased and C loss as BVOCs increased (Ghirardo and others 2020). Results from our study seem to support the opposite trend as warming increased plant growth and CO 2 uptake potential by 50 and 107%, respectively (Silfver and others 2020) and still, no general increase in BVOC emissions appeared. This supports our interpretation that higher herbivory and increasing leaf damages are the crucial mechanisms that are needed to trigger warming effects on BVOC emission rates in Subarctic ecosystem.

Plant Species-specific Responses
As we focused on ecosystem level responses of BVOC emissions to the warming and insect herbivory treatments in our study, we cannot conclude much of the species-specific responses. The monoterpene myrcene and sesquiterpenes palustrol and aromadendrene have earlier been shown to be mainly emitted by R. tomentosum and further being adsorbed and re-released by birches (Himanen and others 2010; Mofikoya and others 2018). In our study, palustrol was one of those compounds that responded to warming and was almost exclusively released from plots that contained R. tomentosum, although not all plots with R. tomentosum emitted palustrol in all samplings. This finding supports the view that plant species composition and the potential changes in composition along with warming likely play a role in BVOC emission responses to climate warming (Faubert and others 2012).

CONCLUSION
The significance of our study lies in the ecosystem level approach and the natural plant and herbivore communities involved in the study. Our results reveal that insect herbivores, despite their low numbers and the seemingly minor leaf damages they currently cause in Subarctic ecosystems, have a dominant role in controlling plant BVOC emissions in these ecosystems. As insect abundances are likely to increase with climate warming, our results forecast increases in plant BVOC emission rates from high-latitude ecosystems in the future warmer climate. Considering the significant role of BVOCs in the mitigation of climate warming, it appears that plant-herbivore interactions and their effects on BVOC emission rates warrant to be considered in climatic models, over and above the instant temperature effects, to understand and foresee Arctic ecosystem-atmosphere interactions and feedbacks.