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.
2.1 Experimental rearing temperatures
Immediately after F2 individuals from each site hatched, we divided them into cold- and warm-temperature 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%.
2.2 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 (Tbs) close to their preferred temperatures (Tprefs). We calculated this using the metric db, after Hertz et al. (1993), by taking the absolute deviation of each individual's (Tb) from its Tpref using protocols described below. We obtained Tpref, db, and critical thermal maximum (CTMax) 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 Tpref, 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 fittings. 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.2, 9.2, 11.2, 15.4, 18.0, 44.3, and 67.2°C, respectively. For each Tpref trial, we placed an individual into the middle partition and started the trial immediately. After a 15-minute period, we recorded Tb of each individual by immediately inserting a hypodermic needle with an embedded thermocouple under the metasternum of the thorax.
To obtain db, we used a protocol identical to that used to determine Tpref 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 field 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 Tb as in the Tpref trials. We used this Tb measurement and each individual's Tpref to calculate db. Being able to obtain this metric in the lab allowed us to obtain an estimate of Tprefs unbiased by field 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 Tprefs (consistent with other studies, e. g. Forsman 2000), we analyzed the potential effects of weight, site, acclimation, and their interaction terms on sex-specific Tprefs using Kruskal-Wallis and Scheirer-Ray-Hare tests. Values of db were also not normally distributed, but conformed to a normal distribution when we applied a log10 transformation. We used a 4-way ANOVA to analyze them by using site, sex, rearing temperature, weight, and their interactions as predictor variables and db as a response variable. We used total residual sums of squares to obtain F values in this and all subsequent ANOVAs.
2.3 Upper Thermal Limits (CTMax)
We measured CTMax with a thermal ramping experiment, following the protocol described in Preston and Johnson (2020). An important limitation to the CTMax trials was that there was only one IL female raised in the cold environment that survived to be tested. As preliminary analyses indicated a significant effect of sex on CTMax, we analyzed each sex separately with 3-way ANOVAs using CTMax as a response variable and site, rearing temperature, weight, and their interactions as predictor variables. Where we detected a significant interaction effect between site and rearing temperature, we visualized differences between mean CTMax in the cold-reared individuals and mean CTMax 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.
2.4 Heat Stress Symptom Principal Component Analysis
While the ramping test to determine CTMax was ongoing, we recorded air temperatures at which symptoms of heat stress manifested, one of each for each individual. These symptoms were 1) restlessness, defined as at least ten seconds of continuous movement after at least five minutes of motionlessness, 2) jumping, 3) onset of back leg spasms, 4) loss of coordinated movement, defined as the continued movement of extremities without ambulation after at least five 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 significant 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.
2.5 Differences in Gene Expression Among Males
For select individuals that we tested for upper thermal limits as described above, we used whole transcriptome profiling to assess differences in gene expression among treatments. For logistical and financial reasons, and to control for sex effects in gene expression, we only used males for this assay. Immediately after testing for CTMax, 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; Thermofisher Scientific, Waltham, United States).
We assessed total purified RNA in each sample using a Qubit fluorometer (Thermofisher Scientific, 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 assessed the quality of raw reads with FastQC version 0.11.9 (Andrews 2010) and subsequently trimmed low-quality bases using Trimmomatic version 0.39 (Bolger et al. 2014). Afterwards, we used the Trinity package version 2.8.6 (Grabherr et al. 2011; Haas et al. 2013) to construct a de novo transcriptome. We then used BUSCO version 3.0.2 (Seppey et al. 2019) to validate the completeness of our assembly against the Insecta database from the BUSCO website (https://busco.ezlab.org/), searching a total of 1658 BUSCO groups. Finally, we used Bowtie2 version 2.3.5 (Langmead & Salzberg 2012, Langmead et al. 2019) to index the assembled transcriptome.
We used the TagSeq samples to generate a count table in order to examine differential gene expression. We used FASTX-Toolkit’s (version 0.0.14; Gordon & Hannon 2010) Fastx_clipper tool in conjunction with the Cutadapt tool (Martin 2011) to remove common TagSeq adapters and primers. We then used Bowtie2 version 2.3.5(Langmead & Salzberg 2012, Langmead et al. 2019) to map the TagSeq reads to the reference transcriptome. We generated read-counts-per gene isoform using a custom perl script (Mikhail Matz, personal comm.), then counted hits per gene using SAMtools (Li et al. 2009).
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 specified site of origin and acclimation treatment as predictor variables, and log2 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 significantly differentially expressed (DE; α = 0.01) genes (392 total). We used the discontiguous-MEGABLAST algorithm to compare significantly differently expressed sequences against the National Center for Biotechnology Information's (Wheeler et al. 2007) standard nucleotide collection (nr) database, applying the Orthoptera taxonomy filter (taxid: 6993). We used a BLAST expectation value of 10− 1, a word size of 11, and a high-scoring 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 significant 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.