Evolutionary trade-offs between heat and cold tolerance limit responses to fluctuating climates.

The evolutionary potential of species to cope with short-term temperature fluctuations during reproduction is critical to predicting responses to future climate change. Despite this, vertebrate research has focused on reproduction under high or low temperatures in relatively stable temperate climates. Here, we characterize the genetic basis of reproductive thermal tolerance to temperature fluctuations in the ostrich, which lives in variable environments in tropical and subtropical Africa. Both heat and cold tolerance were under selection and heritable, indicating the potential for evolutionary responses to mean temperature change. However, we found evidence for a negative, genetic correlation between heat and cold tolerance that should limit the potential for adaptation to fluctuating temperatures. Genetic constraints between heat and cold tolerance appear a crucial, yet underappreciated, factor influencing responses to climate change.

INTRODUCTION

Accelerated climate change is resulting in higher and more variable temperatures, posing new challenges for species (–). Previous research on vertebrates has focused on mean temperatures, rather than examining tolerance to shifts between high and low temperatures (–). While it is important to investigate the effect of mean temperatures, population persistence in the face of climate change depends on whether reproductive success is maintained during temperature fluctuations (–).

Climatic variability can reduce reproductive success jointly selecting for heat and cold tolerance. Evolutionary responses to selection can occur if there is heritable variation in heat and cold tolerance and their genetic association is weak or positive. Direct estimates of selection on heat and cold tolerance are, however, rare (), and estimating its genetic basis is extremely challenging (–). Long-term studies are needed where individuals with known genotypes are repeatedly measured, but this requires a notorious amount of effort (). In vertebrates, long-term studies of temperature effects on reproduction have primarily been carried out on temperate species (–). However, climate models predict that increases in temperature volatility will be greatest in tropical and subtropical areas (, ).

Here, we use a unique study system, the ostrich (Struthio camelus), to quantify the potential to evolve tolerance to fluctuating temperatures. The ostrich is the world’s largest bird and inhabits extreme thermal environments in tropical and subtropical Africa. We used daily records of temperature and reproductive success of 1277 individuals in experimental breeding pairs in the Klein Karoo region of South Africa over a 21-year period. Here, daily temperatures can range from −5° to 45°C during the breeding period (). We focus on temperature effects on reproduction, as survival can underestimate how temperature affects fitness (, ). We analyze selection on female egg-laying rates as our recent work shows this is a key determinant of reproductive success and is influenced by temperature (Fig. 1) (). Female egg laying was monitored daily across entire breeding seasons (~6 months), typically for 3 years per individual, providing detailed estimates of changes in reproductive success with increasing and decreasing temperatures, hereafter referred to as heat and cold tolerance. We treated heat and cold tolerance as different traits in this study because many of the mechanisms and loci governing the efficiency of thermoregulation differ under decreasing or increasing temperatures (–). Using repeated estimates of egg-laying rates of females and a nine-generation pedigree, we examined the genetic basis of reproductive heat and cold tolerance.

Fig. 1.

Ostriches (S. camelus) inhabit a variable climate that influences reproduction.

(A) Ostriches experience highly variable thermal environments both within and across populations (left, ostrich pair in temperate De Hoop Nature Reserve; right, ostrich pair in the arid Karoo National Park, South Africa; photos by C.K.C.). (B) Reproductive success, as measured by egg-laying rates, rapidly declines with deviations from a thermal optimum of around 20°C (). Points are averages across females with SEs binned according to the temperature variable. Point size illustrates relative number of females: smallest point = 56 and largest point = 652. Fitted line and 95% credible interval (CI) (shaded area) were extracted from animal random regression model (table S1).

RESULTS

There was significant stabilizing selection on heat and cold tolerance (Fig. 2). Females with egg-laying rates that were more resistant to increases and decreases in temperature had the highest reproductive success (Fig. 2). Similar patterns of stabilizing selection were also evident at the genetic level, indicating that genotypes that are more robust to temperature change have higher reproductive success {genetic correlation (rg)reproductive success-cold2 [credible interval (CI)] = −0.51 (−0.71, −0.17), pMCMC = 0.004; table S6; rgreproductive success-heat2 (CI) = −0.54 (−0.70, −0.16), pMCMC = 0.014; table S7; see “estimating selection using multi-response models”}. These results suggest that tolerance to temperature shifts during the breeding period is important in thermally volatile environments (Fig. 2), whereas in temperate species the importance of the timing of breeding has been highlighted (, ).

