Herbivores disrupt clinal variation in plant responses to water limitation

Plasticity in plant traits, including secondary metabolites, is critical to plant survival and competitiveness under stressful conditions. The ability of a plant to respond effectively to combined stressors can be impacted by crosstalk in biochemical pathways, resource availability, and evolutionary history, but such responses remain underexplored. In particular, we know little about intraspecific variation in response to combined stressors or whether such variation is associated with the stress history of a given population.

Here, we investigated the consequences of combined water and herbivory stress for plant traits, including relative growth rate, leaf morphology, and various measures of phytochemistry, using a common garden of Asclepias fascicularis milkweeds. To examine how plant trait means and plasticities depend on the history of environmental stress, seeds for the experiment were collected from across a gradient of aridity in the Great Basin, USA. We then conducted a factorial experiment crossing water limitation with herbivory.

Plants responded to water limitation alone by increasing the evenness of UV-absorbent secondary metabolites, and to herbivory alone by increasing the richness of metabolites. However, plants that experienced combined water and herbivory stress exhibited similar phytochemical diversity to well-watered control plants. This lack of plasticity in phytochemical diversity in plants experiencing combined stressors was associated with a reduction in relative growth rates.

Leaf chemistry means and plasticities exhibited clinal variation corresponding to seed-source water deficits. The total concentration of UV-absorbent metabolites decreased with increasing water availability among seed sources, driven by higher concentrations of flavonol glycosides, which are hypothesized to act as antioxidants, among plants from drier sites. Plants sourced from drier sites exhibited higher plasticity in flavonol glycoside concentrations in response to water limitation, which increased phytochemical evenness, but simultaneous herbivory dampened plant responses to water limitation irrespective of seed source.

Synthesis: These results suggest that climatic history can affect intraspecific phytochemical plasticity, which may confer tolerance to water limitation, but that co-occurring herbivory disrupts such patterns. Global change is increasing the frequency and intensity of stress combinations, such that understanding intraspecific responses to combined stressors is critical for predicting the persistence of plant populations.

Methods
Study system

Asclepias fascicularis (narrowleaf milkweed) is one of the most widely distributed milkweed species in the western United States and an important food-plant species for the monarch butterfly (Danaus plexippus), among other specialist herbivores (Woodson 1954, Dilts et al. 2019). Narrowleaf milkweed is found across a wide range of water availabilities, including in very dry locations (down to at least 100 mm of annual precipitation) (Woodson 1954). Milkweeds are well known for their chemical defenses, which are sequestered by monarch butterflies. Unlike many milkweeds, A. fascicularis contains few cardiac glycosides, but the leaves contain high concentrations of flavonols and both leaves and roots contain diverse pregnane glycosides (Rasmann and Agrawal 2011, Mundim and Pringle 2020, Diethelm et al. 2022). Flavonols are hypothesized to play a role in the mitigation of plant water stress by functioning as antioxidants, and the biological role of the pregnane glycosides remains unknown (Kaminska-Rozek and Pukacki 2004, Zehnder and Hunter 2007, Diethelm et al. 2022).

            To explore the role of climatic history in plant trait responses to combined stressors, we collected seeds from five sites spanning a west-to-east transect of 385 km of the Great Basin Desert, USA, in fall 2016 (Fig. S1). Annual precipitation at these five sites encompasses the 5th–95th percentiles of annual precipitation among 453 occurrence records for A. fascicularis in the western United States (Dilts et al. 2019). To estimate the typical drought stress at each of these sites, we calculated the cumulative annual climatic water deficit (CWD) using monthly precipitation and temperature data from the PRISM Data Explorer tool (PRISM Climate Group 2019) for the years 2004-2016 (for more details, see Diethelm et al. 2022). The seed-source sites with high CWD (hereafter, dry) experience more water limitation on an annual basis than the sites with low CWD (hereafter, wet). From dry to wet, these sites were: Fallon, NV (FN; 989.0 ± 19.7 mm), Georgetown, CA (CA; 985.7 ± 39.4 mm), Battle Mountain, NV (BM; 847.7 ± 29.5 mm), Reno, NV (RN; 587.2 ± 20.0), and Verdi, NV (VE; 460.3 ± 16.7) (Fig. S1 and Table S1). In a previous study conducted in a glasshouse, plants from drier source populations had higher constitutive concentrations of leaf flavonols, whereas plants from wetter source populations exhibited higher induction of leaf flavonols under acute water stress (Diethelm et al. 2022).

 Water limitation ´ herbivory experiment 

