Large-scale genome-wide study reveals climate adaptive variability in a cosmopolitan pest

Genomic data

The foundational resource for this study was a dataset of 40,107,925 nuclear SNPs sequenced from a worldwide sample of 532 DBM individuals collected in 114 different sites based on our previous project15. DNA was extracted from each of the 532 individuals using DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol, and eluted from the DNeasy Mini spin column in 200 μl TE buffer. Genomic sequencing was performed with Illumina HiSeq 2000 at BGI, Shenzhen, China, to produce 90 bp paired-end reads for every individual. Using custom scripts, raw reads were processed to filter out poor reads with 10 ambiguous “N” bases, >40% low-quality bases, or identical sequences at the two ends and obtain clean reads. The clean reads were mapped onto the DBM reference genome (v2)21 with Stampy (v1.0.27) using default parameters. SNP calling was then performed using the GATK HaplotypeCaller with parameters –emitRefConfidence GVCF –variant_index_type LINEAR –variant_index_parameter 128,000. The 40,107,925 nuclear SNPs generated present one variant on average in every six bp of the reference genome, which is the densest variant map for any organism, including the recently released data on human58 and Arabidopsis thaliana genome sequences59. The SNP dataset is available at www.ebi.ac.uk/ena with the accession code PRJEB24034.

In the present study, to investigate the genetic variation associated with climate, we excluded samples from the regions that are only seasonally suitable for growth of DBM (with the Ecoclimatic Index EI = 0)16. This was done because in these regions, populations are unlikely to receive perennially unpunctuated selection by local environmental variables and genetic variation cannot be continuously passed on to future generations over years. Specifically, regions that are only seasonally suitable for DBM growth and development (i.e., with an ecoclimatic index, EI = 0) are too harsh to allow survival in low temperature conditions during the winter. Annual recolonization of those regions from areas where DBM can overwinter (with an ecoclimatic index, EI > 0) has been biologically and genetically confirmed14,60,61,62,63. No genetic differentiation was found among different geographical populations spanning from overwintering regions to seasonally inhabited regions15,62,63. If migration from one habitat overwhelms the other, migration from the source introduces new genetic variation that may prevent local adaptation64,65. The retained samples included 372 individuals from 78 sampling sites in the year-round persistence regions of DBM across the world (where EI > 0). These samples were collected from different continents, with 13 samples from Africa, 29 from Asia, 5 from Europe, 13 from North America including Hawaii, 12 from South America, and 6 from Oceania (Supplementary Fig. 1 and Supplementary Data 1).

The retained 372 individuals shared a subset of 34,969,375 SNPs, which accounted for 87.19% of the total SNPs (40,107,925) and represented most of the genomic variation among 532 individuals (Supplementary Data 2)15. Using VCFtools v.0.1.6, we excluded the SNPs with minor allele frequency (MAF) < 5% and missing rate >10%. We then sampled data by examining single SNPs in small, 25 bp DNA window to focus on loci independent of linkage disequilibrium. After quality filtering, a total of 200,055 bi-allelic SNPs across 357 DBM individuals collected from 75 different sites worldwide was retained for further analysis (Supplementary Fig. 1, Supplementary Data 1 and 2).

Climate Data

Nineteen climate variables related to temperature and precipitation (Supplementary Table 1), which are known to have impacts on physiological and ecological traits of insects, were retrieved with high resolution from WorldClim, a public database (www.worldclim.org)66.

Climate associated genomic variation

Identification of SNPs under climate selection

Three models, Samβada v0.5.317, latent factor mixed model (LFMM) v1.418 and Bayenv 219, were used to identify putatively adaptive loci associated with climate variables. Samβada identifies associations between specific genetic markers and environmental variables by logistic regression17. Simple univariate and multivariate logistic regression models for each climate variable were fitted. A single SNP was considered to be a candidate locus when the log-likelihood ratios (G scores) and/or Wald scores were significant with Bonferroni correction at a 99% confidence level. LFMM is based on population genetics, ecological modelling, and statistical analysis to identify the candidate loci that are highly correlated with environmental variables18. SNPs showing an association with climate variables were identified based on z-scores, which was computed using 10,000 cycles and 5000 sweeps for burn-in. We used the R package LEA to estimate the median z-scores of 5 runs and re-adjusting p-values with FDR correction67. SNPs with median z-scores above the absolute value of 4 and corresponding to P value < 10−5 were considered as significant locus. In Bayenv 2, a covariance matrix based on putatively neutral markers is used as a null model to control for demographic effects when testing relationships between the genetic differentiation and a given environmental variable19. We randomly sampled SNPs at 200 SNP intervals from the SNP dataset. A total of 117,887 SNPs with loose linkage disequilibrium were obtained for developing the covariance matrix, which was estimated with 100,000 iterations. We then assessed the correlations between individual SNP and 19 climate variables at 100,000 Markov chain Monte Carlo (MCMC) for Bayes factor analysis. Five independent runs of the Bayenv program were performed with different random seeds68. The results were presented as Bayes factors (BFs). An averaged log10(BF) value of five runs >1.5 is considered as high support for a model where environmental parameters have significant effects on allele frequencies69. A total of 3648 putative adaptive SNPs were identified by at least one of the three models (Samβada, LFMM, and Bayenv 2; Supplementary Data 3, 4 and 5)

