Alternative developmental and transcriptomic responses to host plant water limitation in a butterfly metapopulation

Predicting how climate change affects biotic interactions and their evolution poses a challenge. Plant-insect herbivore interactions are particularly sensitive to climate change, as climate-induced changes in plant quality cascade into the performance of insect herbivores. Whereas the immediate survival of herbivore individuals depends on plastic responses to climate change induced nutritional stress, long-term population persistence via evolutionary adaptation requires genetic variation for these responses. In order to assess the prospects for population persistence under climate change, it is therefore crucial to characterise response mechanisms to climate change induced stressors, and quantify their variability in natural populations. Here, we test developmental and transcriptomic responses to water limitation induced host plant quality change in a Glanville fritillary butterfly (Melitaea cinxia) metapopulation. We combine nuclear magnetic resonance spectroscopy on the plant metabolome, larval developmental assays and an RNA seq analysis of the larval transcriptome. We observed that responses to feeding on water limited plants, in which amino acids and aromatic compounds are enriched, showed marked intrapopulation variation, with individuals of some families performing better on control and others on water limited plants. The transcriptomic responses were concordant with the developmental responses: Families exhibiting opposite developmental responses also produced opposite transcriptomic responses, e.g. in growth associated intracellular signalling. The opposite developmental and transcriptomic responses are associated with between families differences in organic compound catabolism and storage protein production. The results reveal heritable intrapopulation variability in plasticity, suggesting potential for evolutionary responses to drought-induced changes in host plant quality in the Finnish M. cinxia metapopulation.

To examine the mechanisms via which host plant water stress cascades into larval performance, 88 we combined host plant metabolic profiling with development assays and full-transcriptome 89 sequencing of herbivore larvae. First, we profiled metabolic differences between well-watered 90 and water-limited host plants using proton nuclear magnetic resonance spectroscopy ( 1 H- The ribwort plantain 99 The ribwort plantain (Plantago lanceolata, Plantaginaceae) is the natural host of M. cinxia in 100 its studied range. It produces anti-herbivore and anti-fungal chemicals (iridoid glycosides and 101 phenolic compounds) of which the amounts can vary with plant genotype, age, and 102 environmental conditions (Bowers et al., 1992;Bowers & Stamp, 1993). 103 We collected seeds from a natural population in the Åland Islands (60.196° Lat., 20.704° Lon.) 104 and after germination planted 360 plants in 0.75 litre pots (two saplings each). We reared the 105 plants for three months in controlled greenhouse conditions (ca. 40 ml water / pot daily, 15L:9D 106 photoperiod with 26:18 °C temperature cycle) before initiating the water limitation treatment. 107 We exposed 240 plants to a water limitation treatment in which daily watering was reduced by 108 50% compared to controls (20ml per pot). This watering scenario was developed in a pilot study 109 in which we experimented with minimum watering allowing the plants to stay alive. More 110 plants were allocated to the water limitation treatment than to the control treatment, because 111 plants produced substantially less leaf biomass in the water limitation treatment ( Figure S1). In 112 order to minimize temporal trends in plant quality caused by the plants acclimatising to altered 113 7 water availability, we initiated the water limitation treatment well in advance (47 days) to the 114 larval exposure (see below).

115
Leaf metabolomics assays 116 Each morning prior to watering (9-10 am), we randomly harvested P. lanceolata leaves from 117 control and water limited plants and cut them into 2.25cm 2 pieces, discarding the basal and tip 118 parts of the leaves. We used these pieces to feed larvae during the experiment (see below), and 119 selected a random subset of six pieces from both treatments for metabolomics assays. For the 120 assays, we recorded the fresh biomass of each piece of leaf, snap froze the pools in liquid 121 nitrogen and stored them in -80 °C. We then freeze dried the samples for 48 hours after which    Table S1). We allowed the larvae to continue diapause in a 147 climate chamber (+5 °C, 95% air humidity) for five months and after breaking diapause reared 148 them to adults (12L:12D photoperiod, 28:15 °C temperature cycle) with daily ad libitum 149 provision of control reared P. lanceolata leaves. To maintain original population genetic 150 structure while minimizing risk of inbreeding in the larval groups, we mated virgin females 151 with males derived from the same habitat patch but different overwintering nests. We placed 152 the mated females on host plants surrounded by mist nets immediately after they had finished 153 mating and collected egg clutches laid on the plant until the female was found dead. 154 We monitored hatching of egg clutches daily and -in order to have two replicates per full-sib 155 family per treatment -picked two larval groups of a minimum of eighty larvae each from each 156 female to enter the experiment. To minimize potential quality differences caused by clutch rank (Rosa & Saastamoinen, 2017), we selected the larval groups from within the first three larval 158 groups.