Fig. 2.

Heat and cold tolerance are under selection and heritable.

Females whose egg-laying rates were least affected by temperature decreases (A) and increases (B) had the highest reproductive success (quadratic selection gradients estimated using linear models were significant for both increasing and decreasing temperatures: table S2 and S3; see “estimating selection using multiple regression”). The relative cold tolerance of a female was calculated as the change in egg-laying rate when temperatures decreased standardized against her average egg-laying rate over these temperatures. Relative heat tolerance was calculated with the same approach across increasing temperatures from the optimum. Stabilizing selection on heat and cold tolerance was validated using multiresponse mixed models (tables S4 to S7; see “estimating selection using multi-response models”). (C and D) Predicted phenotypic and additive genetic variance in egg-laying rates went up as temperatures increased and decreased from the optimum. (E) Heat and cold tolerance were heritable. Three extreme data points in (A) are not shown (fig. S1). In (C) to (E), estimates and 95% CIs were calculated from the posteriors of an animal random regression model including variance of within-individual slopes (table S8).

For thermal tolerance to evolve, it needs to be heritable. The change in egg-laying rates with decreasing and increasing temperatures showed significant heritability (h2) (Fig. 2; see “modeling thermal tolerance using random regression”). Of the phenotypic variance in egg-laying rates, 17% was explained by genetic differences in heat tolerance, and 10% in cold tolerance [Fig. 2E; h2heat tolerance (CI) = 0.17 (0.08, 0.26); h2cold tolerance (CI) = 0.10 (0.03, 0.23); table S8]. Estimates of heritability of thermal tolerance slopes were robust to alternative modeling approaches: Character state models that estimate the correlation between egg-laying rates at high, benign, and low temperatures produced comparable results (see “modeling thermal tolerance using character-state models” fig. S2, and table S9).

Although heritabilities were relatively constant across temperatures (fig. S3), both phenotypic and additive genetic variances in egg-laying rates were sensitive to temperature change (Fig. 2, C and D, and table S8). Additive genetic and phenotypic variances increased with deviations from the optimum temperature (Fig. 2, C and D). The higher additive genetic variance at extreme temperatures can be important as it suggests that evolvability might increase in thermally stressful environments (–), perhaps due to the expression of cryptic genetic variation that is not detected under more benign environmental conditions ().

The sign and magnitude of the genetic correlation between hot and cold tolerance are predicted to influence the evolutionary response and rate of adaptation to fluctuating temperatures (–). If heat-tolerant genotypes are also more cold tolerant, then fluctuating temperatures should result in a higher rate of evolutionary change and faster adaptation compared with when heat and cold tolerance evolve independently. Alternatively, if heat and cold tolerance are negatively genetically correlated, the evolutionary response to selection will be constrained. Whether heat and cold tolerance are genetically correlated is unknown for most organisms in natural populations, with previous studies focusing on laboratory populations of microbes (, ) and fruit flies (, ).

We found a significant negative correlation between heat and cold tolerance at both the phenotypic and genetic levels [phenotypic correlation (r)heat-cold (CI) = −0.90 (−0.96, −0.79), pMCMC = 0.001; rgheat-cold (CI) = −0.72 (−0.93, −0.35), pMCMC = 0.008; Fig. 3, tables S1 and S10; see “modeling the relationship between heat and cold tolerance”). Consequently, females that were able to maintain egg-laying rates under higher temperatures produced fewer eggs as temperatures decreased, and vice versa. Females that were most tolerant to increasing temperatures (top 50% of predicted values from model) had a 21% reduction in their egg-laying rate as temperatures decreased by 5°C from the optimum, compared to a reduction of just 6% for the least heat-tolerant females (bottom 50%). This suggests that there is a negative and genetically based fitness trade-off between heat and cold tolerance (Fig. 3 and table S1).

Fig. 3.

Cold tolerance is negatively related to heat tolerance at the phenotypic and genetic levels.

