Experimental warming increases the vulnerability of highelevation plant populations to a specialist herbivore

may differ in their interactions with herbivores under climate warming. Without genetic adaptation, existing physiological responses of high-elevation populations to warmer temperatures may leave these populations vulnerable to the increases in herbivore pressure predicted under climate change.

Our understanding of these indirect effects of temperature is currently limited to differences among species, rather than variation within species, yet the latter is more useful for understanding adaptive responses of plants to future changes in the biotic environment (e.g.Anderson et al., 2015;Halbritter et al., 2018).Exploring differences in temperature-mediated herbivore resistance of populations from high-and low-elevations can provide insight into the relative vulnerability of high-elevation populations to ongoing climate change, as well as the need and potential capacity for evolutionary adaptation in high-elevation populations facing climate change.
Warmer temperatures have previously been shown both to increase herbivore activity (Bale et al., 2002;Lemoine et al., 2014) and to alter plant defence trait expression and resistance to herbivory (Jamieson et al., 2012;Kissen et al., 2016;Paudel et al., 2020;Pincebourde et al., 2017), although relevant empirical studies often do not account for the possibility that plant phenotypic responses to temperature variation may affect herbivore feeding rates and performance.Those that explicitly test this effect suggest that by altering plant nutritional quality and defence investment, plant responses to warmer temperatures can reduce larval performance (DeLucia et al., 2012;Jamieson et al., 2015;Pham & Hwang, 2020;Zvereva & Kozlov, 2006;Kharouba & Yang, 2021).Moreover, the magnitude of such

Accepted Article Accepted Article
plant-mediated effects of temperature on herbivore feeding rates and performance can be comparable to the direct effects of temperature on herbivores (Bauerfeind & Fischer, 2013;Nechols et al., 2020;Pham & Hwang, 2020), suggesting they may significantly influence plant vulnerability to herbivory under climate warming.There is also considerable evidence that the effect of warming on herbivore resistance varies among plant species (Descombes et al., 2020;Hamann et al., 2021;Zhang et al., 2019); however, relatively little is known about the specific traits responsible for such variation or about the extent of variation in temperature-mediated herbivore resistance among populations within species.
Adaptation to different thermal environments can yield divergent trait plasticity in response to temperature change (De Kort et al., 2020;Drake et al., 2015;McKay et al., 2001;Münzbergová et al., 2017;Peterson et al., 2018); for example, warm and cold-adapted species can differ in their phenological sensitivity to temperature variation (Gugger et al., 2015), or show divergent gene expression profiles when grown at different locations along an elevation gradient (Walter et al., 2020).Consequently, examining within-species variation in temperature-mediated trait plasticity along thermal gradients can provide insight into potential adaptation to changing climate (e.g.Gugger et al., 2015;Schmid et al., 2017).Populations in cooler parts of species ranges often display elevated rates of growth under warmer temperatures as they escape their lower physiological limits (e.g.Drake et al., 2015;Midolo & Wellstein, 2020); however, high-elevation populations transplanted to lower elevations can also show reduced rates of survival, as they are exposed to novel biotic and abiotic stressors (Midolo & Wellstein, 2020).Adaptation to cold environments may therefore result in maladaptive or attenuated phenotypic responses to climate warming (De Kort et al., 2020;Drake et al., 2015;Münzbergová et al., 2017), the challenges of which might be exacerbated by exposure to more frequent and intense interactions with invertebrate herbivores.
Alpine plants are particularly useful for examining the extent of within-species variation in responses to changing thermal environments as they are exposed to contrasting temperature and herbivore pressure across steep altitudinal gradients over short geographic distances (Halbritter et al., 2018;Moreira et al., 2018).Populations from different elevations frequently diverge in traits relevant to herbivore resistance, including leaf trichome density (Buckley, Widmer, et al., 2019;Løe et al., 2007) and chemical defence production (Buckley, Pashalidou, et al., 2019;Rasmann, Buri, et al., 2014;Wagner & Mitchell-Olds, 2018), but also traits, such as primary metabolism (Alonso-Amelot, 2008;Moreira et al., 2018), that indirectly influence resistance.While recent reviews suggest that resistance to herbivory and rates of herbivory do not consistently vary with the elevational origin of species and populations (Moreira et al., 2018;Rasmann, Pellissier, et al., 2014), several studies support the hypothesis that high-elevation populations are more vulnerable to herbivory than

Accepted Article Accepted Article
populations at low elevations, consistent with adaptation to low herbivore pressure at high altitudes (Buckley, Widmer, et al., 2019;Callis-Duehl et al., 2017;Descombes et al., 2017;Pellissier et al., 2012Pellissier et al., , 2014)).To date, however, only a few experiments have explicitly tested the role of temperature in shaping elevational trends in defence trait variation (Dostálek et al., 2020;Pellissier et al., 2014).Pellissier et al. (2014) found that lower temperatures reduced chemical defence production in Plantago lanceolata, though this did not impact performance of a generalist herbivore; rather, high-elevation plants were more susceptible to herbivory irrespective of growing temperature.By contrast, Dostálek et al. (2020) found higher rates of herbivory on Impatiens species grown under warmer conditions but could only explain part of this variation in herbivore performance through changes in traits such as leaf area and density.
Both these studies suggest that temperature can alter herbivore performance, but neither provides clear evidence for elevational divergence in responses to temperature, potentially due to limited population-level sampling within plant species.
To address this lack of knowledge regarding the extent of within-species variation in temperature-mediated herbivore resistance and therefore the relative vulnerability of highelevation populations, we grew populations of an alpine plant with a broad elevational range under different temperatures and then exposed them to a specialist herbivore.This enabled us to test whether population responses to temperature varied with elevation and whether this altered interactions with an herbivore.The alpine plant Arabis alpina shows elevational trends in growth and defence-related traits, with high-elevation plants displaying smaller rosettes, reduced morphological plasticity and greater defence inducibility (Buckley, Widmer, et al., 2019;de Villemereuil et al., 2018).Field surveys also suggest that herbivore pressure, as assessed by leaf damage, declines with increasing elevation (Buckley et al. 2019), with specialist herbivores notably absent from the highest elevation sites (unpublished data).In a previous study, Buckley et al. 2019 found that high-elevation populations of A. alpina are more palatable to a specialist herbivore, on average, than low-elevation populations; however, the role of abiotic factors, specifically temperature, in shaping the observed trait variation was not explored.The aims of the current study were therefore to: (1) assess elevational trends and population-level variation in phenological, chemical and defence traits under different growing temperatures; (2) relate this temperature-driven trait plasticity to the performance of a specialist herbivore on sampled leaf tissue; and (3) examine whether elevational trends in trait plasticity and herbivore resistance might have implications for responses of highelevation populations to climate change.The results presented below provide insight into the indirect effects of temperature change on plant-herbivore interactions, and also highlight the potential vulnerability of high-elevation populations to the increased levels of herbivory predicted under future climate warming.

