Quantitative genetics of temperature performance curves of Neurospora crassa

Earth’s temperature is increasing due to anthropogenic CO2 emissions; and organisms need either to adapt to higher temperatures, migrate into colder areas, or face extinction. Temperature affects nearly all aspects of an organism’s physiology via its influence on metabolic rate and protein structure, therefore genetic adaptation to increased temperature may be much harder to achieve compared to other abiotic stresses. There is still much to be learned about the evolutionary potential for adaptation to higher temperatures, therefore we studied the quantitative genetics of growth rates in different temperatures that make up the thermal performance curve of the fungal model system Neurospora crassa. We studied the amount of genetic variation for thermal performance curves and examined possible genetic constraints by estimating the G-matrix. We observed a substantial amount of genetic variation for growth in different temperatures, and most genetic variation was for performance curve elevation. Contrary to common theoretical assumptions, we did not find strong evidence for genetic trade-offs for growth between hotter and colder temperatures. We also simulated short term evolution of thermal performance curves of N. crassa, and suggest that they can have versatile responses to selection.

: A) An illustration of a hypothetical temperature performance curve. Temperature is on the horizontal axis and growth rate is on the vertical axis. T opt shows the optimal temperature, where growth rate has its maximum value, µ max . Temperature where growth rate reaches zero as temperature increases is denoted as CT max . B) Change in reaction norm elevation shifts the reaction norm on the vertical axis. C) A horizontal shift. D) Change in reaction norm shape.
peratures with the expense of reduction of performance in cold temperatures resulting in a hot-cold thermodynamic effect (Sørensen et al., 2018).
To calculate the optimum temperature for each genotype without using a specific model that 150 may fit for some genotypes better than others, we fitted natural splines for each genotype. We 151 extracted the maximum growth rate (µ max ) and optimum temperature (T opt ) from the spline fit. We 152 then fit a model 153 ln(y i ) ∼ N(µ i , σ) (1) where y i is the ith maximum growth rate, α is the intercept, β is the slope, and T opt,i is the ith 154 optimum temperature. We used weakly regularizing priors: a normal distribution for α and β, and 155 a half-cauchy distribution for σ with location 0 and scale 2. MCMC estimation was done using two 156 chains, with a warmup of 1000 iterations, followed by 4000 iterations of sampling. For this analysis Estimation of genetic variance and covariance components 167 We were interested in estimating the genetic variance and covariance components for growth rates 168 at different temperatures that together describe different aspects of temperature performance curves. 169 Because Neurospora can be propagated clonally, we can estimate genetic variance components 170 using clonal analysis. Among genotype variance is an estimate of the genetic variance and within 171 genotype variance is an estimate of the environmental variance (Lynch and Walsh, 1998). We used 172 a multivariate model to estimate genetic variance components at each temperature and the genetic 173 correlations of all possible temperature pairs. The advantage of this approach is that we do not have 174 to assume any particular shape for the temperature reaction norm. The multivariate model was where α is a vector of intercepts, α g [i] is a vector of genotypic effects, S G and S E are 6 × 6 diagonal 176 matrices with genetic or environmental standard deviations on the diagonal, and R G and R E are 177 matrices for genetic and environmental correlations respectively. We used weakly informative 178 priors by using the half location-scale version of the student's t distribution with three degrees 179 of freedom and 10 as the scale parameter. Thus, the prior for intercept effects was hT (3,2,10) hT (3,3,10) hT (3,4,10) hT (3,4,10) hT (3,4,10) hT(3, 3, 10) for growth rates from 20 to 40°C. The prior for each standard deviation in the model was 182 σ ∼ hT(3, 0, 10), and we used an lkj prior for the correlation matrices: R E , R G ∼ LKJ(1).

183
For MCMC estimation two chains were run with a warmup period of 1000 iterations, followed by 184 5000 iterations of sampling, with thinning set to 2. By inspecting MCMC traceplots ( Figure S2) 185 and the diagnostic summary statisticR, which was 1 for all parameters, we found no evidence of 186 convergence problems.

