Abstract

   Thermal tolerance of an organism depends on both the ability to
   dynamically adjust to a thermal stress and preparatory developmental
   processes that enhance thermal resistance. However, the extent to which
   standing genetic variation in thermal tolerance alleles influence
   dynamic stress responses vs. preparatory processes is unknown. Here,
   using the model species Drosophila melanogaster, we used a combination
   of Genome Wide Association mapping (GWAS) and transcriptomic profiling
   to characterize whether genes associated with thermal tolerance are
   primarily involved in dynamic stress responses or preparatory processes
   that influence physiological condition at the time of thermal stress.
   To test our hypotheses, we measured the critical thermal minimum
   (CT[min]) and critical thermal maximum (CT[max]) of 100 lines of the
   Drosophila Genetic Reference Panel (DGRP) and used GWAS to identify
   loci that explain variation in thermal limits. We observed greater
   variation in lower thermal limits, with CT[min] ranging from 1.81 to
   8.60°C, while CT[max] ranged from 38.74 to 40.64°C. We identified 151
   and 99 distinct genes associated with CT[min] and CT[max],
   respectively, and there was strong support that these genes are
   involved in both dynamic responses to thermal stress and preparatory
   processes that increase thermal resistance. Many of the genes
   identified by GWAS were involved in the direct transcriptional response
   to thermal stress (72/151 for cold; 59/99 for heat), and overall GWAS
   candidates were more likely to be differentially expressed than other
   genes. Further, several GWAS candidates were regulatory genes that may
   participate in the regulation of stress responses, and gene ontologies
   related to development and morphogenesis were enriched, suggesting many
   of these genes influence thermal tolerance through effects on
   development and physiological status. Overall, our results suggest that
   thermal tolerance alleles can influence both dynamic plastic responses
   to thermal stress and preparatory processes that improve thermal
   resistance. These results also have utility for directly comparing GWAS
   and transcriptomic approaches for identifying candidate genes
   associated with thermal tolerance.

   Keywords: thermal limit, CT[min], CT[max], heat shock, cold shock,
   genomics, transcriptomics

Introduction

   Temperature directly affects performance, survival, fitness, and
   consequently, the geographic distribution of organisms ([51]Angilletta,
   2009; [52]Dowd et al., 2015). Ectotherms are particularly vulnerable to
   changes in temperature, and these organisms have evolved a suite of
   adaptations to cope with thermal variability. An ectotherm’s thermal
   tolerance is determined by both fixed genetic factors and plastic
   changes in behavior, morphology, physiology, and gene expression.
   Genetic variation in thermal tolerance is well-documented (e.g.,
   [53]Sørensen et al., 2001; [54]McMillan et al., 2005; [55]Rako et al.,
   2007) and can occur through changes in basal stress tolerance and/or
   changes in the ability to quickly respond to thermal challenges
   ([56]Ayrinhac et al., 2004). These heritable differences within
   populations permit evolutionary shifts in thermal response as selection
   acts ([57]Hoffmann et al., 2003), and adaptive differences in thermal
   tolerance across latitudinal gradients and thermal environments are
   common ([58]Hoffmann et al., 2002; [59]Fallis et al., 2014).
   Specifically, populations from higher latitudes often are more tolerant
   of low temperatures than populations from lower latitudes, and the same
   pattern is also seen for heat stress, where populations that extend to
   lower latitudes often have improved survival at high temperatures
   (e.g., [60]Calosi et al., 2010; but see [61]Castañeda et al., 2015).
   Thus, thermal tolerance is a trait that is both highly plastic and
   highly adaptable, and understanding the genetic basis of thermal
   tolerance is critical for predicting future responses to environmental
   change.

   Genome Wide Association Studies (GWAS) and other quantitative genetic
   approaches have characterized the genetic architecture of thermal
   tolerance and identified genes that regulate temperature-dependent
   traits (e.g., [62]Morgan and Mackay, 2006; [63]Rako et al., 2007;
   [64]Svetec et al., 2011; [65]Rohde et al., 2016). A series of recent
   studies with the Drosophila Genetic Reference Panel (DGRP; [66]Mackay
   et al., 2012) have identified a number of candidate loci associated
   with thermal tolerance. [67]Rolandi et al. (2018) found 12 SNPs
   associated with variation in critical thermal maximum (CT[max]), and
   most of these SNPs were located within intronic regions, suggesting
   that variation in the heat stress response could be mediated by
   regulatory changes in gene expression or splicing. For cold, distinct
   but related traits often have non-overlapping genetic architectures,
   suggesting these traits have the capacity to evolve independently. For
   example, two plastic responses to cold, rapid cold hardening and
   developmental cold acclimation, have non-overlapping SNPs associated
   with them, although the genes associated with these traits share some
   functional similarities ([68]Gerken et al., 2015). Similarly, [69]Teets
   and Hahn (2018) found minimal overlap in genes associated with cold
   shock response and chill coma recovery, and [70]Freda et al. (2017)
   found no overlap in genes associated with adult and larval cold
   hardiness. The candidate genes identified in [71]Teets and Hahn (2018)
   were functionally tested with RNAi, and knockdown of most genes
   affected cold tolerance, indicating that GWAS is a robust method for
   identifying genes with functional roles in thermal tolerance. Taken
   together, the various GWAS studies of thermal traits indicate that the
   thermal stress response is a highly polygenic trait, but additional
   studies linking these polymorphisms to their functional consequences
   are needed to clarify their role in thermal tolerance.

   One way the genetic makeup of an organism influences thermal tolerance
   is by modifying gene expression changes in response to temperature
   change ([72]Stucki et al., 2017). Transcriptional responses to thermal
   variability have been described at various levels, including whole
   transcriptomic studies of specific stress treatments (e.g., [73]Qin et
   al., 2005; [74]Sørensen et al., 2005, [75]2016; [76]Teets et al.,
   2012), targeted experiments for specific candidate genes (e.g., Frost
   in [77]Goto, 2001; [78]Sinclair et al., 2007; [79]Zhu et al., 2017),
   and comparisons of transcriptomic responses to thermal stressors both
   among (e.g., in damselflies; [80]Lancaster et al., 2016) and within
   populations (e.g., [81]Telonis-Scott et al., 2009). A consistent theme
   from these studies is that changing temperatures can cause substantial
   changes in gene expression. For example, in D. melanogaster,
   acclimation that enhanced the cold response led to nearly one third of
   the transcriptome being differentially regulated (with around 60% of
   these genes being downregulated; [82]MacMillan et al., 2016). This cold
   acclimation included upregulation of genes already known to have an
   association with stress and temperature responses, such as Frost and
   many genes encoding for heat shock proteins. Similar sets of genes are
   also upregulated following brief cold shock in D. melanogaster and the
   flesh fly Sarcophaga bullata ([83]Qin et al., 2005; [84]Teets et al.,
   2012), indicating that anticipatory acclimation responses share some
   mechanisms with dynamic responses that occur during and after stress.
   For heat stress, most genes that are differentially expressed following
   short-term heat hardening ([85]Sørensen et al., 2005) and heat shock
   ([86]Telonis-Scott et al., 2013) in D. melanogaster are downregulated,
   with the exception of heat shock proteins, which are generally
   upregulated. However, despite the large body of literature on
   transcriptional responses to thermal stress, additional work is needed
   to clarify the functional consequences of these transcriptomic changes
   and determine how segregating variation in thermal tolerance relates to
   these transcriptional mechanisms.

   Thermal tolerance is a combination of dynamic plastic changes that
   occur during and after a stress event (i.e., processes that actively
   counter, repair or minimize the consequences of damage) and preparative
   processes that enhance stress resistance (i.e., processes that prevent
   damage; [87]Roy and Kirchner, 2000; [88]Wos and Willi, 2015). Plastic
   processes that occur during and after thermal stress largely involve
   production of stress proteins (e.g., heat shock proteins), often at the
   expense of other biological processes ([89]Feder and Hofmann, 1999).
   Preparative processes that enhance thermal resistance include
   production of protective osmolytes (e.g., cryoprotectants; [90]Yancey,
   2005; [91]Storey and Storey, 2012), changes in membrane composition and
   cell structure that permit membrane function at extreme temperature
   ([92]Sinensky, 1974; [93]Koštál, 2010), and anticipatory production of
   stress proteins during dormancy and/or thermal acclimation
   ([94]Manjunatha et al., 2010; [95]Colinet et al., 2013; [96]MacMillan
   et al., 2016). Thus, an allelic variant may contribute to basal
   tolerance to extreme temperature by altering either of these two
   components: enhancing the plastic ability to adjust to stress by
   participating in or regulating the dynamic temperature response, or by
   better preparing the organism for that stress. However, for genes
   associated with variation in thermal tolerance, whether these genes
   primarily affect dynamic plastic processes or preparative processes is
   unclear.

   Here, we used a combination of GWAS and RNA-seq to address the extent
   to which genes associated with thermal tolerance variation are involved
   in the dynamic stress response and preparative responses. We measured
   critical thermal minimum (CT[min]) and CT[max] ([97]Schou et al., 2017)
   in 100 lines from the DGRP, and the resulting phenotypic data were used
   in conjunction with genome-wide polymorphism data to identify genes
   associated with variation in thermal limits. These candidate genes were
   then compared to differentially expressed genes identified via
   transcriptomic assays of a single genotype exposed to heat and cold
   shock treatments to identify their roles in the stress response. Three
   non-mutually exclusive hypotheses were considered. To identify
   candidates involved in dynamic stress responses, we tested the
   following two specific hypotheses: (H1) Genes associated with thermal
   tolerance are part of the dynamic response, and are directly up- or
   downregulated during thermal stress; (H2) Genes associated with thermal
   tolerance are transcription factors and regulatory genes that regulate
   the dynamic transcriptional response to thermal stress. In support of
   H1, we predict that GWAS candidates will be more likely to be up- or
   downregulated in response to thermal stress, and these candidates will
   include genes directly activated during the stress response (e.g., heat
   shock proteins) as well as genes downregulated because they are
   incompatible with stressful temperatures (e.g., certain metabolic
   processes and reproduction). To identify genes involved in preparative
   processes that enhance thermal stress resistance, we tested the
   following hypothesis: (H3) Genes associated with thermal tolerance are
   involved in preparatory developmental and physiological processes that
   influence the condition of the organism at the time of thermal stress.
   Here we predict that specific GWAS candidates will be involved in the
   stress resistance processes discussed above (e.g., cell membrane
   remodeling, osmolyte production, etc.), but these candidates will not
   necessarily be part of the dynamic stress response. Our results
   indicate that genes associated with the thermal response have diverse
   functional roles that contribute to thermal tolerance in all three of
   these ways. There is considerable overlap between the genes associated
   with quantitative variation in thermal tolerance and those that are
   differentially expressed in response to thermal stress, and our GWAS
   analysis indicated an abundance of genes involved in developmental
   processes and cell morphogenesis that may have a role in enhancing
   stress resistance. Testing these hypotheses will advance our
   understanding of the functional consequences of genes polymorphisms
   associated with thermal tolerance. Furthermore, these results also have
   utility for directly comparing two commonly used methods for
   identifying and characterizing candidate genes associated with thermal
   tolerance.

