The study was conducted at three sites (ca. 100, 500 and 900 masl; Low, Mid and High elevations, respectively) along the southside of the mountain Storehaugfjellet in Kaupanger, inner Sognefjord, Norway (61.2°N, 007.2°E) from 2016 to 2018. The inner Sognefjord region is a mountainous area that offers marked topographical variation. Storehaugfjellet starts at sea level and reaches up to 1169 masl with a mean slope of ca. 30%. The bioclimatic variation along the mountain side progresses from temperate and hemiboreal mixed forest, via boreal coniferous dominated forest, to alpine environments with a timberline largely at 800- 900 masl (Moen, 1999). The specific climate of lowlands in the Sognefjord-region is characterised by an average growing season (≥ 5°C) of ca. 160-170 days, mean January temperatures of ca. 0 to -4°C and mean July temperatures of ca. 12-16°C, annual temperatures of ca. 8-6°C, and annual precipitation of < 1000mm. The vegetation of the Storehaugfjellet is dominated by pine-bilberry forest with common boreal tree species in tree layer and dwarf shrubs in the field layer. Above the timberline the same dwarf-shrubs occur, but here the vegetation was patchier with more bare rock and soil. Supplementary Table S1 , Appendix 1 gives an overview of area and site characteristics.

Vaccinium myrtillus (bilberry or European blueberry) and Vaccinium vitis-idaea (lingonberry or cowberry) are both long-lived dwarf shrubs. Bilberry is relatively fast-growing species has deciduous leaves and evergreen stems, usually 10−60 cm high (Ritchie, 1956), whereas lingonberry is a relatively small evergreen shrub of ca. 10 to 30 cm height with more conservative growth (Ritchie, 1955; Tolvanen, 1995). Both species reproduce clonally in both forest and open habitats by substantial propagation belowground, and sexually mainly by bee-pollinated seeds that are encapsuled in abundant berries (Nuortila et al., 2002) (Jacquemart, 1993). Both Vaccinium species co-occur in forests and heathlands’ in Eurasiatic boreal and alpine areas, of intermediate productivity, acidity, and light availability (Ritchie, 1956). Bilberry may constitute up to 40% of the ground cover whereas lingonberry is a bit less abundant (Kuusipalo, 1988; Økland, 1996). A recent elevational study in the same region revealed that bilberry size, cover and berry production may all peak at mid-elevations (Auestad et al., 2023).

We selected three open woodland sites at different elevations and bioclimatic zones: 1) Low: at ca. 100 masl in the hemiboreal zone, 2) Mid: at ca. 500 masl in the boreal zone, and 3) High: at ca. 900 masl in an alpine environment at the timberline ( Figure 1 ; Supplementary Table S1 , Appendix 1 , for details on sites). All sites were located along the south-facing slope of the Storehaugfjellet to reduce impact of microclimate variation on experimental responses, and were selected to meet the following criteria: a) high cover and viable populations of the study species; b) located in a woodland clearing such that shadow/light effects and variation was negligible; and c) if the clearing was a previous clear-cut (Low and Mid sites), the forestry activities should have been performed >20 years ago to avoid recent disturbance effects and minimise influence of historical site variation.

The aerial distance between Low and High site was ca. 4.5 km ( Figure 1 ). At each sites there are identical anorthosite bedrock and moraine sediments (Geological Survey of Norway, 2023), and the sites have comparable slope, light availability and soil conditions (Benevenuto et al., 2020; Supplementary Table S1 ). The combination of these features indicate that the influence of confounding abiotic or biotic factors as well as biogeographic variation, which are often a problem in latitudinal climate gradient studies (cf. Fukami and Wardle, 2005), were minimised in our study design such that differences among sites were highly likely to be due to the elevational-induced climatic variation across the sites.

