Thermal performance curves are sensitive to genotype and environment

Under all three rearing temperature conditions (16°, 23° and 28°), the TPCs for locomotion are unimodal with declines in performance toward low and high temperatures. For animals reared at cool to moderate temperatures, the TPC shows the classic shape of a steeper drop in performance toward high than low temperatures, but those reared at a high 28° actually tend to show sharper performance declines toward low temperatures (Figure 2A,B; Figure 3B,E; Supplementary Figure S2). Thermal performance curves also shift position and width depending on rearing temperature, with cooler rearing temperatures yielding wider TPCs that show more locomotory activity at lower temperatures (Figure 2A, Figure 3C,F; Supplementary Figure S2).

Despite rearing temperature yielding large and consistent differences in locomotion that give all genotypes a similar characteristic TPC shape (Ayrinhac et al. 2004; Hoffmann et al. 2005), the genetic identity of strains also contributes substantial variation in temperature-dependent locomotion profiles (Figure 2B). While strain-specific genotypic differences appear more pronounced than general differences among phylogeographic groups, visual inspection of TPCs for the genotypes from the Kerala “clade” strains exhibit cold-shifted TPCs compared to most other strains (Figure 2B, Figure 3F). By contrast, the two Montreal “clade” genotypes tend to perform poorly at colder temperatures and the Hubei genotype has the highest and one of the broadest TPCs (Figure 3).

To quantify these trends, we applied a 3-parameter logistic function (Equation 1) to data from each half of the TPC (the ‘cool’ half with subscript “C” and the ‘hot’ half with subscript “H”). From these function fits to data, we extracted the fit parameter values as simple and biologically interpretable summaries of locomotion, treating the reaction norm phenotype as a function-valued trait to compare TPCs across strain and rearing condition treatments.

First, we characterized baseline swimming behavior at benign temperatures with the asymptote parameter (α), which is related to maximum locomotion. We observed α to be highest for most strains when animals were raised at a benign 23° and α lowest when raised at 16° (ANOVA α_C F_2,66 = 18.48, P < 0.0001; α_H F_2,66 = 15.57, P < 0.0001; Tukey’s posthoc tests show differences among all rearing temperatures for both α_C and α_H) (Figure 3A,D; Supplementary Table S4). Animals reared at 23° appear visibly healthier than when reared at the extreme temperatures, and this may translate into them being capable of faster locomotion at initial assay temperatures. Alternatively, the thermal jump from rearing temperature to the initial assay temperature (25°) is largest for animals reared at 16°, which might shock them into slow initial locomotion as captured by α. We also observed significant genetic variation in both α_H and α_C (non-overlapping 95% CIs among genotypes; Figure 3A,D). The Hubei strain VX0034 shows consistently fast α across rearing temperatures, whereas the three strains from Nairobi consistently exhibit especially slow baseline locomotion when reared at 28° (Figure 3A,D). Strain genotypes also differ considerably in their degree of plasticity to rearing conditions, with some strains having nearly identical α_C values regardless of rearing temperature (overlapping 95% CI across rearing temperatures for e.g., VT847, NIC65, GXW23) and others being very sensitive (non-overlapping 95% CI across rearing temperatures for e.g., NIC21, JU1348, VX0034; Figure 3A). This variation in trait canalization also extends to the hot end of the TPC (e.g., compare low vs. high plasticity in α_H for VT847, QR24, NIC20 vs. ED3083, NIC21, VX0034; Figure 3D). These strain differences in plasticity point to genetic variation for the phenotypic robustness of baseline locomotion behavior in response to rearing conditions, in addition to heritable variation for the raw locomotory response.