187
Genetic correlations and temperature differences 188 We were also interested in how the genetic correlation of growth rates changes as temperatures 189 are further apart. In order to examine how correlations change in a statistically rigorous manner, 190 we calculated pairwise temperature differences for each estimated genetic correlation (n = 15), 191 and fitted a Bayesian linear model with genetic correlation as the response, taking into account the 192 uncertainty in the estimated genetic correlations. This is a linear model with measurement error 193 where uncertainty in the estimated genetic correlations is propagated to the intercept and slope 194 estimates of the linear model, see McElreath (2015) for details. We compared models with or 195 without slope effects for temperature and whether genetic correlations involving growth rate at 40 196°C had a different intercept or slope (Table 2). We used leave-one-out cross-validation method for x est,i ∼ N(µ i , σ) (4) where x obs,i is the median of ith observed genetic correlation, x sd,i is the observed standard devi-201 ation of ith genetic correlation, x est,i is the estimated genetic correlation for ith observation, α is 202 the intercept, α 40 is the intercept effect when one of the temperatures is 40°C, c i is an indicator 203 variable whether one of the temperatures is 40°C, β is the slope effect, and d i is the temperature 204 difference for the ith observation. MCMC estimation was done using two chains, with a warmup 205 of 1000 iterations, followed by 4000 iterations of sampling.

206
Quantitative genetics 207 We estimated heritability, the proportion of genetic variance of the total variance, for each temper-208 ature as

214
where σ 2 A is the additive variance and σ 2 AA is the additive × additive epistatic variance, σ 2 AAA is the 215 additive × additive × additive variance, and so on (Lynch and Walsh, 1998). With our experimental 216 design we cannot estimate the epistatic variance terms, as is the case with many other common 217 quantitative genetic designs, and going further we assumed that epistatic variances were small and were ignored. This seems like a strong assumption, but there is some justification for doing so: 219 even if there is plenty of epistasis at the level of gene action, this is not necessarily translated into 220 epistatic variance (Hill et al., 2008;Mäki-Tanila and Hill, 2014). Empirical data also suggest that 221 most genetic variation is additive (Hill et al., 2008). The genetic covariance of traits 1 and 2 is 222 cov G 1,2 = σ G 1 σ G 2 r G 1,2 , where r G 1,2 is the correlation of the standard deviations or the genetic 223 correlation. Thus, genetic correlation for traits 1 and 2 can be defined as 225 In addition to heritabilities, we used coefficients of variation to compare genetic and environ-226 mental variances. Heritability can be influenced by changes in either genetic or environmental 227 variance, and genetic variance by itself is not a unitless variable (Houle, 1992). The genetic coeffi-228 cient of variation was: wherez is the mean phenotype. Accordingly, the environmental coefficient of variation was CV E = 231 100σ E /z.

232
We obtained the G-matrix to describe how the growth rates at different temperatures were cor-233 related and to be able to calculate multivariate response to selection for thermal performance curves.

234
This matrix contains genetic variance components on the diagonal and covariance components on 235 off-diagonals, so for n traits G is an n × n matrix: 237 For environmental variance it is possible to construct an analogous E-matrix that is the environ-238 mental variance-covariance matrix.
We performed eigen decomposition of the G-matrix to gain insight into genetic constraints 240 of reaction norm evolution. The eigenvector corresponding the leading eigenvalue, or the first 241 principle component, gives the direction of multivariate evolution with the least genetic resistance 242 (Schluter, 1996). We obtained these components by principle component analysis of the G-matrix.

243
To assess uncertainty in the eigen decomposition we constructed a G-matrix for each posterior

260
To asses whether evolvability along a certain selection gradient is particularly high or low it is possi-  where S is a vector of selection differentials for each temperature, G and P are the genetic and 281 phenotypic variance-covariance matrices respectively, and R is the response to selection. Response 282 to selection can also be expressed in terms of the selection gradient, β, as where β = P −1 S. The biological interpretation of selection differential and selection gradient are different, as selection differential of 0 for a given trait does not imply selective neutrality but rather 286 stabilizing selection, whereas selection gradient of 0 for a trait implies that the trait is selectively 287 neutral. See figure S3 for an illustration of the differences between these concepts. When we asked 288 how would evolution proceed in a particular direction, we simulated response to selection using 289 selection gradient (Equation 13). And when we asked whether selection could generate a particu- ized to be 0.6 mm/h. We estimated evolvability and conditional evolvability along these gra-308 dients as explained above. Then we used different selection regimes to examine how we could shapes ( Figure 2A), possibly reflecting that these genotypes were poorly suited to lab conditions, 330 due to specific nutritional requirements for example.

331
Thermodynamics of thermal performance curves 332 We examined whether differences between genotypes could be explained by a thermodynamic 333 effect, i. e. does the maximum growth rate increase with optimum temperature. We obtained µ max 334 and T opt from the natural spline fits and plotted ln(µ max ) against 1/(kT opt ) ( Figure 2B). For the and removing outliers agree, it seems that removing the outliers is quite reasonable in this case.

347
In order to analyse the data without forcing the tolerance curves to fit any predetermined shape, or  (Table 1).