Prediction of climate-associated genomic variation

Generalized dissimilarity modelling (GDM)20, a distance-based method, can account for the nonlinear relationship between genetic variation and environmental/geographical factors, and has been recently used to map ecological adaptation from genomic data under current and future climates4. First, we examined a spatially explicit selection process for each of the putative adaptive SNPs using GDM4, with the R package gdm70. We subsampled the genetic dataset to include only populations with sample size ≥5 (60 populations) to obtain accurate allele frequencies. To reduce the number of bioclimatic variables, we preferentially discarded those with multiple correlated variables (Supplementary Table 5). We then ran the GDM model using 3,648 SNPs and the variables with Pearson |r| < 0.8 (bio02, bio03, bio07, bio08, bio09, bio12, bio15, bio18, and bio19), one of bio01 and bio11 (|r| = 0.92), one of bio5 and bio 10 (r| = 0.94), as well as one of bio14 and bio17 (|r| = 1.0). The subset of 12 bioclimatic variables with highest value of explained deviance (28.17%) in the GDM model were retained, including bio01, bio02, bio03, bio07, bio08, bio09, bio10, bio12, bio14, bio15, bio18, bio19, for further analysis. Pairwise FST matrix among 60 populations were calculated for each of the 3648 SNPs using the R package hierfstat71, and rescaled between 0 and 1. Geographical distance in the GDM was based on Euclidean distance as the thirteenth variable to test whether genetic variation across environmental gradients was better explained by climate variables than geographical distance, which effectively acts as a screening for SNPs that may respond predominantly to neutral genetic process including isolation by distance72. The relative importance of the 12 climate variables and geographical distance was ranked based on the fitted I-Splines in GDM (Fig. 1b). The maximum value of each variable in the fitted I-Splines was rescaled between 0 and 1. Those SNPs with geographical distance ranking as one of the 3 most important variables were excluded in the following GDM analysis. In addition, we randomly sampled 200 SNPs as a “reference group” to test its explainable proportion of the GDM deviance. According to our GDM record, the reference SNP group accounted for 11.2% of the GDM deviance for the entire model, so that those SNPs with a < 11.2% contribution to the GDM deviance were also excluded in further GDM analysis. After additional filtering, 517 of 3,648 SNPs were retained. The 517-SNP-based genetic distance matrix was further integrated with geographical distance and climate variables to be used in the entire GDM model, that explained 42.80% of the deviance for the 517 SNPs. To predict the climate adaptation of DBM, we finally retrieved current climate variables at 61,655 gridded points across the world from WorldClim, using ArcGIS 10.2. The gdm.transform function was used to predict and map the pattern of climate-associated genomic variation along environmental gradients across the world (Fig. 1c). The genetic turnover was summarized using a principal component analysis (PCA), with the top three components transformed for visualization in a red-green-blue (RGB) color scale as suggested in the GDM manual70. Loadings based on the principal components indicate the direction and magnitude of association with adaptation to different predictors (Fig. 1c). The genetic variation along environmental gradients in DBM across the world was visualized, with similar pattern of genetic composition at climate-adaptive loci illustrated by similar colors (Fig. 1c).

Genetic vulnerability and adaptive potential

Prediction of genetic vulnerability