Our function fit captures the sensitivity to changes in ambient temperature as the sharpness of the decline in locomotory activity, with this slope parameter (β) also influenced by rearing temperature and genotype (ANOVA β_C F_2,66 = 55.94, P < 0.0001; β_H F_2,66 = 3.75, P = 0.029; non-overlapping 95% CIs among genotypes) (Figure 3B,E). The magnitude of the slope (|β|) is steeper for most strains when the assay increments toward hot temperatures than when it approaches cool temperatures, at least when they were reared under cool or benign conditions (Figure 3B,E; Supplementary Table S4). For animals derived from these same cool or benign rearing conditions, we also observe more genetic variation in β when they experience the warming portion of the TPC compared to cooling half of the TPC (F-test for variances: 16° rearing F_22,22 = 0.137, P < 0.0001; 23° rearing F_22,22 = 0.097, P < 0.0001; 28° rearing F_22,22 = 0.828, P = 0.66). Moreover, the cooler the rearing temperature, the shallower the slope for the decline in locomotion toward low temperatures (Tukey posthoc tests indicates β_{C 16°} < β_{C 23°} < β_{C 28°}; Figure 3); but the cooler the rearing temperature, the steeper the decline in locomotion toward high temperatures (Tukey posthoc tests indicates β_{H 16°} ≤ β_{H 23°} ≤ β_{H 28°} with β_{H 16°} < β_{H 28°}; Figure 3B; Supplementary Table S4). Despite considerable overlap in slope estimates among strains for the three rearing temperatures, especially at the hotter portion of the TPC, we also observe heritable strain differences in β in response to rearing conditions (Figure 3B,E). For example, strain JU1338 exhibits virtually identical values of β_C regardless of the rearing temperature (similar for QX1547 and JU1348), whereas β_C values range threefold depending on rearing temperature for GXW23 (similar for HK104 and QR25). As we also described for α, these findings implicate heritable variation for phenotypic robustness of the dynamic behavioral response encapsulated by β as a function of rearing conditions. Even though both α and β derive from the same locomotory TPC, distinct genotypes can show robustness or plasticity separately for each trait, suggesting that these behavioral responses involve integration of the same environmental inputs with distinct genetic pathways.

The inflection point parameter (τ) identifies the temperature at which the TPC flips from concave to convex, heuristically indicating the center of the transition from normal to stopped locomotion. For all strains, low rearing temperatures led to colder values for the lower inflection point (τ_C; ANOVA F_2,66 = 215.75, P < 0.0001) (Figure 3C). However, the upper τ values (τ_H) are similar regardless of whether worms were reared at 23° or 28° (τ_H; Figure 3F; ANOVA F_2,66 = 39.29, P < 0.0001, Tukey’s posthoc test indicates no significant difference between 23° or 28°), and 5 strains actually have a colder value of τ_H when reared at 28° than 23° (non-overlapping 95% CI for JU1338, HK104, EG4181, JU1348, NIC20). This suggests that limits to tolerance for high ambient temperature might not permit further upward shifts of the TPC with higher rearing temperature, with the notable exception of the Nairobi genotypes due to their unusually low τ_H with 23° rearing (Figure 3F). Plasticity in τ as a function of rearing conditions is pervasive, with very few strains showing “robust” insensitivity to rearing temperatures (e.g., VT847 for τ_C, JU1908 and VX0034 for τ_H; Figure 3C,F). Temperate clade strains consistently show the greatest plasticity for τ_C in response to different rearing temperatures, with nearly a 4° shift among rearing treatments, whereas the Nairobi strains exhibit consistently high plasticity for τ_H (Figure 3C,F).