Materials and Methods

Insect Rearing

   The Drosophila Genetic Reference Panel (DGRP) was established from a
   natural population in Raleigh, North Carolina, and isofemale lines were
   isogenized with 20 generations of full-sib mating ([98]Mackay et al.,
   2012). DGRP lines were obtained from the Bloomington Drosophila Stock
   Center, maintained at 25°C on a 12:12 light–dark cycle, and fed a
   standard cornmeal/soy flour diet consisting of 0.58% agar, 1.73% yeast,
   7.31% cornmeal, 1.00% soy flour, 0.13% Tegosept (w/v), 7.69% light corn
   syrup, and 0.48% propionic acid (v/v) in H[2]O. To generate flies for
   CT[min] and CT[max] assays, 15 females and 10 males were added to vials
   containing food and dry active yeast and were allowed to mate and lay
   eggs for 4 days. Restricting the number of adults in each vial and
   limiting the time to lay eggs prevented vials from becoming
   overcrowded, as extremely high larval densities can impact thermal
   tolerance ([99]Sørensen and Loeschcke, 2001). Ten days after removing
   the parental adults, adults of the resulting progeny were removed and
   held for 24 h to ensure that all flies had an opportunity to mate.
   After 24 h, males and females were sorted and placed into separate
   vials in groups of 20. Flies were held in the vials for 3–4 days prior
   to measuring CT[min]. For CT[max], flies were held for 2–3 days, and 24
   h prior to the experiment, flies were lightly anesthetized with CO[2]
   and individually transferred to small screw-top vials with food. All
   flies were between 4 and 9 days old at the time of the experiment.

   Seven-day-old D. melanogaster Canton-S female flies were used to
   characterize gene expression responses after a cold or heat shock. Only
   females were used to minimize confounding variation due to sex
   differences in gene expression and to include gene expression responses
   associated with protection of egg production, which is strongly related
   to fitness and may be expected to be under selection in nature. We
   selected the Canton-S background for these experiments to (1) address
   the extent to which GWAS candidates predict gene expression in a
   standard, well-characterized genetic background, to increase the
   generalizability of our results, and (2) provide candidate genes for
   future functional experiments, as most mutant and transgenic strains
   are in the Canton-S background. To generate flies of known age for
   RNA-seq assays, all adults were removed from mixed-sex stock vials that
   had been maintained at approximately 50 flies per vial and all newly
   eclosed adults were sampled daily. Same-day cohorts were maintained in
   mixed-sex vials at a density of ∼30 flies/vial on a nutrient-rich
   medium, consisting of 0.88% agar, 8.33% yeast, 10% cornmeal, 0.33%
   Tegosept (w/v), 4.66% molasses, and 0.66% propionic acid (v/v) in
   dH[2]O ([100]Buchanan et al., 2018), at 25°C on a 12:12 light–dark
   cycle. Cohorts were transferred to fresh food vials after 4 days.