To predict the population-level variation in “genetic vulnerability” (GV) under future climate scenarios, we used the “Genetic offset” (GO) method of Fitzpatrick and Keller4. The GO represents the extent of mismatch between current and expected future genetic variation based on genotype-environment relationships modelled by GDM analysis73. The projected future climate variables of 2050 and 2080 for four different greenhouse gas scenarios, Representative Concentration Pathways (RCPs), including RCP2.6, RCP4.5, RCP6.0 and RCP8.5 based on the NorESM1-M Global Climate Model (GCM) across the world were retrieved from WorldClim, using ArcGIS 10.2. Those four RCPs represent different gas emission scenarios, reflecting mild (RCP2.6) to extreme (RCP8.5) conditions74. The GO was predicted by predict.gdm function in R package gdm. A metric of GO for each of the gridded climate points was obtained, which implied that the populations with greater genetic offset would be more vulnerable or less tolerant to the future climate change. The resulting GOs were rescaled between 0.1 and 0.975,76, and then mapped with ArcGIS 10.2 to show the geographical distributions of population-level variation in genetic tolerance to future climate changes (Fig. 2a).

Prediction of habitat suitability

The CLIMEX model, which has been shown to be effective for examination of species distribution under future climate scenarios77, was used to predict the habitat suitability for DBM in 2050 under the RCP8.5 scenario, the only one that is available in CliMond (www.climond.org/). The CLIMEX model for DBM in Zalucki and Furlong16 was developed based on temperature, moisture, and stress indices. Temperature indices include limiting minimum (DV0) and maximum temperature (DV3), lower (DV1), and upper (DV2) optimal temperature. Moisture indices include minimum (SM0) and maximum (SM3) tolerable soil moisture, lower (SM1) and upper (SM2) optimal soil moisture. Stress indices include cold stress, heat stress, dry stress, wet stress, and hot-wet stress. The values of these parameters for prediction of habitat suitability in our study were taken from Zalucki & Furlong16. We altered one of the parameters in CLIMEX16, changing the hot-wet stress temperature threshold from 30 °C to 32 °C based on a study on the relationship between temperature and developmental rate showing that DBM can survive and develop at temperatures <32 °C78 (Supplementary Table 6). The resulting ecoclimatic index (EI) values based on our CLIMEX simulation were used to calculate the difference in ecoclimatic index (DEI) between current and projected future climate scenarios using the equation: DEI = EIF – EIC, where EIF is the ecoclimatic index under the projected future climate scenario and EIC is the ecoclimatic index under the current climate. Functions of Inverse Distance Weighting (IDW) and Overlay in ArcGIS 10.2 were performed to generate maps showing the predicted distribution of DEI values (with DEI > 0 showing the EI-increased regions and DEI < 0 showing the EI-decreased regions) for each of the gridded climate points across the year-round persistence regions worldwide (Fig. 2b).

Prediction of eco-genetic adaptation

The ecoclimatic index (EI) values based on CLIMEX simulations are highly generalized and do not reflect evolutionary adaptation and variation of tolerance levels among populations, because of the absence of information on fitness traits. Based on our recently-generated dataset of genome-wide single nucleotide polymorphisms (SNPs) sequenced from a worldwide DBM sample of different locations (sites) across a diverse range of biogeographical regions15, we developed a new metric, the “eco-genetic index” (EGI), which combines the predictions based on genetic offset (GO) with the difference in ecoclimatic index between current and projected future climate scenarios (DEI). EGI allows us to incorporate our population-specific genomic data into EI for better predicting the habitat suitability of DBM subject to climate in 2050 under RCP8.5 scenario based on NorESM1-M GCM. Because DEI and GO, are correlated (Pearson’s R79 = 0.53, P < 0.001; Supplementary Fig. 5), we used a linear normalization transfer function to make values of deii and goi dimensionless, and then used the weighted geometric averaging (WGA) operator and the artificial bee colony (ABC) algorithm to optimize the value of alpha and improve the algorithm of EGI. Here, only regions of decreasing EI (DEI < 0) were considered because in these regions, DBM populations will be challenged by habitat suitability decline under the RCP8.5 scenario for 2050. Each gridded point in ArcGIS was described as a vector: Pi = {deii,goi}, where i = 1, 2, …, n, and n is the number of gridded points. Let the EGI and GO of each point be egii and goi, calculated as follows:

Step 1: Calculate the absolute values of deii and goi and normalize75,76 them into [0.1, 0.9].

Step 2: Combine the normalized deii and goi with the weighted geometric averaging (WGA) operator80 as:

$${eg}{i}_{i}={de}{{i}_{i}}^{alpha }g{{o}_{i}}^{1-alpha }$$

(1)