To investigate how exposure to combined stressors affects plant traits, we conducted a factorial common-garden experiment manipulating water and herbivory. This experiment was conducted from May–October 2017 in a 30x10-m outdoor plot at the University of Nevada, Reno, Valley Road Agricultural Station. Asclepias fascicularis seeds from each of the five seed-source sites, representing four maternal families per site and 16 plants per family (N = 320), were germinated in May 2017 and transferred to 0.5-L pots of Full Circle® Soar potting soil to grow for two months. In July 2017, the 236 plants that had successfully germinated were transplanted into the plot in a randomized block design (Fig. S2). Some of these plants failed to establish in the plot, reducing the total number of plants to 147 (see Table 1 for n per treatment). All plants were fertilized once per month with 1:1:1 Triple Pro® NPK fertilizer pellets.

            To evaluate plant responses to single versus combined stressors, plants were randomly assigned to one of four possible treatments: well-watered control (well watered + no herbivory); dry (dry + no herbivory); herbivory (well watered + herbivory); or dry + herbivory. Well watered and dry irrigation lines alternated from east to west through the plot. To enhance establishment, all plants were watered daily to soil saturation using drip-irrigation for the first 8 weeks. For 38 days beginning 6 September 2017, well watered plants received 84 mL of water to an ~157 cm2 soil area per day while the dry lines received no irrigation. The well watered treatment was designed to simulate precipitation at the wettest seed source site (VE), whereas the drought treatment was designed to bring plants to their wilting point. After one month, plants in the dry treatment were beginning to show severe wilt, so we supplemented them with 84 mL over two days in early October.  

            Herbivory treatments were implemented over 3 weeks, starting in mid-September. Leaves were mechanically damaged to produce ~15% tissue loss, using a rolling leather punch to simulate herbivory (Baldwin 1990). To test if this mechanical damage effectively simulated damage by herbivores, a randomly selected subset of plants (n = 30) received ~5% leaf removal by monarch (Danaus plexippus) caterpillars in addition to ~10% mechanical damage. Caterpillars were obtained from a laboratory colony at the University of California, Davis, and were maintained on A. fascicularis plants in a glasshouse. In the common garden, larvae were restrained to ~5% of aboveground plant tissue using netted bags. We checked plants periodically and removed caterpillars once all of the leaves in the netted section had been consumed.

Plant responses to the treatments were quantified as follows. To calculate changes in aboveground relative growth rate (RGR), we measured: (i) plant height (cm) and (ii) the number of  >4-cm-long stems branching directly from the main plant axis at the beginning and end of treatments. We then multiplied (i) by (ii) to estimate plant size. Relative growth rate was calculated per day using the natural logarithms of the change in size, divided by the 38 days of the experiment. At harvest, six leaves from the second whorl of the main axis of each plant were separated to assess the following traits: leaf width, leaf mass per area (LMA), foliar water content, trichome density, and secondary metabolites. We measured leaf width (mm) at the widest point of three leaves and averaged those measurements per plant. To determine LMA (mg dry mass /mm2), we cut a rectangle of 1-cm length of the 4th leaf, measured the width, and dried it at 60° C for 48 h before reweighing. Foliar water content was calculated as the difference between the wet and dry masses divided by the area of the same rectangle. To determine trichome density, we counted the number of trichomes on a randomly selected 2-mm wide leaf edge. To assess plant water-use efficiency, we analyzed 13-carbon isotopes (Farquhar et al. 1989) from homogenously ground leaf and stem material at the UC Davis Stable Isotope Facility. These same samples were also analyzed for percent carbon (C) and percent nitrogen (N), which we used to calculate the C:N ratio.

            To investigate how our treatments affected plant secondary metabolites, we performed an untargeted analysis of the UV-absorbent metabolites in the remaining two leaves. Leaves were stored in a –80 °C freezer until they were freeze-dried, ground, and extracted in methanol with a digitoxin internal standard. Extracts were then run on an UPLC-UV system; for detailed methods, see Appendix S1. UV absorbance spectra were recorded from 200 to 330 nm, and peak areas were quantified at 219 nm. Compounds were differentiated based on retention time and the UV spectrum (Appendix S2) and concentrations were estimated as digitoxin equivalents. We also quantified the richness (S), evenness (J), and the Shannon diversity index (H) of all of the UV-absorbent peaks for each sample with a minimum area of 8,000 absorbance units (AU) and a minimum height of 5,000 AU. The chemical identities of key peaks were verified by LC-Q-ToF-MS at the Max Planck Institute for Chemical Ecology (Appendix S1,S3).

Statistical analysis

To analyze plant responses to experimental treatments, we used generalized linear mixed models  (GLMMs) from the glmmTMB package (Brooks et al. 2017) in R version 3.6.1 (R Core Team 2022). To compare effect sizes, we z-transformed all continuous variables before analysis using the BBmisc package (Bischl et al. 2017), and we report beta coefficients (β) from the models with standard errors. We assessed the residuals of each fitted model, and we square-root or log-transformed response variables when these transformations provided a better fit to the gaussian distribution as needed. All plant-trait models used plant maternal family nested within the seed-source site as random intercept effects to account for non-independence within families or sites.

To determine the predictors of plant traits, we used backward selection (Zuur et al. 2009) in the MuMIn package (Barton 2019); see Tables S2-S3 for model selection details. To determine whether interactive effects between stressors depended on seed source, each saturated model began with a three-way interaction among water treatment, herbivory treatment, and seed-source CWD. Neither the concentration nor the richness of UV-absorbent plant metabolites was changed by the addition of monarch herbivory to mechanical damage (respectively, Welch's t = -0.12, df = 56.1, P = 0.9; Welch's t = 0.53, df = 65.8, P = 0.6; Fig. S3), so all herbivory treatments were pooled in subsequent analyses.

Models were selected based on the lowest sample-size corrected Akaike Information Criterion (AICc), and any marginally significant predictors were evaluated using log-likelihood ratio tests (LRT) between models in the lmtest package (Zeileis and Hothorn 2002). Marginal and conditional R2 values were calculated for selected models in the MuMIn package. Data are archived in Dryad (Diethelm et al. 2023).