353
By plotting the model means and genetic correlations it appeared that genetic correlation be-354 tween adjacent temperatures was generally high, and decreased as temperatures were further apart 355 and correlations involving 40°C also seemed lower ( Figure 3A). We tested this idea formally 356 and fitted a model of genetic correlations and temperature differences. We compared the differ-357 ent models, and the best model had different intercepts for correlations involving 40°C and for 358 correlations not involving 40°C, and identical slopes for these two groups (Table 2) growth at 40°C than in lower temperatures.

368
Most of the variation observed in growth rates was due to genetic variation present among the 369 strains. Heritabilities for growth at different temperatures were high, around 0.89 for temperatures 370 from 20 to 35°C ( Figure 3C). As temperature increased further heritability dropped to 0.63 at 40 371°C ( Figure 3C). However, this lower heritability was not due to decreased genetic variation but 372 increased environmental variance at 37.5 and 40°C ( has the most genetic variance (Table 1), but this would have been also misleading as coefficient 377 of genetic variation reveals that growth at 40°C has the most genetic variation followed by the other 378 temperatures in decreasing order ( Figure 3D). The same was true for coefficient of environmental 379 variation ( Figure 3D).  (Table 1).

394
When looking trait specific evolvabilities we also observed that growth rate at 40°C had the 395 highest conditional evolvability and the highest autonomy (Table 3). This indicates that out of all 396 of the growth rates, growth rate at 40°C can evolve by itself most easily. The rest of the traits had 397 very low autonomies reflecting their high genetic correlations (Tables 1 and 3).

398
Evolution of performance curves 399 In order to examine how a performance curve of a population that has the same G-matrix as es-  to maintain the original phenotype at the other temperatures there were also correlated responses 419 but these were less uniform ( Figure S3). Accordingly, selection at a single temperature often lead 420 to correlated responses in nearby temperatures ( Figure S4). Selection at multiple temperatures 421 lead to stronger responses to selection and correlated responses ( Figures S5 and S6). For instance, 422 selection at 25 and 30°C increased growth rate also at 20°C ( Figure S5). and 30°C, response to selection after 5 generations of selection was 3.01 (2.98-3.04, 95% HPD).

429
Thus, it was not possible to change a certain temperature completely independently of the others, 430 but often extreme temperatures could be changed without affecting the growth at the other extreme.

431
We then asked is it possible to create similar evolutionary responses in performance curves as 432 shown in Figure 1. We were able to find a set of selection differentials that were able to generate 433 changes in elevation, horizontal shift, or shape ( Figure 5). This shows that despite strong genetic 434 correlations it is possible for the performance curves to evolve in almost any manner if selection 435 favors such a performance curve. However, selection regimes involving horizontal shifts require 436 selection for increased growth rate in some temperatures and decreased growth rate in others (Fig-437 ure 5). Evolvabilities and conditional evolvabilities were highest for elevation changes. For opti-438 mum shifts and shape changes conditional evolvabilites were lower than the average conditional 439 evolvability over all phenotypic space ( Figure 5C). This indicates that elevation changes are less 440 constrained than changes in optimum temperature or performance curve shape.   Figure 5: Simulated responses to selection using selection differentials, S. A) Selection differentials for each selection regime. Black line is the mean empirical performance curve and dots represent values of selection differentials for each temperature. Blue dots represent selection for increased growth, red for decreased growth and black dots indicate stabilizing selection at this temperature. Selection regime 1 selects for increased elevation, regime 2 selects for decreased elevation, regime 3 selects for lower optimum temperature, regime 4 selects for higher optimum, regime 5 selects for broader shape, and regime 6 selects for narrower shape. B) Simulated responses to selection. C) Evolvability and conditional evolvability for each of the selection gradients implied by the selection differentials. Red horizontal lines show the average conditional evolvability (c) across the entire phenotypic space. Average evolvability is higher than the y-axis scale and is not shown.

442
We have shown that there is substantial genetic variation in thermal performance curves of Neu-443 rospora crassa. Most of this variation is in performance curve elevation and there is very little 444 evidence of strong trade-offs. Genetic variation in growth is strongly correlated among nearby 445 temperatures but there is a threshold before or at 40°C after which this correlation drops, indi-446 cating that physiological processes at 40°C are different than those at lower temperatures. Such