where α ϵ [0, 1] is the weight of normalized deii, i = 1, 2, …, n.

Step 3: The optimal value of α can be determined with the following steps:

Step 3.1: Take the natural logarithm of Eq. (1):

$${ln},e,g{i}_{i}=alpha * {ln},{de}{i}_{i}+left(1-alpha right){ln},g{o}_{i},i=1,,2,,cdots ,,n$$

(2)

Obviously, ({ln},{de}{i}_{i}, < ,0) and ({ln},g{o}_{i}, < ,0)

Step 3.2: Assuming a, b > 0, there exists a theorem such that:

$$frac{a}{b}+frac{b}{a}-2=frac{{left(a-bright)}^{2}}{{ab}}ge 0$$

(3)

The equality in Eq. (3) holds only if a = b.

To balance the indices of deii and goi, we minimized the total deviation between egii and deii, goi. According to Eq. (3) and Zhou et al.81, the deviations di and ti are formulated as follows:

$${d}_{i}=frac{-{ln},eg{i}_{i}}{-{ln},{de}{i}_{i}}+frac{-{ln},{de}{i}_{i}}{-{ln},eg{i}_{i}}-2$$

(4)

$${t}_{i}=frac{-{ln},eg{i}_{i}}{-{ln},g{o}_{i}}+frac{-{ln},g{o}_{i}}{-{ln},eg{i}_{i}}-2$$

(5)

Hence, the minimized model (M − 1) can be expressed as:

$$(M-1):min ,Y,= ,mathop{sum }limits_{i=1}^{n}({d}_{i}+{t}_{i})\ = ,mathop{sum }limits_{i=1}^{n}left[!left(frac{-,ln,eg{i}_{i}}{-,ln,de{i}_{i}}+frac{-,ln,de{i}_{i}}{-,ln,eg{i}_{i}}-2!!right)+left(frac{-,ln,eg{i}_{i}}{-,ln,g{o}_{i}}+frac{-,ln,g{o}_{i}}{-,ln,eg{i}_{i}}-2!!right)!right]$$

(6)

Based on Eq. (3), the minimized model (M − 1) can be equivalently written as:

$$(M-2):,min Y,= ,mathop{sum }limits_{i=1}^{n}left[frac{alpha ,ln,de{i}_{i}+(1-alpha )ln,g{o}_{i}}{ln,de{i}_{i}}+frac{ln,de{i}_{i}}{alpha ,ln,de{i}_{i}+(1-alpha )ln,g{o}_{i}}right.\ +left.frac{alpha ,ln,de{i}_{i}+(1-alpha )ln,g{o}_{i}}{ln,g{o}_{i}}+frac{ln,g{o}_{i}}{alpha ,ln,de{i}_{i}+(1-alpha )ln,g{o}_{i}}-4right]$$

(7)

Step 3.3: The Artificial bee colony (ABC) algorithm82 is used to solve model (M − 2), which is a fractional programming problem where Y is the objective function and α is the independent variable. The parameters are set as: population size = 20, number of iterations = 30, and limit = 5. The convergence curve of the ABC algorithm for solving the above model is plotted in Supplementary Fig. 6. The algorithm converged at the sixth iteration with the optimal estimate α = 0.5046 (Supplementary Fig. 7). Thus, we get ({eg}{i}_{i}={de}{{i}_{i}}^{0.5046}{g}{{o}_{i}}^{0.4954}). The resulting EGI values were then mapped with ArcGIS 10.2 to show its geographical distribution in the EI-decreased regions under the projected future climate scenario (Fig. 2c). It is noteworthy that this analysis used high-quality genomic data with individuals collected from 75 sampling sites across a wide range of eco-climate regions worldwide. This allowed us to estimate the population-level genetic variation in response to climate change. However, we currently do not have the site-specific physiological DBM data available for calculation of the eco-climate index. Future detailed studies on physiological variation in different populations will further improve our prediction of eco-genetic adaptation to climate change.

Genetic basis of climate adaption

Gene expression analysis by RT-qPCR

The wild-type strain of DBM, Geneva (G88), was used in this assay. This G88 strain was collected from the New York State Agricultural Experiment Station in 1988 and has since been maintained on artificial diet without exposure to insecticides83. It was provided by Dr. Antony M. Shelton (Cornell University, USA) to the Institute of Applied Ecology, Fujian Agriculture and Forestry University in 2016. Since then, we maintained this strain on artificial diet without exposure to insecticides at 26 °C that is a favorable temperature to rear and maintain this wild-type strain.