Individual cold tolerance (slope: change in egg laying as temperatures decrease from 20°C), heat tolerance (slope: change in egg laying as temperatures increase from 20°C), and egg-laying rates under benign conditions (population optimum: 20°C). Estimates of correlations were obtained from random regression models (tables S1 and S10; see “modeling the relationship between heat and cold tolerance”).

The negative genetic correlation between heat and cold tolerance may occur via two different mechanisms. First, the optimum temperature for reproduction may differ among genotypes. In this case, genotypes with an optimum at lower temperatures will suffer greater heat stress, and genotypes with a higher optimum will suffer greater cold stress. Alternatively, the optimum temperature for reproduction may be similar across genotypes, but more cold-tolerant genotypes may be less heat tolerant, and vice versa (, ). These two, nonmutually exclusive, possibilities are difficult to disentangle because of the resolution of data required to separate their relative effects (text S1). Nevertheless, a comparison of the individuals estimated to have the highest (top 50%) and lowest (bottom 50%) heat and cold tolerance pointed to differences in reproductive thermal optima being important, suggestive of specialist genotypes being adapted to either high or low temperatures (fig. S4; see also fig. S5 for the top and bottom 30%).

Next, we investigated how the negative genetic correlation between heat and cold tolerance may have evolved. One possibility is that different combinations of heat and cold tolerance promote local adaptation to specific thermal conditions. In environments where heat stress is more pervasive, selection may favor adaptations that confer greater heat tolerance while disfavoring adaptations that increase cold tolerance if they are costly (“correlational selection”) (). Alternatively, genetic correlations between heat and cold tolerance may result from some universal genetic mechanism that pleiotropically links cold and heat tolerance across different populations, regardless of their climatic conditions. These two scenarios can be disentangled by examining whether correlations between heat and cold tolerance are population specific or are present in populations that have historically inhabited different environments.

There are several genetically and phenotypically differentiated populations of ostriches () that inhabit different areas of Africa with different temperature regimes. Three different populations are kept at the study site: Zimbabwean Blues (ZB), South African Blacks (SAB), and Kenyan Reds (KR). These populations are named after their area of origin and color variation in their skin (see “study site and population” and “modeling thermal tolerance across populations”). We investigated whether the phenotypic correlations between heat and cold tolerance differed, or were similar, across these populations. We found negative phenotypic correlations between heat and cold tolerance of similar magnitude across the divergent genetic backgrounds of southern and eastern Africa populations [rheat-cold (CI): ZB = −0.88 (−0.96, −0.36), pMCMC = 0.017; SAB = −0.90 (−0.95, −0.75), pMCMC = 0.001; KR = −0.90 (−0.98, −0.47), pMCMC = 0.002; table S13].

To further test whether the genetic basis of heat and cold tolerance was conserved across populations, we examined hybrid females. If the coregulation of heat and cold tolerance within populations evolved via similar genetic mechanisms with additive effects, hybrids should have thermal tolerances with intermediate values between parental populations. We again found that heat and cold tolerance were negatively correlated in hybrids and that their thermal optima were intermediate between parental populations (Fig. 4, table S14; see “modeling thermal tolerance across populations”). Parental populations did, however, differ in their thermal optima. Heat tolerance was prioritized over cold tolerance in east African compared with southern African populations [Fig. 4; heat tolerance: KR vs. ZB (CI) = 1.12 (2.25, −0.13), pMCMC = 0.088; KR vs. SAB (CI) = 0.76 (1.64, −0.37), pMCMC = 0.227; cold tolerance: KR vs. SAB (CI) = −2.94 (−1.32, −5.35), pMCMC = 0.001; KR vs. ZB (CI) = −2.91 (−0.64, −5.20), pMCMC = 0.010; table S13]. Together, these results suggest that there is an additive genetic basis to the coregulation of heat and cold tolerance and that their negative relationship is not a local phenomenon but the product of more universal genetic mechanisms. Consequently, if temperature fluctuations increase, then the evolution of higher heat tolerance is likely to come at a cost to cold tolerance across ostrich populations. Responses to mean temperature change may, nevertheless, be possible via an adaptive shift in how heat and cold tolerance are prioritized.

Fig. 4.

Heat and cold tolerance among three ostrich populations and their hybrids.