Phenotypic Assays

   To measure CT[min] and CT[max], we used a dynamic ramping approach in
   which flies were gradually cooled or heated until motor function was
   lost. To assess CT[min], we used a vertical jacketed column (modified
   from [101]Huey et al., 1992) connected to a temperature-controlled
   fluid bath, and the temperature was monitored inside the column with a
   type T thermocouple ([102]Supplementary Figure 1A). For detailed
   assembly instructions for the jacketed column, see [103]Awde et al.
   (2020). For each line, ∼20 males and ∼20 females were combined in the
   column and submitted to the following thermal program: 25°C for 5 min,
   25°C to 10°C at 0.5°C min^–1, 10°C for 2 min, then 10°C to −10°C at
   0.25°C min^–1. The ramping rates are in line with other studies of
   CT[min] and were designed to maximize throughput and prevent cold
   hardening during the procedure (e.g., [104]Sinclair et al., 2015). At
   10°C we began collecting flies as they reached their CT[min] and fell
   through the column into collection vials containing 70% ethanol. New
   vials were placed under the column at 0.25°C intervals as the
   temperature decreased. Flies were typically at the top or on the walls
   of the column at the beginning of a trial (since they are negatively
   geotropic), and any flies that remained at the bottom of the column
   were discarded once the temperature reached 10°C. Flies from each vial
   were then sexed and counted, and the CT[min] was recorded as the
   maximum temperature for a given interval (e.g., flies collected between
   10°C and 9.75°C had a CT[min] of 10°C). CT[min] for each line was
   estimated by averaging the CT[min] of all flies tested across two
   independent cohorts. Due to variation in line productivity, escaping
   flies, and discarded flies, the total number of flies measured per line
   ranged from 15 to 44 for males (median = 28 and mode = 28) and from 9
   to 38 for females (median = 26 and mode = 26).

   CT[max] was assessed using the same apparatus as CT[min], except the
   jacketed column was arranged horizontally and flies were contained
   individually in 2 ml screw-top vials to prevent them from voluntarily
   walking out of the column as temperature increased ([105]Supplementary
   Figure 1B). For each line, ∼18 males and ∼18 females were individually
   placed in vials attached to a wooden dowel ([106]Supplementary Figure
   1B). The wooden dowel with the vials was placed inside the column and
   submitted to the following ramping program: 25°C for 5 min, 25°C to
   35°C at 0.5°C min^–1, then 35°C to 45°C at 0.25°C min^–1. Flies were
   checked for movement after the temperature reached 35°C by flicking the
   wooden dowel every 0.2°C. The CT[max] of flies was recorded when flies
   were motionless and no longer responded to stimulus. As with CT[min],
   CT[max] for each line was estimated by averaging the CT[max] of all
   flies tested across two independent cohorts. Due to variation in line
   productivity and escaping flies, the total number of flies per line
   ranged from 23 to 53 for males (median = 33 and mode = 33) and from 24
   to 51 for females (median = 33 and mode = 32).

   We also tested the extent to which CT[min] and CT[max] were correlated
   with other life-history parameters and other measures of thermal
   tolerance using previously collected phenotype data for the DGRP. We
   obtained data for lifespan and fecundity from [107]Durham et al.
   (2014), Wolbachia infection status and chill coma recovery time from
   [108]Mackay et al. (2012), rapid cold hardening and chronic and acute
   survival from cold from [109]Gerken et al. (2015), CT[min] from
   [110]Ørsted et al. (2018), cumulative cold tolerance from [111]Teets
   and Hahn (2018), heat knockdown from [112]Rohde et al. (2016), CT[max]
   from [113]Rolandi et al. (2018), and cold and heat hardness from
   [114]Freda et al. (2019). We used Pearson correlations (cor.test) to
   test for linear correlation between these measures in R (version 3.6.1;
   [115]R Core Team, 2019).

Heritability and Genome Wide Association Study (GWAS)

   Broad sense heritability (H^2), defined as the proportion of the total
   phenotypic variation that is due to all genetic factors, was estimated
   as H^2 = σ^2[L]/(σ^2[L] + σ^2[ε]), where σ^2[L] is among-line and
   σ^2[ε]is within-line variance components ([116]Mackay and Huang, 2018).
   Variance components were estimated using a linear mixed model and
   treating line as a random effect, with the lme4 package ([117]Bates et
   al., 2015) in R.

   Genome wide associative mapping was used to identify genetic
   polymorphisms associated with CT[min] and CT[max] using the GWAS
   platform available on the DGRP website^[118]1 ([119]Mackay et al.,
   2012). This analysis associates the phenotypic variation of DGRP lines
   with single nucleotide polymorphisms (SNPs), insertions, deletions, and
   multiple nucleotide polymorphisms (MNPs). One of the lines tested (208)
   was removed from the GWAS analysis by the DGRP server, thus the GWAS
   analysis included 99 lines. Variants with p-value < 1E-4 (using the
   average mixed p-value of the two sexes) were considered significant and
   were annotated to genes using FlyBase annotation v5.57. The average
   mixed p-value of the two sexes for GWAS analysis was chosen because
   both CT[min] and CT[max] were significantly correlated across sexes
   (see the section “Results”). To identify transcriptional regulators of
   thermal stress in support of H2 (see the section “Introduction”), we
   compared GWAS candidates to annotated transcription factors in FlyBase
   annotation v5.57.

   As an alternative to GWAS on individual variants, we also conducted
   gene-based GWAS to test the aggregated effect of a set of SNPs (e.g.,
   SNPs within a gene) on CT[min] and CT[max] phenotypes. Gene-based
   p-values were calculated by contrasting the observed T value to an
   empirical distribution generated from resampling under the null
   hypothesis with permutations using PLINK (version 1.9; [120]Purcell et
   al., 2007; [121]Liu et al., 2010). We controlled for confounding
   genetic relatedness between the DGRP lines used in this study, and used
   the Tracy-Windom test in the AssocTests package ([122]Wang et al.,
   2015) to evaluate eigenvalues from 20 principal components (PCs) of
   genotype. We retained the first eight PCs as covariates in the PLINK
   model as described above ([123]Patterson et al., 2006). As inversions
   and Wolbachia infection status can also influence the phenotypes of the
   DGRP lines, we used the adjusted phenotypes for these factors outputted
   from the DGRP2 website. Variants with MAF ≥ 5% and a genotype rate of
   20% were used as well as the FlyBase v5.57 gene annotations. A total of
   1,939,313 variants were tested over 8,954 and 8,270 genes with at least
   one significant variant for CT[min] and CT[max], respectively.
   Genome-wide significance was determined by controlling for FDR using
   the q value method ([124]Storey and Tibshirani, 2003).