We used nine temperature treatments and one control at 26 °C. A male or female individual was placed in a 1.5 ml plastic vial (4.0 cm height) with a pinhole in the side wall to allow air exchange. Before treatments, all vials with DBM were placed into the incubator at 26 °C. A group of thirty vials (15 vials with DBM females and 15 with males) was frozen in liquid nitrogen and used as controls. Additional groups of thirty vials (15 containing females, 15 containing males) were exposed to each of nine distinct temperature treatments. These treatments were defined from a previous study on lethal temperature limits of DBM84 and our preliminary experiments: three high-temperature treatments: (1) H1: 40 °C for 30 min, (2) H2: 43 °C for 30 min, (3) H3: 43 °C for 30 min with 24 h of recovery at 26 °C; and six low-temperature treatments: (1) L1: −14 °C for 30 min, (2) L2: −14 °C for 30 min with 24 h of recovery at 26 °C, (3) L3: −17 °C for 30 min, (4) L4: −17 °C for 30 min with 24 h of recovery at 26 °C, (5) L5: −20 °C for 15 min, and (6) L6: −20 °C for 15 min with 24 h recovery at 26 °C (Fig. 3). High and low temperature treatments were conducted in incubators and freezers, respectively. The exposure duration of moths at −20 °C was set for 15 min because moths started to die when exposed to −20 °C for over 20 min. After each of the treatments, moths were immediately frozen in liquid nitrogen. The thirty moths from each treatment were grouped into three replicates of 5 females and 5 males in tubes of 1.5 ml each (in total, six tubes each containing 5 individuals of same sex for each treatment including controls). All tubes with frozen moth samples were stored at −80 °C before RNA extraction.

Total RNA was extracted using Eastep® Super Total RNA Extraction kit (LS1040, Promega, USA), following the manufacturer’s protocol. A NanoDrop Spectrophotometer (ND2000, Thermo Scientific, USA) and 2% agarose gel electrophoresis were used to determine the quality and quantity of RNA. cDNA was synthesized with 1 μg of total RNA using GoScriptTM Reverse Transcription System (A5001, Promega Corporation, USA). RT-qPCR was carried out using the target gene-specific primer (forward: 5′-AACCCCCCCTTCATCCAAG-3′, reverse: 5′-CTGCTGAGGCTGTAGGTCATG-3′), which was designed by Oligo 7 and the reference gene primer for normalization (forward: 5′-CAATCAGGCCAATTTACCGC-3′, reverse: 5′-CTGGGTTTACGCCAGTTACG-3′) according to the previous reports85. qRT-PCR was conducted with 2 µL cDNA, 0.4 µL of each primer, 0.15 µL CXR Reference Dye, 7.05 µL DEPC water and 10 µL SYBR Green Supermix (A6001, Promega Corporation, USA) in a 20 µL total reaction mixture using QuantStudioTM 6 Flex real-time PCR system (ThermoFisher Scientific, USA). PCR conditions were set as follows: 10 min at 95 °C followed by 40 cycles of 15 s at 95 °C, 30 s at 60 °C and then 15 s at 95 °C, 1 min at 60 °C, 15 s at 95 °C for a melt curve. The relative expression level of PxCad in each treatment was normalized to the abundance of samples under control temperature, using the 2−ΔΔCt method86. In the present study, we focused on the comparison of difference in gene expression between each of the temperature treatments and the control rather than the comparison between different non-control temperature treatments. Therefore, we used independent t-tests to perform in R to establish differences in gene expression between each temperature treatment and the control.

CRISPR/Cas9-based genome editing

The wild-type strain of G88 was used for functional validation of PxCad. We selected two sgRNA target sites, 5’-AGGGAGATGTCATCATTGCA-3’ and 5’-ACTCGTCCATGTAGATGTTG-3’, in PxCad exon 3 according to the principle of 5’-N20NGG−3’ (with the PAM sequence underlined). To obtain the templates for in vitro transcription of the two sgRNAs, we performed PCR with two primer pairs (Supplementary Table 7) using the KOD-Plus-Neo Kit (TOYOBO, Osaka, Japan), and then purified the products using the Gel Extraction Kit (Omega, Morgan Hill, GA, USA). In vitro transcription was then performed to generate two sgRNAs using the HiScribe™ T7 Quick High Yield RNA Synthesis Kit (New England Biolabs, Ipswich, MA, USA) and Cas9 mRNA using the HiScribe™ T7 ARCA mRNA Kit (with tailing) (New England Biolabs, Ipswich, MA, USA) using linearized PTD-T7-Cas9 vector as a template87. The synthesized sgRNAs and Cas9 mRNA were further purified using phenol–chloroform extraction and ethanol precipitation and stored at −80 °C until use.