The TPCs from animals reared at cooler temperatures exhibit a ∼1° wider span between the two inflection points (Δτ = τ_H - τ_C), shifting position more on the cool end of the TPC and thus reflecting greater lability in response to cold ambient temperatures (ANOVA F_2,66 = 6.39, P = 0.0029, Tukey’s posthoc tests indicate Δτ_16° = Δτ_23° > Δτ_28°; Figure 4A-C; Supplementary Table S4). Likewise, most strains have a narrower TPC (smaller Δτ) when reared at 28° than when reared under cooler conditions (Figure 4C). The widest TPCs across the wild isolate strains are found consistently in those strains from the Temperate clade when reared at 16° or 23° (Figure 4A-C), which could be interpreted as a feature of a “generalist genotype.” Rearing temperature alters TPC midpoint even more dramatically than does genotype (ANOVA F_2,66 = 157.95, P < 0.0001, Tukey’s posthoc tests indicate significant differences among all rearing temperatures; Figure 4A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Correlations among behavioral reaction norm features as a function of rearing conditions and genotype. (A) C. briggsae thermal performance curve midpoint shifts warmer with higher rearing conditions (locomotion index LI1) (ANOVA F_2,66 = 157.95, P < 0.0001, all conditions differ by Tukey’s post-hoc test). The midpoint of the TPC also correlates positively with its width across wild isolate genotypes, for benign and hot rearing temperatures (Spearman’s ρ = 0.64, P = 0.001 for 28°C; ρ = 0.47, P = 0.024 for 23°C; ρ = 0.22, P = 0.3 for 16°C). (B) The temperature of peak performance (T_opt) correlates positively with TPC width across genotypes reared under benign or high temperature conditions (Spearman’s ρ=0.75, P < 0.0001 for 23°C; ρ=0.41, P = 0.051 for 28°C; ρ=-0.046, P = 0.84 for 16°C). (C) TPC widths (Δτ) are narrower on average for wild isolate strains reared at 28°C than at lower temperatures (ANOVA F_2,66 = 6.39, P = 0.0029, Tukey’s post-hoc test). (D) The temperature of peak T_opt is greatest for wild isolate strains reared at high temperature (ANOVA F_2,66 = 55.33, P < 0.0001, Tukey’s posthoc tests: T_{opt 28°} > T_{opt 23°} > T_{opt 16°}). (E) The peak performance for locomotory activity (LI1) P_max, however, occurs for wild isolate strains reared under benign intermediate rearing conditions (ANOVA F_2,66 = 20.25, P < 0.0001, Tukey’s posthoc tests: P_{max 23°} > P_{max 28°} > P_{max 16°}).

Much of the thermal performance curve literature emphasizes peak performance. Most strains reach their peak performance (P_max) when reared at 23° rather than the other two temperatures (ANOVA F_2,66 = 20.25, P < 0.0001, Tukey’s posthoc tests indicate P_{max 23°} > P_{max 28°} > P_{max 16°}) (Figure 4E); this is similar to the pattern for the asymptote parameter (α) with respect to rearing temperature. Just three strains have faster peak locomotion at 28° instead of 23° (ED3083, VT847, HK104; Supplementary Figure S3). To estimate the temperature for peak locomotion (T_opt), which the logistic function fit does not capture, we fit a Gaussian curve (Equation 2) to locomotion data between 18° and 32.5° where the asymmetry in TPC shape is least pronounced (Figure 2). The temperature of peak performance (T_opt) is warmest when we reared animals at 28° compared to 16° (ANOVA F_2,66 = 55.33, P < 0.0001, Tukey’s posthoc tests indicate T_{opt 28°} > T_{opt 23°} > T_{opt 16°}), although for many strains there is little difference between rearing at 23° vs. 28° (Figure 4D). This is similar to the pattern we see with inflection point values (τ) and lends support to the idea that rearing temperature acts to shift the entire curve. We observed no significant correlation between peak performance P_max and the temperature at which peak performance occurred (T_opt) (all P > 0.53). We also found no significant correlation between peak behavioral activity (P_max) and TPC breadth (Δτ) across the 23 wild isolate genotypes for any rearing temperature, despite a positive trend for animals reared under benign 23° conditions (Spearman’s ρ=0.38, P = 0.073). By contrast, we observed significantly higher temperatures of peak performance (T_opt) for genotypes that have broader temperature breadths (Δτ), provided the animals were not reared under cold conditions (Spearman’s ρ=0.75, P < 0.0001 for 23° rearing; ρ=0.41, P = 0.051 for 28° rearing; ρ=-0.046, P = 0.84 for 16° rearing; Figure 4B).

Robustness of thermal performance among individuals within a genotype

Our plots of TPCs exclude error bars to reduce visual clutter, but we also explicitly considered variation within strains and how this variation responds to assay temperature and to rearing temperature to learn more about individual thermal responses. Of course, the lowest standard deviation (SD) in locomotory activity among individuals occurs when all worms are paralyzed at the temperature extremes (Figure 2C; Supplementary Figure S2). The SD is highest near temperatures that are slightly more benign than the inflection points in the TPC, where locomotion just begins to decline (Figure 2C). However, high temperatures induce greater inter-individual variability for a given genotype as some worms stop moving entirely while others continue relatively normal locomotion, being most pronounced when we reared the animals at 28° (Figure 2C; Supplementary Figure S2), which may associate with a wider range of individual body condition. The lower overall SD toward the cold end of the TPC indicates that the decline in locomotion tends to be more consistent across individuals, with this more stereotyped reaction indicative of a more canalized behavioral response to cooling than to warming ambient temperatures.

Genetically distinct TPCs of Tropical and Temperate strains

We next sought to characterize differences between the “Tropical” AF16 and “Temperate” HK104 strains in detail, as the parental genotypes of a library of recombinant inbred lines with which we could explore the genetic architecture of trait variation. When reared at 23°, the locomotory behavior of these two strains differed most at high temperature near 36°, also reflected by the change in proportion of animals moving, despite similar locomotory activity of the two strains in benign and cool portions of the TPC (Figure 5A-C; Supplementary Figure S4). The LI7 locomotion index metric has greater resolution than LI1 when animals are moving slowly (Stegeman et al. 2019), and LI7 reveals that near 36° many individuals are moving slowly in AF16 whereas HK104 ceases locomotion altogether. This distinction in performance near 36° is consistent across four large repeated experiments (Supplementary Figure S4). An instantaneous jump in temperature from a benign 23° to 35.6° also showed AF16 and HK104 to respond distinctively: the movement for AF16 worms declined more quickly over time, but retained slow activity for the duration of the ∼21 min assay (Figure 5D). By contrast, locomotion of HK104 declined more slowly toward complete immobility after ∼13 min of exposure to 35.6° (Figure 5D), suggesting that AF16 can tolerate 35.6° longer than HK104 can. Because of these consistent and conspicuous differences in high temperature responses for which the Tropical AF16 strain outperforms the Temperate HK104 strain, we devoted further effort to this aspect of high temperature locomotory performance.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Thermal performance for Tropical strain AF16 and Temperate strain HK104. (A) Decrementing and incrementing temperature TPCs show decline of locomotory performance with more extreme temperatures with three distinct metrics: (A) locomotion index with 1/15^th second resolution (LI1), (B) locomotion index with 7/15^th second resolution (LI7), and (C) the proportion of moving individuals. Upper and lower halves of the TPC start at an initial 25°C assay temperature (vertical black line). (D) Thermal responses for AF16 and HK104 with increasing temperature steps (dashed lines, same data as in A) or with an instantaneous increase from 23°C to 35.6°C followed by a constant hold temperature (solid lines); top panel shows TPCs and bottom panel shows experimental thermal profiles). Error bars indicate ± 1 SEM (E) Comparison of locomotion index (LI7) TPC and (F) values at 35.6°C (± 1 SEM) for F1 heterozygotes to their homozygous parental genotypes (Tropical AF16 and Temperate HK104; purple = F1_HA with HK104 as maternal parent; pink = F1_AH with AF16 as maternal parent). (G) Comparison of proportion moving TPC for F1 and parental genotypes, with (H) point estimates at 35.6°C (± 1 SEM). (I-K) Parameter estimates from logistic function fits to LI7 TPCs for F1 heterozygotes and parental genotypes show more extreme values for baseline locomotion (α) in heterozygotes and intermediate F1 values for thermal sensitivity (β). All animals reared at 23°C; AF16 shown in red, HK104 show in blue.