159
Treatments & developmental assays 160 On the day after hatching, we divided the larval groups into four smaller groups of twenty larvae 161 each and placed them on separate petri dishes (9 cm diameter, 1.5 cm deep) lined with filter 162 paper. We then randomly assigned each of the dishes to one of four different treatments

173
During the experiment, we recorded development time, body mass at diapause and mortality 174 during development. Once the last larva on a petri dish had entered diapause, we allowed them 175 to spend another four days in normal rearing temperature and photoperiod, after which we 176 measured their body mass and placed them in climate chambers (+5 °C, 95% air humidity) for 177 diapause. We allowed the larvae to diapause for six months, after which we woke them up and 178 recorded overwintering mortality.

179
With one exception, unexplained mortality during the rearing (i.e. mortality that could not be 180 explained by accidents during handling) was low in all families and in all treatments (mean = 181 1.1 larvae / petri dish, sd=1.4 larvae). Only the control treatment of the first replicate larval 182 group of family F-5 had a mortality of 55% (11 larvae). As the larvae in this group were 183 developing poorly in general, we concluded it to be an outlier case, potentially suffering from 184 a disease or some other unknown agent, and decided to exclude this larval group from any 185 further analyses.

186
Transcriptomics sampling and sequencing 187 When more than fifty percent of the larvae on a petri dish had spent two full days in the 4 th 188 larval instar, we sampled either ten (from the first larval group of each full-sib family) or five 189 (from the second larval group of each full-sib family) larvae for RNA and DNA extraction (see 190   Table S2 for exceptions). Before noon on the day of sampling, we provided the larval groups 191 with leaf tissue matching their treatment and monitored that they fed on the plant before 192 sampling to ensure sampling larvae that are feeding. We weighed and sampled the larvae for 193 RNA and DNA extraction approximately 1.5-2h after feeding, by immersing them in liquid 194 nitrogen and stored the sampled larvae in -80 °C until further processing.

195
To extract RNA and DNA, we placed the larvae in dry ice and homogenized the frozen larvae within the larval family as a group-level intercept. Like with the above described GLM for 247 water content, we did this in Stan via the R packages brms and Rstan.

248
Because the majority of larvae in three of the families (F-7, F-8 and F-9; Figure S4, Table S1)

300
Host plants shift metabolome but not water content upon water limitation treatment 301 We observed that the leaf tissue water content did not differ between the control and water 302 limited plants (Table 1) suggesting that the plants in the two treatments became more similar as the experiment 319 proceeded ( Figure S5, Table S3).  The first three of these axes illustrate the overall transcriptomic differences between families 360 and the fourth highlights both the responses to host water limitation and the response differences 361 between the families (i.e. a family-by-treatment interaction; Figure 3a, Figure S7). An 362 additional RDA model maintaining all temporal water limitation levels produced practically 363 identical results (Table S6).

364
The first RDA axis separates families in which performance was improved in the early  (Table S9). Transcripts belonging to the latter group include those that translate 385 into peripherin-2 like protein, lachesin and gustatory receptor 22 (Table S9) (Table S10). gives further support to the interpretation that the responses are mostly related to larval growth 408 (Table S11). Some transcripts point towards the production or functioning of specific tissues 409 such as the cuticle (i.e. larval cuticle protein F1-like) or the nervous system (i.e. neurexin-1 410 alpha) but most transcripts are related to cellular signalling associated with cell proliferation 411 and growth (Table S11) (Table S11). All the above mentioned transcripts are negatively 417 loaded on RDA4, whereas the opposite is true for transcripts coding for moderately methionine-418 rich storage protein a, which was also associated with RDA1 (see above). Thus, in family F-2,  Table S12).

449
In families F-1 and F-5, the processes were not as clearly associated with endocytosis.  Intrapopulation variability in performance on water stressed host plants 475 We found that whereas the larvae of four M. cinxia families grew larger, developed faster and/or 476 had lower overwintering mortality when feeding on water stressed host plants in which amino 477 acids and phenolic compounds were enriched, the opposite was true for three of the families, 478 and mixed responses were observed in two families (Figures 1 and 2). The observation of  (Figures 2 and 3a, Table S4).

529
The protein products associated with the between families differences (i.e. those most strongly 530 loaded on the first RDA axis) revealed differences regarding metabolising and storing nutrients.

531
Individuals of the family F-2, in which performance was better on control treated plants, had 532 higher baseline expression of transcripts involved in the production of storage proteins 533 (moderately methionine-rich protein a and arylphorin) and fat accumulation (transmembrane 534 protein 135) (Table S9). In holometabolous insects, these proteins are typically associated with 535 preparing for metamorphosis, diapause, egg production, or periods of environmental/nutritional

573
Although we know that regional population growth rates are linked with precipitation across

578
The observed between families variability in responses to host plant water limitation is a 579 promising sign of heritable genetic variability in drought tolerance within the system. However, 580 since the parental generation of the current study was collected from the wild, we cannot