Effects Of Acclimation, Population, And Sex On Behavioral Thermoregulation, CTMax, Symptoms Of Heat Stress, And Gene Expression of Melanoplus Differentialis, A Generalist Grasshopper - Does Temporal Thermal Heterogeneity Prepare Populations For A Warming World?

Insects thermoregulate using both canalized and plastic mechanisms. Populations of insects utilize these mechanisms to different extents, and while it is posited that the degree of thermal uctuation a population experiences can determine the optimal combination of mechanisms to utilize, this is still being elucidated. We used three populations of the generalist grasshopper, Melanoplus differentialis (Thomas, 1856), from sites experiencing different degrees of thermal heterogeneity to test for correlations between thermal heterogeneity and 1) behavioral thermoregulation, 2) upper temperature tolerance, 3) the ability to thermally acclimate, and 4) gene expression. We found that 1) behavioral thermoregulation did not differ among sites, 2) CT Max of males, but not females, was higher at more thermally heterogeneous sites, 3) there was acclimation in some of the tested traits, but thermally heterogeneous sites did not always have the most plastic individuals, and 4) there were differences in gene expression among sites, but these differences were not between the most and least thermally heterogeneous sites. We concluded that thermal heterogeneity may play a selective role in some, but not all, of the measured thermoregulatory traits and their plasticity.


Introduction
All insects experience uctuating temperatures. While uctuation within a range of temperatures around a thermal optimum is tolerated and can be bene cial (Adamo and Lovett 2011), exceeding this range in either direction can also result in detrimental effects (Colinet et al. 2015). One of the mechanisms insects have evolved to enhance thermal tolerance is canalized physiological thermoregulation. This has been demonstrated in the lab, where arti cial selection can produce individuals with increased cold tolerance (Condon et al. 2015), and in the eld, where extreme thermal specialists can utilize canalized thermoregulation to withstand temperatures as high as 62.2 (Bernstein 1979) and as low as -24° C (Slabber et al. 2007). The ability to thermoregulate in this manner can vary among species and populations (Slatyer et al. 2016). While a canalized ability to withstand hot or cold temperatures can alleviate some of the negative effects of temperature stress on populations, evidence indicates that a canalized response alone is not enough to alleviate effects of hazardous temperatures. For example, Mitchell & Hoffmann (2010), who tested Drosophila melanogaster for heritability of thermal resistance, concluded that there was little evolutionary potential to widen their thermal tolerance via canalized thermoregulation. Thus, plastic responses (e.g. Andrew et al. 2013) may be favored as complementary mechanisms to offset deleterious temperatures, as thermotolerant species can also be plastic (Verberk et al. 2018).
Plastic responses have an advantage over canalized responses in that they can be induced when needed and are adaptive in situations of relatively sudden and extreme shifts in temperature (Angilletta 2009), as in the current situation of climate change. Acclimation is a plastic response that has been de ned as an adaptive compensatory change in an individual caused by a shift in an environmental variable or variables (Bullock 1955, Precht 1955, Prosser 1955). Thermal acclimation can enhance performance near the low (Mutamiswa et al. 2018) and high (Sømme 1982, Dietz & Somero 1992) ends of an individual's debate. However, what drives the evolution of increased plasticity in the eld is less clear. Extreme temperatures alone are not expected to drive the evolution of plasticity, as, in some species, selection for enhanced thermal tolerance can reduce plasticity in thermal tolerance (Kelly et al. 2017). Theory posits that populations will be more plastic when they experience environmental heterogeneity (Lynch & Gabriel 1987, Moran 1992, Scheiner & Holt 2012, Scheiner 2013, particularly when it is predictable (Reed et al. 2010). This has been supported by results from modeling experiments (e.g. Scheiner (Semlitsch et al. 1990, Lind et al. 2010) environmental heterogeneity as drivers of the evolution of developmental (i.e. irreversible) plasticity. While spatial variation may be important, a modeling experiment by Moran (1992) indicated that physiological plasticity was more readily maintained under a situation of temporal variation as compared to a situation of spatial variation, and that species that inhabit narrower geographic ranges and experience less thermal uctuation generally do not possess the physiological machinery necessary to maintain homeostasis across more broad temperature ranges (e.g. Strickland et al. 2016). However, Scheiner (2013) conducted an individualbased simulation of the evolution of plasticity that incorporated both spatial and temporal variation at different times and concluded from the results that spatial heterogeneity was more important than temporal. Thus, the most important type of heterogeneity driving phenotypic plasticity is not universally agreed upon and may not be generalizable.
In a previous experiment, we quanti ed thermal heterogeneity at ve sites across the Midwest of the United States over a 101-year period during our experimental species' active season (Preston and Johnson 2020). The animals we used in this experiment were from three of those ve sites; Minneapolis, Kansas (KS), Marshall, Missouri (MO), and Mount Vernon, Illinois (IL; Suppl. Figure 1). In that study, we determined that KS had a higher daily temperature range (Suppl. Figure 2 and more variance in daily maximal temperatures than MO and IL. Missouri and IL were more comparable, but MO did have more variance in daily maxima than IL (Suppl. Figure 3).
We raised an F2 generation of the generalist grasshopper Melanoplus differentialis (Thomas, 1856) bred from F1 individuals (described in Preston and Johnson 2020) in a warm and cold rearing environment to test the hypotheses that populations experiencing more thermal heterogeneity will: 1) be able to more effectively behaviorally thermoregulate, 2) have a greater ability to withstand extremely high temperatures, 3) display more pronounced acclimation effects in one or more physiological and behavioral thermoregulatory traits in response to variable rearing temperatures, and 4) show differences in gene expression in response to stressfully high temperatures when compared to populations with less thermal heterogeneity.