The sites, with their specific elevation, were selected to cover the main altitudinal and climatic distribution of the study species, and as such represent populations living under different climatic conditions. We assumed that the Mid elevation represented close to optimal climatic conditions for bilberry in particular, whereas the low and high would represent sub-optimal climatic conditions towards warm and cold range limits, respectively. This pattern underlines the assertion that the gradient selected mainly reflects temperature variation according to the bioclimatic zones hemiboreal, boreal and alpine ( Supplementary Table S1 , Appendix 1 ), and additional recordings are undertaken to test this assumption (see Data collection). Temperatures were considerably higher (ca. 3-4 °C) in 2018 than other years, which was an extreme warm and dry year across the whole Northern Hemisphere (Vogel et al., 2019).

At all three sites we established a randomised block design with six study blocks of 4 × 4 m, minimum 10 m apart but distributed over areas of ca. 1 ha per site, in areas with relatively high cover of the two study species. We divided blocks into four 2 x 2m grids each with sixteen 50 x 50cm squares and then randomly selected one square from each using a random number generator. Each selected square was established as a permanent plot (0.25m²; marked with sticks and metal corner tubes in the ground for relocation if required) if it met the selection criteria of a relatively high cover of bilberry and presence of lingonberry, did not contain a tree, a tree stump or large bare rock, and was a minimum distance of 3 m from other permanent plots. The square block design was not possible at the High elevation due to the more pronounced patchiness of the vegetation. Therefore, plots in blocks were arranged linearly where necessary, ensuring that each plot was at least 3 m away from its nearest neighbour. The study plots within the six blocks at all sites were randomly assigned to one of the following treatments:

The warming treatment intended to simulate a mean climate warming predicted for the following decades. The elevated temperature was achieved using open top chambers (OTCs), hexagonal miniature greenhouses constructed of transparent polycarbonate panels to a standard design often used in alpine and arctic areas to experimentally warm small areas (Marion et al., 1997; Halbritter et al., 2020). The OTC had a base diameter of ca. 170 cm, and an opening diameter of ca. 110 cm leaving enough space between the plot edge and panel surface to avoid edge effects (i.e., permanent plots was located approximately in the centre of OTCs). To quantify the degree to which the temperature was elevated, soil surface temperatures were measured every hour inside OTCs in three blocks at each elevation, with temperature loggers (2016 and 2017: LogTag HAX0-8, New Zealand, 2018-and onwards: Tinytag TGP-4017, Gemini). Temperature loggers were paired with ambient data loggers located 1 m to the east or west of each measured OTC. We report mean temperature data at site level from 2016 to 2021 based on the ambient control plot recordings. We also measured soil moisture as an indication of potential drought effects both at the start and end of the project. In 2016 (25 August) and 2018 (7 June) we used an AT Delta-T moisture meter (type HH2 SM300 v 4.0, by Delta-T Devices Ltd.), probing the soil outside each of the four corners of the permanent plots.

The OTCs (and data loggers) were placed on the plots on the same days as the first methyl jasmonate (MeJA) applications in 2016 (see more below). OTCs were removed at the end of the growing season (28 October 2016 and 11 October 2017) to avoid the accumulation of snow during the winter and replaced when plots were snow free in the spring of 2017 and 2018. Thus, our warming treatment only elevated growing season temperatures. Warming treatment using OTCs may have some side-effects and thereby affect some other abiotic conditions in addition to temperature (Dunne et al., 2004; Hollister et al., 2023).

The treatment to induce plant resistance to herbivory (hereafter Resistance) intended to simulate the physiological impact that invertebrate (insect) herbivore attacks may have on plants (i.e., herbivory resistance and associated responses). To induce plant defence responses we applied MeJa to plots or plants of bilberry (Hegland et al., 2016; Seldal et al., 2017). MeJA is a volatile organic compound that is central to the activation of induced plant defence responses, particularly in the jasmonate dependent defensive pathway (Rodriguez-Saona et al., 2013), and previous studies on bilberry had shown that treatment with MeJA caused ca. 2000 genes to be upregulated and ca. 1500 genes to be downregulated (Benevenuto et al., 2019). The herbivory resistance responses induced by MeJA-application generally implies resource allocation from growth and reproduction to defence that can be observed via ecological responses on both plants and their insect herbivores (Hegland et al., 2016; Seldal et al., 2017; Benevenuto et al., 2018). The MeJA was applied on the whole plot with a 5L knapsack sprayer as in Hegland et al. (2016). Each application involved steadily spraying the liquid over the plot on three separate occasions (Low: 26 May, 3 and 14 June, Mid: 2, 10 and 19 June, High: 10, 19 and 30 June) in 2016 to ensure complete coverage and to simulate the plant response to repeated herbivory during the growing season.