RNA-Sequencing and Differential Gene Expression

   To characterize gene expression responses to thermal shock, three
   replicates of three females each were exposed to cold or heat shock
   conditions by placing flies in sealed 15 × 150 mm glass test tubes and
   submerging in a circulating water bath programmed to cool or heat at a
   rate of 0.25°C min^–1 until the temperature reached 4°C or 37°C,
   respectively. Flies were held at the final temperature for 5 min and
   then collected into pre-filled bead homogenization tubes (Benchmark
   Scientific) under CO[2] anesthesia, immediately snap-frozen in liquid
   nitrogen and held at −80°C. Control flies were similarly handled but
   remained at 25°C until collection and flash-freezing. Whole flies were
   homogenized using a Bullet Blender Bead Homogenizer (Next Advance) in
   300 μL TRIzol Reagent (Life Technologies) followed by purification with
   a Direct-zol RNA MicroPrep Kit (Zymo Research #R2060) according to the
   manufacturer’s instructions. DNA was removed using DNAse I on the
   column followed by two washes with RNA Wash Buffer. Total RNA was
   eluted with 15 μL DNAse/RNAse-free water, quantified using a Qubit 2.0
   Fluorometer (Thermo Fisher), and RNA quality was assessed using an
   Agilent 2100 Bioanalyzer (Agilent Technologies). RNA-seq libraries were
   generated from 1 ug of rRNA-depleted total RNA using the NuQuant
   Universal RNA-Seq Library Preparation Kit (Nugen #M01506) according to
   the manufacturer’s protocol with 12 cycles of PCR. A total of nine
   libraries were pooled and sequenced on a single lane of an Illumina
   HiSeq1500 single-read flow cell. The sequence quality of the resulting
   raw Illumina reads was assessed using FastQC (version 0.11.4) and reads
   were aligned to the D. melanogaster reference genome (Release 6) using
   STAR aligner ([125]Dobin et al., 2013). Genes were quantified using
   featureCounts (part of the Rsubread package, version 2.0.0, [126]Liao
   et al., 2019) against the DM6 build. Differential expression was
   performed using the DESeq2 package (version 1.24.0; [127]Love et al.,
   2014) in R. All data have been deposited into the NCBI SRA database
   with accession Bioproject [128]PRJNA612361.

Pathway Enrichment Analysis

   Overrepresentation analysis (ORA) of significantly differentially
   expressed genes (Benjamini–Hochberg corrected p-value < 0.01,
   fold-change > 2) and the GWAS genes (p-value < 1E-4) was performed
   using WebGestalt ([129]Wang et al., 2017; minimum five genes per
   category, maximum 2,000 genes) and a false discovery rate cut-off of
   0.05. The results of the overrepresentation analysis were the primary
   means by which we identified GWAS candidates in support of H3 (see the
   section “Introduction”).

Integration of Differential Gene Expression With GWAS

   Genes associated with SNPs identified using GWAS were matched with
   corresponding genes in the expression data set. For cases where
   multiple genes were associated with a single SNP, each gene was
   included. The distribution of log-fold changes of the expression of all
   genes in the heat shock and cold shock versus control was compared to
   the fold-change distribution of the genes significantly associated with
   the corresponding thermal performance limit with Kolmogorov–Smirnov
   tests implemented in R. These analyses were the primary means by which
   we identified inducible genes in support of H1 (see the section
   “Introduction”).

Results

Genetic Variation in Thermal Tolerance

   In this study we measured thermal limits (CT[min] and CT[max]) in a
   subset of lines from the DGRP. CT[min] values across the DRGP lines
   varied considerably more than CT[max] ([130]Figures 1A,B). CT[min]
   ranged from 0.81 to 8.55°C for males and from 2.29 to 8.64°C for
   females, while CT[max] ranged from 38.63 to 40.72°C for males and from
   38.51 to 40.80°C for females ([131]Supplementary Table 1). The sex of
   the flies did not affect CT[min] (p = 0.39), but it did affect CT[max]
   (p < 0.001), with the males being slightly more heat tolerant than the
   females. However, the effect size was small (effect size: 0.13°C). The
   interaction of sex and line was also significant for both CT[min] and
   CT[max] (p < 0.01 and p < 0.001, respectively). Within each phenotype,
   values were significantly correlated across sexes (CT[min], r = 0.85,
   p-value < 0.001; CT[max], r = 0.52, p-value < 0.001; [132]Supplementary
   Figures 2A,B). The sex-averaged CT[min] and CT[max] values were not
   significantly correlated across lines (r = 0.06, p-value = 0.54;
   [133]Figure 1C).

FIGURE 1.

   FIGURE 1
   [134]Open in a new tab

   Distribution of thermal limits in the DGRP. Histograms of mean
   phenotypes across 100 DGRP lines for (A) CT[min] and (B) CT[max]. (C)
   Correlation between line-means of CT[min] and CT[max].

   We tested for trade-offs associated with thermal tolerance by comparing
   CT[min] and CT[max] with previously collected lifespan and fecundity
   data ([135]Durham et al., 2014); however, we found no evidence of
   trade-offs among these traits, as neither thermal tolerance measurement
   was correlated with longevity or fecundity ([136]Supplementary Table
   2). DGRP lines have variable infection status by Wolbachia pipientis, a
   ubiquitous endosymbiont in insects that can significantly modulate
   physiology ([137]Werren, 1997). Within the lines studied, we found no
   evidence that Wolbachia infection impacts CT[min] or CT[max] (p = 0.07
   and p = 0.19, respectively). We also performed correlations between our
   data and other measures of thermal tolerance, and we found no
   correlation among our results and chill coma recovery time ([138]Mackay
   et al., 2012), measures of thermal plasticity (i.e., rapid cold
   hardening and survival from cold; [139]Gerken et al., 2015), cumulative
   cold tolerance ([140]Teets and Hahn, 2018), and heat knockdown time
   ([141]Rohde et al., 2016; [142]Supplementary Table 2). [143]Ørsted et
   al. (2018) measured CT[min] in males reared at 26°C (using a ramping
   method of 0.1°C min^–1), and we found significant correlations among
   their CT[min] values with our CT[min] values for males and females (r =
   0.60, p-value < 0.01 and r = 0.61, p-value < 0.01, respectively;
   [144]Supplementary Table 2). As in our study, [145]Rolandi et al.
   (2018) also measured CT[max] using a ramping method (0.25°C min^–1). We
   found significant correlations between their CT[max] values for males
   with our CT[max] values for males and females (r = 0.59, p-value < 0.01
   and r = 0.55, p-value < 0.01, respectively; [146]Supplementary Table
   2).

Genetic Architecture of Thermal Limits

   The estimated broad sense heritability (H^2) for CT[max] was 0.29 and
   for CT[min] was 0.25.

   Using available genomic data for the DGRP, we identified genetic
   polymorphisms associated with CT[min] and CT[max]. For the 99 lines
   measured, more than 1.9 million variants were analyzed (mostly SNPs),
   and we found ∼550 unique allelic variants associated with these traits
   (p-value threshold of 1E-4; [147]Supplementary Table 3). We identified
   348 allelic variants (319 SNPs, 15 deletions, 13 insertions, and 1 MNP)
   significantly associated with CT[min], and 193 allelic variants (173
   SNPs, 9 deletions, 8 insertions, and 3 MNPs) with CT[max] ([148]Table
   1). Polymorphisms associated with CT[min] and CT[max] were identified
   on all chromosomes ([149]Figures 2A,B). CT[min] and CT[max] did not
   share any allelic variants ([150]Figure 2C). Among these allelic
   variants, 262 mapped to 151 unique genes for CT[min], and 169 mapped to
   99 unique genes for CT[max]. Three genes (iab8, Btk29A and Sp1) were
   common between both traits ([151]Figure 2D). From all the genes
   associated with the allelic variants, 8% of the CT[min] and 9% of the
   CT[max] genes encode transcription factors ([152]Table 1), including
   one of the genes common to both traits (Sp1). Most of the allelic
   variants significantly associated with both traits were located in
   introns (55% for CT[min] and 71% for CT[max]; [153]Supplementary Table
   4). The distribution of the direction of effect sizes differed between
   SNPs that underlined CT[min] vs. SNPs that underlined CT[max]
   ([154]Figure 3A; Kolmogorov–Smirnov test, D = 0.87, p-value < 1E-9),
   such that the CT[min]-associated alleles that were most common in the
   DGRP population caused individuals to have higher CT[min] (i.e., worse
   cold tolerance; [155]Figure 3B), whereas the CT[max]-associated alleles
   that were most common in the population overwhelmingly caused
   individuals to have higher CT[max] (i.e., better heat tolerance;
   [156]Figure 3B). Additionally, there was a negative exponential
   relationship between effect size and minor allele frequency for both
   CT[min]- and CT[max]-associated SNPs ([157]Figure 3B).

TABLE 1.

   Overlap between the CT[min] and CT[max] SNPs and expression data.
   GWAS    Relaxed SNPs Strict SNPs Unique genes DEGs TFs DE TFs
   CT[min]     348          53          151       72  12    3
   CT[max]     193          21           99       59   9    9
   [158]Open in a new tab

   Counts of Relaxed SNPs (p-value < 1E-4), Strict SNPs (p-value < 1E-5),
   Unique genes (based on relaxed SNPs), differentially expressed genes
   (DEGs; FDR < 0.01), transcription factors (TFs), and differentially
   expressed transcription factors (DE TFs).

FIGURE 2.

   FIGURE 2
   [159]Open in a new tab

   Results of GWAS to identify polymorphisms associated with thermal
   tolerance. Manhattan plots of the results from the GWAS for (A) CT[min]
   and (B) CT[max]. The blue line corresponds to p-value < 1E-4 and the
   red line corresponds to p-value < 1E-5. The overlap of (C) allelic
   variants and (D) unique genes between CT[min] and CT[max] is also
   shown.

FIGURE 3.

   FIGURE 3
   [160]Open in a new tab

   Distribution of effect sizes for SNPs associated with CT[min] and
   CT[max]. (A) Histogram of the tolerance effect size scaled by the
   phenotypic standard deviation (σ[p]) for the same allelic variants. (B)
   Scatter plot of the minor allele frequency versus the tolerance effect
   size scaled by the phenotypic standard deviation (σ[p]) for the allelic
   variants associated with CT[min] and CT[max] (p-value < 1E-5).

   There is some evidence that low abundance alleles can be underpowered
   in the DGRP ([161]Ivanov et al., 2015), so as an alternative to the
   variant-based GWAS performed above, we also conducted gene-based GWAS.
   However, in the gene-based GWAS analysis, almost no genes were detected
   as significant using the p-value threshold of 1E-4 (three for CT[min]
   and one for CT[max]; [162]Supplementary Table 5). Thus, for the
   remainder of the paper, we will focus on the results of the
   variant-based GWAS described above.

Differential Gene Expression Following Acute Thermal Exposures

   RNA-seq and differential expression analysis were used to determine the
   gene expression responses of the Canton-S strain of D. melanogaster
   under ramped cold shock and heat shock conditions. In total, 15,844
   genes were expressed across all treatment groups. Principal component
   analysis (PCA) of expressed genes showed clustering of replicates for
   each condition ([163]Supplementary Figure 3). Pairwise comparison of
   cold shock and heat shock to controls revealed a large number of
   significantly differentially expressed genes (p-adj. < 0.05,
   fold-change > 2; [164]Figure 4). Among the differentially expressed
   genes, more were downregulated (5,126 in cold shock and 6,241 in heat
   shock) and fewer were upregulated relative to controls (1,826 in cold
   shock and 2,314 in heat shock). The direction of regulation for most
   differentially expressed genes was consistent across treatments, with
   6,081 genes changing similarly in both magnitude and direction in
   response to cold and heat shock conditions ([165]Figure 4).

FIGURE 4.

   FIGURE 4
   [166]Open in a new tab

   Scatter plot of the gene expression log[2] fold-change of heat-shock
   (37°C) and cold-shock (4°C) relative to the 25°C control. “Shared” DEGs
   indicate a significant expression relative fold-change in the same
   direction in both heat and cold shock. “Sig1 Shared” genes indicate
   that the relative fold-change is the same under heat shock and cold
   shock, but only significantly differentially expressed in one
   condition. “Unique” DEGs indicate that the genes were differentially
   expressed in each condition but in opposite directions. “Sig1 Unique”
   genes indicate that the genes are expressed in opposite directions
   relative to the control, but only significantly so in one condition.
   “NS” indicates genes that are not differentially expressed in both heat
   shock and cold shock. The density plots surrounding the figure indicate
   the density of the expression of genes that are associated with CT[min]
   (top; blue) or CT[max] (right; red) relative to the log[2] fold-change
   expression of all other genes (gray).

Integration of Transcriptomics and GWAS

   Of the 151 unique genes found for CT[min], 72 (47.7%) were
   differentially expressed under cold shock. Of the 99 unique genes
   associated with CT[max], 59 (59.6%) were differentially expressed under
   heat shock. The distribution of log-fold change expression of heat
   shock versus control for the GWAS candidates associated with CT[max]
   significantly differed from the distribution of all other genes under
   heat shock ([167]Figure 4; p-value < 0.01). For the genes associated
   with CT[min], there was a similar trend, but the log-fold change
   distribution of those genes only marginally differed from the
   background expression of all other genes ([168]Figure 4; p-value =
   0.062).

Overrepresentation Analysis

   Among the GWAS candidates, we identified six overrepresented GO
   biological process categories for CT[min] and 17 categories for CT[max]
   that met the FDR cut-off of 0.05 ([169]Figure 5). For both traits we
   found many enriched GO terms related to development, differentiation,
   and morphogenesis. Among the differentially expressed genes from the
   RNA-seq experiments, we identified 99 enriched GO biological process
   categories for cold shock and 113 enriched categories for heat shock
   ([170]Supplementary Table 6). Many of these categories were also
   related to development and differentiation. There was overlap between
   the enriched categories for GWAS and gene expression data. Five of the
   six categories enriched among genes associated with CT[min] were also
   significantly enriched among genes differentially expressed under cold
   shock, and 13 of the 17 categories identified for CT[max] were also
   enriched under heat shock ([171]Figure 5).

FIGURE 5.

   FIGURE 5
   [172]Open in a new tab

   Expression patterns of GWAS genes within significantly overrepresented
   GO biological process categories. All expression patterns are expressed
   relative to the 25°C control. Bolded categories indicate those also
   over-represented in the full DEG dataset for the corresponding
   temperature extreme. Bolded genes indicate significant differential
   gene expression at an FDR < 0.01.

Discussion

   Here, we characterized the genetic architecture of thermal tolerance
   and identified candidate genes that contribute to both dynamic
   responses to thermal stress and preparative processes that enhance
   stress resistance. Our study suggests that GWAS candidates are involved
   in both dynamic stress responses and preparative processes that
   influence the condition of the insect at the time of thermal stress.
   Together, our results indicate diverse functions for genes involved in
   thermal tolerance and allow us to generate new hypotheses for the
   genetic basis of thermal tolerance. Below, we discuss the genetic
   architecture of thermal tolerance in general, followed by a discussion
   of our three specific hypotheses to test the relative contribution of
   dynamic and preparative processes in shaping thermal tolerance.

Genetic Architecture of Thermal Tolerance

   While several studies have separately assessed the genetic basis of
   cold and heat tolerance, here we measured both CT[min] and CT[max]
   across 100 lines of the DGRP. We found variation in both measures,
   although CT[min] varied considerably more than CT[max] ([173]Figure 1).
   This pattern of variation in upper and lower thermal limits is also
   seen across species and populations with distinct geographic ranges,
   for both latitudinal and altitudinal gradients (e.g., [174]Gaston and
   Chown, 1999; [175]Addo-Bediako et al., 2000; [176]Chown, 2001;
   [177]Hoffmann et al., 2002; [178]Nyamukondiwa et al., 2011; [179]Sunday
   et al., 2011, [180]2019; [181]Kellermann et al., 2012a, [182]b). These
   patterns of variation in thermal limits, both within and across
   species, likely reflect stronger latitudinal and interannual variation
   in winter conditions relative to summer conditions ([183]Williams et
   al., 2015). In addition, our results are consistent with previous
   studies indicating that upper and lower limits have distinct underlying
   mechanisms (e.g., [184]Chown, 2001; [185]Nyamukondiwa et al., 2011), as
   we found no phenotypic correlation between CT[min] and CT[max] across
   lines.

   The variation in both CT[min] and CT[max] had a strong genetic
   component. Broad sense heritability was high for both CT[max] (H^2 =
   0.29) and CT[min] (H^2 = 0.25), which is consistent with previous
   heritability estimates for thermal responses in the DGRP. Our
   heritability estimate for CT[max] was higher than a previous estimate
   for a smaller subset of the DGRP population (H^2 = 0.14, [186]Rolandi
   et al., 2018). For CT[min], heritability was within the range observed
   by [187]Gerken et al. (2015) for acute and chronic cold exposure (H^2 =
   0.15 and 0.44, respectively). The strong heritability for both traits
   suggests high evolutionary capacity for thermal tolerance in the
   mid-latitude population from which the DGRP was derived.

   While variation in CT[min] and CT[max] was explained by distinct
   allelic variants, some variants mapped to the same genes, suggesting
   that some genes can affect both heat and cold tolerance. Of the 247
   total unique genes, three genes were common to both traits: iab-8 (a
   non-coding regulatory RNA), Btk29A (a tyrosine kinase involved in
   cellularization and morphogenesis), and Sp1 (a zinc finger
   transcription factor involved in ventral thoracic appendage
   specification, leg growth and in the development of type-II
   neuroblasts). At the gene level, there were also some candidate genes
   in common between our study and previous work. The gene CG42673 was
   associated with CT[max] in this study and also with chill coma recovery
   time in [188]Mackay et al. (2012). CG42673 is a putative nitric-oxide
   synthase binding protein, and while nitric oxide has not been linked to
   thermal tolerance in insects, nitric oxide is an important mediator of
   both heat and cold tolerance in plants ([189]Parankusam et al., 2017;
   [190]Costa-Broseta et al., 2018). In both our study and in [191]Ørsted
   et al. (2018), Mur89F was associated with CT[min], and this gene is
   involved in chitin metabolic process and extracellular matrix.

   Several patterns in our data suggest that thermal tolerance, especially
   heat tolerance, is an important component of organismal fitness in
   nature. Alleles that enhance heat tolerance (i.e., raise CT[max]) were
   more common in the DGRP population than alleles that impair heat
   tolerance ([192]Figure 3A). However, the pattern is opposite for cold
   tolerance; most alleles that improve cold tolerance (i.e., lower
   CT[min]) were relatively infrequent in the population ([193]Figure 3A).
   These results were counter to our expectations, since there is
   generally stronger latitudinal variation in cold tolerance than heat
   tolerance (see above). However, intraspecific latitudinal clines for
   heat tolerance have been observed in D. melanogaster ([194]Hoffmann et
   al., 2002; [195]Sgrò et al., 2010), indicating that there is selection
   for heat tolerance in this species. In the case of the DGRP, which
   originates in mid-latitude North Carolina, selection for cold tolerance
   may be lower than for heat tolerance, which could explain the relative
   rarity of alleles that improve cold tolerance. Alternatively, alleles
   that improve cold tolerance may have negative effects on other
   fitness-related traits, which would prevent these alleles from
   increasing in frequency in the population. Further, some polymorphisms
   in D. melanogaster oscillate in allele frequency across seasons
   ([196]Bergland et al., 2014), so depending on when the DGRP was
   collected (presumably summer, although exact details are not provided
   in [197]Mackay et al., 2012), cold tolerance alleles might be less
   common in the panel. Previously, a similar pattern was shown for
   oxidative stress resistance in the DGRP – i.e., most alleles that
   improve oxidative stress response are present at low frequency in the
   population ([198]Weber et al., 2012). Finally, our study did not
   address plasticity in thermal tolerance, which may be an important
   means of response to cold challenges. Thus, the rarity of alleles
   improving basal cold tolerance may not be relevant in this population
   if the population possesses alleles for thermal plasticity. However,
   regardless of the reason for these results, these allele frequency
   patterns between heat and cold tolerance are not random and suggest
   that different forces may be affecting patterns of standing genetic
   variation for the two traits despite their similar overall
   heritabilities.

   We also observed a negative exponential relationship between effect
   sizes and minor allele frequencies for SNPs that underlie both CT[min]
   and CT[max] ([199]Figure 3B). Overwhelmingly, the SNPs that have the
   greatest effect sizes in both directions are present at the lowest
   frequencies in the population. Many other studies have reported this
   same pattern for genetic polymorphisms that underlie a wide range of
   different traits in the DGRP, including oxidative stress resistance,
   startle response, starvation resistance, chill coma recovery time, and
   position effect variegation ([200]Mackay et al., 2012; [201]Weber et
   al., 2012; [202]Kelsey and Clark, 2017). Overall, this pattern suggests
   that large-effect alleles that underlie CT[min] or CT[max] have
   undergone selection, either to increase the frequency of large-effect
   alleles in the population by positive selection ([203]Barton and
   Keightley, 2002), thus driving the alternative allele to low frequency,
   or to remove large-effect alleles from the population by purifying
   selection ([204]Keightley and Lynch, 2003).

   Some of the genes underlying variation in CT[min] and CT[max] in the
   DGRP are also likely important for temperature adaptation in natural
   populations. Previous work by [205]Fabian et al. (2012) and
   [206]Bergland et al. (2016) used F[ST] outlier analyses to identify
   loci across the genome that are likely to be undergoing adaptive
   divergence among populations of D. melanogaster that inhabit different
   thermal environments in eastern North America and eastern Australia.
   Among the 109 candidate genes that showed convergent patterns of clinal
   latitudinal differentiation in North America and Australia ([207]Fabian
   et al., 2012; [208]Bergland et al., 2016), seven of these genes were
   associated with CT[max] and one was associated with CT[min] in the
   DGRP. The clinal genes associated with CT[max] were beat-VII, dpr8,
   CG33970, CG42322, Tsp66E, A2bp1, and Moe, and the sole clinal gene
   associated with CT[min] was fru. Most of these clinal CT[max] genes
   also changed expression following heat stress: three genes were
   downregulated (beat-VII, dpr8, and CG42322), two genes were upregulated
   (Tsp66E and Moe), and two genes (CG33970 and A2bp1) were not
   differentially expressed ([209]Supplementary Table 7). These genes are
   involved in a myriad of cellular and developmental processes, but a
   general theme is the potential role of neuronal processes in the
   thermal adaptation of heat tolerance (see below for additional
   discussion on the nervous system). It is also interesting to note the
   potential role of regulatory genes in thermal adaptation. Specifically,
   the CT[max] clinal gene A2bp1 is an RNA-binding Fox protein that
   regulates transcription and mRNA translation ([210]Usha and
   Shashidhara, 2010; [211]Carreira-Rosario et al., 2016), and the sole
   CT[min] clinal gene fru is a zinc finger transcription factor known to
   be a regulator of transcriptional activity of many genes across various
   tissues ([212]Sato and Yamamoto, 2020). While these genes did not
   dynamically respond to heat or cold stress ([213]Supplementary Table
   7), they may be important for setting up the developmental and/or
   cellular contexts in which stress responses operate.

H1: Genes Associated With Thermal Tolerance Are Involved in Dynamic Stress
Responses

   Thermal tolerance is shaped by a combination of preparative processes
   that improve stress resistance and dynamic changes in gene expression
   and activity that occur during and after stress. Dynamic changes in
   gene expression are well-established responses to thermal stress, and
   here we observed sweeping changes in gene expression in response to
   temperature change. Ramping at 0.25°C min^–1 toward both temperature
   extremes elicited transcriptome-wide gene expression changes, with
   43.9% and 54.0% of detected genes differentially expressed under cold
   and heat shock, respectively. These values are substantially larger
   than those reported in other studies, which range from minimal
   transcriptional response to ∼15% of the transcriptome, depending on the
   methodology used ([214]Zhou et al., 2012; [215]Telonis-Scott et al.,
   2013; [216]von Heckel et al., 2016; [217]Königer and Grath, 2018).

   By pairing GWAS and RNA-seq using similar ramping methodologies, we can
   assess the extent to which GWAS candidates are directly involved in
   dynamic stress responses. GWAS-associated CT[max] genes were
   significantly more temperature-responsive than the transcriptome at
   large and CT[min] genes were marginally so ([218]Figure 4), suggesting
   that genes associated with thermal tolerance are involved in the
   dynamic response to thermal stress. This congruence was also mirrored
   in the intersection of overrepresented GO terms in the two datasets,
   with five of six categories overrepresented in the CT[min] GWAS also
   enriched among genes of cold shock response, and 13 of 17 CT[max]
   categories shared with the heat shock response ([219]Figure 5).

   Despite a lower total number of genes identified that underlie CT[max,]
   the total set of overrepresented biological process categories was more
   diverse than for CT[min] and included cell signaling, muscle structure
   and function, and the sensory system ([220]Figure 5). This may indicate
   a stronger role of active defensive responses in setting CT[max]. At
   the individual gene level, the majority of top GWAS hits for CT[max]
   were thermally responsive ([221]Table 1) and were similarly diverse in
   function, including genes involved in neuropeptide signaling, mRNA
   processing, and protein dephosphorylation and catabolism. In contrast
   to the cold response, which included a mix of downregulated and
   upregulated genes, the majority of CT[max] GWAS hits found within
   thermally responsive categories were upregulated under heat ramp
   conditions, suggesting that they are involved in active protection from
   or in response to heat damage ([222]Figure 4). Thus, although the
   magnitude of phenotypic variation in CT[max] is substantially lower
   than that of CT[min], standing genetic variation may be mediated via a
   wider range of defensive mechanisms, each of small effect.

   Interestingly, the GWAS analysis did not indicate a role for the genes
   most commonly associated with thermal tolerance in experimental work.
   Much of the early literature on thermal limits focused on the effects
   of copy number and regulatory control of the heat shock protein (hsp)
   genes ([223]Welte et al., 1993; [224]Feder et al., 1996), and natural
   selection may also affect hsp allele frequencies (e.g., hsp70;
   [225]Bettencourt et al., 2002). However, neither hsp genes nor their
   regulatory factors (e.g., hsf-1) were identified from the GWAS analysis
   as causal drivers of variation for either heat or cold limits within
   the DGRP. These hsp genes and their regulatory factors (hsf) did
   increase in expression in response to both cold and heat shock, so
   while these canonical stress genes had clear roles in dynamic stress
   responses, polymorphisms in these genes were not associated with
   thermal limits in the DGRP. For cold tolerance, Frost (Fst) is commonly
   upregulated in response to cold acclimation and during recovery from
   cold stress ([226]Goto, 2001; [227]Qin et al., 2005; [228]Sinclair et
   al., 2007; [229]Colinet et al., 2010), and it is also located within a
   QTL for chill coma recovery time ([230]Morgan and Mackay, 2006). In our
   study, Frost was not associated with variation for thermal limits.
   Further, it was not differentially expressed following cold shock but
   was upregulated after the heat shock. The lack of upregulation during
   cold stress is likely due to flies being sampled at 4°C, as Frost
   expression typically only increases during recovery from cold stress
   ([231]Bing et al., 2012). The unexpected upregulation during a heat
   ramp suggests that Frost may be involved in heat stress, in addition to
   its role for cold and desiccation stress ([232]Sinclair et al., 2007).
   Likewise, Starvin (stv), a poorly studied gene that is strongly
   upregulated during recovery from cold stress ([233]Colinet and
   Hoffmann, 2010), was not associated with thermal tolerance but was
   upregulated after heat shock in this study.

H2: Genes Associated With Thermal Tolerance Have Regulatory Functions

   Changes in gene transcription are one of the primary cellular responses
   to cold and heat stress in this and other studies ([234]Gasch et al.,
   2000; [235]Leemans et al., 2000; [236]Gracey, 2007; [237]Lockwood et
   al., 2010; [238]Brown et al., 2014; [239]Sørensen et al., 2016).
   Moreover, there is a direct connection between whole-organism stress
   tolerance and the ability to transcriptionally respond to stress, as
   organisms with limited transcriptional stress responses have lower
   survival following exposure to stress ([240]Welte et al., 1993;
   [241]Feder et al., 1996; [242]Hofmann et al., 2000). Therefore, we
   expected to find candidate genes for thermal tolerance that have gene
   regulatory functions, such as transcription factors. Polymorphisms in
   regulatory genes or genomic regions may modify transcriptional
   responses to thermal stress, and thereby confer phenotypic differences
   in whole-organism thermal tolerance ([243]Zatsepina et al., 2001;
   [244]Bettencourt et al., 2002; [245]Lerman et al., 2003). We report
   evidence for this potential regulatory effect among the SNPs that
   underlie both heat and cold tolerance, but consistent with the stronger
   pattern of upregulation in the dynamic response to heat, our results
   suggest a larger role of transcription factors in driving genetically
   based variation in CT[max] than in CT[min].

   Overall, there were nine CT[max]-associated transcription factor genes
   ([246]Table 1), and all were differentially expressed in response to
   heat stress ([247]Supplementary Table 7). Indeed, the top four SNPs
   that were associated with CT[max] (lowest p-value; [248]Supplementary
   Table 3) were in two genes that encode transcription factors, Oaz and
   lola. Oaz encodes a transcription factor known to be involved in
   spiracle development ([249]Krattinger et al., 2007), and thus may
   mediate developmental mechanisms that impact heat tolerance, especially
   since spiracles facilitate gas exchange and failings of aerobic
   respiration may set upper thermal limits ([250]Dahlhoff and Somero,
   1993; [251]Pörtner, 2002). Oaz may also be involved in regulating the
   dynamic transcriptional response to heat stress, as it was the most
   strongly downregulated transcription factor following heat stress
   ([252]Supplementary Table 7). lola is involved in a diverse array of
   cellular and developmental processes ([253]Thurmond et al., 2019), and
   is represented among several GO categories in [254]Figure 5. All four
   of the top CT[max] SNPs lie in introns of the coding sequences of Oaz
   or lola, suggesting that these polymorphisms influence gene regulation
   ([255]Bicknell et al., 2012). Indeed, in the case of lola both
   mutations lie in a region that is a putative transcription factor
   binding site. While previous work also showed that these two genes
   respond to heat stress in D. melanogaster ([256]Brown et al., 2014), to
   our knowledge no previous studies have identified a functional role for
   these genes in heat tolerance. Another notable CT[max]-associated SNP
   lies in two overlapping genes that encode the transcription factors
   HmgD and HmgZ; the CT[max] SNP lies in the 5′ UTR intron of HmgZ and in
   the putative upstream regulatory region of HmgD. These genes encode
   proteins that belong to the family of high mobility group box
   transcription factors that are known to facilitate gene transcription
   by promoting DNA structural flexibility via chromatin remodeling
   ([257]Štros, 2010), and both of these genes were differentially
   expressed following heat stress ([258]Supplementary Table 7).
   Interestingly, high mobility group proteins have previously been
   reported to show expression patterns that track environmental
   temperature in killifish, Austrofundulus limnaeus ([259]Podrabsky and
   Somero, 2004). Thus, the regulation of gene transcription may be a key
   aspect of heat tolerance in D. melanogaster.

   The genetic architecture of cold tolerance also included genetic
   variation in transcription factor genes, but most of the
   CT[min]-associated transcription factor genes did not change expression
   following cold stress, suggesting a qualitative difference in the role
   of transcription factors in heat vs. cold tolerance. While there were
   12 transcription factor genes with significant associations with
   CT[min]([260]Table 1), only three of these genes changed expression
   following cold stress ([261]Supplementary Table 7). Importantly, one of
   the top SNPs in association with CT[min] lies in the gene blistered
   (bs) ([262]Supplementary Table 3). bs encodes a transcription factor
   known to be involved in a variety of developmental processes, including
   wing morphogenesis ([263]Dworkin and Gibson, 2006), tracheal
   development ([264]Affolter et al., 1994), and neural system development
   ([265]Donlea et al., 2009; [266]Thran et al., 2013). Similar to lola,
   the thermal tolerance SNP in bs lies in an intron with a transcription
   factor binding site; however, in contrast to lola and the other
   CT[max]-associated transcription factor genes, bs did not change
   expression in response to cold shock ([267]Supplementary Table 7).

   One of the three genes associated with both CT[min] and CT[max], the
   long non-coding RNA (lncRNA) iab-8, was downregulated in both cold and
   heat shock. LncRNAs have been implicated in a range of biological
   processes and are emerging as key regulators of gene expression at
   transcriptional and post-transcriptional levels ([268]Li et al., 2019).
   In vivo studies of lncRNAs revealed that dysregulated expression of
   lncRNAs in Drosophila may result in poor stress resistance
   ([269]Lakhotia et al., 2012). The iab-8 lncRNA, expressed in the
   embryonic abdominal segment 8, represses the expression of the abd-A
   gene in the posterior central nervous system ([270]Li et al., 2019).
   The abd-A gene is linked with neural system development ([271]Bello et
   al., 2003; [272]Cenci and Gould, 2005), and as discussed below, the
   nervous system likely plays an important role in thermal tolerance.

H3: Genes Involved in Thermal Tolerance Affect the Developmental and
Structural Context

   Thermal tolerance occurs within a developmental and structural context,
   making physiological systems more or less resistant to temperature
   challenges ([273]González-Tokman et al., 2020). Thus, genes associated
   with thermal tolerance may do so by altering the physiological
   condition of the organism at the time of thermal stress. Because these
   genes alter baseline preparedness prior to application of cold or heat,
   the expression of these genes may not directly respond to temperature
   changes. Moreover, we would expect their biological function to be
   concentrated in processes underlying thermal stability of physiological
   functions, such as the central and peripheral nervous system, cell
   membrane composition, and proteome stability ([274]Cossins and Prosser,
   1978; [275]Gu and Hilser, 2009; [276]Cooper et al., 2012; [277]Fields
   et al., 2015; [278]MacMillan et al., 2015a; [279]Willot et al., 2017).

   Our results suggest that although segregating variation in both heat
   and cold tolerance is likely to include some structural effects,
   preparative processes that enhance thermal stress resistance appear to
   play a stronger role for CT[min]. Fully half of the CT[min] GWAS genes
   did not change significantly in expression in response to cold
   exposure, regardless of the significance cutoff used ([280]Table 1).
   These included a cluster of functionally related genes (dar1, fru,
   NetA, RhoGEF64C, trio, twf) involved in nervous system development,
   which was reflected in overrepresentation of the nervous system and
   cell morphogenesis and differentiation GO categories in the CT[min]
   GWAS gene set. Neuronal failure operationally defines both CT[min] and
   CT[max] ([281]Andersen et al., 2018; [282]Andersen and Overgaard, 2019;
   [283]Jørgensen et al., 2019), and dynamic stabilization of the
   neuromuscular circuit under temperature stress is a likely mechanism
   for altering thermal limits. Indeed, previous investigation of the
   genetic architecture of cold hardiness and electrophysiological
   analyses of the rapid hardening response both suggest an important role
   for stabilization of ion channels and cytoskeletal structures
   supporting the synapse and neuromuscular junction ([284]Klose and
   Robertson, 2004; [285]Robertson and Money, 2012; [286]Freda et al.,
   2017). Aside from genes involved in neural morphogenesis, the GO term
   extracellular structure organization was overrepresented among cold
   tolerance genes, and this was the only overrepresented category that
   did not overlap with the differential expression categories. The
   glial-derived extracellular matrix is integrally involved in
   development, stabilization and plasticity of neuronal synapses, and is
   involved in promoting cell survival, facilitating repair, and
   maintaining synaptic current amplitude under stress conditions
   ([287]Dityatev et al., 2010; [288]Faissner et al., 2010; [289]Wang et
   al., 2018). Together, these results suggest that stabilization of
   nervous function is an important component of cold tolerance, which is
   consistent with recent physiological literature (reviewed by
   [290]Overgaard and MacMillan, 2017).

   For CT[max], more of the GWAS candidate genes were differentially
   expressed, especially when considering the strongest candidates
   ([291]Supplementary Table 7). Relatively few of the non-responsive
   genes were found within overrepresented categories, with a range of
   only 0–3 included in each of the 17 enriched CT[max] GO categories
   ([292]Figure 5). The few genes that consistently appeared in
   overrepresented categories, fra, fz2 and ptc, are also functionally
   associated with the nervous system, including axon and dendrite
   guidance and synapse organization. Spreading depolarization of the
   central nervous system (triggered by failure to maintain ion gradients
   between the intra- and extracellular compartments) is linked with heat
   tolerance across Drosophila species ([293]Jørgensen et al., 2019),
   indicating that neuronal failure is an important component of heat
   tolerance in addition to cold tolerance. An additional set of three
   genes, Pax, sls and Zasp66, co-occurred in several associated
   categories of muscle structure and development ([294]Figure 5).
   Depolarization of muscles has been associated with lower thermal limits
   (e.g., [295]MacMillan et al., 2014, [296]2015b) but to our knowledge
   has not been linked to upper limits.

Conclusion

   Several studies have separately assessed the genetic architecture and
   plastic transcriptional responses to thermal stress, but the extent to
   which genes associated with thermal tolerance are involved in
   preparative and dynamic stress responses has not been assessed. Here,
   we show that genes associated with variation in thermal tolerance,
   identified via GWAS, included differentially expressed genes directly
   involved in the dynamic stress response, as well as a number of
   transcription factors that likely regulate these processes. However,
   while GWAS candidate genes were more likely to be differentially
   expressed, genes commonly associated with thermal stress, such as heat
   shock proteins (hsp), were not identified among GWAS genes. These core
   stress genes tend to be highly conserved, so it is likely that these
   genes have little genetic variation, especially in functional regions.
   In addition, consistent with previous studies (e.g., [297]Sørensen et
   al., 2005; [298]Telonis-Scott et al., 2013; [299]MacMillan et al.,
   2016), most of the differentially expressed genes were downregulated
   for both hot and cold stresses. This result suggests an important role
   for shutting down certain biological processes during stress, and
   future studies should address these processes that are incompatible
   with stress tolerance, in addition to the well-studied protective
   pathways that are activated by stress.

   A noteworthy finding of our study is a prominent role for genes
   involved in preparatory physiological processes that influence the
   condition of the organism at the time of thermal stress. While we did
   not observe an abundance of genes involved in processes commonly
   associated with preparation for thermal stress, such as cell membrane,
   circadian function, and immune response ([300]Sinensky, 1974;
   [301]Koštál, 2010; [302]Teets and Hahn, 2018), we found a strong
   representation of nervous system processes for both CT[min] and CT[max]
   GWAS genes. Many of these genes have defined roles in development or
   morphogenesis, suggesting that developmental processes can influence
   thermal tolerance later in life, or that some of these genes are
   co-opted for thermal tolerance later in life. Future validation with
   other tools, such as RNA interference and transgenic overexpression,
   can clarify the precise role of these genes in thermal tolerance.

   While all three of our hypotheses were supported for both cold and heat
   tolerance, genes associated with upper thermal limits tended to be more
   involved in dynamic stress responses than those associated with cold
   tolerance. Furthermore, some genes associated with thermal tolerance
   appear to play multiple roles, highlighting that our three hypotheses
   are not mutually exclusive and that some genes likely have pleiotropic
   roles to shape thermal tolerance. Both upper and lower thermal limits
   had a strong genetic component, but the genetic signatures for these
   traits were largely distinct, with no overlapping SNPs and only three
   overlapping genes. Thus, heat and cold tolerance involve distinct
   molecular processes and can likely independently involve in response to
   changing environmental conditions.

Data Availability Statement

   All datasets generated for this study are included in the
   article/[303]Supplementary Material. Sequencing data have been
   deposited into the NCBI SRA database with accession Bioproject
   [304]PRJNA612361.

Author Contributions

   NT, BL, SC, JW, and HA acquired the financial support for the project
   leading to this publication. NT, BL, and SC conceived and designed the
   study. ML, DA, LU, NJ, MW, EM, BP, and KB performed the experiments and
   data collection. ML, DA, TO’L, SF, BL, NT, and SC analyzed the data.
   ML, DA, TO’L, SF, BL, NT, and SC wrote the manuscript. All authors
   approved the final version of the manuscript.

Conflict of Interest

   The authors declare that the research was conducted in the absence of
   any commercial or financial relationships that could be construed as a
   potential conflict of interest.

Acknowledgments