The reciprocal F1 progeny of AF16 and HK104 exhibited largely intermediate TPC phenotypes at those high temperatures that most distinguish the parental strains (35.6-36.7°), with the F1s appearing slightly more similar to their maternal parent to indicate some parent-of-origin asymmetry (Figure 5E,F). These F1 data cannot distinguish mitochondrial vs. epigenetic parental effects, although we can exclude sex-linkage as a factor because all individuals were hermaphrodites (XX karyotype); QTL mapping with RILs found no evidence for a mitochondrial contribution (see below). The largely intermediate F1 heterozygote phenotypes imply that trait variation cannot be explained by a single dominant gene of large effect, although the inflection point parameter values are more ambiguous on this point (Figure 5F-K). These results also tell us that heterozygous genotypes do not exhibit outbreeding depression for locomotion in response to increasing temperature, in contrast to some other phenotypes in C. briggsae (Dolgin et al. 2007), with the asymptote parameter actually showing overdominant trait values that exceed both parents (Figure 5I). We next applied quantitative trait locus mapping using recombinant inbred lines to further dissect the genetic architecture of these TPC differences.

Locomotory thermal performance variation among recombinant genotypes for QTL mapping

We quantified locomotory thermal performance for 153 RILs, plus their two parental lines, with the liquid microdroplet assay for the ‘hot’ portion of the thermal performance curve. We then tested for associations between genotype and phenotype, yielding multiple QTL across the genome that each explain up to H² = 27.6% of the variation in behavioral thermal performance traits (Supplementary Table S5). We summarized trait variation in temperature-dependent locomotory responses using a diversity of phenotypic metrics and reduced their dimensionality through function-valued trait analysis of the reaction norms. Most RILs produced phenotypes that spanned the trait space between parental values, but those that fell beyond parental performance indicate transgressive segregation and imply multiple loci controlling trait variation (Figure 6A,B). For example, 52 RIL strains have steeper point estimates for slope (|β_H|) than the HK104 parent whereas only 7 RILs had a shallower estimated slope than AF16 (Figure 6D). Similarly, only 6 RILs showed higher locomotory activity than AF16 at 35.6° (LI7) and only 9 RILs have inflection point values (τ_H) higher than AF16 (Figure 6E); β_H and τ_H correlated only weakly across RIL strains (r=0.127 P = 0.0132), suggesting that the genetic basis to variation in these traits may be distinct. Many strains (83) have asymptote estimates lower than either parent, indicating that they have slower baseline locomotion (α_H; Figure 6C). Overall, the range of variation among the RILs is similar to (or greater than) that seen among wild isolate genotypes of C. briggsae: mixing up the genomes of just two genotypes to create the RILs produced a range of phenotypic variation similar to the natural range of variation seen in the entire species (Figure 6D).