At each site we recorded several environmental variables to control the assumptions of our study design, i.e. that the main variation along the elevational gradient were climatic, in addition to the general site characteristics on elevation, bedrock, soil type etc (see Supplementary Table S1 , Appendix 1 ). For this purpose we recorded variation in a) Snow cover as the snow free-dates in the permanent plots in 2017 and 2018; b) Canopy cover, i.e. shading effects, as number of 24 squares of a densiometer at each permanent plot in 2016; c) Slope as degrees with a compass at each permanent plot in 2016; d) Aspect with compass to obtain main direction in each permanent plot in 2016; e) Soil acidity as pH of soil sample from each permanent plot in 2016 diluted in water and recorded with pH-meter in lab afterwards; f) Temperatures by using the control plot recordings at site level from 2016-2018 (see Study design for details); g) Soil moisture by using the plot level recordings (see Study design for details) from the start of the study in 2016; and h) Abundance of both focal species as cover estimates in each permanent plot in 2017 (also used as a covariate in the data analysis for reproductive output).

Within each plot, four ramets of bilberry and lingonberry were selected and tagged for repeated recordings. The plants were selected systematically by placing a 0.5 × 0.5 m 16 square grid quadrat over the plot, and then inserting marker sticks at the grid cross-points closest to the four corners of the plots. The nearest ramets of each plant species to each stick was then marked with a cable tie and small coloured bead for recognition.

To investigate the site and treatment effects (Warming and Resistance) on growth of the two study species in the study plots we recorded several variables for each ramet twice during each growing season. For bilberry, we measured 1) ramet height, 2) number of leaves, and 3) number of shoots. For lingonberry, we measured 1) ramet height, 2) crown diameter at the widest point, 3) crown diameter perpendicular to 2), and 4) number of shoots. The first measurements were conducted before the treatments were applied on each elevation (Late May/early June from Low to High) and they were repeated approximately 8-9 weeks later (Mid July to early August, respectively), when plants were assumed to have reached peak growth. The recordings were repeated in the same way in 2017 and 2018, although in 2017 sampling at high was a little later because of a deeper layer of snow and a colder spring. The various size variables were used to estimate ramet size (i.e. dry mass) using the approach for non-destructive sampling of growth and biomass developed for bilberry by Hegland et al. (2010) required for such population studies. We adapted the model to the size variables recorded in this study, and developed a novel dry mass equation for lingonberry (see Appendix 2 for the model-based equations).

To examine the reproductive effects on the plant species we also recorded the number of flowers (May-June) and berries (July-August) at the plot level. In 2017 we were not able to count berries of lingonberry due to logistical problems. To study the impact on plant ramet mortality we used the recordings of whether a marked ramet survived or died (found withered or not retrieved) from autumn to spring or during the growing season to represent a total mortality at plot level. To estimate the herbivory impact on the plants, we counted number of leaves chewed by insects (not 2016 for lingonberry), at the same dates as for the plant size recordings.

Prior to statistical analysis we converted the plant size (i.e. dry mass) variables to growth data representing the change between the two sampling occasions of each year. We processed the growth and mortality data by averaging the variables at the plot level as the distributions of ramet level data for seasonal growth were extremely leptokurtic violating the assumption of normality for linear mixed modelling. We converted the insect herbivory data to a binomial variable (1 = leaves chewed, 0 = no leaves chewed) as the data contained many zeros and efforts to model them as count data lead to non-normal residuals (see below) and modelled herbivory at the ramet level using only data from the second sampling in each year (first sample was mostly too underdeveloped).