465
Another possibility is that genetic differences between the strains in how well they are able to 466 grow in lab conditions are thermodynamically amplified, as increasing temperature also increases 467 metabolic rate (Schulte, 2015). However, our estimates of activation energy were much lower than  (Whitlock, 1996;Kassen, 2002 one temperature did not limit plasticity (Fragata et al., 2016;Manenti et al., 2015), most genetic 494 variation has been observed for overall performance (Klepsatel et al., 2013;Latimer et al., 2015), or 495 that adaptation was largely temperature specific with no apparent trade-offs (Bennett et al., 1992).
Genetic correlations between growth rates at nearby temperatures were strong, which is to be 497 expected, as a difference of a few°C is likely to be a very similar physical environment for an 498 organism. However, growth rate at 40°C had a lower genetic correlation to growth rates at other 499 temperatures. This suggests that at 40°C there was some physiological process activated, which 500 has genetic variation, but that was not active or was at much lower level of activity in lower tem-501 peratures. The most obvious candidate for such a process is the heat shock response (Piper, 1993;502 Feder and Hofmann, 1999;Sørensen et al., 2003). Previously the heat shock response of N. crassa 503 has been studied at 42°C or higher (Mohsenzadeh et al., 1998;Plesofsky-Vig and Brambl, 1985;504 Guy et al., 1986) but it probably occurs already at lower temperatures, as we observed significant 505 slow down of growth at 40°C. The canonical heat shock proteins are important for the physio-506 logical heat shock response, but there can be additional mechanisms involved: there is evidence 507 that the sugar trehalose plays some role in N. crassa heat shock response (Bonini et al., 1995).   Table 1: Genetic variances, covariances, correlations, and environmental variances for growth rates in different temperatures estimated from the multivariate model. The diagonal (in bold) contains genetic variances (σ 2 G ), upper triangle contains genetic covariances (σ Gx σ Gy r Gx,y ), and lower triangle contains genetic correlations (r Gx,y ). The last column contains environmental variances (σ 2 E ). Estimates are posterior means with 95% HPD intervals shown in parenthesis.
(      Figure S1: A) Phenotypic means for each genotype, those genotypes that were removed from the thermodynamic analysis as outliers are coloured red. B) Logarithm of maximum growth rate, µ max , plotted against inverse of kT opt . Datapoints that were removed as outliers are coloured red. Black regression line is ordinary regression fitted to the data without outliers (red points removed), slope (±SE) is −0.16(±0.03). Red line is ordinary regression fitted to all of the data, slope is −0.34(±0.06). Blue line is robust regression with bisquare weights fitted to all of the data, slope is −0.17(±0.03).  Figure S3: Illustration on how selection differential and selection gradient generate different responses to selection. Top row: selection based on selection differentials. A) There is selection for increased growth at at 20°C and to maintain the original phenotype for the other traits. Black line is the empirical performance curve and blue dots represent selection differentials for each temperature, S = {0.2, 0, 0, 0, 0, 0}. B) The realized selection gradient (β = P −1 S) implied by this selection differential. C) Blue line shows the mean phenotype after 5 generations of selection, shaded area contains 95% of the simulations. Bottom row: selection based on selection gradient. D) There is selection for increased growth at 20°C as in A) but selection gradient is β = {0.2, 0, 0, 0, 0, 0}. E) Now realized selection gradient is the same as in (D), there is selection for increased growth in one temperature but phenotypes of other temperatures are selectively neutral. E) Response to selection after 5 generations of selection as in (C) for this selection gradient. Selection using similar selection differentials and selection gradient leads to different phenotypic responses.  Figure S5: Selection for increased growth rate in two temperatures. Selection differential is 0.2 mm/h at each generation for each temperature, so 0.4 in total for each selection regime. A) Selection differentials for each selection regime. B) Simulated responses to selection.  Figure S6: Selection for increased growth rate in three temperatures. Selection differential is 0.2 mm/h at each generation for each temperature, so 0.6 in total for each selection regime. A) Selection differentials for each selection regime. B) Simulated responses to selection.  Family C This study -C3

Supplementary Tables
Family C This study -C4 Family C This study -C5 Family C This study -C6 Family C This study -C7 Family C This study -C8 Family C This study -C9 Family C This study -C10 Family C This study -C11 Family C This study -C12 Family C This study -C13 Family C This study -C14 Family C This study -C15 Family C This study -C16 Family C This study -C17 Family C This study -C18 Family C This study -C19 Family C This study -C20 Family C This study -C21 Family C This study -C22 Family C This study -C23 Family C This study -C24 Family C This study -C25 Family C This study -C26 Family C This study -C27 Family C This study -C28 Family C This study -C29 Family C This study -C30 Family C This study -C31 Family C This study -C32 Family C This study -C33 Family C This study -C34 Family C This study -C35 Family C This study -C36 Family C This study -C37 Family C This study -C38 Family C This study -C39 Family C This study -C40 Family C This study -C41 Family C This study -C42 Family C This study -C43 Family C This study -C44 Family C This study -C45 Family C This study -C46 Family C This study -Continued on next page. . .