Materials And Methods
All animals used in this experiment were treated in accordance with the U.S. National Research Council's Guide for the Care and Use of Laboratory Animals, the U.S. Public Health Service's Policy on Humane Care and Use of Laboratory Animals, and the National Research Council's Guide for the Care and Use of Laboratory Animals.
The experimental population in this study was an F2 generation produced from F1 individuals from three states in the United States: Kansas (KS), Missouri (MO), and Illinois (IL; described in Preston and Johnson 2020). We allowed these individuals to mate randomly with others whose parents were from the same site. For convenience, we refer to these groups of F2 individuals (KS, MO, and IL) and their sites of origin interchangeably as populations hereafter. We subsequently followed one of two rearing protocols described in Preston and Johnson (2020), using the diapause treatment that subjected eggs to 0°C (as opposed to 22°C), as it yielded more hatches in that experiment. We only altered this protocol by raising the F2 generation in either a high or low temperature treatment (see Experimental rearing temperatures below) until they reached adulthood, after which we measured them for a variety of traits.

Experimental rearing temperatures
Immediately after F2 individuals from each site hatched, we divided them into cold-and warmtemperature rearing treatments and raised them to adulthood. For the cold rearing environment, we originally had planned to use 15.6°C. We chose this temperature because it 1) is close to what has been used for rearing Melanoplus species in the past (e.g. Harrison 1988), 2) approximates the 101-year mean May (when M. differentialis nymphs are developing) daily minimum temperature (DTMin) that our experimental populations have experienced (11.2°C), and 3) was the coolest temperature we were able to maintain with a reasonable amount of error. However, in a pilot assay, we determined that mortality was too high during rearing at 15.6°C for these populations. Thus, we used 21°C instead.
For the warm rearing environment, we used 30°C. Similar to the temperature originally chosen for the cold-acclimated individuals, this temperature 1) is close to what has been used for rearing Melanoplus species in the past (Jonas & Joern 2013), 2) approximates these populations' previously experienced maximum 101-year mean June and July daily maximum temperature (31.4°C), and 3) was the warmest temperature we were able to maintain within a reasonable amount of error. For both acclimation treatments, we kept relative humidity at 40%.