To address the hypotheses, we applied Linear Mixed-Effects modelling (LMM) to a) the growth variables as they were normally distributed and contained negative values; and b) Generalised Linear Mixed-Effects modelling (GLMM) to the mortality data (proportion of dead/alive ramets modelled with binomial error distribution), the reproduction data (number of flowers and berries; negative binomial error distribution), and the leaf herbivory (chewed or not; binomial error distribution). In all cases, we first parameterised models with the fixed effects of Year and Site, and the Warming and Resistance treatments. We used Block nested within Site as random effects, except that the herbivory model used Plot nested within Block nested within Site as this variable was evaluated at the ramet level. Interactions were included starting with three-ways interactions as these and simpler interactions were relevant to our research questions, but removed interactions if non-significant to simplify models until only significant interactions were included. We inspected the significance of model terms using analysis of deviance with type III sums of squares.

All models were built in the R programming environment (v4.1.2, RCoreTeam, 2021) (v4.1.2, using the nlme package (Pinheiro, 2010) for the linear mixed-effects models, and the glmmTMB package (Brooks et al., 2017) for the GLMMs. We subsequently performed graphical checks for the assumptions of linear mixed modelling, including variance homogeneity and normality of residuals, and normality of random effects. For the LMM we performed this step using the built-in residual functions of the nlme package, and for the GLMMs this was performed using the DhARMA package (Hartig, 2020). In some cases, the residuals showed heteroskedasticity, particularly when plotted against Site. We therefore included a variance function with fixed variances using the varIdent function of the nlme package, using Site as the variance covariate. For the models of flower and berry numbers, we included a covariate of (log transformed) percent cover of the modelled species to account for the variation in flowers and berries associated with species distribution, and we included dry mass as a covariate in models for herbivory to account for the impact that the size of the ramet may have on herbivore activity. Finally, we present the results as point estimates with 95% confidence intervals and average Least Significant Difference (LSD) intervals calculated using the predictmeans package (Luo et al., 2023). LSDs are intervals calculated using the standard error of the difference between group means and allow for visual comparison between predicted means based on scientific importance, avoiding overuse of probability values (Hector, 2021). However, LSDs were calculated for every pairwise comparison posing problems for simple graphical visualisation. We therefore plot the average LSDs produced by predictmeans for each significant interaction, and LSDs that do not overlap were considered to be significantly different at the 95% confidence level.

Temperatures were consistently several degrees higher during the growing season (April-September) at the Low (hemiboreal) site in all study years in accordance with study design ( Table 1 ), but in warm periods of the summer Mid (boreal) and High (alpine) sites exhibited similar temperature profiles ( Supplementary Figure S1 ; Appendix 3 ). The OTCs were generally effective at raising the ground temperature at all sites in 2016 (mean difference 0.57 to 1.10°C) and 2018 (0.33 – 0.74°C), but not in 2017 (year with late spring; 0.01 – 0.07°C). Generally, the manipulated temperature difference was larger at the High site ( Supplementary Table S2 , Appendix 3 ). Plots were much earlier snow free at Low, than Mid site which again were earlier than the High site ( Table 1 ). Soil moisture on the other hand was highest at Mid site, but relatively similar at Low and High sites ( Table 1 ). The OTC-plots had ca. 20% lower soil moisture than ambient plots, and the driest year (2018) was ca. 60% drier than the wettest year (2016; Supplementary Table S3 , Appendix 3 ). Canopy cover was overall very low in accordance with study design, but the mean values of densiometer was slightly higher at the Low site than the two other sites ( Table 1 ). Soil pH recordings showed that soils at all sites were acidic and showed little across site-variation which also hold true for slope and aspect.

The data from the ambient control plots showed that bilberry was most abundant at Mid (ca. 50% cover vs. <20% at Low and High in 2017), and exhibited largest reproductive output (numbers of flowers and berries) in two of three study years at the Mid site ( Appendix 4 ). However, lingonberry, had highest cover and reproductive output at Low. These results verify our expectations that Mid would be optimal for the most common Vaccinium species, bilberry, and Low relatively more optimal for the more drought tolerant lingonberry (see introduction).

Growth of bilberry populations (i.e., seasonal changes in dry mass) was weakly affected by Warming but the effect was dependent on Year and Site as indicated by two significant two-way interactions ( Table 2 ). Warming had a negative impact at Low, neutral at Mid and positive at High ( Figure 2A ). Also, the bilberry growth was highest in 2016 and Warming had a slight positive effect in 2016 and 2017, whereas there was a relatively marked negative effect on growth in the warm and dry year of 2018 ( Supplementary Figure S2A , Appendix 5 ). The resistance treatment with MeJa-application (Resistance) had no statistically significant effect on bilberry growth during any of the study years ( Table 1 ).

Lingonberry growth was only marginally affected by warming in a site-dependent way ( Table 2 , Warming x Site: p = 0.065), and Figure 2B indicates a small positive effect of warming at the High, neutral impact at Mid and negative effect at the Low in a similar pattern to that found for bilberry. There was also a significant two-way interaction between Site and Year, such that lingonberry growth at Mid fluctuated across years with highest growth in 2017, and lowest in 2016 ( Supplementary Figure S2B , Appendix 5 ).

The warming treatment had an impact on mortality of both focal species ( Table 3 ). For bilberry, there was a significant three-way interaction between Year, Site and Warming, suggesting that in 2018 the warmed plots experienced higher mortality at Low, but lower at Mid and very low at High (6 ramets in total), indicating decreasing mortality with elevation in warmed plots ( Figure 3A ). For lingonberry ramets, there was an interaction between Site and Warming ( Table 3 ). The warmed plots had an overall declining mortality with elevation (site), whereas the control (ambient) at Mid stands out with very low mortality. As such lingonberry ramets in the warmed plots at Mid experienced higher mortality than in controls, but lower mortality than in controls at High ( Figure 3B ).

Reproductive data showed several significant two-way interactions that signalled strong treatment effects ( Table 4 ). Bilberry flower numbers were strongly reduced by the Resistance treatment, but this effect depended on Warming, Site, and Year (all two-way interactions in Table 4 ). In other words, flower numbers of bilberry were most reduced by MeJa-application in warmed plots ( Figure 4A ), with increasing elevation from Low to High sites ( Figure 4B ), the year after treatment, i.e., in 2017 ( Supplementary Figure S3B , Appendix 5 ), and the site-dependent effect was largely due to the variation in flowers across sites in 2018 ( Supplementary Figure S3A, Appendix 5 ), Lingonberry flower production was also depressed by the Resistance treatment and here the effect also increased along the elevational gradient to the extent that at High no flowers were produced in Resistance plots ( Table 4 , Figure 4C ). There was also a Year x Resistance interaction, with the strongest negative effect of MeJA-application on lingonberry flower numbers in 2017, the year after application ( Supplementary Figure S5 , Appendix 5 ).

Berry numbers in bilberry were negatively impacted by both Resistance and Warming, but the negative impact of MeJa-application appeared strongest in the warmed plots ( Figure 5A , interaction in Table 3 ). The Site x Warming interaction suggested that berry numbers were depressed in warmed plots at both Mid and High, but not at the Low site ( Figure 5B ). There was a Year effect of Resistance treatment as berry numbers decreased in 2017, the year after MeJa was applied ( Supplementary Figure S4A , Appendix 5 ) whereas for the other years the effects varied more across sites with no consistent pattern ( Supplementary Figure S4B , Appendix 5 ).

The Resistance and Warming treatments did not appear to interact on the berry production of lingonberry ( Table 4 ). A Year by Site interaction showed that berry output increased along the climatic gradient in 2016, but by 2018 the sites had similar berry numbers ( Table 4 ; figure not shown). However, the lack of lingonberry berry data in 2017 makes interpreting these results difficult (see methods).

There were several two-way interactions for the probability of plants to experience leaf herbivory involving all fixed factors, but only for bilberry ( Table 5 ). Bilberry ramets were more likely to experience leaf herbivory in Warming plots at Low and Mid, but decreased probability at High ( Figure 6 ). Also, there was a decreased likelihood for leaf herbivory in Resistance plots the year after treatment in 2017, but then herbivory was more likely in these plots in 2018 ( Supplementary Figure S6A , Appendix 5 ). Lastly, herbivory was most likely at Low and Mid in the first two years of study, but more likely at High in the warm year of 2018 ( Supplementary Figure S6B , Appendix 5 ). The probability of herbivory of lingonberry leaves was not affected by any of the experimental treatments ( Table 5 ).

The interaction between site and experimental warming, suggested a likely climate- dependency of the warming treatment, supporting our predictions, with reduced aboveground growth of both dwarf-shrub species in lowland and increased growth at the alpine site ( Table 5 ). Although both species experienced reduced growth in lowland OTCs, the warming induced decline in bilberry growth across the three year-study period may also hint towards a possible species-dependent warming response that deserves further attention. The increased growth of dwarf shrubs in general as well as for Vaccinium species is common in tundra or timberline sites (Shevtsova et al., 1997; Walker et al., 2006; Anadon-Rosell et al., 2017; Fazlioglu and Wan, 2021). Most warming studies, however, do not offer contrasts with lowland or southern populations of same plant species, but Lindborg et al. (2021) found that northern-most populations of a common wetland species (Caltha palustris) showed the strongest growth responses to experimental warming along a latitudinal climate gradient. The contrast in warming-induced growth responses between high- and lowland populations ( Table 5 ) may support theory that predicts increased growth with climate warming if other resources (i.e., moisture, nutrients etc) are not limiting. Populations at the alpine site are likely to be limited by temperature and hence a short growing season. The additional warmth afforded by the OTCs may therefore have improved conditions at this alpine site, such as extending the growing season via earlier vegetative growth and enhancing nitrogen availability (Rinnan et al., 2009; Anadon-Rosell et al., 2017). Also, the decreased mortality of both bilberry (in warm year 2018) and lingonberry (overall) in particular with warming at the alpine site highlights how warmer climate may benefit these shrub species in cold conditions and suggest possibilities of future population growth and expansion.

At the lowland site, closer to the warm-range limit of the dwarf shrubs, growth may be limited by water availability, and we reported ca. 20% less soil moisture in lowland than at mid elevational site and. The decreased growth, and higher mortality of bilberry, at the lowland site therefore lends support to predictions of stress-induced responses in plant populations due to elevated temperatures and/or drought (Jamieson et al., 2012). The larger effect of warming on bilberry growth in the extremely dry year of 2018 underlines this potential warming lowland effect, comparable to growth reductions and increased mortality in populations of Norway spruce (Picea abies) under warm and dry conditions and opposite responses towards colder distribution limits (Andreassen et al., 2006). It has already been shown that variable snow conditions and reduced snow cover are phenomena that together with winter warming events in winter increase mortality and reduce cover of Vaccinium species, particularly bilberry (Odland and Munkejord, 2008; Bokhorst et al., 2010, Bokhorst et al., 2011; Taulavuori et al., 2011; Kreyling et al., 2012). The decreased lowland plant growth, and for bilberry also a distinct negative-temporal warming effect and increased shoot mortality, in response to growing-season warming in our study may add to the list of responses and factors that may shape the common phenomenon of “trailing edge” populations (Parmesan, 2006). Such demographic trends may eventually result in reduced abundance of dwarf shrubs in temperate vs. boreal areas (Hedwall and Brunet, 2016), and perhaps even in boreal vs. alpine areas in future (Puchałka et al., 2023).

Contrary to our expectations, flowering showed little site dependency and variation across populations, and was thus negatively impacted by experimental warming in general. For fruit numbers only bilberry showed a significant negative effect to OTC-warming ( Table 5 ; controlled for abundance). Most OTC experiments in tundra locations report increased flowering numbers for shrubs under warming although these may be species- or trait-dependent (Arft et al., 1999; Fazlioglu and Wan, 2021). In an alpine study Alatalo et al. (2021) discovered that warming effects generally decreased fruit production in deciduous shrubs but increased it in evergreen shrubs. Our species-dependent fruiting effects were comparable to these findings and in line with our prediction that the faster-growing deciduous bilberry would experience more negative warming effects on reproduction than the more slow-growing evergreen lingonberry. In addition, the more pronounced negative effect of warming on berry production in bilberry at the alpine site, coupled with an increase in growth, may indicate a trade-off resource allocation strategy by bilberry, similar to that shown by Cassiope tetragona (Johnstone and Henry, 1997). However, it is noteworthy that while we have manipulated the conditions faced by plants at our sites, we have not been able to alter conditions for their insect pollinators and the results on reproductive output may thus not represent the full population response to climate warming in these outcrossed species. Our warming treatment may, for example, have altered flowering phenology, which can impact flower visitation (Gillespie and Cooper, 2022) and we have not accounted for interannual variations in local site-level pollinator communities. Similarly, we cannot rule out the possibility that the OTCs interfered with pollinators or pollination, although the evidence about this effect is equivocal (Hollister et al., 2023).

The lack of a site-dependent effect of warming on flowering in bilberry (and lingonberry) populations, in contrast to fruiting, may indicate that the climate impact on flower buds was mainly determined by climatic influence in earlier years or by other factors. Indeed, correlation studies on long-term harvest index data have shown that a warm autumn, which is the main period of flower bud formation in bilberry, may reduce reproduction in the subsequent year (Selås, 2000; Selås et al., 2015). Flowering and fruiting in bilberry could also be affected by a range of other factors that could not be assessed in our short-term study with only one site per elevation, such as previous summer temperatures (+), winter warming events (-), late winter snow cover (+), winter herbivory (-), summer herbivory (-), spring frost (-), pollination (+) and drought during ripening (-) (Jacquemart, 1997; Tolvanen, 1997; Selås, 2000; Bokhorst et al., 2011; Boulanger-Lapointe et al., 2017; Lilleeng et al., 2021).

Leaf herbivory did not vary with site, but in bilberry it increased with warming in the mildest climatic sites (Low and Mid). The difference among species may reflect the different palatability and attractiveness of the focal dwarf shrub species due to their deciduous vs. evergreen plant traits (Tolvanen, 1994). However, the decreased leaf herbivory on bilberry when experimentally warmed at the high site was rather unexpected as insect herbivory was hypothesised to increase with increasing temperature (Rasmann et al., 2014; Birkemoe et al., 2016). The most likely explanation for the climate-dependent effect is that warming, in particular in colder environments, could affect the plant species ability to defend itself or that leaf nutrient content changes more in response to warming at higher elevations (Anadon-Rosell et al., 2017).

The decreased reproduction (both species; Table 6 ) and leaf herbivory (bilberry) particular the year after treatment ( Supplementary Figure S3–S5 , Appendix 5 ) that we observed in response to MeJa-application occurred as hypothesised and followed previous experience with this chemical treatment (Hegland et al., 2016; Seldal et al., 2017). The stronger resistance effects in the deciduous bilberry than the evergreen lingonberry ( Table 6 ) was also as expected (Tolvanen, 1994; Taulavuori et al., 2013). The climate dependency in the herbivory resistance responses were evident through the overall more pronounced reduced reproduction at the alpine than in the lowland site. According to optimal defence theory plants may often have relatively lower levels of constitutive defences and higher inducibility when attacks are relatively rare, and there is lower risk but higher costs involved in saving resources (McKey, 1974; Zangerl and Rutledge, 1996; Rasmann et al., 2014). In our study system, leaf herbivory decreased in likelihood along the climatic gradient, and the corresponding increased response to MeJa-application under colder conditions is in line with this theory of higher inducibility when herbivory risk is relatively low, or attacks vary in time or space.

The effect of MeJA-application was thus strongest for bilberry as expected, and the combined treatment effect of warming and herbivory resistance only had a synergistic effect on reproduction in bilberry. This type of combined treatment in the field has, as far as we know, never been explored. We would argue that the effect itself can be understood as an extension of the findings for the herbivory resistance treatment along the climate gradient, and thus explained by optimal defence theory. That is, only when the risk of herbivory increases (as with warming) the plants should theoretically invest more in defence. Moreover, for species living in cold climates, warmth is often a limiting resource, and as such the combined effect could also be a result of improved herbivory resistance at the cost of reproduction. In several cases it has been observed that plant species of lower elevations often have the highest inducibility (Pellissier et al., 2016; Benevenuto et al., 2020) but when warmed, as in our study, the growing conditions and resource availability at higher elevations are likely to approach those of lowlands and plants respond with higher resistance. Several authors have stressed the importance of increased herbivory and even insect outbreaks when the climate is warming (Tylianakis et al., 2008; Jamieson et al., 2012; Jepsen et al., 2013; McDowell et al., 2020). This combined effect is of special importance in boreal and tundra ecosystems, as outbreaks and dynamic population fluctuations are considered more frequent here (Ims et al., 2008; Kausrud et al., 2008). Although herbivory affects plant performance directly, indirect effects caused by herbivory resistance and immune responses are also likely to influence plant responses (Jamieson et al., 2012; Rasmann et al., 2014), which we attempted to simulate in our combined treatment. The synergistic negative effects of drought and herbivory on plant growth and mortality are well documented (Zandalinas and Mittler, 2022), and although combined effects have been predicted to affect plant reproduction there are still few recorded examples. The impact of multiple stressors on plant reproduction, both flowering and fruiting (as in our study), in a warming world may thus be considerable and deserves more attention.

Although bilberry responded more clearly to treatment and treatment combinations in this study, lingonberry still showed climate-dependent effects of herbivory resistance on its flowering. The zero flowers at High in 2017, one year after MeJa-application, may suggest an inducible defence system that redirects resources away from reproduction in lingonberry, albeit weaker than in the deciduous and faster-growing bilberry. As with bilberry, flowering buds are formed the previous autumn for lingonberry (Tolvanen and Laine, 1997) and the reduction in flowers in the year following treatment of both study species suggests that the resistance treatment may have impacted this resource allocation immediately. The resource costs involved in reproduction are likely to be greater in lingonberry than for bilberry, and bilberry can probably recover more rapidly from herbivory than can lingonberry (Tolvanen and Laine, 1997). The finding of a potential inducible defence in lingonberry contrasts earlier reporting of the plant species as having mainly constitutive defence (Tolvanen, 1994; Taulavuori et al., 2013), although its evergreen American relative Vaccinium macrocarpon (cranberry) has been reported to have such an inducible defence system (Rodriguez-Saona et al., 2013). The lack of reduced herbivory, as well as growth and berry production) in lingonberry after MeJa-application, may be an argument against an induced defence system in this species, but further investigations are warranted.

Our study sites were selected to represent a range of bioclimatic conditions for our focal plant species, from a hemiboreal lowland, via a boreal to an alpine site, to investigate whether experimental warming and herbivory resistance may have species- and context- dependent effects. The variation in temperatures across sites were consistent with the general predicted decline of ca 0.6-0.7°C per 100 m altitude (Wallace and Hobbs, 2006) and the snow cover data, that span ca. 40 days difference between the lowland and the alpine site, also indicates that the variation across study sites clearly represent differences in climatic factors. Although our study design was specifically targeted to avoid confounding effects, we recorded many other non-climatic environmental variables that overall were similar across sites ( Table 1 , Supplementary Table S1 in Appendix 1 ). While we cannot fully rule out that other environmental variables also may be influential, we argue that our three sites primarily represent a climatic gradient, similar to that used and interpreted as such in other studies (Rasmann et al., Descombes et al., Benevenuto et al., Fukami and Wardle, 2005). Moreover, we did not replicate the elevational sites due to practical (site availability) and resource restrictions (labour and funding). As such we may generalise our findings to the population level of the studied plant species, but our results should be inferred with care when considering the effects of climate warming and herbivory resistance beyond the site level.