• Download figure
• Open in new tab
• Download powerpoint

Figure 6

Variation among recombinant inbred lines for function-valued traits and synthetic “PC traits” of temperature-dependent locomotory behavior. (A) Thermal performance curves for 153 recombinant inbred lines (RILs) (gray) derived from C. briggsae parental genotypes (Tropical AF16, red; Temperate HK104, blue). Locomotion indices (LI7) for each strain at 35.6°C (dashed vertical line) is shown as a histogram of phenotypic values in (B). (C-E) Distribution of parameter estimates from function fits of locomotion index scores to each strain (LI7; see Equation 1 and Figure 1). Distribution of slope estimates (β_H) for 23 wild isolate strains (orange) in (D) is superimposed on the distribution for the 153 RILs (gray) for function fits to the LI7 metric from TPCs of worms reared at 23°C. Red and blue lines show values for strains AF16 and HK104. (F) Distributions of ‘synthetic’ locomotion behavior trait values for 153 RIL strains for the first eight principal components (PCs). Red lines indicate AF16, blue lines indicate HK104.

A function-valued trait approach provides one powerful way to distill a dynamic multivariate phenotype down to a simpler representation (Stinchcombe and Kirkpatrick 2012). Principal Component Analysis (PCA) represents another common approach in QTL analysis to reduce dimensionality, providing an interesting comparison to function valued traits alone. Therefore, we applied PCA to reduce the dimensionality of the 167 metrics for the 155 strains to yield 8 principal components that are statistically independent and together explain two-thirds of the total variation (Figure 6F; Supplementary Figure S5). Correlated phenotypes get weighted strongly together in a given principal component (PC) axis to represent an integrated “synthetic” phenotype, and we analyzed this subset of 8 PCs as composite traits for subsequent analysis. The first principal component (PC1) explains ∼20% of the phenotypic variability and relates primarily to locomotion at low to moderate temperatures up to 32.5°, where locomotion is normal or just starts to decline for many strains (Supplementary Figure S6). Subsequent orthogonal components explain smaller portions of phenotypic variation and become more dominated by a few phenotypes with high weights (Supplementary Figure S6). The variation in PC2 captures ∼16% of the variation among the RIL phenotypes, generally being associated with responses to high temperatures above 35.6° (Supplementary Figure S6). We also included measures of within-strain variation itself as “penetrance” phenotypes: the standard deviations for locomotion metrics and standard error for parameter estimates. These penetrance metrics tend to cluster together in PCs 3, 4, 6 and 7 that are distinct from the PCs loaded primarily with simple locomotion metrics (PCs 1, 2, 5; Supplementary Figure S6).