Accuracy of Behavioral Thermoregulation
We conducted behavioral thermoregulation trials to determine individuals' accuracy of behavioral thermoregulation, i.e. the ability to keep their body temperatures (T b s) close to their preferred temperatures (T pref s). We calculated this using the metric d b , after Hertz et al. (1993), by taking the absolute deviation of each individual's (T b ) from its T pref using protocols described below. We obtained T pref , d b , and critical thermal maximum (CT Max ) from each individual, allowing a 24-hour period to pass between each type of trial to allow individuals to recover. We have listed the sample size for this and all other non-gene expression trials in Suppl. Table 1).
To determine T pref , we placed each individual into a custom-made 10.16 x 91.44 x 10.16 cm wooden shuttle box (Suppl. Figures 4 and 5). We covered the shuttle box with a clear acrylic sheet secured onto the main body of the shuttle box with adhesive hook and loop fasteners. We divided the shuttle box into seven 10.16 x 13.06 x 10.16 cm partitions (7.62 x 9.8 x 4.92 cm of open space in each partition; Fig. 1). We placed a space heater on one side and a window air conditioning unit on the other and connected each to the shuttle box with aluminum duct and detachable ttings. We cut 2.54 cm square notches into the wooden dividers between each partition to allow individuals free passage. The partitions had a temperature range of 6.9 (closest to the air conditioning unit) to 72.9°C (closest to the space heater). The mean steady state temperatures of the partitions were, from the partition closest to the air conditioning unit to that closest to the space heater, 8. To obtain d b , we used a protocol identical to that used to determine T pref described above, with the exception that, every three minutes, we rotated the shuttle box 180 ° such that the space heater was heating the end that the air conditioning unit was previously cooling and vice versa. This presented individuals with a microenvironment that varied drastically in temperature over time, similar to quick shifts in temperature an organism might experience in the eld from a change in cloud cover, wind speed, an incoming front, etc. After four rotations over 15 minutes, we removed each individual from the shuttle box and recorded its thoracic T b as in the T pref trials. We used this T b measurement and each individual's T pref to calculate d b . Being able to obtain this metric in the lab allowed us to obtain an estimate of T pref s unbiased by eld conditions. We carried out all behavioral trials between 10 AM and 2 PM and took care to keep lighting constant from one trial to the next.
Preferred temperatures were not normally distributed and were not easily transformed. Thus, we analyzed them using non-parametric analyses. As a preliminary analysis determined that there were sex differences in T pref s (consistent with other studies, e. g. Forsman 2000), we analyzed the potential effects of weight, site, acclimation, and their interaction terms on sex-speci c T pref s using Kruskal-Wallis and Scheirer-Ray-Hare tests. Values of d b were also not normally distributed, but conformed to a normal distribution when we applied a log 10 transformation. We used a 4-way ANOVA to analyze them by using site, sex, rearing temperature, weight, and their interactions as predictor variables and d b as a response variable. We used total residual sums of squares to obtain F values in this and all subsequent ANOVAs.

Upper Thermal Limits (CT Max )
We measured CT Max with a thermal ramping experiment, following the protocol described in Preston and Johnson (2020). An important limitation to the CT Max trials was that there was only one IL female raised in the cold environment that survived to be tested. As preliminary analyses indicated a signi cant effect of sex on CT Max , we analyzed each sex separately with 3-way ANOVAs using CT Max as a response variable and site, rearing temperature, weight, and their interactions as predictor variables. Where we detected a signi cant interaction effect between site and rearing temperature, we visualized differences between mean CT Max in the cold-reared individuals and mean CT Max in the warm-reared individuals for each population with generalized reaction norms calculated among individuals originating from each population. We did not include IL females when visualizing generalized reaction norms due to having only one cold-acclimated IL female for that analysis.