Fresh eggs from G88 moths laid within 15–20 min were injected with a mixture of Cas9 mRNA (380 ng/μl) and sgRNAs (100 ng/μl) using an IM 300 Microinjector (Narishige, Tokyo, Japan). Microinjected eggs were then incubated at 26 °C and allowed to develop to adult (G0). Virgin G0 adults were crossed with the wild-type G88 strain in single pairs to produce the G1 progeny. After oviposition, we sacrificed G0 adults and extracted the genomic DNA using the Tissue DNA Kit (Omega, Morgan Hill, GA, USA). gDNA was then used to amplify the fragment containing the sgRNA target sites using PCR with the Phanta Max Super-Fidelity DNA Polymerase (Vazyme, Nanjing, China) and a primer pair (Supplementary Table 7). The generated PCR products were Sanger sequenced by the Biosune Biotech Company (Fuzhou, China) to examine the mutation in the target sequence region. After mutation identification, we focused on G1 progeny derived from mutated G0 parents. Retained G1 siblings were crossed in single pairs to generate G2 progeny and progeny of heterozygous G1 parents were kept after molecular identification. Similarly, we performed single-pair crosses between G2 siblings and kept only G3 from homozygous mutant G2 parents. Retained G3 siblings were pooled to establish the PxCad knockout strain (MU, G88-Cad).

To verify the absence of PxCad protein in the PxCad knockout strain, we isolated midgut brush border membrane vesicle (BBMV) proteins from the G88 strain and knockout strain using the differential magnesium precipitation method88. Total extracted BBMV proteins (30 μg) were separated on 7.5% SDS-PAGE, and the regions (~120–250 kDa) predicted to contain PxCad proteins were excised from the gel stained with Coomassie blue. These gel slices were subjected to tryptic digestion and nano-LC-MS/MS analysis at Huada Protein Research Center (Shenzhen, China).

Bioassay: behavioral responses to different temperatures

DBM strains used in this assay were wild-type (G88 or WT) and mutant (G88-Cad or MU) adults. Female or male individuals emerging within 24 h were placed into a 1.5 ml clear plastic vial (4.0 cm height) with a pinhole in the side wall. Twenty vials with females and twenty with males were then put into a plastic container (16.5 cm length, 11.5 cm width, 3.5 cm height).

Four plastic containers with DBM adults were placed into a climate incubator (DRX-400-DG, Ningbo Saifu Experimental Instrument CO., LTD, China) where DBM adults were exposed to one of the following heat shock treatments: 40 °C, 41 °C, 42 °C, 43 °C, and 44 °C, for 2 h. After each heat shock treatment, the plastic containers were removed from the climate incubator and placed in a temperature-controlled room at 26 °C for 24 h, after which survival was recorded.

A similar procedure was used for cold stress, with four plastic containers of DBM adults being placed into a freezer. Three temperature treatments were used: −14 °C, −17 °C, and −20 °C. When the temperature was set at −14 °C, the adults were exposed for periods of 150, 210, 270, and 330 min. At −17 °C, durations were 40, 60, 80, 100, and 120 min. At −20 °C, durations were 20, 25, 30 and 35 min. After cold treatments, adults were transferred to a temperature-controlled room at 26 °C for 24 h, after which survival was recorded.

Statistical analysis

Logistic regression models were fitted to the survival data to test behavioral responses to different temperatures. Results of the cold-exposure experiment were analyzed first using the model: log[p/(1−p)] = aijk+bijk t where p is the probability of survival, i is the temperature index (−14, −17, or −20 °C), j is the strain index (WT or MU), k is sex (male or female), and t is exposure duration (min) used as a continuous predictor. The responses were complex, with significant interactions to the third order. Therefore, to allow clearer interpretation of the results, a separate analysis was performed for each temperature using the simpler model log[p/(1 − p)] = ajk+bjk t.

The results of the heat-shock experiment were analyzed using the model log[p/(1 − p)] = ajk + bjk T where j is the strain (WT or MU), k is sex (male or female), and T is temperature (°C) used as a continuous predictor.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Read more here: Source link