(A) The response of three ostrich populations with the popularized names SAB (nfemales = 494), ZB (nfemales = 68), and KR (nfemales = 26) to increasing and decreasing temperatures from the 20°C optimum (see “modeling thermal tolerance across populations”). Points are averages with SEs binned according to temperature. Point size illustrates relative number of females: smallest point = 2 and largest point = 494. Fitted line and 95% CI (shaded area) were extracted from a random regression model (table S13). (B) Populations were hybridized, and the thermal tolerance of female offspring was examined at sexual maturity (aged 3+) (SAB × KR, nfemales = 30; and SAB × ZB, nfemales = 97). Posterior means of thermal optima and their CIs were extracted from a random regression model (table S14). (C and D) Individual estimates of phenotypic thermal tolerance (slope of egg laying with change in temperature) of the two populations of hybrids (table S14).

DISCUSSION

The evolution of thermal tolerance is a central component of adaptation to climatic change (, , ). For ostriches, reproductive thermal tolerance is heritable and under selection, two prerequisites for adaptive evolutionary change. However, evolutionary responses to fluctuating temperatures may be constrained by the negative genetic correlation between heat and cold tolerance. These results suggest that predictions about responses to climate change based solely on either heat or cold tolerance will likely be inaccurate.

Whether the genetic constraints between heat and cold tolerance observed in this study are widespread across endotherms remains to be investigated. In microbes and fruit flies, there is evidence that heat and cold tolerance are genetically linked (–), suggesting that constraints on thermal adaptations may be universal across different groups of organisms. Whether certain taxa suffer more from such constraints than others and how heat versus cold tolerance is prioritized need to be established. As climatic conditions become more volatile in the future, understanding the genetic basis of traits involved in thermoregulation under decreasing and increasing temperature will become ever more important in forecasting the influence of climate change on biodiversity.

MATERIALS AND METHODS

Study site and population

The study site is situated at the Oudtshoorn Research Farm in the arid Klein Karoo of South Africa (GPS: 33°38′21.5"S, 22°15’17.4"E). Here, we used 197 enclosures of natural Karoo habitat (~0.25 ha) to monitor the reproductive success of ostrich breeding pairs. Previously, we used reproductive data from this site to show that egg-laying rate is highly sensitive to temperature fluctuations and that this reproductive parameter varies independently of other gametic traits at extreme temperatures (). The ostriches used for quantitative genetic analyses in this study were derived from 139 founding individuals, consisting of individuals classified into one of two populations with the popularized names SAB (S. camelus) or ZB (Struthio camelus australis). ZB is named by its origin in Namibia and Zimbabwe and is also referred to as the South African ostrich. The ancestry of the smaller, sometimes blacker necked, SAB is uncertain, but they potentially originated from North African (Struthio camelus camelus) and ZB ostriches hybridizing. SAB are also referred to as Struthio camelus domesticus ().

Ostriches received a diet designed for breeding individuals (120 g of protein, 10.5 MJ of metabolizable energy, 26 g of calcium, and 6 g of phosphorus per kilogram feed) and water ad libitum. A reduction in diet content across years to 90 g of protein and 7.5 MJ of metabolizable energy was done to reduce feeding costs as effects on fertility were negligible (, ). Since 1990, offspring from breeding pairs have been recruited back into the population each year, enabling relatedness between all individuals at the study site to be estimated through a nine-generation pedigree. Maximum daily temperature records were obtained from a local weather station 600 m from the study site. Ethical clearance was obtained from the Western Cape Department of Agriculture (DECRA R12/48).

Reproductive data

One adult male and one adult female ostrich were assigned to an enclosure in May/June of each year and kept together until the end of the breeding season in December/January. Enclosures were checked for eggs twice a day, yielding 652 female records for egg-laying rates and a total of 1830 year-by-female records (on average, 2.8 years per female from 1998 to 2018). Removal of eggs each day results in females laying throughout the breeding season. This has the advantage of allowing temperature effects on reproduction to be consistently monitored across long periods (~6 months), resulting in very large sample sizes: Ostriches are able lay one egg every other day, giving on average 94 opportunities for egg laying per female per year. While the experimental removal of eggs is artificial, it is not entirely different from their natural biology. Females can lay in multiple nests, foregoing incubation, that together with a high rate of nest failure (, ) results in females laying across breeding seasons.