Heat Stress Symptom Principal Component Analysis
While the ramping test to determine CT Max was ongoing, we recorded air temperatures at which symptoms of heat stress manifested, one of each for each individual. These symptoms were 1) restlessness, de ned as at least ten seconds of continuous movement after at least ve minutes of motionlessness, 2) jumping, 3) onset of back leg spasms, 4) loss of coordinated movement, de ned as the continued movement of extremities without ambulation after at least ve minutes of continuous ambulation, and 5) heat coma. We analyzed air temperatures at which these symptoms manifested with a principal component analysis (PCA), using the air temperatures at which the symptoms above manifested to construct principal components (PCs). As a preliminary analysis indicated that there was a signi cant effect of sex, we analyzed the sexes separately and subsequently used PCs in 3-way ANOVAs with site, rearing temperature, and weight as predictor variables.

Differences in Gene Expression Among Males
For select individuals that we tested for upper thermal limits as described above, we used whole transcriptome pro ling to assess differences in gene expression among treatments. For logistical and nancial reasons, and to control for sex effects in gene expression, we only used males for this assay. Immediately after testing for CT Max , we took whole-head samples and macerated them in Trizol. We stored each sample at -80°C before extracting total RNA using a standard Trizol protocol (Invitrogen; Thermo sher Scienti c, Waltham, United States).
We assessed total puri ed RNA in each sample using a Qubit uorometer (Thermo sher Scienti c, Waltham, United States). Afterwards, we assessed RNA fragment size and integrity using a Bioanalyzer automated electrophoresis tool (Agilent, Santa Clara, United States). We outsourced cDNA library construction and transcriptome sequencing to the University of Texas at Austin, where we had the transcriptomes of two MO individuals sequenced using the NextSeq 500 system (Illumina, San Diego, United States) according to standard RNASeq protocols expressly for the purpose of building a reference transcriptome. The remainder of samples were processed by the HiSeq 2500 system (Illumina, San Diego, United States) using TagSeq (Meyer et al. 2011), a sequencing protocol that costs much less than RNASeq but can be just as or more accurate and informative (Lohman et al. 2016).
We performed the transcriptome assembly in a Linux environment using the RNASeq samples. We We used the DESeq2 package (Love et al. 2014) in R (version 3.6; R Core Team 2019) to normalize and analyze reads from the count table for potential differences in gene expression. One of the reasons we chose this package is because of our low sample size; DESeq2 biases shrinkage by applying it more strongly where information for a gene is low (Love et al. 2014), thus reducing the probability of false positives. This package uses a generalized linear model, in which we speci ed site of origin and acclimation treatment as predictor variables, and log 2 fold change (FC) for each sequence as response variables. Prior to analysis, we removed all rows (genes) that had < 10 reads. We used the package DEVis (Price 2019) to visualize differences in gene expression among sites and treatments by calculating dissimilarity measures for each pairwise comparison of samples to create a multidimensional scaling hull plot. After visualization, it was apparent that one sample (KS, warm-acclimated) was of poor quality and had to be removed from the analysis. Afterwards, we re-ran the analysis and re-visualized the results without the sample.
We used Blast2GO (https://www.blast2go.com/) to infer functionality of signi cantly differentially expressed (DE; α = 0.01) genes (392 total). We used the discontiguous-MEGABLAST algorithm to compare signi cantly differently expressed sequences against the National Center for Biotechnology Information's (Wheeler et al. 2007) standard nucleotide collection (nr) database, applying the Orthoptera taxonomy lter (taxid: 6993). We used a BLAST expectation value of 10 − 1 , a word size of 11, and a highscoring segment pair length cutoff of 33. We used Blast2GO and InterProScan (Quevillon et al. 2005) to obtain two separate lists of gene ontology (GO) terms. Afterwards, we merged the GO terms into a single list. Finally, we used Blast2GO to generate multi-level pie charts to categorize DE sequences by biological, cellular, and molecular GO terms.
After we performed the above steps, we failed to obtain a signi cant match for one of the genes with the most difference in expression between two treatments. Thus, we subjected it to another discontiguous MEGABLAST, but we extended the query to the Insecta database.

Accuracy of Behavioral Thermoregulation
Weight did not affect male (Kruskal-Wallis test, df = 27, c 2 = 28.659, p = 0.38) or female (Kruskal-Wallis test, df = 36, c 2 = 37.892, p = 0.383) T pref . Similarly, neither site, rearing temperature, nor the interaction between the two affected male or female T pref (Scheirer-Ray-Hare test, all p values > 0.05). The only signi cant factor affecting log 10 d b was rearing temperature; cold-reared individuals had a lower mean value of log 10 d b (Table 1, Fig. 2). Table 1 Results from a 4-way ANOVA run on select individuals of Melanoplus differentialis from three populations, reared at warm (30°C) and cold (21°C) temperatures. We used the log 10 of d b (the deviation of an individual's preferred from its actual body temperature) as a response variable and population (Site), sex, rearing temperature (Acclimation), and weight as predictor variables. Signi cant result highlighted in yellow.  Fig. 4). loss of movement occurred. PC2 loaded very strongly and positively for air temperatures at which spasms began, and strongly and negatively for air temperatures at which coma set in (Table 3). The only factor signi cantly affecting male PC1 was weight, whereas only acclimation signi cantly affected PC2 (Table 4). Heavier males tended to have higher values of PC1 (r = 0.509), and warmacclimated males had a higher mean value of PC2 (Fig. 5). Thus, heavier males withstood higher temperatures before manifesting most symptoms, and warm-acclimated males withstood higher temperatures before spasming. In the female PCA, PC1 and PC2 had eigenvalues of 2.640 and 1.221, respectively. The remaining PCs had eigenvalues < 1 and were not considered further. The proportion of variance explained by PC1 and PC2 were 0.528 and 0.244, respectively. PC1 strongly and positively loaded for air temperatures at which restlessness, jumping, and spasm set in. PC2 loaded strongly and positively for air temperatures at which loss of movement and coma occurred (Table 3).
Warm-acclimated females had a signi cantly higher mean value of PC1 (Table 5; Fig. 6), and warmacclimated MO females had a higher mean value of PC2 than cold-acclimated MO females, whereas warm-acclimated KS females had a lower mean value of PC2 than cold-acclimated KS females (Fig. 7).
While this interaction effect was signi cant, its effect size (η p 2 ) was only half a percent (Table 5). PC2 values of IL females are not visualized in Fig. 7, as there were no cold-acclimated females in this trial.  MO had the greatest number of DE genes between acclimation treatments, with a greater number of genes downregulated than upregulated in the warm acclimation treatment when compared to the cold acclimation treatment. (Fig. 8). When we compared differences in expression among sites while keeping acclimation treatment constant, the largest difference in expression existed in the KSWarm-MOWarm comparison, with more genes upregulated in the MOWarm treatment ( Fig. 9), with the total number of DE genes being approximately one and a half times that of the MOWarm-MOCold comparison (Fig. 8).
The Gene expression pairwise distances as determined by multidimensional scaling did not have any discernible pattern when grouped by acclimation temperature (Fig. 12). However, they appeared to group somewhat strongly by site of origin (Fig. 13).
Of the ve sequences that had an absolute log 2  We obtained signi cant BLAST hits for 299 out of 392 sequences when running the discontiguous MEGABLAST. Of these, we were able to obtain GO annotation for 182 sequences. The three most common categories of GO terms assigned to DE sequences, when grouped by biological process, were 1) ATP synthesis, 2) mitochondrial electron transport, and 3) aerobic electron transport (Fig. 14). When grouped by cellular component, they were 1) integral component of membrane, 2) respiratory chain complex IV, and 3) mitochondrial proton-transporting ATP synthase complex, coupling factor (Fig. 16).