Plant material and specialist herbivore used in study
Seeds were collected from nine field sites across mountain regions of Switzerland, from low (~1000m above sea level) to high elvations (~3000m a.s.l.) (Figure 1a; Table S1 in supporting information).At each site, ripe fruits were sampled (in either 2015 or 2016) from individual plants separated by at least 2m, as described in Buckley et al. (2019), then stored in darkness at room temperature.As the paternal contribution to these seeds is unknown, we refer to seeds collected from each mother plant as a maternal family.
Pieris brassicae is a specialist on host plants in the Brassicaceae family and larvae were sourced from a lab colony at ETH Zurich reared on cultivated Brassica oleracea.Little is known about the interaction of A. alpina and specialist Brassicaceae herbivores, though other Pieridae (e.g.P. napi) have been observed feeding on A. alpina in the field (JB personal observation).First instar larvae soon after hatching were used in the larval performance assay.

Experimental growing conditions
For all experiments, seeds were stratified for one week at 4°C (8hrs light; 16hrs dark) in 54-cell (6 x 9) trays filled with a low nutrient alpine wildflower soil mix (Ricoter, Switzerland; ingredients given in Supporting Information Appendix S1) before being moved to different temperature treatments for germination.We used "warm" and "cold" growing conditions approximating low-elevation and intermediate-elevation summer temperatures (based on weather station data from MeteoSchweiz; https://www.meteoschweiz.admin.ch).
Average summer air temperature from 1 st June-30 th August during the day (8am-8pm) and night (8pm-8am) was calculated using data from four different weather stations at 800-1200m above sea-level (low-elevation) and 1800-2200m a.s.l (intermediate-elevation).Warm conditions were thereby defined as 19°C:14°C (day:night) and cold conditions as 13°C:10°C (Figure 1b).Temperature data was available from one high-elevation weather station at 2691m (8°C: 6°C day:night), but such low temperatures were difficult to maintain in the available growth chambers.Due to space limitations, we used one large growth chamber for each set of conditions in each experiment, though these chambers were of identical specification (Kälte 3000, Landquart, Switzerland) and located adjacent to one another in the same room, minimising potential confounding effects of cabinet location.Temperature, light intensity (22Klux) and humidity (50% day, 60% night) were continuously monitored by internal sensors.When possible, we moved plants between chambers and reprogrammed

Accepted Article Accepted Article
temperature conditions (switching which chamber was warm and cold) to minimise chamber effects.
Once seeds had germinated (1-2 weeks), excess seedlings were thinned, removing the smallest seedlings to reach the number required.Tray position within each growth cabinet was rotated every 2-3 days until plants were repotted to minimise any confounding effects of position in chamber.In an initial experiment, we grew plants under long daylength regimes comparable to a mid-summer time point in Switzerland (15 th July: ~15.5hrs light, 5.5hrs dark, 1.5hr transition either side); however, this regime resulted in widespread stress-induced flowering, particularly within certain populations (Appendix S1; Table S2; Figure S2; data in Appendix S2).To enable consistent vegetative tissue sampling, we therefore shifted to an equal 12hr:12hr light:dark photoperiod for experiments 1 and 2 (as used in Buckley et al., 2019).

Experiment 1: Constitutive variation in the metabolome among populations and under different temperature regimes
Plants from the nine study populations were grown from seed in a randomised arrangement in growth chambers set to different temperature regimes.We germinated seeds from six maternal families from each population in a random arrangement across a 54-cell tray, aiming for one individual per family per temperature treatment.The same randomised germination layout was used for both growing temperatures (Figure S1b).After 44d (warm) or 65d (cold), when plants were at similar stages of development (visually assessed by rosette size and root development), plants were transferred from trays to 7cm clay pots.Then at 84d (warm) and 99d (cold) after stratification two fully expanded leaves from one plant per family were sampled into liquid Nitrogen and transferred to a -80°C freezer.The leaves were then freezedried and finely ground to a powder in a Genogrinder (Spex Sample Prep, Metuchen, NJ, USA) with steel grinding beads, and 9.8-10.2mgground dry leaf tissue was weighed into 4ml glass vials for metabolomics extraction.We used leaf material from 6 individuals per population per temperature treatment (N=54 in warm conditions; but N=52 due to lower germination rates in cold conditions).
We performed a non-targeted metabolomics extraction adapted from a published protocol (Broeckling et al., 2005) to identify a broad diversity of polar and non-polar metabolites (detailed protocol in Appendix S1).Briefly, non-polar and polar compounds were extracted by respectively adding 1.5mL chloroform with a docosanol standard at 10μg/mL or 1mL HPLC grade water (MilliQ: Merck, Darmstadt, Germany) with 25μg/mL ribitol standard to an initial extraction in 70% methanol.The non-polar and polar phases were then separated, dried down and derivatised with MSTFA + 1% TMCS (Sigma-Aldrich, St. Louis, MI), as

Accepted Article Accepted Article
described in Appendix S1.Samples were extracted in three randomised batches and then 1μL of each derivatised and concentrated sample was injected in a randomised order with an automatic Agilent injector 7693 autosampler (Santa Clara, CA, USA) into an Agilent 7890B (Santa Clara, CA, USA) GC.The sample was injected splitless, but after separation split and analysed by both the connected MS Agilent 5977A (Santa Clara, CA, USA) and an FID, allowing for quantification.We used an initial temperature of 70C and a 5C/minute ramp to 325C at a constant flow rate of 0.7mL/min and an Agilent HP-5MS UI column (30m length, 0.25mm diameter, film 0.25m).During the GC/MS run, four extraction blanks and six pyridine blanks were run to check for contamination.Unfortunately, a problem with the nonpolar extractions resulted in poor quality chromatograms, so we focussed subsequent analyses on only the polar metabolites.
Samples were processed using Agilent MassHunter Qualitative and Quantitative software (Agilent Technologies, Inc. 2008) to remove unreliable peaks (difficult to consistently integrate across samples), leaving 52 robustly integrated polar metabolites for further analysis.Where possible, metabolites were annotated using matching scores and visual comparisons of acquired fragmentation patterns to those available in the NIST2014 MS library, otherwise they were labelled 'unidentified_1' to 'unidentified_X'.Based on these identifications, compounds were also allocated a KEGG Brite classification (for metabolism), to examine changes in groups of similar metabolites (Table S3).Final quantifications were in ng/L relative to the amount of the internal standard ribitol.For each metabolite, the average peak area for the extraction blanks was subtracted from the peak area for each sample.This corrected data matrix was used for all subsequent analyses (and is provided in Appendix S3).
Experiment 2: The effects of growing temperature on anti-herbivore defences and performance of a specialist herbivore Using the same nine populations and temperature conditions as in experiment 1 (with a 12hr:12hr photoperiod), we next tested whether growing temperature and plant elevation of origin might impact the performance of a specialist herbivore (larval mass) and investment in defence-related traits (constitutive glucosinolates, trichome density and specific leaf area).We germinated seeds from 7-8 maternal families per population, where possible using the same families and experimental design as in experiment 1 (Figure S1c).After 27d (warm treatment) or 41d (cold treatment) growth, the seedlings were transplanted to small clay pots.We haphazardly sampled one plant per family for measuring glucosinolates, trichome density and leaf density (SLA), and used the remaining plants from each population for the larval performance assays.Due to limited seed availability and variable germination rates,

Accepted Article Accepted Article
particularly under the cold conditions, the number of plant families used for the larval performance assay varied from that used for trait measurements (see Table S4).Trait measurement data from experiment 2 is provided in Appendix S4.
Plant leaf material for trait measurements was sampled in two sets (each containing half the plants from each population grown under the two temperature treatments) 74d and 76d post-stratification over a 3hr period from 11.00-14.30.Three large similar-sized leaves per rosette were sampled, frozen in liquid nitrogen and then freeze-dried for 24hrs.Rosette width and phenological status was recorded and two leaf discs were sampled from two leaves per plant using a 6mm borer (a total of 4 per plant).The leaf discs were frozen and freezedried before weighing on a MT5 balance (Mettler Toledo, Greifensee, Switzerland) to the nearest microgram.Two freeze-dried leaf discs per plant were then photographed using a Leica MC170 HD camera attached to a LEICA M420 microscope (Leica Biosystems, Nussloch, Germany) and the number of trichomes counted on both the lower and upper leaf surfaces.
Glucosinolate extractions were performed as described in a recent HPLC protocol (Grosser & van Dam, 2017), with minor modifications as described in Appendix S1.
Glucosinolates were annotated using MS spectra and UV profiles (see Table S5), and independently confirmed through MS spectra by Michael Reichelt (Max Planck Institute).
Coelution of two glucosinolates at an early retention time made it difficult to quantify these separately.Extracted ion chromatograms for their major desulfo ions (progoitrin m/z = 308.08;glucoraphanin m/z = 356.09) in the Agilent MassHunter quantitative software were therefore used to estimate relative amounts of each compound at this retention time.
Automatic peak integrations by the MassHunter quantitative software were manually inspected and corrected where necessary.
For the larval performance assay, we used plants at 74 days post stratification and L1 instar Pieris brassicae larvae from a lab colony at ETH Zurich.The assay was performed in a greenhouse chamber set to warmer temperatures than those used for growing plants (22°C: 19°C light:dark).Three microboxes of 11cm diameter (Sac O2, Deinze, Belgium) were used per population per temperature treatment, with each containing a piece of moistened filter paper and 10 L1 larvae.In total, there were 54 pots randomly arranged across a shaded greenhouse bench.Every 1-2d, leaf tissue was sampled from 1-2 haphazardly selected plants from each of the nine populations grown under the two temperatures (each represented by plants from multiple maternal families) and then immediately fed to larvae ad libitum in the microboxes.At least one individual of each maternal family (per population) was sampled over the course of the performance assay.After harvesting leaf material to feed to larvae, plants were discarded.On day 1, we provided leaf discs to larvae, but due to limited plant

Accepted Article Accepted Article
availability, we then fed whole leaves to larvae from day 2 onwards.The performance assay lasted for 8 days (25 th June to 3 rd July), although three pots in the cold treatment (SFH3-C, SFH1-C, 12-1-C) were frozen after 6 days due to insufficient leaf tissue from the available plants; these microboxes were excluded from statistical analyses.Growth rates of P. brassicae are relatively slow on A. alpina compared to B. oleraceae (JB personal observation); most larvae did not develop beyond the L2 stage during the experiment.At the end of the experiment, larvae were frozen at -20°C and freeze-dried before being weighed to the nearest microgram on a MT5 balance (Mettler Toledo, Greifensee, Switzerland).

Statistical analysis
Model selection: For all models below containing multiple explanatory factors, we identified minimal adequate models using backwards model selection, removing non-significant factors and comparing nested models using likelihood ratio tests until only significant variables remained (Crawley, 2013).Adjusted R 2 values represent the proportion of variance explained by the minimal models.All data analyses were done using the base stats package of R statistical software (R Core Team, 2020) unless otherwise specified.Linear mixed effects models were implemented using lme4 (Bates et al., 2015) and R 2 estimates for mixed models obtained using the R package piecewiseSEM (Lefcheck, 2016).
We tested for interactions between growing temperature and elevation to assess whether populations from different elevations differed in their phenotypic response to temperature variation (suggesting variation in trait plasticity across populations).We implemented models both with population elevation as a continuous variable (in metres above sea level) or population as a fixed factor with nine levels.Previous work identified significant population-level variation independent of elevation (Buckley et al. 2019), so we wanted to again assess the relative importance of population and elevation as explanatory factors.Model fit was assessed by examining plots of residuals against fitted values and quantile-quantile plots to identify deviations from that expected under normality.If necessary, squareroot or log transformations of the response variable were used to improve model fit.Details about the models used and relevant covariates are described for each experiment below.
Experiment 1: We first used multivariate analysis to test for population divergence in the polar metabolome under the two growing temperatures.A small number of plants were bolting or flowering at time of sampling (13 out of 106 plants), so we controlled for this factor in analyses (defining plants as vegetative or reproductive).Matrices of amounts of polar metabolites were first analysed by PCA using the function prcomp (variables were scaled to have unit variance).Sample groups identified via PCA were similar to those using

Accepted Article Accepted Article
NMDS analysis (Figure S3; metaMDS function in the vegan R package; Oskanen et al., 2019), so we focused on the PCA output for further analysis.We then tested whether different experimental factors explained multivariate physiological responses using both PERMANOVA (testing for the marginal effects of each factor; vegan R package) and then linear models with the first three principal components (PCs) as response variables.Two different models were fitted for each response variable: Identifying traits explaining variation in larval mass: To assess which of the above traits might best explain observed variation in larval performance, we tested whether larval performance could be predicted by the morphological and chemical traits analysed above.We included SLA and trichome density on the lower leaf surface, as well as the first three principal components of polar metabolite variation and first three principal components of glucosinolate variation.As larvae were fed leaf material from multiple maternal families per population, we could not use maternal family as the unit of replication.Instead, we estimated mean values of all above traits and larval performance at the level of population and temperature (N = 18, except for morphological traits where data was missing from population FC leaving N = 16).Using mean values also accounted for some differences in the maternal families used between experiments 1 and 2 (due to insufficient seeds for some high-elevation families).The following full model was fitted and simplified using backwards selection as described at the beginning of the statistical analysis section.

Experiment 1: Temperature and plant population explain observed variation in the composition of polar leaf metabolites
Plots of PC1 vs PC2 showed clear clustering of plants grown under different temperatures (Figure 2a).Consistent with these patterns, PERMANOVA revealed no significant interaction between growing temperature and elevation or population (marginal effect of elevation * temperature: R 2 = 0.01, F = 1.36,P = 0.241; population * temperature: R 2 = 0.04, F = 1.06,P = 0.373), but did find significant main effects of temperature and elevation (or population) on variation in polar metabolite composition (Table 1a,b).
The first principal component, explaining 19.2% polar metabolite variance, reflected a significant effect of temperature, elevation and phenological status (R 2 = 0.58; Table 2).
Similarly, when population replaced elevation in the model, variation in PC1 was best

Accepted Article Accepted Article
explained by significant effects of temperature, population and phenological status (R 2 = 0.72).The interaction between temperature and elevation (or population) did not explain variation in any of the first three principal components.Rather, variance in PC2 was explained by significant effects of temperature and PC3 explained by significant effects of elevation or population together with extraction set (Table 2).When the metabolites were grouped by class, it was clear that changes in total disaccharides dominated changes in total metabolite abundance (Figure 2b), with the change dominated by a 2.1-fold increase in one disaccharide (sucrose) under cold temperatures (Table S3).The annotated metabolites showed a wide range of responses to temperature, with 30 metabolites increasing under warmer conditions and 21 decreasing under warmer conditions (one showed no noticeable change).Principal component loadings suggest widespread changes in most annotated metabolites in PC1 (Figure S4a), but PC2 and PC3 are shaped by changes in specific metabolite groups.Specifically, PC2 shows strong positive (and negative loadings) of 9 of 10 monosaccharides, but also strong negative loadings of one disaccharide (maltose) and 4 of 7 carboxylic acids (Figure S4b).PC3 on the other hand shows strong negative loadings for many metabolites across metabolite groups, except a strong positive loading for one disaccharide (sucrose; Figure S4c).

Experiment 2: Low-elevation plants grown under warm temperatures are more resistant to a specialist herbivore than when grown under cold temperatures
A similarly high proportion of larvae survived when fed on plants grown under cold conditions (99.6% survival rate) or warm conditions (97.8% survival rate), but larvae gained significantly less mass when fed on plant material grown under warm than cold temperatures (ꭓ 2 = 44.3,df = 1, P< 0.0001; R 2 = 0.20; Figure 3a; summary statistics in Table S6).Although the interaction of elevation and temperature was not significant, there were significant effects of elevation and temperature in an additive model explaining 20% of variance in larval mass (Table 3a; parameter estimates in Table S7a).When samples were separated by temperature, there was a significant effect of elevation on larval mass for plants grown under warm temperatures (ꭓ 2 = 4.89, df = 1, P = 0.027; R 2 = 0.09; Figure 3b; Table S7b), but not under cold conditions (ꭓ 2 = 0.61, df = 1, P = 0.691; Figure 3c; Table S7c).Parameter estimates for these models are given in Table S7b and S7c.Notably, the interaction between population and temperature was significant and explained 48% of variation in larval mass (Table 3a), suggesting differences among populations in plant responses to temperature (Figure S5; Table S7d).
Table 3:  Trichome densities on the upper and lower leaf surface were significantly correlated (Pearson correlation coefficient = 0.674, P <0.0001; Figure S6a).Growing temperature did not significantly impact trichome density on either leaf surface (Table 3b; Figure S6b,d), but trichome density showed a significant decline with population elevation (Table 3b; Figure S6c,e).More variance in trichome density was explained by population than elevation (e.g.lower trichomes: elevation R 2 = 0.19; population R 2 = 0.45), reflecting significant variation independent of elevation; for example, trichome densities differed between two highelevation populations (DM trichome densities higher than Aal29; Figure S6c,e; Table S6).
In contrast to trichome density, plants growing at colder temperatures showed significantly higher SLA values indicating reduced leaf density (Table 3b; Figure S7a) and there was a significant population-by-temperature interaction effect (Table 3b; Figure S7b).
This interaction was driven by the intermediate population SFH, which showed a 0.7-fold SLA reduction under cold growing temperatures relative to warm temperatures, whereas most populations showed a 1.3-1.6-foldincrease in SLA under colder temperatures.
Total glucosinolate levels were significantly higher under cold treatment conditions (Table 4a; Figure 4a), with eight individual glucosinolates showing a 1.1 to 9.3-fold increase under cold conditions and four glucosinolates showing a 0.4 to 1.0-fold decrease under cold temperatures (Table S5).The first principal component explained 38.9% variation in glucosinolates, with increasing values of PC1 representing decreasing total glucosinolate content (all 12 glucosinolates showed negative loadings on PC1; Figure S8a; Figure S9a).
Phenological status explained significant variation in both PC1 and PC2 (Table 4a), as a result of a high proportion of flowering plants in population SFH.Plots of PC1 against PC2 or PC3 do not show samples grouping by temperature, but highlight the effect of phenological status (Figure S9b,c).PC2 plotted against PC3 reveals some separation of high and low-elevation samples, though both overlap with samples from intermediate elevation (Figure 4b).
Supporting this, variation in PC3 (13.0% of glucosinolate variation) was explained by significant effects of both temperature and elevation (Table 4a).PERMANOVA also reveals a significant effect of the interaction between elevation and temperature explaining 19.0% of glucosinolate variation, though the interaction between population and temperature explains 1.9x more variation.As with other defence traits, population consistently explains a greater proportion of glucosinolate variation than elevation, suggesting significant population-level effects independent of elevation (Table 4b vs 4c).
Two glucosinolates coeluted at 1.12mins and could only be quantified using their major desulfo-ion.The relative amounts of these compounds varied among and within populations (Table S8).High-elevation populations tended to be dominated by the presence of 4-methylsulfinylbutyl glucosinolate (glucoraphanin), whereas populations 04 and 12 (low

Regression analysis reveals association between larval performance and both morphological and chemical traits
We first fitted a full model containing all morphological and chemical explanatory variables, which showed significant negative effect of trichome density and PC2 (polar) and significant positive effects of PC1 (polar) and PC1 and PC2 (glucosinolates).However, the backwards selection process was impacted by correlated chemical and morphological variables (PC1_polar and SLA; Table S9a), as well as different sample sizes for chemical traits (N=18) and morphological traits (N=16 due to missing population FC).To avoid these issues, we modelled the effects of chemical and morphological traits separately.The chemical trait model revealed significant positive effects of both PC1 (polar), PC3 (polar) and PC2 (glucosinolates) on larval performance (Table 5a, Table S9b).PC1 (polar) explained 47.0% of the variance in mean larval mass (Figure 5a), 1.7x more variance than the next best predictor PC2 (glucosinolates)(Figure 5b).The morphological trait model showed a small significant negative effect of trichome density and significant positive effect of SLA (Table 5b, Table S9c).SLA explained more variation in larval mass in this model than trichome density (34% vs 19%; Table 5b).Notably though the significant chemical traits together explained 2.1x more variance in larval performance than the significant morphological traits.
Table 5: Statistical significance of the fixed effects of (a) chemical traits and (b) leaf morphological traits on variation in larval performance.Separate models for chemical and leaf morphological variables were fitted as explained in the main text.The F-statistics, degrees of freedom (df) and P-values indicate whether a factor had a significant effect following removal from the model (a backward selection process).The adjusted R 2 value for the simplified (minimal adequate) model is given, as well as the R 2 estimates for each factor in the minimal adequate model.The terms are given in sequential order they were removed from the model (i.e., the top variable is first one removed).The first three principal components for the polar metabolome ('_polar') and glucosinolates ('_gsls') are used.Note that morphological measurements from one population are missing (see sample sizes, N).
Parameter estimates for these full models are given in Table S9b and S9c

DISCUSSION
Our results show reduced performance of a specialist herbivore feeding on an alpine plant grown under warm than cold temperatures, though this difference was greater for plants originating from low-elevations than high-elevations.The reduced performance of larvae on plants grown under warmer conditions was correlated with several morphological and chemical defence traits, with changes in polar metabolites explained most variance in larval performance.The elevational trend in herbivore performance on leaf tissue grown under warm not cold temperatures is consistent with low-elevation populations showing adaptation to abiotic and/or biotic factors that increase their resistance to the higher levels of herbivory they naturally encounter.In addition to providing further evidence that plant responses to temperature variation can alter subsequent interactions with invertebrate herbivores, these findings suggest that trait plasticity under different temperatures may need to evolve in highelevation plant populations to reduce their vulnerability to biotic stressors in a warming climate.
Building on previous work demonstrating significant indirect effects of temperature on plant resistance and herbivore performance (Bauerfeind & Fischer, 2013;Hamann et al., 2021;Nechols et al., 2020;Kharouba & Yang, 2021), as well as effects of experimental warming on interactions of high-elevation plant communities with herbivores (Birkemoe et al., 2016;Descombes, Kergunteuil, et al., 2020;Descombes, Pitteloud, et al., 2020), the current findings demonstrate that the response of plant populations to warming can vary across broad environmental gradients, with implications for their interactions with herbivores.
Growing A. alpina under warm conditions resulted in reduced performance of herbivores, consistent with the warming-induced resistance seen in other low-elevation Brassicaceae (e.g.Bauerfeind & Fischer, 2013;and Rorippa dubia;Pham & Hwang, 2020), yet the reduction in herbivore performance was greater on low-elevation populations compared to that on high-elevation populations.This result is consistent with a previous experiment using whole, intact plants and similar warm growing conditions (Buckley et al. 2019).Relatively few studies have tested for variation in larval performance on multiple populations from different elevations (Buckley, Pashalidou, et al., 2019;Rasmann, Pellissier, et al., 2014), and even fewer have examined how the elevational origin of populations might influence plant resistance to herbivores under different growing temperatures (Dostálek et al., 2020;Pellissier et al., 2014).Although in the present study we do not directly measure herbivore resistance (amount of leaf material eaten), our results on larval performance together with plant morphological and chemical trait data suggest plant populations vary in their resistance to herbivores or in their quality as food for herbivores.Repeating our larval performance assays by exposing larvae to different temperatures (similar to those under which the plants were grown) would further help disentangle the direct and indirect effects of temperature on herbivore performance (e.g.Bauerfeind & Fischer, 2013).Nevertheless, our results provide novel insight into the extent of variation among populations in the indirect effects of temperature change on plant-herbivore interactions, and also suggest that highelevation populations might be more vulnerable to the combined effects of warming and the upward spread of herbivores predicted under future climate change.

Sinapsis arvensis;
Understanding the traits responsible for changes in herbivore performance is challenging given that experimental warming induces change in plant traits both directly and indirectly linked to herbivore resistance (Pincebourde et al., 2017).Indeed, in Arabis alpina we observed effects on nearly all measured traits, including several phenological traits; A. alpina grew faster and flowered more rapidly under warm than cold conditions (Appendix 1).
Comparing temperature-driven changes in relevant traits to the performance of the specialist herbivore Pieris brassicae in feeding assays revealed multiple significant morphological and chemical traits, but the first principal component of polar metabolite variation was the best predictor of larval performance.Total glucosinolates (reflected by PC1 of glucosinolate variation) was not associated with herbivore performance, though another glucosinolate PC axes was a significant predictor.The lesser effect of defensive compounds relative to the polar metabolites is perhaps expected given that Pieris brassicae is adapted for feeding on glucosinolate-containing plants (Jeschke et al., 2016).Rather, our data suggests that leaf nutritional content is an important mediator of herbivore performance, a result supported by several other experimental warming studies (Jamieson et al., 2015;Pham & Hwang, 2020).
Given the broad range of metabolites measured in our study, we were unable to isolate the specific metabolites driving the changes in herbivore performance (plants grown under warm temperatures showed higher levels of 22 polar metabolites and reduced levels of another 14 polar metabolites).Furthermore, the polar metabolites were not measured on the same plant material as that fed to the caterpillars, limiting our ability to detect the role of individual metabolites.However, the changes we observed in non-defensive metabolites were consistent with physiological responses to cold exposure in other studies (Fürtauer et al., 2019); for example, the accumulation of disaccharides, such as sucrose and fructose, in leaves of plants exposed to 4/5C chilling conditions has been observed in A. alpina accessions from different elevations (Wingler et al., 2015), and in the Brassicaceae relatives Arabidopsis lyrata (Davey et al., 2009) and Arabidopsis thaliana (Cook et al., 2004).Sucrose has been identified as a likely phagostimulant in P. brassicae (Schoonhoven & van Loon, 2002), which might explain increased larval growth rates on cold treatment plants.The specific changes impacting herbivore performance might therefore involve such sugars or correlated shifts in other groups of metabolites (e.g.N-rich compounds such as amino acids; Lemoine et al. 2013;Jamieson et al. 2015).Although not directly measured in our study, the observed changes in primary metabolites might also reflect shifts in C:N ratios, with reduced nitrogen levels associated with reduced larval performance in other experimental warming studies (e.g.Pham & Hwang, 2020;Jamieson et al. 2015;Couture et al. 2015).Further experimental work is needed to move beyond these correlative associations and identify the specific metabolic drivers.Nevertheless, it seems likely that these consistent physiological responses to warming alters nutrient availability (and thereby food plant quality) such that herbivores show reduced performance on this alpine plant.
As our plants were grown in a common environment, our findings further suggest that the observed variation in leaf physiology and associated impacts on herbivores among populations has a genetic basis that could support adaptive evolution in response to varying herbivore pressures.As the extent of leaf damage by herbivores on field populations of A. alpina also declines with increasing elevation (Buckley, Widmer, et al., 2019), our results are consistent with adaptation of low-elevation populations to more intense herbivory, a pattern seen in other systems (e.g.Anderson et al., 2015).However, the relatively weak associations between herbivore performance and known herbivore resistance traits (e.g., glucosinolates), suggests that responses to herbivores may reflect adaptation to abiotic factors that vary with elevation, which in turn impact the quality of plant tissue as food for herbivores (Moreira et al., 2018).Leaf morphological traits, such as trichomes, can be under selection from herbivores and abiotic stressors (e.g., aridity; Kessler et al. 2007).Additional data is therefore needed to conclusively determine whether herbivory is indeed a key selection pressure shaping the observed variation in plant traits.In particular, our study employed a model specialist herbivore, rather than one known to inflict significant damage on natural populations of A. alpina (for example Athalia rosae or Plutella xylostella) or relevant generalist herbivores (e.g.snails; JB personal observation).Future studies could assess the observed vulnerability of high-elevation A. alpina populations to different herbivores under warmer temperatures over the full developmental period of both specialist and generalist herbivores to confirm the ecological relevance of this pattern.
The observed elevational variation in herbivore performance under warm but not cold temperatures is consistent with reduced plasticity of traits that alter herbivore resistance in high-elevation plants.It is difficult to identify specific traits driving this pattern based on the current data, but it seems likely that changes in physiological plasticity are playing a significant role.By contrast, trichome density, a structural defence against herbivores, showed limited plasticity with respect to temperature, with consistent differences among populations under warm and cold conditions.Greater physiological plasticity in low-elevation populations might reflect the more variable abiotic and biotic conditions that this alpine plant experiences (e.g., as a result of limited snow protection over winter).An increasing number of studies are identifying divergent plastic responses between high-and low-elevation species and populations (de Villemereuil et al., 2018;Pellissier et al., 2016;Schmid et al., 2017), highlighting the potential importance of adaptive variation in plasticity for understanding the responses of populations from different elevations to climate change (Midolo & Wellstein, 2020).Furthermore, our finding of elevational divergence in plant responses to temperature variation suggests that without genetic adaptation plants at high elevations are more vulnerable to herbivores under warm temperature than their low elevation counterparts.
Colonisation of high-elevation sites by low-elevation plants is possible in the long-term, but unlikely to occur at rates sufficient to outpace the upward spread of herbivores.Rather, the significant variation we find in temperature-mediated trait plasticity across populations independent of elevations might facilitate future evolutionary responses to changing abiotic and biotic factors at high-elevations.Future research into the responses of alpine species to future climate warming should therefore account for both the indirect effects of temperature and population-level divergence in thermal plasticity when trying to predict the outcomes of changing plant-herbivore interactions.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article.
Appendix S1-additional materials and methods for chemical analysis and collecting phenological data from experiment 1, and supporting tables and figures:  Appendix S3data collected from experiment 1 on variation in 52 polar metabolites.
(1a) PC ~ temperature * elevation + extraction set + phenological status (1c) PC~ temperature * population + extraction set + phenological status To assess how the 52 metabolites contribute to the observed multivariate patterns, we extracted the individual metabolite loadings for the three principal components and plotted these by the KEGG-brite classifications to look for significant changes in key metabolite groups.Experiment 2: Mean upper and lower leaf disc trichome density (based on two discs) and mean SLA (based on up to four leaf discs) was calculated per plant and used as the response variable.Due to an error in sample labelling, one of the high-elevation populations FC was excluded from this analysis.Leaf tissue sampling set (A or B) and phenological status (vegetative vs reproductive) was included to account for the small number of plants entering the reproductive phase (bolting/flowering) during this experiment (8 out of 112).We fitted the following linear models separately for mean trichome density and mean SLA: (2a) trait mean ~ temperature * elevation + phenological status + sampling set.(2b) trait mean ~ temperature * population + phenological status + sampling set.Variation in individual larval mass was then modelled using a linear mixed effects model including a random effect of microbox (three replicate containers for each population-bytemperature combination): (2c) Larval mass ~ temperature * elevation + (1|microbox) (2d) Larval mass ~ temperature * population + (1|microbox) To explicitly test whether larvae perform worse on low-than high-elevation plants, we additionally tested the effect of elevation (in metres a.sl.) for the warm and cold treatment data separately: (2e) Larval mass ~ elevation + (1|microbox) Variation in both total glucosinolates and the first three principal components of glucosinolate variation were modelled, including the covariates of leaf tissue sampling set, glucosinolate extraction set (set 1, 2 or 3) and phenological status.Similar to analyses of the polar Accepted Article Accepted Article metabolome, we used PERMANOVA to test for the marginal effects of each factor on glucosinolate variation in the following models: (2f) glucosinolate ~ temperature * elevation + extraction set + sampling set + phenological status (2g) glucosinolate ~ temperature * population + extraction set + sampling set + phenological status of fixed factor effects on the first three principal components of polar metabolome variation (PC1-PC3) following simplification of the following full models.The upper five rows represent results for the model including elevation: PC ~ temperature * elevation + extraction set; and the lower five rows give results for the same model but with population replacing elevation as a fixed factor in the model.The degrees of freedom for the likelihood ratio test are given in the first column next to the corresponding factor, and F-statistics and P-values given under the column for each principal component axis of variation.The variance explained by the minimal adequate model after simplification is given by the adjusted R 2 Statistical significance of the fixed effects of elevation (or population), temperature and relevant covariates on (a) larval performance, and (b) morphological traits (trichome density and SLA) following model simplification.For both (a) and (b) separate models were used to test the interaction between elevation and temperature (upper rows), and population and temperature (lower rows).For modelling variation in trichome density and SLA we included covariates of sampling set (leaves were sampled in two sets) and plant phenological status (vegetative or reproductive).Transformations were applied to response variables where necessary to improve model fit.For all analyses, the degrees of freedom for the likelihood ratio test are given next to the corresponding factor, and F-statistics and P-values given under 'significance of effect'.The adjusted R 2 explained by the significant fixed factors in the minimal adequate model (highlighted in bold) is given.
elevation) and SFH (intermediate-elevation) were dominated by 2(S)-hydroxy-3-butenyl (progoitrin); other populations varied in the dominant glucosinolates at this retention time.Table 4: (a) Upper five rows: statistical significance of fixed effects of elevation, temperature and relevant covariates on total glucosinolates and the first three PCs of glucosinolate variation after model simplification; Lower five rows: same starting model but with population replacing elevation.The degrees of freedom for the likelihood ratio test are given next to the corresponding factor, and F-statistics and P-values given under the relevant column.The variance (adjusted R 2 ) explained by the fixed factors in the minimal adequate model after simplification is given below the results of each model, with included factors highlighted in bold.(b, c) PERMANOVA output testing the marginal effects of factors on glucosinolate variation using separate models for (b) an elevation * temperature interaction, or (c) population * temperature interaction.The Sum of squares ('SS') and variance explained by each term are also provided in this output.As the interaction between elevation and temperature was significant, PERMANOVA does not provide the main effects of elevation and temperature.

Figure 1 :
Figure 1: Populations sampled for the experiments and the experimental design used.(a) A digital elevation map of Switzerland marked with the nine sampled populations (coloured by broad elevation class).(b) Overview of experimental design used for the three experiments: The nine sampled populations were grown under warm and cold temperature conditions based on mean summer temperatures at approximately 1000m and 2000m elevation.Further details about experimental design are given in the methods and Figure S1.

Figure 2 :
Figure 2: (a) Principal Components Analysis reveals multivariate divergence in the polar metabolome with respect to temperature.(b) The mean abundance in nanograms per microlitre of 52 polar metabolites grouped into KEGG Brite metabolite classes (see inset key) under the two temperature treatments.

Figure 3 :
Figure3: Variation in larval mass after feeding on leaf material from nine different populations grown under warm and cold growing conditions.(a) Violin plot showing larvae gained a significantly greater mass feeding on plant material from cold rather than warm conditions (P<0.0001;R 2 = 0.20).(b,c) Larval mass significantly increases on plants from higher elevations under warm (red) conditions (P= 0.027; R 2 = 0.09), but under cold (blue) conditions larval mass shows no significant trend with elevation (P = 0.691).The relevant statistics are detailed in the main text (experiment 2 results) and in TableS7b and Table S7c.

Figure 4 :
Figure 4: (a) Variation in total glucosinolate production across populations grown under different temperatures.Glucosinolate amounts are expressed in micromoles per gram dry tissue weight (molg -1 DW), with amounts estimated relative to the internal standard sinigrin.(b) Multivariate glucosinolate variation: PC2 plotted against PC3 showing some clustering of samples by elevation (samples coloured by broad elevation classes; open symbols indicate the plant was flowering, filled indicate vegetative).The variance explained and statistical significance of the temperature effect on total glucosinolates and temperature/elevation effect on PC3 is given in Table4a.

Figure 5 :
Figure 5: Plots of two significant chemical predictors of variation in larval performance.(a) Regression of mean PC1 (polar metabolites) per population per treatment on mean larval mass; (b) Regression of mean PC2 (glucosinolates) per population per treatment on mean larval mass.The regression lines illustrate the positive effect of both factors on larval mass.

Figure S1 :
Figure S1: Overview of experimental design for experiments Figure S2: Overview of phenological data from initial experiment Figure S3: NMDS plot of polar metabolite variation Figure S4: PCA loadings for PCs1-3 of polar metabolite variation Figure S5: Boxplot of larval mass variation across populations and temperatures Figure S6: Plots summarising variation in number of trichomes Figure S7: Plots summarising variation in Specific Leaf Area (SLA) Figure S8: PCA loadings for PCs1-3 of glucosinolate variation Figure S9: Additional plots of principal component variation for glucosinolate data

Table 1 :
Marginal effects of explanatory factors in a PERMANOVA analysis of the polar metabolome dataset.(a) Main effects of elevation, temperature and covariates in the model: temperature + elevation + extraction set + phenological status; (b) Main effects of an identical model but with population replacing elevation.These statistics reflect the effect of a factor in a model containing all other explanatory variables.Interaction effects were not significant as described in main text (tested in separate models).The degrees of freedom (df), sum of squares (SS), variance explained (R 2 ) and relevant F-statistics/P-values are given.

Table S1 :
Description of populations used in the experiments.

Table S2 :
Summary of phenological data from initial experiment

Table S3 :
Mean difference in 52 polar metabolites between growing temperatures

Table S4 :
Sample sizes for traits measured in experiment 2

Table S5 :
Mean amounts of 12 glucosinolates under different growing temperatures

Table S6 :
Summary statistics for traits measured in experiment 2

Table S7 :
Parameter estimates for statistical models (larval performance data)

Table S8 :
Relative amounts of two co-eluting glucosinates

Table S9 :
Parameter estimates for statistical models (multiple regression analysis)