To perform QTL mapping, we integrated the RIL phenotype data with the RIL genotypes for 1031 SNP markers (Ross et al. 2011), which comprised 430 distinct genetic blocks in the genomes of the 153 RILs used here. We primarily focus our QTL mapping analysis and interpretation on the three TPC function fit parameters for the ‘hot’ portion of the thermal performance curve rather than PCs (slope β_H, asymptote α_H, inflection point τ_H). These analyses uncovered at least one QTL on chromosome V that explains H² =12–28% of the variation in locomotion at high temperature (depending on the specific phenotype) and another major QTL on chromosome II associated with locomotion at lower temperature (H² = 12–19%; Supplementary Table S5). The function fit parameters as phenotypes are especially interesting because they give a summary of characteristics of the entire thermal performance curve and correspond strongly to some of the less intuitive summaries of phenotype captured with PCs. Nevertheless, we also mapped the eight principal component axes as “synthetic traits,” using 2500 permutations to determine genome-wide significance thresholds based on a Bonferroni correction for multiple phenotypes. This reduction in data dimensionality allowed the maintenance of power to detect QTL while avoiding excessive multiple hypothesis testing and overly redundant phenotypes among the 167 metrics of temperature-dependent locomotory behavior. We then validated the biological interpretations of the “PC-QTL” with genomic associations involving a key univariate trait (LI7 at 35.6°).

Locomotion at benign temperature maps to QTL on Chromosome II

The asymptote parameter (α_H), which summarizes baseline movement behavior at benign temperatures, maps to the center of Chromosome II (Figure 7A). Both of the synthetic traits corresponding to locomotion at benign initial temperatures in the liquid micro-droplet assay, PC1 and PC6, similarly map to this region of the genome (Figure 8A,B). Other univariate phenotypes for locomotion at benign temperatures also had LOD peaks on Chromosome II in exploratory univariate analyses (Supplementary Figure S7; Supplementary Figure S8), many of which weigh heavily in PC1 and/or PC6 (Supplementary Figure S6). These QTL on Chromosome II may have a common genetic basis that affects all of these phenotypes. The HK104 allele at the marker nearest to the QTL peak is correlated with higher values of locomotion metrics (Figure 7C; Figure 8C), consistent with the parental strain differences (Figure 5I).

• Download figure
• Open in new tab
• Download powerpoint

Figure 7

QTL mapping for function-valued traits describing behavioral reaction norms. (A) QTL mapping for three function fit parameter phenotypes for locomotion index (normalized LI7). Dashed lines indicate threshold for significance from 2500 permutations of the data (Bonferroni corrected P = 0.05/12). (B) LOD peaks (solid lines) and 95% Bayes credible intervals (vertical dashed lines) for chromosomes with significant QTL for each function parameter phenotype; significance thresholds as in (A). (C) Phenotype values of RIL strains segregated by parental genotype at the marker closest to the QTL peak in (B) indicated above each panel. (D) Multiple imputation QTL mapping analysis of the Locomotion index (LI7) phenotype at 35.6°C, for all chromosomes and an expanded plot for the QTL on Chromosome V. LOD peak on Chromosome V at 37.35cM indicated by solid vertical line, with dashed vertical lines indicating 95% Bayes credible interval. Horizontal lines indicates the 5% genome-wide significance LOD threshold after 2500 permutations and Bonferroni correction as in Figure 8. Panel on right shows phenotype values of RIL strains segregated by parental genotype at the markers closest to the QTL peaks on Chromosome V (mean ± 1 SE for each marker allele shown adjacent to the distribution of points); values for parental strains shown as dashed horizontal lines (red AF16, blue HK104). Mean ± 1 SE across RILs for each marker allele shown adjacent to the distribution of points; dashed horizontal lines indicate parental strain phenotypes (red AF16, blue HK104). See Supplementary Table S5 for physical map positions of QTL.

• Download figure
• Open in new tab
• Download powerpoint

Figure 8

QTL mapping for synthetic “PC traits.” (A) QTL maps for the 8 ‘synthetic’ principal component locomotion behavior traits show logarithm of the odds (LOD) scores plotted against marker position by chromosome (“mt” = mitochondrial genome). Dashed lines indicate the 5% genome wide significance threshold from 2500 permutations and Bonferroni correction (P = 0.05/12). PCs with non-significant QTL indicated with gray colors in (A) and Figure 6F. (B) Significant QTL peaks shown for individual chromosomes (colored as for PC traits shown in A; significance thresholds as in A). Solid vertical lines indicates LOD peak for the marker loci shown in (C); dashed vertical lines indicate the 95% Bayes credible interval. (C) Phenotype values of RIL strains segregated by parental genotype at the marker closest to the QTL peak in (B) indicated above each panel. Parental strain values indicated by horizontal dashed lines (red AF16, blue HK104; PC colors as in Figure 6F).