Accuracy of Behavioral Thermoregulation
No tested factors had an effect on T pref s. It is possible that this was an artifact of our protocol; T prefs of warm-acclimated individuals may actually been higher than that of cold-acclimated individuals, but steeper temperature gradients near the warm end of the shuttle box may have made it harder for individuals selecting warmer temperatures to maintain a constant temperature, whereas individuals selecting colder temperatures may have been able to more systematically choose and maintain T prefs in the relatively gradual temperature gradient present near the colder end of the shuttle box. Contrary to our results, Forsman et al. (2002) found that acclimation affected orthopteran T pref . Their study used an 800 x 350 x 3 mm copper plate with a heat retaining hot tray at one end and a polystyrene foam box with ice at the other. The plate had a temperature range of 14-50 ˚C, but they do not disclose the steepness of the temperature gradient. If their gradient was much smoother than the one we utilized, then it is likely that individuals used in their experiment were more easily able to maintain T pref s regardless of rearing temperature. As T pref was a critical component of our calculation of d b , this may have been responsible for detection of a lower d b in cold-acclimated individuals.

Upper Thermal Limits (CT Max )
Male CT Max was higher in KS males than in IL males. Preston and Johnson (2020) found that the KS site had more variable temperatures and more days per month with extreme heat events than the IL site, supporting our hypothesis of a positive correlation between thermal tolerance and heterogeneity.
However, that study did not detect signi cant effects of population on male CT Max . It's possible that maternal effects had more of an effect on weight (which affected CT Max ) in that study, as those males were an F1 generation raised from eggs laid by eld-collected individuals, and the males used in this study were F2 individuals raised from eggs laid by F1 individuals.
Illinois females had a higher mean CT Max than MO females, which was surprising given that IL males had the lowest mean CT Max . A study with another Melanoplus species found that males and females in a single population differed in their epicuticular lipid composition (Gibbs & Mousseau 1994), which affects water loss at high temperatures (Noble-Nesbitt 1991), which in turn could affect CT Max . If M. differentialis also has sex differences in epicuticular lipid composition, and these differences are more pronounced in different populations, these differences may be responsible for the mismatch in relative CT Max among males and females at each site.
There were no signi cant acclimation effects in males, agreeing with the conclusion by Gunderson and protein transcription than cold-acclimated (0.4 ˚C) lobsters when exposed to heat stress. While identifying the function of every differentially expressed gene is beyond the scope of this study, it is reasonable to conclude that, as warm-acclimated males responding to stressfully high temperatures had more differences in gene expression, these genes may be relevant to thermoregulation and thermotolerance.
The largest (by far) number of DE genes were between KS and MO. Gene expression grouped more strongly by population than by acclimation treatment, with the notable exception of sample 15C, a coldacclimated IL male. This agrees with the relative pooled mean values of the Wald tests as well as the relative amount of log 2 FCs. This may be due to substantial genetic variation among these populations, minimal ability to acclimate in these populations, or a combination of the two. Dunning et al. (2014) observed similar results with stick insects and gene expression in response to low temperatures; while all populations showed a transcriptional response to cold, the majority of the unigenes identi ed were population speci c. Thus, while populations (or groups of species) may mount a transcriptional response to a common stressor, the genes involved and the strength of the response can differ by virtue of genetic differentiation among populations. As populations in this study grouped more strongly by site than acclimation treatment, it may be that the transcriptomic responses of the populations we sampled were more in uenced by population-speci c evolutionary in uences than they were by different rearing temperatures.
When categorized by biological process, the three most common roles of differentially expressed genes with GO terms were ATP synthesis coupled proton transport, mitochondrial electron transport, and facilitation of the aerobic electron transport chain. This indicates that males may have differed in their rate of cellular respiration. As temperature increases, insect respiration and metabolism are both expected to increase up to a critical thermal limit (Neven 2000). The oxygen limitation model posits that thermal limits of performance (in this case, CT Max ) are set by the point at which aerobic respiration fails to meet energetic needs. In the case of high temperatures, this would occur when ventilation and circulation fall below the level required to supply mitochondria with su cient oxygen (Pörtner 2001   Generalized reaction norms for critical thermal maximum (CTMax) of Melanoplus differentialis females from two populations reared at warm (30 °C) and cold (21 °C) temperatures; ± 1 SE.        Distribution of node scores derived from differentially expressed sequences obtained from head samples of male Melanoplus differentialis from three populations reared in warm (30 °C) and cold (21 °C) acclimation treatments, grouped by biological function. These scores were determined from data generated by discontiguous-MEGABLAST and subsequent gene ontology annotation performed in Blast2GO. Figure 15