We removed data from pairs where the male or female was replaced during the breeding season, which occurred occasionally when individuals were injured or died. Data on the egg-laying rates of these replacement pairs indicated that acclimation to enclosures and/or partners takes approximately 45 days (). On the basis of this information, we removed data from the first 45 days of each season. Two-year-old females have substantially lower reproductive success than older breeders (); so, these were also removed from the data. In 2 years, the breeding season was extended to February/April, but we removed data from these months to ensure consistency across years. Similarly, we removed data from pairs that spent fewer than 200 days in their enclosure in a given year. Pairs that laid fewer than 10 eggs per year were removed to avoid including incompatible pairs and individuals that are not in breeding condition.

Statistical analyses

Measuring responses to temperature change

Before modeling the effect of temperature on egg-laying rates, the following needed to be established: time-lag effects of sensitivity to temperature change, the optimum temperature for egg laying, and the best way to measure temperature change. A previous investigation showed that females are most sensitive to fluctuations in ambient temperature 2 to 4 days before egg laying and that egg-laying rates were optimized when the daily temperature maximum was ~20°C (). To analyze the effects of temperature change on egg laying, we measured temperature differences from the 20°C optimum by grouping egg-laying records into seven temperature bins (continuous variable): four increasing (20° to 23.2°C, 23.3° to 26.1°C, 26.2° to 29.5°C, and 29.5° to 43.6°C) and three decreasing (20° to 18.5°C, 18.5° to 16.8°C, and 16.8° to 8.9°C). Temperature bins were created so that egg-laying rate could be modeled as a binomial response variable: the proportion of potential laying days (days/2, hereafter referred to as 2 days) within each temperature bin where an egg was laid for each female during a given year (#eggs/2 days). This avoided treating egg-laying rate as a binary variable (egg laid or not), which often have convergence problems in complex quantitative genetic models. To make the intercept of statistical models represent the most benign temperature, we set 20°C to 0. The variance of slopes (see below) depends on the scale of the environmental parameter. We therefore standardized data by dividing by 24°C, the maximum temperature change from the optimum (20° to 44°C), resulting in 1 being the maximum temperature change.

To allow heat and cold tolerance to vary independently, we defined a factor, temperature type, to denote whether the change in temperature was due to increases from the optimum (20° to 44°C, referred to as heat tolerance; nrecords = 128,451) or decreases (20° to 9°C, referred to as cold tolerance; nrecords = 43,989). Heat and cold tolerance were analyzed as separate traits because this allowed us to test for “general temperature tolerance” (positive correlations between heat and cold tolerance), independent heat and cold tolerance (no correlation), and antagonistic relationships between heat and cold tolerance (negative correlations). Although the lower absolute temperatures recorded in our study are relatively mild, we follow the current convention and refer to reductions in individual performance from the optimum temperature as cold tolerance (, , ).

General modeling strategy

We performed three main sets of analyses: (i) selection for heat and cold tolerance using multiple regression, (ii) selection for heat and cold tolerance using multiresponse models, and (iii) the genetic basis to heat and cold tolerance.

Data were analyzed using Bayesian generalized linear mixed models (GLMMs) implemented in the R package MCMCglmm v.2.2 () in R v.3.6.0 (). All models had the following basic structure unless otherwise stated. The fixed effects included were “female age” (continuous, mean centered, and scaled to unit variance), “female population” (factorial: SAB, ZB, or hybrid), and “population of pair male” (factorial: SAB, ZB, or hybrid). The random effects included were “enclosure,” which varied in size and vegetation cover and were repeatedly used across years, and “female ID” to account for the repeated sampling of each female. For random terms, we used the weakly informative inverse-Gamma distribution [scale = 0.001, shape = 0.001, i.e., V = diag(n), nu = n − 1 + 0.002, with n being the dimension of the matrix] as priors. We examined the sensitivity of a model to the prior by using a parameter expanded prior with a lower pull toward zero [V = diag(n), nu = n, alpha.mu = rep(0, n), alpha.V = diag(n)*252, with n being the dimension of the matrix] and found similar results in all cases. Each model was run for 9,100,000 iterations, of which the initial 100,000 were discarded, and only every 5000th iteration was used for estimating posterior probabilities. The number of iterations was based on inspection of autocorrelation among posterior samples in preliminary runs. Convergence of the estimates was checked by running the model three times and inspecting the overlap of estimates in trace plots and the level of autocorrelation among posterior samples. We report the posterior mode of random effects and posterior mean of fixed effects, in addition to 95% CIs.

Estimating selection using multiple regression

Selection gradients for heat and cold tolerance were estimated using GLMMs. Reproductive success was our proxy for fitness, measured as the egg-laying rate (total #eggs/2 days) of a female across a year. To create a measure of relative fitness, this proportion was divided by the mean within each year and therefore transformed to a Gaussian response variable. We then entered our measure of relative fitness as a response variable in one model for cold tolerance and one model for heat tolerance. These models were implemented using the basic model structure (see “general modeling strategy”) but with additional fixed effects of the linear and quadratic terms of relative cold or heat tolerance (continuous, mean centered, and scaled to unit variance). Nonlinear selection gradients were estimated by multiplying the quadratic regression coefficient by two (). Heat and cold tolerance were measured as the relative change in egg-laying rate between adjacent temperature bins (Eq. 1)

Relative changes were averaged across all adjacent pairs of bins to produce one estimate of heat and cold tolerance per individual per year and to capture the change in egg-laying rate across the full scale of increasing or decreasing temperatures. The change in egg-laying rate between two adjacent temperature bins was scaled by the mean egg-laying rate for that temperature type (increasing versus decreasing temperature, Eq. 1; see text S2 for additional statistical support) because individuals with high egg-laying rates by definition show the largest absolute reduction in laying. The estimate of relative heat tolerance is therefore the average change in egg-laying rate across four temperature bin comparisons (see “measuring responses to temperature change”), and relative cold tolerance the average across three temperature bin comparisons. Records where a female did not produce any eggs under a given temperature type in that year were omitted (nrecords decreasing temperatures = 38, nrecords increasing temperatures = 0).

We verified that measuring hot and cold tolerance by taking an average across temperature bins was appropriate by confirming that the relative change in egg laying between adjacent pairs of temperature bins had similar relationships with our fitness proxy (text S2). We took the average instead of the difference between the lowest and the highest bins of increases in temperatures, as this would be inaccurate for individuals with very high or low egg-laying rates at intermediate bins of temperature increase. Estimates of heat and cold tolerance are also produced by models outlined in “modeling thermal tolerance using random regression”. However, we avoided using these model estimates as these compound error across analyses, causing anticonservative results ().

Estimating selection using multiresponse models

Estimating selection gradients (see “estimating selection using multiple regression”) is the standard approach to investigate the presence of selection (). However, it does not allow the strength of selection on phenotypes and genotypes to be separated and does not allow the error in fitness and thermal tolerance to be modeled simultaneously. Consequently, we estimated stabilizing and directional selection at the phenotypic and genotypic levels using multiresponse models, one for heat tolerance and one for cold tolerance. Each model included three response variables: (i) reproductive success, measured as egg-laying rate (total #eggs/2 days) over the year, modeled as a binomial response variable using a logit link-function (“multinomial2”); (ii) a linear term of relative heat or cold tolerance modeled as a Gaussian response variable; and (iii) a quadratic term of relative heat or cold tolerance modeled as Gaussian response variable. Covariance matrices (see below) between reproductive success and the linear term of relative heat tolerance provide information about directional selection, while covariance matrices between reproductive success and the quadratic term of relative heat tolerance provide information about stabilizing selection. The same approach applies to relative cold tolerance. Both the linear and quadratic terms of heat or cold tolerance were mean centered and scaled to unit variance before modeling.

Multiresponse models of selection were implemented using the basic model structure (see “general modeling strategy”) with the following modifications. We estimated the fixed effect of female population (factorial: SAB, ZB, or hybrid) separately for each response variable. We accounted for environmental effects that varied across years, such as diet, by having year as a random effect. As the variance in one response variable caused by year effects may differ from another response variable, we estimated year variances separately for each response variable. This approach was also applied to enclosure.

To measure female variance in each of our response variables, and estimate covariances among responses variables, we included a random effect of female ID as a 3 × 3 unstructured variance-covariance matrix. To capture uncertainty in each of our response variables within each female, a second 3 × 3 unstructured variance-covariance matrix of residual variance was included. This is possible as females were monitored for several years, and each estimate of our response variables was based on data from only 1 year. Such partitioning of variance ensures that the trait covariances among females needed to measure selection are not biased by covariance between traits within a year. To estimate genetic variances and covariance across response variables, we ran a set of models where an “animal” term linked to the pedigree was included as a third 3 × 3 unstructured variance-covariance matrix. Selection was quantified by calculating the correlations between reproductive success and either the linear (directional selection) or quadratic term (stabilizing selection) of thermal tolerance through the variance-covariance matrix [rtrait1-trait2 = covariancetrait1,trait2/sqrt(vartrait1*vartrait2)]. Selection at the phenotypic level (r) was quantified through the female ID term (model without the animal term) and at the genetic level (rg) through the animal term. Each model was run for 3,100,000 iterations, the initial 100,000 were discarded, and every 3000th iteration was used for estimating posterior probabilities.

Modeling thermal tolerance using random regression

The genetic basis of heat and cold tolerance in egg-laying rates was examined using random regression animal models (RRAMs) () in a mixed model framework (, ). Egg-laying rate (#eggs/2 days) was fitted as a binomial response variable using a logit link-function (multinomial2) and modeled using the basic model structure (see “general modeling strategy”), with the following modifications. Models included the fixed effects of temperature change across bins (continuous, ranging from 0 to 1) and temperature type (factorial: decreases or increases) (see “measuring responses to temperature change”). The interaction between temperature change and temperature type was modeled with a common intercept for decreases and increases, as the way temperature change was calculated dictated identical intercepts. We included interactions between female population, temperature change, and temperature type. Temperature change and temperature type were interacted with female ID, to estimate permanent environment variance (pe), and with an animal term linked to the pedigree, to estimate additive genetic variance (a), in both intercepts and slopes (). These were modeled as two 3 × 3 unstructured variance-covariance matrices composed of the intercept (peint or aint), slope during temperature decreases (pesl-cold tolerance or asl-cold tolerance, i.e., the cold tolerance slope), and slope during temperature increases (pesl-heat tolerance or asl-heat tolerance, i.e., the heat tolerance slope). Last, we included year as a random effect. Narrow sense heritability (h2) at the optimum temperature (temperature change = 0, corresponding to 20°C) was then estimated as the proportion of intercept variance explained by the additive genetic variance

Variation in heat and cold tolerance slopes represents the genetic and phenotypic variation in responses to temperature change across females. To express the proportion of individual slope variance that is heritable (i.e., the heritability of thermal plasticity), we constructed a second set of models. In these models, we added a third 3 × 3 unstructured variance-covariance matrix of individual-by-year (id-yr) combinations, capturing the within-individual variance in slopes. Variance in individual slopes is at a different scale to that of intercepts and is also dependent on the scaling of temperature change. For these reasons, we followed a recently introduced practice (, ) that enabled us to estimate the heritability of thermal plasticity as the proportion of slope variance attributable to additive genetic variance as follows

Using this second set of models, we also estimated the environment-dependent additive genetic variance and heritability for each temperature type x following (, )

There is ongoing debate whether the fixed-effect variance (σ2f) should be included in the denominator when estimating heritabilities (). It has been argued that σ2f should be included if the fixed-effect variance captures natural variation, or excluded if it represents experimental variance (), but the distinction is not always clear. For full transparency, we provide estimates of σ2f excluding variance from the temperature change (σ2f) as this parameter has already been accounted for by the interaction with the random terms. We estimated fixed-effect variance of all terms (σ2f) and of temperature change alone (σ2f) following (), such that σ2f = σ2f − σ2f (tables S1, S8, and S9).

As egg-laying rate is modeled using a logit link-function, these estimates are calculated on the latent scale. While this scale has the benefit of fulfilling the typical assumptions of parametric analyses, it may not reflect the scale at which selection is working, and methods therefore have been developed to make inferences on the observed scale (). Methods are currently lacking for performing this transformation for models using logit link-functions and where the number of trials varies between data points, in our case the number of days with and without eggs per female. Instead, it is possible to calculate estimates of heritability at the expected scale (corresponding to the liability scale in a threshold model) according to equations in () using the R package QGglmm (). Similar methods are not available for the slope variance parameters presented above, and all estimates in the main document are therefore on the latent scale for consistency. When possible, we also provide estimates at the expected scale in the Supplementary Materials (tables S1, S8, and S9).

Modeling thermal tolerance using character-state models

To verify results from random regression models (see “modeling thermal tolerance using random regression”), we examined the genetic basis to egg-laying rate under changing temperatures using character-state models. Here, we model the changes in egg-laying rate across three thermal states (cold, benign, and hot). These types of models produce estimates of the means and variances at each state, as opposed to estimates of the rates of change (slopes) obtained from the random regression models. The ranges for these states were limited by the lower number of cold days compared to hot, according to the thermal optimum cutoff used in the random regression analysis (20°C). To avoid low replication in the cold state relative to the hot state, we classified the lowest 50% of days as “cold” (<17.7°C, nrecords = 21,488), and the highest 30% of days as “hot” (>28.8°C, nrecords = 38,557), with the remainder being classified as “benign” (nrecords = 110,307). The models followed the same general approach as the random regression models with egg-laying rate (#eggs/2 days) fitted as a binomial response variable (see “modeling thermal tolerance using random regression”). The major difference was that thermal state (factorial: cold, benign, and hot) was included as a fixed effect (instead of temperature change and temperature type). Thermal state was interacted with the female ID (pe), and animal term [linked to the pedigree (a)], in two 3 × 3 unstructured variance-covariance matrices composed of the cold, benign, and hot thermal states for pe or a. We also estimated the residual variance separately for each thermal state.

Modeling the relationship between heat and cold tolerance

Phenotypic and genetic correlations between heat and cold tolerance were estimated using covariance estimates from the random regression models described in “modeling thermal tolerance using random regression”. To test the relationship between egg-laying rates at benign temperatures and heat or cold tolerance, we calculated correlations between the intercept and heat or cold slopes. To test the relationship between heat and cold tolerance, the correlations between heat and cold slopes were calculated. Phenotypic correlations were calculated using the variance-covariance matrices of female ID interacted with temperature change and temperature type in models run without the animal term (tables S10 to S12). To calculate genetic correlations, we used the variance-covariance matrices of the animal term linked to the pedigree (table S1 and S9).

Modeling thermal tolerance across populations

To investigate whether the observed relationship between heat and cold tolerance is a population-specific phenomenon or present across populations of ostriches, we examined the relationship between heat and cold tolerance within different ostrich populations. Three genetically different populations are kept at the study site: SAB, ZB, and KR (Struthio camelus massaicus), also referred to as Masai ostrich. KR females were excluded from the quantitative genetic analyses because of their limited number of individuals and fewer years of data collection (nfemales = 26, nyears = 12). The ZB populations also had relatively few individuals (nfemales = 68, nyears = 21), which prevented estimation of genetic correlations separately for each population. However, phenotypic correlations between heat and cold tolerance could be estimated for each population, using a dataset that included all individuals with at least 85% expected relatedness to one of the three populations (SAB, ZB, and KR). With this dataset, we constructed a random regression model following the same general approach as in “modeling thermal tolerance using random regression”, with egg-laying rate (#eggs/2 days) fitted as a binomial response variable. In this model, the interaction between temperature change, temperature type, and the random effect female ID was modeled with three separate 3 × 3 variance-covariance matrices, one for each population. This enabled us to estimate phenotypic correlations specific to each population. This model was run without the pedigree.

To further understand whether the genetic basis to thermal tolerance is conserved across populations, we compared individuals that were hybrids (reciprocal crosses) from two populations. At the field site, the SAB population had been crossed with the ZB population (SAB × ZB: nfemales = 97, nyears = 19) and the KR population (SAB × KR: nfemales = 30, nyears = 9). We performed an analysis identical to the population analyses described above, but with five female population categories (three parental populations and two hybrids). This model allowed us to compare heat and cold tolerance of hybrids with the parental populations, both using the population means and by extracting individual slopes of heat and cold tolerance. Individual slopes were estimated as the sum of the population slope and the slope of the female ID term to capture individual differences in thermal tolerance.