High temperature locomotory activity maps to QTL on Chromosome V

The function-valued trait parameters for slope (β_H) and inflection point (τ_H) both mapped to QTL in the center of Chromosome V, indicating that this genomic region contains important genetic variation explaining locomotory differences at high temperatures (Figure 7A,B). Supporting this result, two synthetic PC traits mapped as QTL to the center of Chromosome V (PC2 and PC3), and PC3 also mapped to a QTL on the right end of the X Chromosome (Figure 8A,B). The QTL peaks for PC2 and PC3 on Chromosome V are 3.9Mb apart, with partially-overlapping 95% Bayes credible intervals, reminiscent of the “twin peaks” pattern of LOD scores for β_H (Figure 7A cf. Figure 8A). Univariate phenotypes related to locomotion at high temperatures (>35°) load most heavily on PC2, whereas PC3 is most closely tied to “penetrance” phenotypes related to variability among individuals at moderately high temperatures (<33°; Figure 9A, Supplementary Figure S6). Consequently, we mapped the univariate 35.6° locomotion index (LI7), also yielding a QTL peak on Chromosome V (Figure 7D), in fact the highest LOD value for any phenotype (Supplementary Table S5), and consistent with the biological interpretation of at least one QTL on Chromosome V associated with high temperature locomotory activity. Additional exploratory univariate QTL analyses reinforce this interpretation, given that 32 univariate traits related to high temperature locomotion associate with this region of Chromosome V even after Bonferroni multiple-test correction (Supplementary Figure S6; Supplementary Figure S9). Strains with Tropical AF16 genotypes at the peak locations of these QTL had values related to greater performance at high temperature, consistent with the difference between parental strains (Figure 8C; Figure 7C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 9

Univariate trait loadings for PC2 and PC3 that map to Chromosome V QTL and two-dimensional QTL scans. (A) Biplot of 167 individual phenotype weights measured for the 153 RILs, for the orthogonal PC2 and PC3 ‘synthetic’ trait axes that map to Chromosome V. Length and direction of vectors indicate weight of the phenotype in each PC with the text indicating individual phenotype metrics (Supplementary Figure S6). For example, the Locomotion index (LI7) at 35.6°C is weighted strongly in both these PC’s and can be found near the top-right corner (“LI7_35”). Circles indicate each RIL strain, colored according to their genotype for markers at QTL peaks for PC2 and PC3 on Chromosome V (marker cbv12146 for PC2 and marker cb49450 for PC3; red = both AF16 alleles, blue = both HK104 alleles, purple indicates different alleles between the peak of PC2 and PC3; peaks were 6cM and ∼3.9Mb apart, Supplementary Table S5). (B-F) Two-dimensional QTL scans show pairwise LOD scores for a given phenotype, shown only for chromosomes with significant QTL. The upper-left triangle above the 1:1 line in each plot contains the additive-effect LOD values after subtracting the LOD values from the single-QTL model. The lower-right triangle below the 1:1 line in each plot shows the LOD values for interactions. Note different color scales in each plot and each off-diagonal triangle zone. Magenta circles highlight the significant QTL locations (Supplementary Table S5).

We also found a QTL on the X chromosome for PC3 (Figure 8A,B). Exploratory analysis with univariate metrics suggest that this QTL mostly corresponds to phenotypes related to within-strain standard deviation for locomotion at relatively benign temperatures (Supplementary Figure S6; Supplementary Figure S7). One possibility is that genetic variation on the X chromosome contributes to the robustness of locomotory responses among individuals following exposure to moderately increasing temperatures.

Given the important role implicated by the Chromosome V QTL for behavioral performance at high temperatures, we sought to affirm a causative role for genetic differences between Tropical AF16 and Temperate HK104 in shaping locomotory thermal performance. Therefore, we created 13 NILs that had regions of HK104’s Chromosome V introgressed into the genomic background of AF16 by backcrossing RILs to AF16 for 7-12 generations using marker-assisted selection. We observed significant differences in locomotion among the NILs at 35.6°, most of which exhibited phenotypes more similar to the HK104 parent (Supplementary Figure S10). Consistent with QTL on Chromosome V, therefore, it can be sufficient for HK104 genotypes in the center of Chromosome V of the NILs to shift the phenotype toward HK104-like trait values (e.g., A1124-gs123 and A1124-gs124). Trait values for some NILs, however, suggest that the genetic architecture on Chromosome V might be more complex than the QTL analysis implies or that our NIL markers miss some introgressed segments of DNA (e.g., A1124-gs124 vs. A1124-gs66). Despite failing to narrow the QTL region, these introgression lines are consistent with the QTL region on Chromosome V containing genetic differences that cause phenotypic differences in behavioral thermal performance at high temperatures.

The QTL region for β and τ spans nearly 11Mb of Chromosome V and contains 2837 putative genes; the Bayes interval for locomotion (LI7) at 35.6° spans almost 7Mb and contains 1784 putative genes. Approximately 70% of these genes have associated gene ontology (GO) terms, implicating potential structural or molecular functions that we cross-referenced with genes for which we detected coding changes (missense, nonsense, or indel alleles). We identified 115 of the 2837 genes (63 of the 1784 genes) to have such coding changes as fixed differences between 20 Tropical and 4 Temperate group strain genomes (Thomas et al. 2015) (Supplementary Figure S11, Supplementary Table S7). Fixed differences between these phylogeographic strain groups should represent changes that are most ecologically relevant to temperature-dependent phenotypic divergence. Sensory and neurological system genes, including membrane proteins, receptors, and signaling pathway genes, were significantly overrepresented relative to genome-wide expectations (PANTHER overrepresentation test; Supplementary Figure S11) (Thomas et al. 2003). Despite the enrichment for these genes with potential influence on behavior, many other kinds of genes could conceivably contribute, in addition to the potential influence of non-coding regulatory or RNA genes on thermal performance.

Testing for multiple QTL on Chromosome V

The LOD peaks on Chromosome V visually appear bimodal for β_H and LI7 at 35.6°, and the broad Bayes credible interval suggests that multiple causal loci might occur under the QTL region (Figure 7A-C). The primary peak is located between 37.4 and 38.0 cM (near the QTL peak for PC3) and the secondary peak across the “valley” is between 25.0 and 31.5cM (near the QTL peak for PC2) (Supplementary Figure S9). Presuming that the two PCs capture distinct orthogonal aspects of the phenotypes, then the ‘twin peak’ pattern on Chromosome V might hint at distinct genetic contributions to temperature-sensitive traits. Can we find statistical evidence for twin peaks? It was not possible to map them separately with the multiple imputation method used for our primary analysis, which assumes only one peak per chromosome, prompting us to consider alternate approaches. Only 15 of the 153 RILs have genotypes that differ between the 37.4cM and 31.5cM locations, however, suggesting limited power to disentangle separate loci. Nevertheless, we explored the possibility of twin peaks using conditional QTL analysis, which did not reveal evidence of a significant pair of QTL peaks on Chromosome V (Supplementary Figure S12).

As an alternate approach, we performed two-dimensional QTL scans, which found evidence for possible new additive and interacting QTL that the single-QTL model could not detect (Figure 9B-D; Supplementary Table S6). These two-dimensional scans test all pairwise combinations of intervals at least 5cM apart for interactions and epistatic effects for 10 key phenotypes, including those with significant QTL from our initial analysis (those in Supplementary Table S5). However, this procedure did not reveal any new QTL for the phenotypes slope (β), inflection point (τ), PC6 or locomotion (LI7) at 35.6°. Instead, the two-dimensional QTL analysis revealed only one pair of loci with an epistatic interaction and eight pairs of loci with additive effects across different chromosomes for the remaining phenotypes (Figure 7; Supplementary Table S6). For asymptote (α) and PC1, this analysis recapitulates QTL at the center of Chromosome II. At least three QTL were also found on Chromosome V associated with the ‘synthetic’ phenotypes PC1, PC2 and PC3. Interestingly, PC2 shows evidence for an additive interaction between QTL near the centers of Chromosomes II and V, which correspond to the two distinct QTL from the single-QTL analysis for the two separate phenotypes.