# Partitioning the phenotypic and genetic variances of reaction norms

Pierre de Villemereuil, Luis-Miguel Chevin

► To cite this version:

Pierre de Villemereuil, Luis-Miguel Chevin. Partitioning the phenotypic and genetic variances of reaction norms. Peer Community Journal, 2025, 5, pp.e15. 10.24072/pcjournal.515. hal-05392851

HAL Id: hal-05392851

https://hal.science/hal-05392851v1

Submitted on 2 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research article

Published
2025-01-27

**Cite as**
Pierre de Villemereuil and
Luis-Miguel Chevin (2025)
Partitioning the phenotypic and
genetic variances of reaction
norms, Peer Community
Journal, 5: e15.

**Correspondence**
pierre.de-villemereuil@mnhn.fr

**Peer-review**
Peer reviewed and
recommended by
PCI Evolutionary Biology,
https://doi.org/10.24072/pci.evolbiol.100708

<figure><img src="image_2.png" /></figure>

This article is licensed
under the Creative Commons
Attribution 4.0 License.

# Partitioning the phenotypic and genetic variances of reaction norms

Pierre de Villemereuil<sup>1,2</sup> and Luis-Miguel Chevin<sup>3</sup>

Volume 5 (2025), article e15

https://doi.org/10.24072/pcjournal.515

## Abstract

Many traits show plastic phenotypic variation across environments, captured by their norms of reaction. These reaction norms may be discrete or continuous, and can substantially vary in shape across organisms and traits, making it difficult to compare amounts and types of plasticity among (or even within) studies. In addition, the evolutionary potential of phenotypic traits and their plasticity in heterogeneous environments critically depends on how reaction norms vary genetically, but there is no consensus on how this should be quantified. Here, we propose a partitioning of phenotypic variance across genotypes and environments that jointly address these challenges. We start by distinguishing the components of phenotypic variance arising from the average reaction norm across genotypes, genetic variation in reaction norms (with additive and non-additive components), and a residual that cannot be predicted from the genotype and the environment. We then further partition the genetic variance of the trait (additive or not) into an environment-blind component and a component arising from genetic variance in plasticity. We show that the additive components can be expressed, and further decomposed according to the relative contributions from each parameter, using what we describe as the reaction norm gradient. This allows for a very general framework applicable from the character-state to curve-parameter approaches, including polynomial functions, or arbitrary non-linear models. To facilitate the use of this variance decomposition, we provide the Reacnorm R package, including a practical tutorial. Overall the toolbox we develop should serve as a basis for a unifying and deeper understanding of the variation and genetics of reaction norms and plasticity, as well as more robust comparative studies of plasticity across organisms and traits.

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Institut de Systématique, Évolution, Biodiversité (ISYEB), École Pratique des Hautes Études PSL, MNHN, CNRS, SU, UA, Paris, France, ²Institut Universitaire de France (IUF), ³CEFE, CNRS, Université de Montpellier, Université Paul Valéry Montpellier 3, EPHE, IRD, Montpellier, France</footnote>
Introduction

The phenotype of a given genotype can vary in response to its environment of development or expression, through a phenomenon broadly described as phenotypic plasticity (Bradshaw 1965; Schlichting and Pigliucci 1998). Phenotypic plasticity is currently attracting considerable interest in the context of rapidly changing natural environments (Chevin et al. 2010; Gienapp et al. 2008; Merilä and Hendry 2014). While the mere existence (and even prevalence) of phenotypic plasticity is uncontroversial, its relative contribution to observed or predicted phenotypic change in the wild (Bonamour et al. 2019; Gienapp et al. 2008; Merilä and Hendry 2014; Teplitsky et al. 2008), as well as the extent of its interplay with population-level processes such as natural selection and population dynamics (de Villemereuil et al. 2020; Reed et al. 2010; Schaum and Collins 2014; Vedder et al. 2013), are very active research areas. Answering these questions requires biologists to be able to dissect and compare phenotypic plasticity in detail in a wide range of traits, environmental contexts and species. This requires a methodology that is appropriate for each context, while being general enough to be comparable across contexts.

The relationship between the phenotype and the environment is captured by the reaction norm (or norm of reaction), which is defined at the level of genotypes (Schlichting and Pigliucci 1998; Wolterek 1909). Reaction norms encompass phenotypic responses to both continuous environments (such as temperature, salinity, etc.) and categorical/discrete ones (such as host plant for a phytophagous insect). Within a simple model of reaction norm, quantifying plasticity may be straightforward. For instance, both empirical (Charmantier et al. 2008; Nussey et al. 2005) and theoretical (Gavrilets and Scheiner 1993a; Lande 2009) work have extensively relied on the assumption of a linear reaction norm, whose slope is used as a metric of plasticity, since it quantifies how much phenotypic change is induced per unit environmental change. However, regression slopes are signed and have units of trait per environment, so even in this simple case some standardisation is needed in order to compare the magnitude of plasticity among studies. Beyond this simple scenario, drawing robust conclusions about phenotypic plasticity requires being able to quantify and compare its magnitude across organisms, traits and environments, in a way that is applicable across the statistical frameworks used to study plasticity.

Beyond *how much* phenotypes change with the environment, *how* they change can also be of importance. First, different reaction norm shapes may come with different biological interpretations. For instance, a bell-shaped (eg quadratic, Gaussian) reaction norm may indicate that some mechanism underlying a measured trait is maximized at an intermediate value of the environment. This is often expected for traits that are direct components of fitness, or that can be interpreted as proxies for performance, for which the reaction norms are generally termed tolerance or performance curves (Angilletta 2009; Deutsch et al. 2008; Lynch and Gabriel 1987). A sigmoid shape, on the other hand, may indicate that plasticity is directional but that the range of possible phenotypes is constrained, or that selection favors discrete-like variation (Chevin et al. 2013; Hammill et al. 2008; Moczek and Emlen 1999; Suzuki and Nijhout 2006). Second, most theoretical models on the evolution of plasticity, especially those based on quantitative genetics which are most directly comparable to empirical data, assume a given reaction norm shape - often linear for simplicity (Lande 2009; Scheiner 1993b; Tufto 2000). The extent to which theoretical predictions on the evolution of plasticity apply to any particular empirical system thus depends on how well the reaction norm shape assumed in the models conforms to observations in this system. In other words, we need some metric for whether a reaction norm is "mostly linear" or "mostly curved", for instance. In addition, when fitting a particular model of reaction norm shape to an empirical dataset, we would like to know how well this model captures the overall plastic variation of the trait across environments.

A third crucial question regarding reaction norms is how (and how much) they vary genetically. It has long been recognized that plasticity can evolve if reaction norms vary genetically (Bradshaw 1965), and theory has predicted how different aspects of reaction norm shape are expected to respond to selection in a variable environment (de Jong 1990; Gavrilets and Scheiner 1993a; Gomulkiewicz and Kirkpatrick 1992). However this theory has been little applied empirically, except for predictions about the slope of linear reaction norms (or phenotypic differences
between two environments). But beyond this, it should also be of interest to identify which aspects of reaction norm shape are more likely to evolve, based on how they vary genetically. For instance, a reaction norm may be highly curved (e.g. quadratic) but have little genetic variability in curvature, instead mostly varying in position, height, or local slope. Distinguishing between the genetic variance of the trait, marginalised across environments, and the genetic variance of plasticity itself, can also be a conceptual and methodological challenge. There is thus a need to compare genetic variation in different components of reaction norm, but previous attempts to do so (in a meta-analysis) were limited by methodological obstacles (Murren et al. 2014, see the Appendix). In fact, comparing genetic variation in the slope versus curvature of a reaction norm, for instance, is not straightforward, as these parameters have different scales and even units (trait per environment, vs trait per squared environment). Moreover, even the notion of average slope and curvature can have different meanings depending on the assumed distribution for the environment. Genetic variation in reaction norm shape can be analyzed by estimating variation in the parameters of a continuous function of the environment, as done by the flexible framework of function-valued traits (Gomulkiewicz and Kirkpatrick 1992; Kirkpatrick and Heckman 1989; Stinchcombe et al. 2012). In addition, it would be useful to be able to compare the relative contributions of variation in different aspects of reaction norm shape to the overall variance arising from plasticity of a trait.

We herein propose a theoretically justified and generally applicable framework to estimate and partition the phenotypic variance of reaction norms, towards three main goals: (i) quantify the contribution of plasticity to the total phenotypic variance in reaction norms; (ii) evaluate the contribution of different aspects of reaction norm shape, and of the full assumed reaction norm model, to overall plastic phenotypic variation; and (iii) quantify heritable variation in the trait and its plastic component, due to the different aspects of the reaction norm. We provide this framework as a new R package Reacnorm, including a tutorial to guide users in applying it. Our hope is that this will stimulate more quantitative investigations of the ways in which phenotypic plasticity contributes to phenotypic variation and evolutionary change.

## Reaction norm models

In the broadest sense, a reaction norm is a decomposition of phenotypic variation among known (often controlled) versus unknown sources of environmental variation. In this sense, we can start by decomposing the phenotypic trait z into two components:

$$ (1) \qquad z = \hat{z} + \tilde{z}. $$

The first term $\hat{z}$ is the reaction norm, that is, the component of phenotypic variation that can be predicted (hence the hat notation) from knowing both the genotype (which we will note g throughout) of an individual and the environment (which we will note ε throughout) in which it developed. Note that by “environment”, we mean either an experimentally controlled environmental variable, or a focal variable (e.g. temperature) within a naturally occurring environmental context. The second term $\tilde{z}$ is the component of the measured phenotype that cannot be predicted from genotype and environment, and arises from unknown environmental factors (usually described as micro-environmental variation), developmental noise, and measurement error.

Types of reaction norms $\hat{z}$ can be further categorised according to the type of environmental variation. The environment may be inherently categorical and unordered, such as host plant for a herbivore insect. It may be ordered but with no (or unknown) quantitative value, such as low, medium, and high treatments. Or it may be ordered quantitatively, with values that are either intrinsically discrete, such as habitat quality, or continuous, such as temperature or salinity.

When environments are categorical, the reaction norm can be studied by treating phenotypic values in different environments as alternative ‘character states’, considered as different traits in a multivariate framework (Falconer 1952; Via and Lande 1985). The mean character state may differ among environments if the trait is plastic; phenotypic and genetic variation may be larger in some environments; and phenotypes may be more or less correlated across environments (Falconer 1952; Via and Lande 1985). Such a modelling framework is readily described by Equation 1 for a genotype g and environment $\epsilon_k$ (where the index k is used to reflect the discrete aspect
<table><caption>**Table 1** – List of the main notations, as well as their source of variation. We here distinguish the “focal” environment, which only concerns the environmental variable used to parametrise the reaction norm, from other putative sources of environmental variation that may influence the phenotypic trait (sometimes described as micro-environmental variation). “Everything” in the table thus includes all (focal and other) sources of environmental and genetic variation, developmental noise and measurement error.</caption><thead><tr><th>Notation</th><th>Explanation</th><th>Varies over</th></tr></thead><tbody><tr><td>z</td><td>Phenotypic value for the trait</td><td>Everything</td></tr><tr><td>ẑ</td><td>Phenotype as predicted from the environment and the genotype</td><td>Focal environment, genotypes</td></tr><tr><td>ε</td><td>Environmental variable</td><td>—</td></tr><tr><td>μ</td><td>Vector of the average value of the phenotypic in each environment</td><td>Focal environment</td></tr><tr><td>G<sub>z</sub></td><td>Additive genetic variance-covariance matrix of trait values across environments (character states)</td><td>—</td></tr><tr><td>θ<sub>g</sub></td><td>Vector of parameter values of the reaction norm for genotype g</td><td>Genotypes</td></tr><tr><td>θ̄</td><td>Vector of mean values of the reaction parameters over the genotypes</td><td>—</td></tr><tr><td>G<sub>θ</sub></td><td>Additive genetic variance-covariance matrix of the reaction norm parameters</td><td>—</td></tr><tr><td>ψ<sub>ε</sub></td><td>Reaction norm gradient, the vector of partial derivatives of the phenotype z against reaction norm parameters θ<sub>g</sub>, averaged over the genotypes at environment ε</td><td>Focal environment</td></tr><tr><td>Ψ</td><td>Variance-covariance matrix of ψ<sub>ε</sub> across environments</td><td>—</td></tr><tr><td>V<sub>P</sub></td><td>Total phenotypic variance in the trait z</td><td>—</td></tr><tr><td>V<sub>Res</sub></td><td>Residual variance, not explained by the reaction norm</td><td>—</td></tr><tr><td>V<sub>Plas</sub>, P<sup>2</sup><sub>RN</sub></td><td>Phenotypic variance arising from changes in the mean reaction norm across environments; divided by V<sub>P</sub> for P<sup>2</sup><sub>RN</sub></td><td>—</td></tr><tr><td>V<sub>Gen</sub>, H<sup>2</sup><sub>RN</sub></td><td>Total genetic variance in the trait across environments; divided by V<sub>P</sub> for H<sup>2</sup><sub>RN</sub></td><td>—</td></tr><tr><td>V<sub>Add</sub>, h<sup>2</sup><sub>RN</sub></td><td>Total additive genetic variance in the trait across environments; divided by V<sub>P</sub> for h<sup>2</sup><sub>RN</sub></td><td>—</td></tr><tr><td>V<sub>A</sub>, h<sup>2</sup></td><td>Environment-blind additive genetic variance of the trait, i.e. based on the mean breeding values across environments, divided by V<sub>P</sub> for h<sup>2</sup></td><td>—</td></tr><tr><td>V<sub>A×E</sub>, h<sup>2</sup><sub>I</sub></td><td>Additive genetic variance arising from plasticity, i.e variance of the mean-centred breeding values, divided by V<sub>P</sub> for h<sup>2</sup><sub>I</sub></td><td>—</td></tr><tr><td>π<sub>Sl</sub>, π<sub>Cv</sub></td><td>Proportion of V<sub>Plas</sub> explained by the average slope (π<sub>Sl</sub>) or curvature (π<sub>Cv</sub>) of the average reaction norm</td><td>—</td></tr><tr><td>φ<sub>i</sub>, φ<sub>ij</sub></td><td>Proportion of V<sub>Plas</sub> explained by parameter i, or by covariation between parameter i and j for a polynomial reaction norm</td><td>—</td></tr><tr><td>γ<sub>i</sub>, γ<sub>ij</sub></td><td>Proportion of V<sub>Add</sub> explained by the additive genetic (co)variation in parameter i (and j)</td><td>—</td></tr><tr><td>ι<sub>i</sub>, ι<sub>ij</sub></td><td>Proportion of V<sub>A×E</sub> explained by the additive genetic (co)variation in parameter i (and j)</td><td>—</td></tr></tbody></table>

of the environmental variable). In practice, such an approach would correspond to an ANOVA (or a mixed model) with discrete environment and genotype-within-environment as (random)
effects of the model. In its most compact form, such a statistical model can be framed as a multivariate Gaussian distribution, with the number of dimensions corresponding to the number of categories in the environment,

$$ (2) \qquad \hat{z} \sim \mathcal{N}(\mu, G_z), $$

where $\mu$ is the vector of expected phenotypic values (across genotypes) within each environment, and $G_z$ is the genetic variance-covariance matrix of trait values within and across environments.

For quantitative environments (both discrete and continuous), the most common approach is to model the reaction norm as a function of environment and genotype:

$$ (3) \qquad \hat{z} = f(\varepsilon, \theta_g), $$

where $\varepsilon$ is the environmental value, and $\theta_g$ is a vector that contains the parameters of the function (e.g. coefficients associated to each exponent for a polynomial) for each genotype $g$; these parameters are thus genetically variable. The parameters $\theta_g$ are generally assumed to be polygenic and thus follow a multivariate Gaussian distribution,

$$ (4) \qquad \theta_g \sim \mathcal{N}(\bar{\theta}, G_\theta), $$

where $\bar{\theta}$ is the vector of average parameter values across genotypes and $G_\theta$ is the additive genetic variance-covariance matrix of the parameters $\theta_g$. This approach has been described alternatively as the “reaction norm” approach, the “polynomial approach”, or a parametric version of function-valued traits. To keep it general here and avoid confusion with the general concept of reaction norm as defined in Equation 1 (which applies even to categorical environments), we will describe it as the “curve-parameter” approach. Note that Equation 4 assumes that the only source of variation in reaction norm parameters $\theta$ is genetic. In cases where reaction norms can be measured in individuals using repeated measurements across environments (individual plasticity sensu Nussey et al. 2007) it can be necessary, or useful, to include other sources of variation in $\theta$, including confounding environmental effects, or permanent environmental effects. For the sake of simplicity, we will assume throughout that all variation in $\theta$ is genetic, but we show in subsection C5 that relaxing this assumption only affects how non-genetic variances are computed.

It can be shown that the character-state and curve-parameter approaches are equivalent, following the spirit of de Jong (1995), who showed that a polynomial curve of sufficient order is exactly equivalent to a character-state model. In particular, the character-state in Equation 2 can be expressed using Equation 3 and Equation 4 by letting $\bar{\theta} = \mu$, $G_\theta = G_z$ and $f$ a function that outputs the $k$th value of $\theta_g$ when evaluated at $\varepsilon_k$ environment (see Appendix A). In the following, we will derive general results using the more general formalism of Equation 3 and Equation 4, and then express them for the particular case of the character-state approach when relevant.

## Partitioning variation in reaction norms

### Complete partition of the variation in reaction norms

The total phenotypic variance in the reaction norm can be partitioned by isolating independent components of variation. The main reasoning will be summarised here, with more mathematical details provided in the Appendix A to Appendix D. For a start, the terms in Equation 1 are assumed to be independent, such that the total phenotypic variance $V(z)$ (usually noted $V_P$) is the sum of the variance predicted by the genotype and the environment $V(\hat{z})$, plus a residual component of variance $V(\tilde{z}_i)$, which we will note $V_{\text{Res}}$. Then, a second distinction can be made between the general, average shape of the reaction norm, and the genotype-specific variation surrounding such an average, as illustrated in Figure 1 using a quadratic reaction norm. The component of phenotypic variance arising from plastic responses to the environment by the mean reaction norm, i.e. after averaging across all genotypes (Figure 1), will be denoted $V_{\text{Plas}}$. This variance can be considered as fully ascribed to the environmental component of phenotypic variation. The component of phenotypic variation attributable to genetic variation in the reaction
norm Figure 1 will be denoted $V_{Gen}$. As these two components are independent by construction, denoting as $E_{g|ε}(\hat{z})$ the expected value of the reaction norm across genotypes at a given environmental value $ε$, we have

$$ (5) \qquad V(\hat{z}) = V(E_{g|ε}(\hat{z})) + V(\hat{z} - E_{g|ε}(\hat{z})) = V_{Plas} + V_{Gen}, $$

such that

$$ (6) \qquad V_P = V_{Plas} + V_{Gen} + V_{Res}. $$

Compared to the classical equation $V_P = V_G + V_E + V_{G×E}$ (Des Marais et al. 2013; Falconer and Mackay 1996; Lynch and Walsh 1998), the correspondence is that $V_E = V_{Plas} + V_{Res}$ and $V_{Gen} = V_G + V_{G×E}$. Also note that both decompositions make the same common assumption that genotypes and environments are not correlated. We have thus decomposed the environmental variance into a component due to phenotypic plasticity in response to $ε$ ($V_{Plas}$) on the one hand, and any other residual source of phenotypic variation ($V_{Res}$) on the other hand, as commonly done in theory (Gavrilets and Scheiner 1993a; Via and Lande 1985) as well as in practice.

<figure><img src="image_4.png" /><figcaption>**Figure 1** – Illustration of the full variance decomposition using quadratic reaction norms. We start from the reaction norms (left graph, grey lines, the residual variance is not illustrated) and compute their average shape across all genotypes (left graph, red line). The phenotypic variance arising from this average shape is $V_{Plas}$. Centering the reaction norms along this average shape directly yields the distribution of the breeding values along environments (middle graph, purple lines), because in this quadratic case, the non-additive genetic variance is $V_{NonAdd} = 0$. The total variance of the breeding values along the environment is $V_{Add}$. The classical, environment-blind additive genetic variance $V_A$ is the variance of the breeding values averaged across environments for each genotype (middle graph, green dots). The $V_{A×E}$ is the variance of the remainder of the breeding values after mean-centring (right graph, blue lines).</figcaption></figure>

The genotypic variance $V_{Gen}$ accounts for all sources of genetic variation, including the genotype-by-environment interaction. Note that this contrasts with a view where the genotype-by-environment interaction is instead associated with the environmental component, e.g. as plastic variance (Falconer and Mackay 1996; Lynch and Walsh 1998; Scheiner 1993a; Scheiner and Lyman 1989). As seen above, $V_{Gen}$ can be decomposed into the genetic variance of the trait, measured using its average genotypic value across environments ($V_G$), and the variance arising from genotype-by-environment interaction ($V_{G×E}$). Here, we will apply such decomposition at the level of the additive genetic variance ($V_{Add}$), relegating all the non-additive parts of $V_G$ and $V_{G×E}$ into a common $V_{NonAdd}$ component (Figure 1), arising from dominance and epistasis (Falconer and Mackay 1996; Lynch and Walsh 1998). Usually, models like Equation 2 or Equation 4 are defined using additive genetic variance-covariance matrices for their basic parameters, meaning that $V_{Add}$ can be directly estimated from the models. As such, we will discard explicit inclusion of dominance or epistasis variance components in a theoretical or statistical model throughout, for
the sake of simplicity. However, non-additive genetic variance can still arise from non-linearity in the (assumed) developmental system (de Villemereuil 2018; de Villemereuil et al. 2016; Morrissey 2015; Rice 2004), meaning that non-additive variance can be generated by the reaction norm itself. Looking at Equation 3 and Equation 4, the ultimate source of any additive genetic variation in the trait z comes from the additive genetic variation in the parameters θ. As a result, non-additivity in the trait arises when the function f(ε, θ) in Equation 3 is non-linear with regard to θ, a situation we will refer to as “non-linearity in the parameters”. Importantly, this means that polynomial (e.g. quadratic) functions, which are linear in their parameters, are such that $V_{NonAdd} = 0$ and $V_{Gen} = V_{Add}$.

When studying the evolution of plasticity, it proves useful to further decompose $V_{Add}$ into two components. The first is the environment-blind additive genetic variance of the trait, arising from differences in average breeding values between genotypes, and typically equal to the classical $V_A$. In other words, $V_A$ is the variance of the breeding values after averaging them across environments (Figure 1), as would be obtained if the genotype-by-environment interaction was ignored altogether. For example, it would be the output of a simple animal model analysis of repeated measurements of a plastic trait in a wild population. The second component of $V_{Add}$ is the additive genetic variance arising from plasticity, which we will note $V_{A×E}$ (for additive genetic component due to genotype-by-environment interactions). $V_{A×E}$ is the remaining additive genetic variance in the reaction norm after removing the mean breeding value for each genotype (Figure 1). This definition is akin to the one used by Albecker et al. (2022), but here more directly expressed in terms of variance of breeding values, i.e. additive genetic variance. It measures the potential for evolution of plasticity in the trait. Notably, if $V_{A×E} = 0$ but $V_{Add} > 0$, then the additive genetic variation in the reaction norms is only due to average differences between genotypes, i.e. the reaction norms of different genotypes are parallel. The variances $V_A$ and $V_{A×E}$ are exactly equivalent to the classical decomposition using $V_G$ and $V_{G×E}$, only applied to the heritable part of the genetic variance. We show below that it is possible to express $V_{Add}$, $V_A$ and $V_{A×E}$ in a way that encompasses all approaches of reaction norm, from a character-state to a curve that is non-linear in its parameters, by computing reaction norm gradients of the trait z with respect to its reaction norm parameters θ, in line with previous theoretical results for the quantitative genetics of non-linear developmental systems and non-Gaussian traits (de Villemereuil et al. 2016; Morrissey 2015).

The complete partition of the phenotypic variance is thus:

$$ (7) \qquad V_P = V_{Plas} + V_A + V_{A×E} + V_{NonAdd} + V_{Res}. $$

From this, it is possible to derive unitless quantities of interest, for instance by standardising by the phenotypic variance, which is more widely applicable and appropriate than mean-standardisation in the context of reaction norms (Pélabon et al. 2020). In particular:

$$ (8) \qquad P_{RN}^2 = \frac{V_{Plas}}{V_P}, $$

is the proportion of the phenotypic variance arising from average plastic responses to environments (depending on the average reaction norm shape). Variance-standardised additive genetic variances are heritabilities. In our case, we can use $V_{Add}$, $V_A$ or $V_{A×E}$ as the numerator, yielding the following relationship:

$$ (9) \qquad h_{RN}^2 = \frac{V_{Add}}{V_P} = \frac{V_A}{V_P} + \frac{V_{A×E}}{V_P} = h^2 + h_I^2. $$

In other words, the heritability of the trait when fully accounting for its reaction norm ($h_{RN}^2$) is equal to the environment-blind heritability of the trait ($h^2$), based on the breeding values averaged across environments) plus the heritability from plasticity ($h_I^2$, based on the breeding values by environment interaction). If it is not possible to measure additive genetic variances due to limitations in the experimental design (e.g. when “genotypes” correspond to populations, accessions or clones), it is possible to perform the same decomposition using “broad-sense heritabilities”,

$$ (10) \qquad H_{RN}^2 = \frac{V_{Gen}}{V_P} = \frac{V_G}{V_P} + \frac{V_{G×E}}{V_P} = H^2 + H_I^2. $$
In all cases, the quantity:

$$ (11) \qquad T_{RN}^2 = \frac{V_{Plas} + V_{Gen}}{V_P} = P_{RN}^2 + H_{RN}^2 $$

would measure the proportion of the phenotypic variance explained by the (possibly plastic and genetically variable) reaction norm, and thus our ability to predict the individual phenotype from the genotype and the environment. In a linear context with respect to the parameters, when the environment is considered a fixed quantity, the quantities $P_{RN}^2$ and $T_{RN}^2$ are analogous to the (resp. marginal and conditional) coefficient of determination of the reaction norm (Johnson 2014; Nakagawa and Schielzeth 2013), but their definition here is given beyond that simple context. Relaxing the assumption that the only source of variation in $\theta$ is of genetic origin (e.g. individual plasticity, Nussey et al. 2007), we show in subsection C5 that only the computation of $V_P$ and $T_{RN}^2$ are slightly affected.

Importantly, so far we are not making any statement about the actual reaction norm shape: $P_{RN}^2$ captures the contribution of the average reaction norm regardless of its shape, and the broad- or narrow-sense heritabilities the contribution of various aspects the genetic variation to the phenotypic variance. The contribution of detailed aspects of reaction norms shape to phenotypic variation are obtained by further partitioning $V_{Plas}$ and the additive genetic variances, as we do below.

## Contributions of reaction norm shape and parameters to $V_{Plas}$

As stated in Equation 5, the general definition of the variance arising from the average reaction norm is $V_{Plas} = V(E_{g|\varepsilon}(\hat{z}))$. Important simplifications arise in more particular cases. For example, when the assumed curve is linear in its parameters, $E_{g|\varepsilon}(\hat{z}) = f(\varepsilon, \bar{\theta})$, where $\bar{\theta}$ is the average value of the parameters across genotypes. In particular, in the case of a quadratic reaction norm (Gavrilets and Scheiner 1993b; Morrissey and Liefting 2016; Scheiner 1993a):

$$ (12) \qquad f(\varepsilon, \theta_g) = (\bar{a} + a_g) + (\bar{b} + b_g)\varepsilon + (\bar{c} + c_g)\varepsilon^2, $$

where $\bar{a}, \bar{b}, \bar{c}$ are the average intercept, first- and second-order parameters of the model, and $a_g, b_g$ and $c_g$ are genotype-specific deviation from these average values for the same parameters, we can express $V_{Plas}$ simply as:

$$ (13) \qquad V_{Plas} = \bar{b}^2V(\varepsilon) + \bar{c}^2V(\varepsilon^2) + 2\bar{b}\bar{c}\text{cov}(\varepsilon, \varepsilon^2). $$

If the environmental variable $\varepsilon$ has been mean-centred and is symmetrical, then $\text{cov}(\varepsilon, \varepsilon^2) = 0$ and the third term vanishes. Finally, in the case of a character-state model, the average phenotype in each environment $\varepsilon_k$ is readily provided by the $\mu_k$ in Equation 2, so that $V_{Plas} = V(\mu)$. Once $V_{Plas}$ is computed, its standardised version $P_{RN}^2$ follows by dividing by the total phenotypic variance.

Pushing the analysis further, we aim to compute the contributions of different aspect of reaction norm shape to the overall environmental plastic variance of the trait, notably the contribution of its slope and curvature, which we will denote as $\pi_{Sl}$ and $\pi_{Cv}$, respectively. For this, at least one of two of the following assumptions must be valid: (i) $\varepsilon$ follows a normal distribution, or (ii) the true reaction norm is quadratic. In all cases, it also require that the environmental variable has been mean-centered. A last requirement is for $f$ to be at least twice differentiable with respect to $\varepsilon$ (which excludes e.g. the character-state approach). In this case, these terms simply depend on the average first- and second-order derivative of $E_{g|\varepsilon}(\hat{z})$ and the variance of $\varepsilon$ and $\varepsilon^2$ (see subsection D1):

$$ (14) \qquad \pi_{Sl} = \frac{E\left(\frac{dE_{g|\varepsilon}}{d\varepsilon}(\hat{z})\right)^2 V(\varepsilon)}{V_{Plas}}, \qquad \pi_{Cv} = \frac{\frac{1}{4}E\left(\frac{d^2E_{g|\varepsilon}}{d\varepsilon^2}(\hat{z})\right)^2 V(\varepsilon^2)}{V_{Plas}}. $$

An important point arising from Equation 14 is that the relative importance of variation in the slope and curvature components of reaction norm depend on variation in the environment, respectively $V(\varepsilon)$ and $V(\varepsilon^2)$ (note that $V(\varepsilon^2) = 2V(\varepsilon)^2$ if the environment is normally distributed).
Crucially, we chose to express this partitioning using the mean environment as the reference environment (as commonly practiced, e.g. Morrissey and Liefting 2016), but any other choice of a reference environment would result in a different π-partition, notably due to a non-null value for Cov(ε, ε²). Fortunately, neither V<sub>Plas</sub> nor P<sub>RN</sub><sup>2</sup> are impacted by this choice in the reference environment. Furthermore, if the reaction norm is linear in the parameters, the derivatives of E<sub>g|ε</sub>(ẑ) can be directly taken as the derivatives of f. In particular, for a quadratic reaction norm as in Equation 12, for a mean-centred environment, those quantities simply are:

$$ (15) \qquad \pi_{\text{Sl}} = \frac{\bar{b}^2 V(\varepsilon)}{V_{\text{Plas}}}, \qquad \pi_{\text{Cv}} = \frac{\bar{c}^2 V(\varepsilon^2)}{V_{\text{Plas}}}, $$

consistent with the fact the first and second order coefficients of a quadratic polynomial correspond to its average slope and curvature, respectively. Only in this configuration do we have π<sub>Sl</sub> + π<sub>Cv</sub> = 1. Unfortunately, this simple, geometric interpretation of the polynomial coefficients is lost above the second-order case (see Appendix D).

Figure 2 shows the values of π<sub>Sl</sub> and π<sub>Cv</sub> for various quadratic reaction norms, assuming ε follows either a normal or uniform distribution, with same mean 0 and variance 1. The values for π<sub>Sl</sub> and π<sub>Cv</sub> translate well the perceived “trendiness” (for large π<sub>Sl</sub>) or “curviness” (for large π<sub>Cv</sub>) of reaction norms, but they may also strongly depend on the statistical distribution of the environmental variable ε, as shown especially in the third example of Figure 2. In this example, the difference arises because the assumed environmental distributions have different kurtosis (the scaled fourth central moment, related to V(ε²) in Equation 15). Because V(ε²) is larger for the Gaussian, this distribution leads to larger π<sub>Cv</sub> than the uniform.

<figure><img src="image_5.png" /><figcaption>Figure 2 - Computation of π<sub>Sl</sub> = π<sub>b</sub> and π<sub>Cv</sub> = π<sub>c</sub>, the relative contributions of linear and quadratic terms to phenotypic variation caused by the mean reaction norm, for different shapes of reaction norms, and two distributions of the environmental variable ε: a standard Gaussian (of mean 0 and variance 1), and a uniform distribution between -√3 and √3 (of mean 0 and variance 1).</figcaption></figure>

When it is not possible to assume that ε is normally distributed (because it is discrete, or experimentally constrained) and a quadratic assumption is not a good fit to the reaction norm, it is always possible to use a higher-order polynomial model to approximate the true reaction norm, in line with theoretical work by de Jong (1990, 1995) and Gavrilets and Scheiner (1993b). In this case, we can conduct an alternative decomposition based on the parameters of the polynomial (rather than the mean slope and curvature of the function), using the fact that a polynomial curve is linear in its parameters. To distinguish this parameter-based decomposition from the specific decomposition in terms of slope and curvature, we use a different notation. The relative contribution of a given exponent m in the polynomial to the variance caused by the mean plasticity becomes (see subsection D2)

$$ (16) \qquad \varphi_m = \frac{\bar{\theta}_m^2 V(\varepsilon^m)}{V_{\text{Plas}}}, $$
and the contribution of the covariance between exponents $l$ and $m$ is

$$ (17) \qquad \varphi_{lm} = \frac{2\bar{\theta}_l\bar{\theta}_m \mathrm{Cov}(\varepsilon^l, \varepsilon^m)}{V_{\mathrm{Plas}}} $$

Note that even with a symmetrical and mean-centred environment, the covariance between higher-order exponents will not be zero in general, contrary to $\varepsilon$ and $\varepsilon^2$ in the quadratic case. Using orthogonal polynomials would solve this issue of covariances, but at the cost of a more complex interpretation of the coefficients. More generally, this $\varphi$-decomposition only relies on the assumption that the reaction norm is linear on its parameters, which includes polynomials as a particularly useful special case. We summarise the requirements and applications for the $\pi$- and $\varphi$-decomposition depending on the context in Figure 3.

<figure><img src="image_6.png" /><figcaption>Figure 3 – Decision tree summarising our suggested workflow for the computation and decomposition of $V_{\mathrm{Plas}}$, depending on the nature of the environmental variable, its normality and the validity of a quadratic approximation of the reaction norm shape.</figcaption></figure>

## Contributions of reaction norm parameters to the genetic variance

We can express the variance of the genotypic values of the reaction norms in Equation 5 in a slightly different, but more operational, manner:

$$ (18) \qquad V_{\mathrm{Gen}} = V(\hat{z} - E_{g|\varepsilon}(\hat{z})) = E(V_{g|\varepsilon}(\hat{z})), $$

i.e. the total genotypic variance of the reaction norms is equal to the environment-specific genotypic variance averaged across environments. As explained above, this total genetic variance can be further decomposed into the genetic variance and the genotype-by-environment variance, i.e. $V_{\mathrm{Gen}} = V_G + V_{G \times E}$ (Des Marais et al. 2013; Falconer and Mackay 1996; Lynch and Walsh 1998). From an evolutionary perspective, the component of main interest is rather the total additive genetic variance of the reaction norm $V_{\mathrm{Add}}$, which will be the main focus of this section. As a reminder, we here assume, that the experimental design allows for the inference of the additive genetic variance of the parameters of the reaction norm ($G_z$ or $G_\theta$ above), and that non-additive variance in the trait $V_{\mathrm{NonAdd}}$ only arises when the reaction norm is non-linear in the parameters (i.e. dominance and/or epistasis were not fitted in the statistical model). This assumption is for the sake of simplicity, as our framework can include such effects into $V_{\mathrm{Gen}}$ if needed.
A general way to relate the additive genetic variance of the trait to the additive genetic variances of the reaction norm parameters is through a vector that we describe as the reaction norm gradient, which we will note $\psi_{\varepsilon}$ (following notations in de Villemereuil et al. 2016),

$$ (19) \qquad \psi_{\varepsilon} = E_g \left( \frac{\partial z}{\partial \theta} \right)_{\varepsilon}, $$

where the subscript $\varepsilon$ makes it clear that $\psi_{\varepsilon}$ will generally be a function of the environment. In the case of a quadratic curve, $\psi_{\varepsilon}$ is the $(1, \varepsilon, \varepsilon^2)^T$ vector (see subsection C3 for a polynomial of arbitrary order). In the case of a character-state model, $\psi_{\varepsilon_k}$ is a vector with 1 for the kth environmental level (or character state), and zero elsewhere. Whether or not the reaction norm is linear in its parameters, the additive genetic variance of the trait in a given environment $\varepsilon$ is (de Villemereuil et al. 2016; Morrissey 2015, and see Appendix B),

$$ (20) \qquad V_{A|\varepsilon} = \psi_{\varepsilon}^T G_{\theta} \psi_{\varepsilon}, $$

where superscript T denotes matrix transposition, $G_{\theta}$ the genetic covariance matrix of reaction norm parameters as defined in Equation 4 for the curve-parameter approach, and $G_{\theta}$ is $G_z$ from Equation 2 for the character-state approach. The total additive genetic variance in the reaction norm, $V_{Add}$, is the average of $V_{A|\varepsilon}$ across environments (see subsection C1):

$$ (21) \qquad V_{Add} = E \left( \psi_{\varepsilon}^T G_{\theta} \psi_{\varepsilon} \right). $$

The environment-blind additive genetic variance of the trait $V_A$, based on breeding values averaged across environments, is (see subsection C2)

$$ (22) \qquad V_A = E(\psi_{\varepsilon})^T G_{\theta} E(\psi_{\varepsilon}). $$

Although some elements of $E(\psi_{\varepsilon})$ and $G_{\theta}$ can be negative, the fact that $G_{\theta}$ is a variance-covariance matrix ensures that $V_A \ge 0$ (see subsection C2). The additive genetic variance arising from plasticity is thus (see subsection C2):

$$ (23) \qquad V_{A \times E} = V_{Add} - V_A = E \left( \psi_{\varepsilon}^T G_{\theta} \psi_{\varepsilon} \right) - E(\psi_{\varepsilon})^T G_{\theta} E(\psi_{\varepsilon}). $$

If we define $\Psi = E(\psi_{\varepsilon} \psi_{\varepsilon}^T) - E(\psi_{\varepsilon}) E(\psi_{\varepsilon})^T$, the variance-covariance matrix of the reaction norm gradients across environments, then a more intuitive way to express $V_{A \times E}$ is as a sum, for all pairs of parameters, of the (co)variance of their reaction norm gradient across environments (in $\Psi$) and their additive genetic (co)variance (in $G_{\theta}$):

$$ (24) \qquad V_{A \times E} = \sum_{i,j} \Psi_{(i,j)} G_{\theta(i,j)} = \mathrm{Tr}(\Psi G_{\theta}), $$

where Tr is the trace of a matrix. All of the quantities above can be divided by $V_P$ to get the corresponding heritabilities.

To illustrate with an example, for a quadratic reaction norm with mean-centred environment as shown in Figure 1, $\psi_{\varepsilon} = (1, \varepsilon, \varepsilon^2)$ and thus we have (see subsection C3)

$$ (25) \qquad \begin{aligned} V_{Add} &= V_a + (V_b + 2C_{ac})E(\varepsilon^2) + V_cE(\varepsilon^4), \\ V_A &= V_a + 2C_{ac}E(\varepsilon^2) + V_cE(\varepsilon^2)^2, \\ V_{A \times E} &= V_bV(\varepsilon) + V_cV(\varepsilon^2), \end{aligned} $$

where $V_a$, $V_b$ and $V_c$ are the additive genetic variances in the parameters $a_g$, $b_g$ and $c_g$, and $C_{ac}$ is the additive genetic covariance between the intercept $a_g$ and the second-order effect $c_g$. Those expressions are reminiscent of classical results from the theory of evolution of plasticity (e.g. de Jong 1990; Gavrilets and Scheiner 1993b), especially regarding the crucial role of $C_{ac}$ in the evolution of quadratic reaction norms, but here distinguishing three important components of the additive genetic variance of reaction norms. In particular, we see how the additive genetic variance arising from plasticity, $V_{A \times E}$, can be simply expressed as the sum of the products of the variances in the reaction norm gradients (here the environment and its squared value) and the corresponding additive genetic variance in the parameters (here $b_g$ and $c_g$ in Equation 12). This means that, in the quadratic case, genetic variances in slope and curvature directly translate into variance arising from plasticity, as they should. By contrast, $V_A$ does not solely depend on the
variance in the intercept $V_a$, but also on the quadratic coefficient, more specifically its covariance with the intercept.

The expressions for these variance components in the character-state approach are best described directly from the $G_z$ matrix. The total additive genetic variance along the reaction norm, $V_{Add}$, is the average of the additive genetic variance in each environment, i.e. the average of the diagonal elements of the $G_z$. The environment-blind additive genetic variance of the trait, $V_A$, is the average of all the elements of the $G_z$ matrix. Finally, the variance $V_{A×E}$ is the sum of the products of the (co)variances in the frequency of each environment and the additive genetic (co)variances in $G_z$. We illustrate in subsection C4 the relationship between the structure in the $G_z$ matrix and the additive genetic variances, but a simplified statement is that $V_{A×E} > 0$ as soon as the correlation between environments are different from 1 and/or variances in the diagonal are not all equal.

To further decompose genetic variation in the reaction norms, we first note that here, the reaction norm parameters are the focus of the decomposition, rather than shape characteristics like the slope or curvature (with the exception of a quadratic reaction norm, the only case were they are formally linked). Because Equation 21 is a sum of products, and since $G_θ$ is a constant, we can isolate each term of the resulting sum as:

$$ (26) \qquad \gamma_i = \frac{E_{\varepsilon}(\psi_{\varepsilon,i}^2) V_g(\theta_i)}{V_{Add}}, \qquad \gamma_{ij} = \frac{2E_{\varepsilon}(\psi_{\varepsilon,i}\psi_{\varepsilon,j}) \text{Cov}_g(\theta_i, \theta_j)}{V_{Add}}, \qquad \sum_i \gamma_i + \sum_{i<j} \gamma_{ij} = 1. $$

Here, $\gamma_i$ provides the contribution of the $i$th parameter in the model to the total additive genetic variance $V_{Add}$, while $\gamma_{ij}$ provides the contribution of the covariation between parameters $i$ and $j$ to $V_{Add}$. As such, this “$\gamma$-decomposition” (where gamma refers to $g$ for Genetics) measures the relative importance of genetic variances and covariances of the parameters to the evolvability of the plastic trait. Large values of $\gamma_i$ indicate that genetic variation in the $i$th parameter translate into a large proportion of the genetic variation in the trait. Also, large positive or negative values for $\gamma_{ij}$ indicate that covariation between parameters $i$ and $j$ can have a large impact in increasing or reducing genetic variation in the trait.

It is also possible to focus on the additive genetic variation arising from plasticity, $V_{A×E}$, which yields:

$$ (27) \qquad \iota_i = \frac{V(\psi_{\varepsilon,i}) V_g(\theta_i)}{V_{A×E}}, \qquad \iota_{ij} = \frac{2\text{Cov}_{\varepsilon}(\psi_{\varepsilon,i}, \psi_{\varepsilon,j}) \text{Cov}_g(\theta_i, \theta_j)}{V_{A×E}}, \qquad \sum_i \iota_i + \sum_{i<j} \iota_{ij} = 1. $$

This “$\iota$-decomposition” (where iota refers to $i$ for Interaction) highlights the fact that $V_{A×E}$ is the sum of the products of (co)variances in elements of the reaction norm gradient $\psi_\varepsilon$ and the additive genetic (co)variances in the parameters.

For a quadratic reaction norm as in Equation 12 with a mean-centred environment, this yields:

$$ (28) \qquad \gamma_a = \frac{V_a}{V_{Add}}, \quad \gamma_b = \frac{V_b E(\varepsilon^2)}{V_{Add}}, \quad \gamma_c = \frac{V_c E(\varepsilon^2)^2}{V_{Add}}, \quad \gamma_{ac} = \frac{2C_{ac} E(\varepsilon^2)}{V_{Add}}, \quad \iota_b = \frac{V_b V(\varepsilon)}{V_{A×E}}, \quad \iota_c = \frac{V_c V(\varepsilon^2)}{V_{A×E}}. $$

Note that since the environment has been mean-centred, we have $V(\varepsilon) = E(\varepsilon^2)$ since $E(\varepsilon)^2 = 0$, and thus $\gamma_b = \iota_b$, i.e. in the quadratic case, all of the genetic variation in the slope contributes to the genetic variance arising from plasticity. Note also that genetic variance in reaction norm intercept $a$ does not contribute to the heritability from plasticity ($\iota_a = 0$).

For the character-state approach, such decomposition would be less informative about the potential for (and constraints on) reaction norm evolution. Instead, we can define an effective number of character states (as proposed for general multivariate phenotypes by Kirkpatrick 2009) as

$$ (29) \qquad n_e = \sum_i \frac{\lambda_i}{\lambda_1}, $$

where $\lambda_i$ is the $i$th eigenvalue of $G_z$ ranked by size (i.e., $\lambda_1$ is the largest eigenvalue). Strong genetic correlations of phenotypes across environments lead to small $n_e$, whereby reaction norm
evolution is highly constrained (with the limit of $n_e = 1$ corresponding to the strongest constraint). Conversely, weak genetic correlations across environments leave more degrees of freedom for reaction norms to evolve, causing a large $n_e$, close to the actual number of assayed environments. This $n_e$ metric does not capture all aspects of reaction norm evolvability, and is best combined with the ratio $V_{A×E}/V_{Add}$ of the proportion of total genetic variance due to genetic variance in plasticity). Unfortunately, $n_e$ is estimated with a strong bias due to the overestimation of the leading eigenvalue of $G_z$ (Lawley 1956), making it less useful in practice than in theory. We thus do not develop this metric further.

## Parameter estimation and variance partitioning in practice

### Estimating the parameters

All the parameters mentioned in the previous section can be estimated through commonly used statistical frameworks. For the character-state approach (Equation 2), a random-parameter model can be used, or alternatively a “multi-trait” model (Mitchell and Houslay 2021; Rovelli et al. 2020). We will focus here on the former, which is more easily implemented while seemingly scarcely used in the literature on plasticity. In the random-parameter model, the environment is considered as a categorical variable, to which a random effect is added using the genotype as the grouping factor. In the curve-parameter approach, the appropriate models will be random-parameter models for a polynomial approach (as mentioned in Morrissey and Liefting 2016), or non-linear mixed models, fitting the reaction norm function $f(\varepsilon, \theta)$ to the data. Genotype-specific parameters, such as the intercept, slope, and any higher-order effects of a polynomial function, are treated as random'

Since the parameters are estimated with noise, it is important to account for the impact of estimation uncertainty when computing variance components. In particular, while variances directly obtained using random effects (e.g. genetic variances) are expected to be unbiased, the variances arising from fixed effects (e.g. variances related to $V_{Plas}$) should be corrected for biases due to uncertainty (as the adjusted $R^2$ does for example). Details are provided in Appendix E.

To compute the total phenotypic variance required to get the estimates $\hat{P}_{RN}^2$, $\hat{H}_{RN}^2$ and $\hat{h}_{RN}^2$, we advise using the sum of all estimated components rather the raw sample variance. The former is common practice in most quantitative genetics inference to account for potential imbalance in the experimental or sampling design (de Villemereuil et al. 2018; Wilson et al. 2010).

We provide an R package, named Reacnorm github.com/devillemereuil/Reacnorm, providing functions implementing the variance decomposition based on raw outputs of statistical models. A tutorial is shipped with the package, as an R vignette, showing how to implement such models using the Bayesian brms R packages (Bürkner 2018), along with Reacnorm.

### Perfect modelling of quadratic curves

We simulated phenotypic data conforming to a quadratic reaction norm, to evaluate the performance of the proposed approach when the reaction norm truly is quadratic. We considered both a discrete and continuous environment. For the discrete environment, we considered $N_{Gen} = 20$ or 5 different genotypes and an environmental gradient of $N_{Env} = 10$ or 4 values, equally spaced from -2 to 2. We sampled $N_{Rep} = N_{Gen}$ individual measures for each genotype within an environment. For the continuous environment, we drew $N_{Env} = 10$ or 4 values from a normal distribution for each of the $N_{Gen} = 200$ or 50 genotypes, without repeats contrary to the discrete case. In both cases, a residual noise was applied around each measure with a residual variance $V_{Res} = 0.25$. In all cases, we defined a quadratic curve with average parameters $\bar{\theta} = (1.5, 0.5, -0.5)$ for intercept, slope and curvature. We then drew $N_{Gen}$ different genotype-specific vectors of curve-parameter $\theta$ from a multivariate normal distribution with mean $\bar{\theta}$ and (genotypic) variance-covariance matrix

$$G_{\theta} = \begin{pmatrix} 0.090 & -0.024 & -0.012 \\ -0.024 & 0.160 & 0.008 \\ -0.012 & 0.008 & 0.040 \end{pmatrix}.$$
<figure><img src="image_7.png" /><figcaption>**Figure 4** – Distribution of the error (difference between the inferred and true value) for each of the three scenarios: two discrete ($N_{Env}$: number of environments, $N_{Gen}$: number of different genotypes, $N_{Rep}$: number of replicates per genotype) and one continuous ($N_{Env}$: number of environment tested per genotype, $N_{Gen}$: number of different genotypes). The grey dots correspond to the average over the 1000 simulations. The character-state approach was impossible for the continuous environment scenario. The yellow boxes on the right show the estimates for $\hat{P}_{RN}^2$ (proportion of variance generated by the plasticity in the mean reaction norm), $\hat{h}_{RN}^2$ (total heritability of the reaction norm), $\hat{h}^2$ (environment-blind heritability) and $\hat{h}_1^2$ (heritability from plasticity) for both the curve-parameter and character-state approaches. For the curve-parameter, the $\pi$-decomposition of $\hat{P}_{RN}^2$ into $\pi_{SI}$ (contribution of the slope) and $\pi_{CV}$ (contribution of the curvature); the $\gamma$-decomposition of $\hat{h}_{RN}^2$ into $\gamma_a$ (genetic contribution of the intercept), $\gamma_b$ (genetic contribution of the slope), $\gamma_c$ (genetic contribution of the curvature) and $\gamma_{ac}$ (genetic contribution of the covariance between the intercept and the curvature) and the $\iota$-decomposition of $\hat{h}_1^2$ into $\iota_b$ (slope) and $\iota_c$ (curvature) are shown.</figcaption></figure>
Figure 1 displays examples of curves resulting from these parameters. The simulation process was repeated 1000 times for each scenario, and for each simulated dataset, we ran estimations using the lme4 R package (Bates et al. 2015) under the curve-parameter (for discrete and continuous environment) and character-state (only for discrete environment) approaches, in order to check how these approaches compare in practice.

From the curve-parameter models, we computed $\hat{V}_{Plas}$ (accounting for the uncertainty in fixed effects), then $\hat{P}_{RN}^2$. We also computed the $\pi$-decomposition ($\hat{\pi}_{Sl}$ and $\hat{\pi}_{Cv}$, Equation 14), since the true reaction norm is quadratic here, as well as $\hat{h}_{RN}^2$, $\hat{h}^2$ and $\hat{h}_I^2$ as in Equation 9. We then applied the $\gamma$-decomposition to $\hat{h}_{RN}^2$ (Equation 26): $\hat{\gamma}_a$ (impact of the genetic variation of the intercept), $\hat{\gamma}_b$ (for the slope), $\hat{\gamma}_c$ (for of the curvature) and $\hat{\gamma}_{ac}$ (for the covariance between the intercept and curvature). Similarly, we applied the $\iota$-decomposition to $\hat{h}_I^2$ (Equation 27): $\hat{\iota}_b$ (for the slope) and $\hat{\iota}_c$ (for the curvature). From the character-state model, we computed only $\hat{P}_{RN}^2$, $\hat{h}_{RN}^2$, $\hat{h}^2$ and $\hat{h}_I^2$.

The yellow boxes in Figure 4 display the theoretical expected values for the different parameters for three scenarios of environmental variation (two discrete, one continuous; other scenarios are shown in Appendix F). Using the first discrete scenario as a reference for now, most of the total phenotypic variance comes from the average plasticity ($P_{RN}^2 = 0.55$). This, in turns, includes a large contribution from the curvature ($\pi_{Cv} = 0.56$) of the average reaction norm, more than from its slope ($\pi_{Sl} = 0.44$). The total heritability of the reaction norm is substantial ($h_{RN}^2 = 0.3$), but interestingly most of it is due to the heritability from plasticity ($h_I^2 = 0.21$), while the environment-blind heritability of the trait is only $h^2 = 0.08$. Contrary to the average shape, most of the additive genetic variation comes from the slope, both when considering the total reaction norm ($\gamma_b = 0.52$), or plasticity alone ($\hat{\iota}_b = 0.76$). All scenarios share the same underlying parameters $\theta$ and $G_\theta$, resulting in very comparable values for our variance decomposition (i.e. $P_{RN}^2$ and the heritabilities) across the different environmental sampling scheme. By contrast, the environmental sampling scheme (especially discrete v. continuous distribution) can substantially impact the expected values of the $\pi$-, $\gamma$- and $\iota$-decompositions. This is especially true when switching from the discrete to the continuous scenarios (e.g. $\pi_{Sl} = 0.44$ for the first discrete scenario while $\pi_{Sl} = 0.33$ for the continuous scenario).

Switching to the error in the estimation of the parameters (left panels of Figure 4), we see first that both the character-state and curve-parameter approaches allow for unbiased inference (Wilcoxon's rank test, $p > 0.05$), apart from a slight bias in the heritabilities ($\hat{h}_{RN}^2$, $\hat{h}^2$ and $\hat{h}_I^2$) and some of their $\gamma$ and $\iota$ components in the discrete scenarios ($< 5\%$ relative bias, Wilcoxon's rank test, $p < 0.05$), notably due to a slight overestimation of the genetic variance of the intercept (visible in the top row of Figure 4). For the discrete case, the precision of the estimates was not much influenced by the number of environments and depended more on the number of genotypes (see Figure S1). For the continuous case, both the number of environments and genotypes influenced the precision of estimates (see Figure S2). As a sanity check, we also verified that $\hat{V}_{Tot}$ (not shown in Figure 4) reflected the raw phenotypic variance with extreme precision (correlation > 99%) in the discrete case and very good precision (correlation > 87%) in the continuous case. The difference between these two types of scenarios is explained by how the stochasticity in environmental values differs among them. Importantly, the results in Figure 4) also illustrate the exact equivalence, in the discrete case, between the curve-parameter and character-state approaches, as the distributions of $\hat{P}_{RN}^2$ and $\hat{h}_{RN}^2$ were nearly identical (Figure 4, correlation > 99%) between the two approaches. This means that our variance partitioning is not impacted by which approach is chosen to study plasticity, as long as the curve-parameter approach captures the true reaction norm shape. When this does not hold, the differences between estimates from these alternative approaches can be exploited efficiently, as we describe below.

## Imperfect modelling of a non-polynomial reaction norm

The true shapes of reaction norms are generally unknown and may be complex, such that any curve-parameter model is likely to be mis-specified to some extent. In the case of a discrete
environment, the character-state approach is arguably more general, as it does not assume anything about the “true” shape of the reaction norm (as pointed out previously by de Jong 1995). Nonetheless, having access to curve-parameters is often very interesting and more actionable (even in cases where the linear and quadratic components cannot be interpreted as the average slope and curvature), especially to predict evolution of phenotypic plasticity (see also de Jong 1995). To get the best of both worlds, we rely on the ability of the character-state approach to recover $P_{RN}^2$, using it as an “anchor”, to assess the performance of a given curve. Note that, under these circumstances, it is not possible to obtain the most natural $\pi$-decomposition in Equation 14, so we instead rely on the $\varphi$-decomposition in Equation 16 (here taken at the second order). Because of this, we need to assess how “bad” our simplification using an imperfect curve is. To do so, we compute the ratio of the variance modelled by the polynomial curve to the total variance due to phenotypic plasticity:

$$ (30) \qquad M_{Plas}^2 = \frac{\hat{V}_{\text{mod}}}{\hat{V}_{\text{Plas}}}, $$

where both $\hat{V}_{\text{mod}}$ and $\hat{V}_{\text{Plas}}$ are bias-corrected. It is important to note here that $M_{Plas}^2$ is just a convenient way to quantify the amount of $\hat{V}_{\text{Plas}}$ explained by the chosen parametric curve, and should not be used to perform model selection. Model selection is a complex matter and we refer the readers to published reviews on this subject (e.g. Johnson and Omland 2004; Tredennick et al. 2021).

In order to demonstrate the soundness and usefulness of this approach, we simulated datasets following relatively common curves that are not well-captured by a second order polynomial: a logistic sigmoid (hereafter sigmoid scenario), or a Gompertz-Gaussian thermal performance curve (hereafter TPC scenario, see Figure 5). We assumed that the environment is sampled at either 10 or 4 values. For each of these conditions, we simulated 1000 datasets, with 10 measures per environment (for the sake of simplicity, and given the focus on $\hat{P}_{RN}^2$ here, we did not include different genotypes in these simulations). We estimated the parameters of a polynomial model, and computed the relative contributions of the first- and second-order parameters using Equation 16. In addition, we computed the unbiased estimates of the variance explained by our polynomial or character-state models to obtain $M_{Plas}^2$.

Our results show that, as expected, the polynomial function is an imperfect proxy of our complex shapes (Figure 5, $M_{Plas}^2 = 0.89$ for the sigmoid and $M_{Plas}^2 = 0.65$ for the TPC), but using the character-state approach allows retrieving the total plastic variance without bias. The approach described here is thus useful to compare a given reaction norm model (e.g. a polynomial function) to an unknown true shape of the reaction norm, in a case where environment is discretised. In more detail, the linear component was the most important component to explain the phenotypic variation for the sigmoid scenario ($\varphi_1 = 0.89$, same as the total model). This was because the quadratic component was always estimated close to zero ($< 10^{-3}$), thus no variance was explained by the quadratic component ($\varphi_2 = 0$). Of course, the sigmoid is not a straight line either, and some remaining variance unexplained by the polynomial curve ($1 - 0.89 = 0.11$) could have been explained by higher-order effects (e.g. cubic effect and higher). By contrast, for the TPC scenario, while the linear component was an important factor ($\varphi_1 = 0.47$), the quadratic component also explained quite a lot of the variance as well ($\varphi_2 = 0.2$). Again, higher-order effect, including at least a cubic effect, would have explained more of the variance arising from the average shape of plasticity.

This example illustrates the usefulness of a combined curve-parameter and character-state approach to study the shape of reaction norms of a discretely sampled environment. While the character-state approach provides a widely applicable estimation of $\hat{P}_{RN}^2$ (if the environment is discretised), the curve-parameter approach provides interpretable information about (at least) first- and second-order parameters of the reaction norm (although they might depart more or less strongly from its average slope and curvature), which helps describing where most phenotypic variance lies. Our ratio $M_{Plas}^2$ can then be used to evaluate how well a chosen polynomial function models an actual reaction norm.
### Sigmoid

$$f(\varepsilon) = \frac{L}{1 + e^{-r\varepsilon}} \quad \begin{array}{l} L = 1 \\ r = 4 \end{array}$$

<figure><img src="image_8.png" /></figure>

### Thermal Perf. Curve

$$f(\varepsilon) = C e^{-e^{\rho(\varepsilon-\varepsilon_0)^{-6}}-\sigma(\varepsilon-\varepsilon_0)^2} \quad \begin{array}{l} C = 1, \varepsilon_0 = 1.2, \rho = 8, \sigma = 0.4 \end{array}$$

<figure><img src="image_9.png" /><figcaption>Figure 5 – Estimation of the variance of the reaction norm when the true shape (sigmoid on the left, Gompertz-Gaussian performance curve on the right, red lines on top graphs) is unknown and approximated from a polynomial function. The estimated reaction norms using a polynomial function (black line, top graphs) only account for a part of the reaction norm shape, while the ANOVA estimation (pink dots, top graphs) fit the true shape more accurately. As a result, the model is expected to explain only a part $M_{Plas}^2$ of phenotypic variance due to plasticity. On the bottom rows, the error distribution are shown for $M_{Plas}^2$, $P_{Plas}^2$, $\varphi_1$ and $\varphi_2$ (grey dots are the average estimated values, black crosses are the expected true values).</figcaption></figure>
Estimation of non-linear models

Although we have focused so far on models that are linear in its parameters, the main strength of our approach is its generality: it can be applied to any arbitrary functions (provided it is differentiable). This requires numerically computing integrals for $V_{Plas}$ (for $\hat{P}_{RN}^2$), $\pi_{Sl}$, $\pi_{Cv}$ and $\psi_{\epsilon}$ (for the heritabilities), but this can be solved with efficient algorithms. We illustrate this by introducing genetic variation in the parameters of the sigmoid and TPC reaction norms illustrated in Figure 5 (top panels). We used a non-zero, but small, residual variance ($V_R = 0.0001$) to avoid numerical issues typical when running thousands of non-linear models. We focused on a continuous environment, and estimated the actual functions used to generate the datasets, using the non-linear modelling function of nlme package (Pinheiro et al. 2009). We used the cubature package (Narasimhan et al. 2023), as in the QGglmm package (de Villemereuil et al. 2016), to compute parameters linked to the variance decomposition, and, further, the $\pi^-$, $\gamma^-$ and $\iota^-$ decomposition. We simulated 1000 datasets for each scenario, consisting of 200 genotypes measured each in 10 different environments, randomly sampled from a normal distribution.

We retrieved our simulated parameters without bias using the nlme function, except for a slight bias (Wilcoxon's rank test, $p < 0.05$) in the variance of $r$ (latent slope) in the sigmoid model and in $C$ (height of the peak) in the TPC model. This translated into significant (Wilcoxon's rank test, $p < 0.05$), but very limited bias (relative bias < 5%) in our derived parameters (Figure 6, bottom panels). Moreover, the sum of variance components ($\hat{V}_{Tot}$) successfully reflects the total phenotypic variance, with a correlation between the two quantities > 91%.

First focusing on the average shape of the reaction norm (Figure 6, top panel), one unfortunate aspect of running a non-linear model is that our bias correction described in Appendix E can no longer be applied. However, this bias is generally small provided the standard error is small for most parameters, and the resulting bias in $\hat{P}_{RN}^2$ is extremely small, and even non-significant for the sigmoid model. This could of course be partly explained by a favourable context here, especially since the residual variance is relatively small. An important distinction here is the difference between the curve defined by the average parameters $f(\epsilon, \bar{\theta})$ (Figure 6, top panel, black curve) and the one defined by the local average phenotype $E_{g|\epsilon}(\hat{z})$ (Figure 6, top panel, red curve), recalling that $\hat{P}_{RN}^2$ is linked to the latter. While the two are very close for the sigmoid case, they differ quite visibly for the TPC one, due to a more pronounced non-linearity in the parameters in the latter. The average slope contributed the most to the overall plastic variance of the mean reaction norm for the sigmoid shape ($\pi_{Sl} = 0.88$), with no impact of average curvature ($\pi_{Cv} = 0$), close to the $\phi$-decomposition in Figure 5. For the TPC scenario, the contribution of the average slope ($\pi_{Sl} = 0.31$) and curvature ($\pi_{Cv} = 0.35$) are similar. In this case, the values are very different from the $\phi$-decomposition in Figure 5 (although note that the distribution of the environment is different between these two scenarios). It might appear as counter-intuitive that the slope contributes so much to variance, since the curve increases from 0 and then decreases toward 0, but this is linked to the fact that the environment is normally distributed, so most values are near $\epsilon = 0$, an area where the slope of the curve is close to being maximised.

Although the variation between genotypes in the top panel of Figure 6 seems quite large, the contribution from the average plasticity $\hat{P}_{RN}^2$ is 1.7 to 3.4 times higher than the one of the genetic variance $\hat{H}_{RN}^2$ (Figure 6, yellow box in first- and second-row panels). This occurs because the genetic variance is actually very low in most environments (Figure 6, brown and purple lines of the second-row panels), and scarcely as high as $V_{Plas}$. As mentioned above, non-linearity in the parameters is less strong for the sigmoid case than for the TPC case, resulting in almost exactly equal values for $\hat{H}_{RN}^2$ and $\hat{h}_{RN}^2$ for the former, while they are slightly different for the latter. In both cases, the small difference between $\hat{H}_{RN}^2$ and $\hat{h}_{RN}^2$ can be explained by the disproportionate importance in the $\gamma$-decomposition of parameters that are actually linearly related to the trait ($\gamma_L = 0.98$ for the sigmoid and $\gamma_C = 0.81$ for the TPC scenarios). In terms of heritability from plasticity, it is substantial in both cases ($h_1^2 = 0.081$ for the sigmoid and $h_1^2 = 0.133$ for the TPC scenario), as can be expected from the non-parallel reaction norms (Figure 6). However, it remains smaller than the environment-blind heritability of the trait in both cases ($h^2 = 0.143$ for the sigmoid and $h^2 = 0.216$ for the TPC scenarios). Interestingly, for the TPC scenario, and
Sigmoid

$$f(\varepsilon) = \frac{L}{1 + e^{-r\varepsilon}}$$

<figure><img src="image_10.png" /></figure>

Thermal Perf. Curve

$$f(\varepsilon) = C e^{-e^{\rho(\varepsilon-\varepsilon_0)-6-\sigma(\varepsilon-\varepsilon_0)^2}}$$

<figure><img src="image_11.png" /><figcaption>Figure 6</figcaption></figure>
**Figure 6 - Caption** Scenarios and results of non-linear modelling of phenotypic plasticity in a continuous environment. On the left: results corresponding to a sigmoid curve scenario; on the right: results corresponding to a TPC scenario. First row: example of the individual curves (each curve corresponds to one individual) simulated in each scenario; yellow box: true parameters for the model and average shape; black curve: $f(\epsilon, \theta)$; red curve: $E_{g|\epsilon}(\hat{z})$. Second row: distribution of the estimations of $V_{G,\epsilon}$ (brown) and $V_{A,\epsilon}$ (purple), along the environment; solid line: average value across simulations; pale ribbon: 95% CI across simulations; yellow box: true values for the genetic variance partition. Third row: $\gamma$-decomposition of $V_{A,\epsilon}$ along the environment, for each parameter and their covariation. Fourth row: distribution of the error for each component of our variance partition, grey dot is the average of estimates over all simulations.

contrary to what happens with the $\gamma$-decomposition, a majority of the additive genetic variance arising from plasticity comes from the variation in the location of the optimum ($\ell_{\epsilon_0} = 0.525$). This is because variation in the location of the optimum shifts the reaction norm along the environment axis (i.e. on the “x-axis”), meaning that even a small shift can generate considerable variation that is non-parallel along the phenotype axis (i.e. along the “y-axis”).

An interesting aspect of our framework is that we can explore the variation of $V_{Gen,\epsilon}$, $V_{A,\epsilon}$ and the $\gamma$-decomposition of $V_{A,\epsilon}$ along the environmental gradient, which can be very informative from an evolutionary perspective. In the case of the sigmoid curve (Figure 6, second and third rows, left panels), the analysis is relatively simple: as the value of the environment increases, the parameter $L$ is multiplied by an increased value (going from 0 to 1 due to the sigmoid function) and thus its genetic variance plays a stronger role. This translates into $V_{Gen,\epsilon}$ and $V_{A,\epsilon}$ increasing with the environment, and $\gamma_L$ accounting for almost all of the genetic variance after the sigmoid inflexion point in 0. The TPC scenario is even more interesting. First, we can see that both $V_{Gen,\epsilon}$ and $V_{A,\epsilon}$ (Figure 6, second row, right panels) are close to zero in the extreme environments and maximised in a region between the optimum and critical maximal temperature, where the reaction norm suddenly drops after the optimum. This maximum also corresponds to the region where $V_{Gen,\epsilon}$ and $V_{A,\epsilon}$ are the most different (and where the red and black departs the most in Figure 6, top row, right panel). Regarding the $\gamma$-decomposition (Figure 6, third row, right panels), the influence of the location of the optimum ($\gamma_{\epsilon_0}$) is maximised at extreme environments, while the influence of the maximum value at the peak ($\gamma_C$) is exactly maximised at the average location of the peak. The influence of the covariation between both ($\gamma_{C\epsilon_0}$) is negative before the peak and positive after.

As these simulations illustrate, our framework allows very finely describing the characteristics of reaction norms, such as how its average shape (slope/curvature) and genetic variation in the parameters influence the phenotypic variance in the trait, while discriminating between total genetic variation of the trait and genetic variation exclusively linked with plasticity itself.

## Discussion

The variance decomposition in Equation 7 is very general, and applicable to any approach used to estimate a reaction norm. In particular, it applies equally well to both the character-state and curve-parameter approaches. Each component and its variance-standardisation provide a different information on the reaction norms: $P_{RN}^2$ quantifies the proportion of phenotypic variance due to the average plastic response across genotypes, while $H_{RN}^2$ or $h_{RN}^2$ quantify the contributions from (broad or additive) genetic variance in the reaction norms. Further, these genetic components can be separated into the environment-blind heritability of the trait ($h^2$) based on the average breeding values across environments, and the heritability from plasticity ($h_1^2$) which is solely based on the gene-by-environment interactions at the level of breeding values. Finally, the sum $T_{RN}^2 = P_{RN}^2 + H_{RN}^2$ quantifies how well we can predict the individual phenotypes based on their genotypes and environments (i.e. genetically variable reaction norms). Those components are efficient summary statistics yielding important information regarding the evolutionary potential of both the trait and its plasticity. Importantly, they are very generally applicable, with
a strict equivalence between e.g. a character-state or a curve-parameter approach. However, they do not provide information regarding the actual shape of the reaction norms. To that end, we further decomposed some of these components in terms of characteristics of the shape or parameters of reaction norms.

The most difficult problem is to decompose the average plastic variance $P_{RN}^2$ into terms arising either from the linear trend ($\pi_{SI}$) or from the curvature ($\pi_{CV}$) of the reaction norm, which we called $\pi$-decomposition. Unfortunately, our estimates for $\pi_{SI}$ and $\pi_{CV}$ are only valid if the environment is normally distributed, or the true reaction norm is quadratic. In other cases, mean slope and curvature loose their simple interpretation, preventing a meaningful $\pi$-decomposition. Nonetheless, for polynomial reaction norms of higher order, we described an alternative decomposition, based on the polynomial coefficients rather than actual slope and curvature, which we called $\varphi$-decomposition. While not as interpretable as the $\pi$-decomposition, this decomposition can serve as a way to compare polynomial shapes across contexts. Based on the equivalence between the curve-parameter and character-state, we introduced $M_{Plas}^2$ as a way to quantify the ability of a polynomial model to recover $V_{Plas}$ compared to an "agnostic" model such as the character-state. Our proposed framework is summarised in Figure 3.

Decomposing $h_{RN}^2$ and $h_I^2$ is comparatively easier, because the model assumed in Equation 3 and Equation 4 ensures that we can always translate additive genetic variance in the parameters $\theta$ into additive genetic variance in the trait $z$, even if the function $f$ is not linear in its parameters. Decomposition of the total heritability of the reaction norm $h_{RN}^2$ into the impact of the parameters $\theta$ leads to the $\gamma$-decomposition. It quantifies the relative importance of genetic variance in different reaction norm parameters to the evolvability of the trait. For instance if a given selection episode concerns individuals that all experienced the same plasticity-inducing environment (i.e. when spatial environmental variation is negligible relative to temporal variation), using the multivariate breeder's equation (Lande 1979), the relative contribution of genetic variation in parameter $\theta_i$ to the response to selection for the trait $z$ is

$$ (31) \qquad \frac{\Delta\theta_i\bar{z}}{\Delta\bar{z}} = \gamma_i + \frac{1}{2}\sum_{i \neq j} \gamma_{ij}, $$

where the $\gamma_i$ and $\gamma_{ij}$ are defined in Equation 26. In other words, the contributions of responses to selection by different reaction norm parameters to overall response to selection by the plastic trait $z$ is directly proportional to their contribution to its genetic variance. Importantly, these contributions will depend on the reaction norm gradient $\psi_{\varepsilon}$ defined in Equation 19, and thus on the environment, as illustrated in Equation 26. In fact, the environment-specific additive genetic variance $V_{\Lambda,\varepsilon}$ is a critical piece of information regarding evolutionary potential, and we can apply the $\gamma$-decomposition within each environment as well. For example, in the TPC scenario investigated above (Figure 6, right panels), the contribution of the peak height parameter $C$ is maximised at the average location of the optimum, where it accounts for 100% of the additive genetic variance. On the contrary, the influence of additive genetic variation in the location of the optimum $\varepsilon_0$ is more important in extreme environments. The complex interaction between the role of $C$ and $\varepsilon_0$ generates a peak for $V_{\Lambda,\varepsilon}$ in the area between the peak and critical maximal value for the environment (where the performance curve reaches zero). In the context of predicting eco-evolutionary response to warming, this would mean that a slight temperature rise above the optimum would provide a very short window of higher evolvability, but followed by a sharp decrease thereof if warming persists. Beyond these simple scenarios, how selection acts on reaction norms and plasticity depends on how the environment varies in space and/or time (de Jong 1999; King and Hadfield 2019; Scheiner 1993b; Tufto 2015), and how the reaction norm gradient $\psi_{\varepsilon}$ and direction selection on the expressed trait $z$ covary across environments. However, an in-depth exploration of how to estimate these selection responses is beyond the scope of the present work.

While the $\gamma$-decomposition is key to understanding and predicting evolution of the trait, it is based on the total heritability of the reaction norm $h_{RN}^2$, which combines additive genetic variation in the trait and its plasticity. To study plasticity in isolation from the environment-blind additive genetic variance in the trait, we decomposed $h_I^2$ in a similar fashion as $h_{RN}^2$, which we
called the *i*-decomposition. The components of the *i*-decomposition measure the contribution of each parameter to the evolutionary potential of plasticity, i.e. to the evolvability of reaction norm shape. In our thermal performance case (TPC) example, the *i* associated to *C* and ε₀ were close to 0.5, meaning that evolution can roughly equally impact the peak height *C* or the location of the optimum ε₀, should selection on the shape of reaction norms occur.

The detailed decomposition that we propose open the door to better comparability across studies, which can be a challenge in meta-analyses of plasticity. Murren et al. (2014) performed such a meta-analysis, comparing genetic variation in different parameters of reaction norm shape across published datasets. However they (i) computed these parameters using only extreme environmental values, instead of the whole range of environments; (ii) did not account for uneven spacing between environments where relevant; (iii) did not account for uncertainty in estimations of reaction norms (as previously highlighted by Morrissey and Liefting 2016); and (iv) assumed the modeled reaction norm shape is true. More details about the analyses in that study are provided in Appendix G. Our approach overcomes all these issues (some of which had been dealt with already by Morrissey and Liefting 2016; Pélabon et al. 2020). Unfortunately the dataset compiled by Murren et al. (2014) does not provide information on uncertainty of phenotypic estimates (related to V<sub>Res</sub>), precluding proper meta-analysis of reaction norm shape variation.

Importantly, our variance partitioning can be implemented through commonly used statistical models, notably (non-)linear mixed models. We showed that even complex non-linear modelling can perform well, only at the cost of using dedicated libraries to compute integrals numerically. This means that biologists can readily seize all the modelling tools introduced here. In particular, although a character-state approach can be performed using a simple random-intercept model, studies of genetic variance in plasticity seem to rather use a multi-trait model, which offers more control, but is more difficult to implement (but see Stirling and Roff 2000). In order to make the variance partitioning introduced here more accessible, we have implemented the computation of all the decomposition mentioned here as an R package named Reacnorm github.com/devillemereuil/Reacnorm, including cases where more than the genetic effect is assumed affecting variation in *θ*. The package also provides a tutorial as a vignette, showing how to implement the models in the Bayesian package brms and use functions from Reacnorm to study the properties of reaction norms. We hope that this will further stimulate interest in investigating variation and evolutionary potential of reaction norms.
# Appendix A. A unified formalism for the curve-parameters and character-state approaches

Despite having different mechanics, the curve-parameter and character-state approaches can be shown to be mathematically equivalent de Jong (1995). We can use this to express both approaches under the same, unified formalism. More precisely, we can express the character-state approach as being a special case of the curve-parameters approach. Under a curve-parameters approach, the reaction norm is seen as a function f of the environment ε and a vector of parameters θg:

$$ (S1) \qquad \hat{z} = f(\varepsilon, \theta_g). $$

The θg's covary across genotypes with a variance-covariance matrix Gθ:

$$ (S2) \qquad \theta_g \sim \mathcal{N}(\bar{\theta}, G_{\theta}). $$

By contrast, in a character-state approach, the reaction norm values of different genotypes across environments are directly provided by sampling from a multivariate normal distribution:

$$ (S3) \qquad \hat{z} \sim \mathcal{N}(\mu, G_z). $$

One way to express the character-state using the same formalism as the curve-parameter is to recognise that Equation S3 can be written as

$$ (S4) \qquad \begin{aligned} \hat{z} &= \mu_g^T u_k, \\ \mu_g &\sim \mathcal{N}(\mu, G_z), \end{aligned} $$

where u_k is the unit vector with 1 at the kth value (corresponding to environment ε_k) and 0 elsewhere. Thus, the character-state model can be expressed using the formalism of Equation S1 and Equation S2, where μ_g in Equation S4 plays the role of θ_g, and thus G_z plays the role of G_θ. In this case, the function f is a function taking the level k of the environment and the parameters μ_g of the genotype g as input, and yielding the evaluated reaction norm ˆz as the output. Evidently, this function f is not continuous and not differentiable along the (categorical) environment. However, it is a continuous, differentiable and even linear function along the (continuous) parameters μ_g. As such, all properties mentioned in the main text and the Appendices pertaining to reaction norms that are “linear in its parameters” also apply to the character-state approach.

# Appendix B. Computation of the additive genetic variance holding environment constant

## B1. Preliminary results

Multiple regression slopes expressed using a variance-covariance matrix. Let us assume a multiple regression between a random variable y and a set of random variables x = (x₁, ..., xₙ)ᵀ such that:

$$ (S5) \qquad y = \mu + x^T \beta + e, $$

where μ is the intercept and e is the residual of the model. Note that in practical regression, the realised sampling of x will be contained in the design matrix of the model. If it exists and is unique, the solution for the vector of multiple regression slopes β can be formulated in terms variance-covariance matrices (see e.g. p.179, Lynch and Walsh 1998):

$$ (S6) \qquad \beta = V(x)^{-1} \text{cov}(x, y), $$

where V(x) is the variance-covariance matrix of x, V(x)⁻¹ is its inverse matrix and cov(x, y) is the column-vector of covariances between the xᵢ and y.
Multivariate version of Stein's lemma. Let us assume that $x = (x_1, \dots, x_{p_x})$ and $y = (y_1, \dots, y_{p_y})$ follow multivariate normal distributions, and that $g$ is a differentiable, $R^{p_x} \to R$ function such that $E(\nabla g)$, where $\nabla g$ is the gradient of $g$ (the vector of partial derivatives), is a vector with finite values, then it can be shown (Landsman and Nešlehová 2008; Landsman et al. 2013) that:

$$ (S7) \qquad \operatorname{cov}(g(\mathbf{x}), \mathbf{y}) = \operatorname{cov}(\mathbf{x}, \mathbf{y})E(\nabla g). $$

Note that covariance matrices of vectors (also known as cross-covariance matrices) are not com-mutative, but are such that $\operatorname{cov}(\mathbf{x}, \mathbf{y}) = \operatorname{cov}(\mathbf{y}, \mathbf{x})^T$. In the case where $p_y = 1$, then $\mathbf{y} = \mathbf{y}$ follows a normal distribution and:

$$ (S8) \qquad \operatorname{cov}(g(\mathbf{x}), \mathbf{y}) = \operatorname{cov}(\mathbf{y}, \mathbf{x})E(\nabla g). $$

Note that $\operatorname{cov}(\mathbf{y}, \mathbf{x})$ is a row-vector and $\operatorname{cov}(\mathbf{x}, \mathbf{y})$ is a column-vector by convention.

## B2. Breeding values in a given environment

Genetics of reaction norms. As mentioned in the main text, a general formalism (including the character-state as a special case) for the reaction norm $\hat{z}$ is given by Equation 3 in the main text, i.e.

$$ (S9) \qquad \hat{z} = f(\varepsilon, \theta_g). $$

The phenotype predicted by the reaction norm $\hat{z}$ thus depends on the environmental value $\varepsilon$, and the reaction norm parameters $\theta_g$ specific to the genotype $g$. When holding the environment $\varepsilon$ constant, the genetic variance is simply the variance of reaction norms across genotypes:

$$ (S10) \qquad V_{G|\varepsilon} = V_{g|\varepsilon}(f(\varepsilon, \theta_g)) $$

If the reaction norms are estimated in such a way that non-additive genetic variance can be separated out from additive genetic variance (e.g. if “genotype” refers to individuals) or are known to be negligible on the one hand; and if the reaction norm is linear in its parameters (i.e. $f$ is a linear function of $\theta_g$, as for a polynomial function) on the other hand, then the additive genetic variance conditional on the environment is readily given by Equation S10, i.e. $V_{A|\varepsilon} = V_{G|\varepsilon}$. In the case where $f$ is not linear in its parameters, it is necessary to rely on the theory in non-linear quantitative genetics (de Villemereuil et al. 2016; Morrissey 2015), as we do below.

Linear relationship between breeding values. The relationship between the breeding value of the trait $A_z$ and the breeding values of the reaction norm parameters $\theta_g$ is the key towards developing a framework that works for any reaction norm, linear in its parameters or not. Let us note $A_\theta$ the vector of breeding values of all the parameters in $\theta$. We will follow the same demonstration as in de Villemereuil et al. (2016), which starts from the point that, by definition, breeding values are all linked through linear relationships (see also Robertson 1966), since they are all linearly linked to the genotype (Lynch and Walsh 1998). More precisely, the breeding value $A_z$ of the phenotypic trait $z$ of an individual linearly depends on a linear combination of its breeding values for the reaction norm parameters $A_\theta$, so that:

$$ (S11) \qquad A_z = \mu_A + A_\theta^T \psi $$

where $\mu_A$ is a constant chosen such that $E(A_z) = 0$, $\psi$ is a vector of slopes that we will shortly describe as the reaction norm gradient.

Derivation of $\psi$. To derive an expression of $\psi$, we can apply the results in Equation S6 to Equation S11, yielding

$$ (S12) \qquad \psi = G_\theta^{-1} \operatorname{cov}(A_\theta, \hat{z}). $$

This assumes that $\operatorname{cov}(A_\theta, A_z) = \operatorname{cov}(A_\theta, \hat{z})$, i.e. that there is no covariance between the environmental values of the phenotype as predicted by the reaction norm and the breeding values of the parameters. This results also assumes that $G_\theta$ is invertible. However, such assumption is already necessary to most statistical algorithms available to infer $G_\theta$ in practice, so that this assumption is not limiting here. Noting that $\hat{z} = f(\varepsilon, \theta)$, we can apply the multivariate version of Stein's lemma (Equation S7):

$$ (S13) \qquad \psi = G_\theta^{-1} \operatorname{cov}(A_\theta, \theta_g) E(\nabla_\theta f) = G_\theta^{-1} G_\theta E(\nabla_\theta f) = E(\nabla_\theta f), $$
where we have used the fact that the covariance of breeding values of reaction norm parameters with their breeding values is their additive genetic covariance matrix $G_θ$. Again, note that this assumes that $f$ is partially differentiable with respect to all elements of $θ_g$. Given that this demonstration was applied when holding the environment constant, the values in $ψ$ generally depend on the environment $ε$, so below and in the main text, we use the notation $ψ_ε$.

Values of $ψ_ε$ in specific contexts. When the reaction norm is linear in its parameters, the values in $ψ_ε$ are (trivially) the linear coefficients of such relation. For a quadratic reaction norm, where $\hat{z} = (\bar{A} + a_g) + (\bar{b} + b_g)ε + (\bar{c} + c_g)ε^2$, such linear coefficients are respectively 1, $ε$ and $ε^2$ for $a_g$, $b_g$ and $c_g$. It results that $ψ_ε = (1, ε, ε^2)^T$ as mentioned in the main text. More generally, if $f$ is a polynomial of order $N$, then $ψ_ε = (1, ε, ..., ε^N)^T$. In the context of a character-state, it can be seen from Equation S4 that the gradient $ψ_ε$ in the parameters will be equal to $u_k$, i.e. a vector of 1 for the $k$th value (corresponding to the environment chosen to be hold constant) and 0 elsewhere.

### B3. Additive genetic variance

By definition, the additive genetic variance of the trait conditional on the environment $V_{A|ε}$ is the variance of the breeding values defined in Equation S11. We can thus express it from the breeding values of the reaction norm parameters (right hand side of Equation S11) as

$$ (S14) \qquad V_{A|ε} = V_{g|ε}(\mathcal{A}_{\theta}^{T} \psi_{ε}) = \psi_{ε}^{T} G_{\theta} \psi_{ε}. $$

This formula holds whether the reaction norm is linear on its parameters or not, and also holds for the character-state approach (although in this case, this formula merely selects the $k$th element of the diagonal of $G_z$).

## Appendix C. Derivation of the general decomposition of variance

### C1. Distinguishing between $V_{Plas}$, $V_{Gen}$ and $V_{Add}$

The phenotype predicted by the reaction norm $\hat{z}$ depends on the environment, and the reaction norm parameters $θ_g$ specific to the genotype $g$. The impacts of environment and genotype are intricately related via the reaction norm shape, but in a given environment, one can still isolate the average impact of the environment from variation among genotypes by computing the average value of the reaction norm across genotypes conditional on the environment, i.e. $E_{g|ε}(\hat{z})$. The variance of $E_{g|ε}(\hat{z})$, taken across environments, is the component $V_{Plas} = V(E_{g|ε}(\hat{z}))$ in the main text, i.e. the phenotypic variance arising from plasticity after averaging across genotypes. The genotypic value $G_z$ of genotype $g$ within the environment $ε$ is then given by

$$ (S15) \qquad G_z = \hat{z} - E_{g|ε}(\hat{z}). $$

Note that, although we removed the average effect of the environment, the genotypic value $G_z$ still depends on both the genotype $g$ and the environment $ε$, because genotypes can vary in their response to the environment. The total genetic variance in the reaction norm is thus $V_{Gen} = V(G_z)$. It is possible to get to the breeding values of the trait in each environment $A_z$ following the process described in Appendix B, i.e. $A_z = μ_a + \mathcal{A}_{\theta}^T ψ_ε$. The total additive genetic variance in the reaction norm is then

$$ (S16) \qquad V_{Add} = V(A_z) = E(V_{g|ε}(A_z)) + V(E_{g|ε}(A_z)) = E(ψ_{ε}^T G_{θ} ψ_{ε}), $$

using the law to total variance and noting that $E_{g|ε}(A_z) = 0$ by construction. In Figure 1 in the main text, the average $E_{g|ε}(\hat{z})$ corresponds to the red line in the left panel of Figure 1 in the main text, while $A_z$ corresponds to the purple lines in the middle panel.

### C2. Distinguishing between $V_{Add}$, $V_A$ and $V_{A×E}$

We can separate the total additive genetic variance of the reaction norm, $V_{Add}$, into two components: the environment-blind additive genetic variance of the trait $V_A$ and the additive
genetic variance arising from plasticity $V_{A×E}$. The first component is given by considering, for a given genotype, its average breeding value across environment:

$$ (S17) \qquad \bar{A} = E_{\varepsilon|g}(A_z). $$

This average corresponds to the breeding value that would be predicted for the same genotype present in all environments (or moving across them, being measured several times), ignoring the impact of the environment. In other words, this average is the predicted breeding value after the impact of the environment has been marginalised. Graphically, it depicts the average shift in the y-axis of the reaction norm, as can be seen in the middle panel of Figure 1 in the main text. The environment-blind additive genetic variance of the trait is

$$ (S18) \qquad V_A = V(\bar{A}) = E(\psi_\varepsilon)^T G_\theta E(\psi_\varepsilon) $$

$V_A$ is here defined as a variance, but there are negative elements in $E(\psi_\varepsilon)$ and $G_\theta$, so in the theory, their product could happen to be a negative scalar. This is not so here, because $G_\theta$ being a variance-covariance matrix, it must be positive semi-definite. By definition of positive semi-definiteness, the product $E(\psi_\varepsilon)^T G_\theta E(\psi_\varepsilon)$ will be positive (or null) for any real vector $E(\psi_\varepsilon)$.

The remaining additive genetic variation after accounting for the marginal breeding value is linked to the impact of genetic variation arising from plasticity, i.e. genotype-by-environment interactions. We can define the part of the breeding values strictly linked to that genotype-by-environment interaction by mean-centring the breeding values, for each genotype:

$$ (S19) \qquad A_I = A_z - \bar{A}. $$

The right panel of Figure 1 depicts these interaction breeding values. The additive genetic variance linked to genotype-by-environment, and thus to variation arising from plasticity, is:

$$ (S20) \qquad V_{A×E} = V(A_I) = V(A_z) + V(\bar{A}) - 2\operatorname{cov}(A_z, \bar{A}) = V(A_z) - V(\bar{A}) = V_{\text{Add}} - V_A, $$

noting that, by construction, $\operatorname{cov}(A_z, \bar{A}) = \operatorname{cov}(\bar{A}, \bar{A}) = V(\bar{A})$. By substituting $V_{\text{Add}}$ and $V_A$ with their values in Equation S16 and Equation S18, we obtain

$$ (S21) \qquad V_{A×E} = E(\psi_{\varepsilon}^T G_{\theta} \psi_{\varepsilon}) - E(\psi_{\varepsilon})^T G_{\theta} E(\psi_{\varepsilon}) = \operatorname{tr}(\Psi G_{\theta}) = \sum_{l,k} \Psi_{l,k} G_{\theta(l,k)}, $$

where $\Psi$ is the variance-covariance matrix of the reaction norm gradient $\psi_\varepsilon$ across the environment. In other words, $V_{A×E}$ is the sum of the products, for all pairs of parameters, of the (co)variance in the reaction norm gradient and the additive genetic (co)variance. The $\gamma$- and $\iota$-decomposition directly comes from dividing each elements of the sums in Equation S16 and Equation S21 respectively by $V_{\text{Add}}$ and $V_{A×E}$, so that the total sums to 1.

### C3. Variance decomposition for a polynomial model

In this section, we will assume a polynomial reaction norm:

$$ (S22) \qquad \hat{z} = \sum_{n=0}^{N} (\bar{\theta}_n + \theta_{n,g}) \varepsilon^n $$

where $\theta_n = \bar{\theta}_n + \theta_{n,g}$ is the nth order coefficient of the polynomial. In this form, it is easy to remark that polynomial reaction norms are linear in their parameters, i.e. there is a linear relationship between the $\theta_n$'s and $\hat{z}$, so that $G_z = A_z$. It results that:

$$ (S23) \qquad G_z = A_z = \hat{z} - E_{g|\varepsilon}(\hat{z}) = \sum_{n=0}^{N} (\bar{\theta}_n + \theta_{n,g}) \varepsilon^n - \sum_{n=0}^{N} \bar{\theta}_n \varepsilon^n = \sum_{n=0}^{N} \theta_{n,g} \varepsilon^n. $$

Taking the derivative of this expression with respect to each of $\theta_{n,g}$ in a given environment $\varepsilon$ would yield a reaction norm gradient equal to the value of each exponent of $\varepsilon$, i.e. $\psi_\varepsilon = (1, \varepsilon, \dots, \varepsilon^N)^T$. The total (additive) genetic variance is thus:

$$ (S24) \qquad V_{\text{Gen}} = V_{\text{Add}} = E(\psi_{\varepsilon}^T G_{\theta} \psi_{\varepsilon}) = \sum_{n} V_n E(\varepsilon^{2n}) + 2 \sum_{n<m} C_{nm} E(\varepsilon^{n+m}), $$
where $V_n$ is the additive genetic variance for $\theta_{n,g}$ and $C_{nm}$ is the additive genetic covariance between $\theta_{m,g}$ and $\theta_{n,g}$. For the quadratic case, if $\varepsilon$ has been mean-centred and is symmetrical, we have $E(\varepsilon) = E(\varepsilon^3) = 0$ and the expression reduces to

$$ (S25) \qquad V_{\text{Gen}} = V_{\text{Add}} = V_0 + (V_1 + C_{03})E(\varepsilon^2) + V_3E(\varepsilon^4). $$

For a given genotype, its average breeding value across environments is

$$ (S26) \qquad \bar{A} = E_{\varepsilon|g}(A_z) = E_{\varepsilon|g} \left( \sum_{n=0}^{N} \theta_{n,g} \varepsilon^n \right) = \sum_{n=0}^{N} \theta_{n,g} E(\varepsilon^n) $$

The environment-blind (additive) genetic variance of the trait is

$$ (S27) \qquad V_G = V_A = E(\psi_{\varepsilon})^T G_{\theta} E(\psi_{\varepsilon}) = \sum_n V_n E(\varepsilon^n)^2 + 2 \sum_{n<m} C_{nm} E(\varepsilon^n) E(\varepsilon^m) $$

For the quadratic case with mean-centred and symmetrical $\varepsilon$, this yields:

$$ (S28) \qquad V_A = V_0 + 2C_{02}E(\varepsilon^2) + V_2E(\varepsilon^4)^2 $$

Finally, the additive genetic variance arising from plasticity itself is

$$ (S29) \qquad V_{A \times E} = V_{\text{Add}} - V_A = \sum_n V_n E(\varepsilon^{2n}) + 2 \sum_{n<m} C_{nm} E(\varepsilon^{n+m}) - \sum_n V_n E(\varepsilon^n)^2 + 2 \sum_{n<m} C_{nm} E(\varepsilon^n) E(\varepsilon^m). $$

By recognising that $V(\varepsilon^n) = E(\varepsilon^{2n}) - E(\varepsilon^n)^2$ and $\text{cov}(\varepsilon^n, \varepsilon^m) = E(\varepsilon^{n+m}) - E(\varepsilon^n)E(\varepsilon^m)$, we can further simplify this expression as:

$$ (S30) \qquad V_{A \times E} = \sum_n V_n V(\varepsilon^n) + 2 \sum_{lk} C_{nm} \text{cov}(\varepsilon^n, \varepsilon^m). $$

For the quadratic case, for a mean-centred and symmetrical $\varepsilon$, all the covariances between the different exponents of $\varepsilon$ are 0, yielding

$$ (S31) \qquad V_{A \times E} = V_1 V(\varepsilon) + V_2 V(\varepsilon^2). $$

### C4. Variance decomposition for the character-state approach

As mentioned in Appendix A, the character-state can be written using a function $f$ such that in environment $\varepsilon_k$ and for genotype $g$, we have

$$ (S32) \qquad \hat{z} = f(\mu_g, \varepsilon_k) = \mu_g^T u_k. $$

In a given environment $\varepsilon_k$, the unit vector $u_k$ is equal to 1 at the $k$th index and 0 elsewhere. The reaction norm gradient is equal to this unit vector, i.e. $\psi_{\varepsilon_k} = u_k$. In the first environment, for example, we have $\psi_{\varepsilon_1} = u_1 = (1, 0, \dots)^T$. As mentioned in Appendix A, the character-state approach is linear in its parameters. We can thus compute the genotypic/breeding values in a given environment $\varepsilon_k$ as

$$ (S33) \qquad G_z = A_z = \hat{z} - E_{g|\varepsilon}(\hat{z}) = \mu_g^T u_k - \mu^T u_k = \mu_{g,k} - \mu_j, $$

where $\mu_{g,k}$ and $\mu_j$ are the $k$th values of the vectors $\mu_g$ and $\mu$. The total (additive) genetic variance is the variance of the breeding values across environments:

$$ (S34) \qquad V_{\text{Gen}} = V_{\text{Add}} = V(A_z) = V(\mu_{g,k}). $$

Since the variance-covariance matrix of $\mu_g$ is the $G_z$ matrix, the variance of all elements $\mu_{g,k}$ taken together is the average of the diagonal elements of $G_z$, which we will note $V_k$. Assuming that all environments are equiprobable for the sake of simplicity (releasing this assumption merely requires to use weighted average), we have

$$ (S35) \qquad V_{\text{Add}} = \frac{1}{K} \sum_{k=1}^{K} V_k. $$

In other words, $V_{\text{Add}}$ is the average of the diagonal elements of the $G_z$ matrix.
The environment-blind (additive) genetic variance of the trait depends on the average of the breeding values across environment for a given genotype:

$$ (S36) \qquad \bar{A} = \frac{1}{K} \sum_k A_{z,k}, $$

where $A_{z,k}$ is the breeding value evaluated at the $k$th environment for a given genotype, still assuming equiprobable environments. It results that the environment-blind (additive) genetic variance of the trait is

$$ (S37) \qquad V_G = V_A = \frac{1}{K^2} \left( \sum_k V_k + 2 \sum_{k<l} C_{kl} \right), $$

where $C_{kl}$ is the genetic covariance between the environment $k$ and $l$. In other words, $V_A$ is the average of all the elements of the $G_z$ matrix.

Finally, the (additive) genetic variance arising from plasticity can be computed as the difference between $V_{Add}$ and $V_A$:

$$ (S38) \qquad V_{G \times E} = V_{A \times E} = V_{Add} - V_A = \frac{1}{K^2} \left( (K-1) \sum_k V_k - 2 \sum_{k<l} C_{kl} \right) $$

A few particular cases are important to note here. The first case is when all environments harbour the same additive genetic variance, say $V$, and are all perfectly correlated with one another. This is a situation generally described as a total absence of genetic variation in plasticity. In our framework, this situation would indeed result in $V_{Add} = V_A = V$ and, indeed, no genetic variation arising from plasticity with $V_{A \times E} = 0$. Note that uneven additive genetic variances across environments, even if genetic correlation are kept perfect across environments, would result in slightly positive genetic variance arising from plasticity with $V_{A \times E} > 0$. This is because, in such context, the trait can still evolve faster in some environments compared to other, hence plasticity can evolve. The second extreme case, is when the environment-blind additive genetic variance of the trait is null, i.e. $V_A = 0$, while all the additive genetic variance in reaction norm is composed of the additive genetic variance arising from plasticity, i.e. $V_{Add} = V_{A \times E}$. This happens when the sum of covariances (the total of which must be negative) exactly compensates the sum of diagonal variances in the $G_z$, meaning that negative genetic correlation between environments are maximised. In this case, it is impossible for directional selection to act on average value of the trait across all environments, but the evolvability of plasticity is maximal. A third, interesting case is when there is absolutely no genetic correlation between environments, i.e. the off-diagonal elements of $G_z$ are all equal to 0. In such case, it is important to note that, because evolution can freely operate across environments, then both $V_A = \frac{1}{K^2} \sum_k V_k$ and $V_{A \times E} = \frac{K-1}{K^2} \sum_k V_k$ are non-zero.

## C5. Decomposition of variance for individual-based reaction norms

In Equation 4, we assumed that the only source of variation in $\theta$ is of genetic origin. This is a classical assumption both in the empirical and theoretical literature (de Jong 1990; Gavrilets and Scheiner 1993a; Via and Lande 1985), but in many cases, it can be useful or needed to include further sources of variation in $\theta$. This is for example the case when studying reaction norms using repeated measurements of the same individual in different environments. In particular, this may require including a further “permanent environment” effect to account for multiple repeats (Wilson et al. 2010) on the same individual, and also allows for the modelling of the reaction norm at the individual level (individual plasticity, Nussey et al. 2007). When other random effects are assumed in the model, we can write the full variation of $\theta$ as:

$$ (S39) \qquad \theta \sim \mathcal{N}(\bar{\theta}, \mathbf{V}_{\theta}), $$
where $V_θ$ is the total variance-covariance matrix of $θ$. Note that Equation 4 is still valid to model the genetic component of $θ$ which we named $θ_g$. In such case, the heritability of the kth component of $θ$ can be computed as the ratio of the kth diagonal element of $G_θ$ to the kth element of $V_θ$, i.e. $h_{θ,k}^2 = \frac{G_{θ,k,k}}{V_{θ,k,k}}$. Because the modelling of $θ_g$ remains unchanged, all our computations of (additive) genetic variances and their decomposition remains completely identical. However, there are two important changes. The first change is that the definition of $V_{Plas}$ does not only depend on averaging over $g$ any more, but on other sources of variations in $θ$ as well, i.e. $V_{Plas} = V(E_{θ|ε}(\hat{z}))$. This means that the marginalisation step conditional to the environment now implies the full $V_θ$ rather only its subcomponent $G_θ$. The second change is that it is not possible to write the total variance of the reaction norm as the sum of $V_{Plas}$ and $V_{Gen}$ anymore, because the latter is only a partial reflection of the full variation in $θ$. Instead, we need to introduce the phenotypic variation in the trait arising from the full sources of variation in $θ$, which we denote here $V_{Param}$:

$$ (S40) \qquad V_{Param} = V(\hat{z} - E_{θ|ε}(\hat{z})) = E(V_{θ|ε}(\hat{z})). $$

Then, we can write the correct formulae for $V_P$ and $T_{RN}^2$:

$$ (S41) \qquad V_P = V_{Plas} + V_{Param} + V_{Res}, \qquad T_{RN}^2 = \frac{V_{Plas} + V_{Param}}{V_P}. $$

The Reacnorm package was designed to be able to input $V_θ$ to compute those quantities if needed.

# Appendix D. Derivation of π- and ϕ-partition of $V_{Plas}$

## D1. The π-decomposition

We have seen in Appendix C how to compute the variance arising from the average shape of reaction norm $V_{Plas}$. In order to go further, we now separate this into a component linked to the average slope of the reaction norm and another linked to the average curvature. For this, we need one or two of the following assumptions to hold true: (i) the environment $ε$ follows a normal distribution; or (ii) the function $f$ is quadratic. In such context, we can isolate the contribution of the slope, $V_{Sl}$, from the contribution of the curvature, $V_{Cv}$ to $V_{Plas}$, based on the best quadratic approximation of $E_{g|ε}(\hat{z})$ (akin to the reasoning in Lande and Arnold 1983, for estimates of selection gradients), as:

$$ (S42) \qquad V_{Sl} = E\left(\frac{dE_{g|ε}}{dε}(\hat{z})\right)^2 V(ε), \qquad V_{Cv} = \frac{1}{4}E\left(\frac{d^2E_{g|ε}}{dε^2}(\hat{z})\right)^2 V(ε^2). $$

As an illustration of why the assumptions above are needed, if $ε$ follows a uniform distribution between -2 and 2; and the average shape of plasticity is the following cubic function, $f(ε) = 2ε - 0.5ε^2 - ε^3$, then the average slope is -2, while the slope from the best quadratic approximation of $E_{g|ε}(\hat{z})$ is -0.4. In such cases, the decomposition in Equation S42 is not valid anymore, due to (i) the impossibility to apply Stein's lemma to a non-normal distribution and (ii) strong covariation between the slope and curvature. This means that whenever the environment is non-normal and the reaction norm is non-quadratic, the π-decomposition can bear little meaning (in the cubic example above, $V_{Sl}$ would be 5.4, while $V_{Plas} = 2.0$, so that $π_{Sl}$ would be largely above 1). A truly quadratic reaction norm is the only case where $π_{Sl} + π_{Cv} = 1$.

## D2. The ϕ-decomposition

In such cases where the environment is non-normal and the reaction norm is non-quadratic, it is always possible to approximate the true shape of the reaction norm using a polynomial function:

$$ (S43) \qquad \hat{z} = \sum_{n=0}^{N} (\bar{\theta}_n + \theta_{n,g})ε^n $$
In the context of decomposing $V_{Plas}$, such polynomial approximation provides a possibility to isolate the (co-)contribution of the (pairs of) coefficients in $E_{g|\epsilon}(\hat{z}) = \sum_{n=0}^{N} \hat{\theta}_n \epsilon^n$:

$$ (S44) \qquad V_{Plas} = V(E_{g|\epsilon}(\hat{z})) = \sum_{n} \bar{\theta}_n^2 V(\epsilon^n) + 2 \sum_{n<m} \bar{\theta}_n \bar{\theta}_m \text{cov}(\epsilon^n, \epsilon^m) $$

From this, we suggest the alternative $\varphi$-decomposition of $V_{Plas}$, with $\varphi_n = \frac{\bar{\theta}_n^2 V(\epsilon^n)}{V_{Plas}}$ and $\varphi_{nm} = \frac{2\bar{\theta}_n \bar{\theta}_m \text{cov}(\epsilon^n, \epsilon^m)}{V_{Plas}}$. It is important to note that this decomposition is based on the coefficients of the polynomial function and, thus, it is unfortunately impossible to simply interpret the $\varphi_n$ in terms of slope (for $\varphi_1$), curvature (for $\varphi_2$), and so on. The only exception is when the reaction norm shape is quadratic, in which case $\pi_{S1} = \varphi_1$ and $\pi_{Cv} = \varphi_2$.

## Appendix E. Correcting for uncertainty in the estimation of fixed effects

Character-state approach. It is easier to start with the character-state approach based on the ANOVA model. We want to compute $V_{Plas}$ as the variance of the group-level effects $\mu$:

$$ (S45) \qquad V_{Plas} = V(\mu) $$

However, we do not have access to the real-world values for $\mu$, but only to the estimated $\hat{\mu}$ from the model. Such estimates, if unbiased, have an expected value of $\mu_k$ in environment $k$ and a standard-error (i.e. the estimation of the sampling standard deviation) $s_k$. In other words, we can state that $\hat{\mu}_k$ is equal to $\mu_k$ up to an additive error:

$$ (S46) \qquad \hat{\mu}_k = \mu_k + \tilde{\mu}_k $$

where $\tilde{\mu}$ is of mean 0 and variance $s_k^2$. Considering each virtual repeat $r$ of the experiment, we can apply the law of total variance:

$$ (S47) \qquad V(\hat{\mu}) = V_{\epsilon}(E_{r|\epsilon}(\hat{\mu})) + E_{\epsilon}(V_{r|\epsilon}(\hat{\mu})) = V_{\epsilon}(\mu) + E_{\epsilon}(s^2). $$

We thus have:

$$ (S48) \qquad V_{Plas} = V_{\epsilon}(\mu) = V_{\epsilon}(\hat{\mu}) - E_{\epsilon}(s^2) $$

This result is equivalent to e.g. the classical computation of the “sire variance” in sire models in quantitative genetics (Lynch and Walsh 1998), although the latter is generally expressed using sums-of-squares.

Curve-parameter approach. There is unfortunately no simple solution to the problem of accounting for the uncertainty of fixed effects in the general context of non-linear modelling. However, for the particular case where the model can be framed as a linear model, as is the case for the polynomial function, then $\hat{z} = X\theta$, where $X$ is the design matrix containing the values for the environment. Noting $\Sigma_X$ the variance-covariance matrix of $X$, we can define $V_{Plas}$ as:

$$ (S49) \qquad V_{Plas} = \theta^T \Sigma_X \theta. $$

Again, the problem is that $\theta$ is unknown, we only have access to the estimated values of the parameters, $\hat{\theta}$, that are inferred with an error provided by the variance-covariance matrix of standard errors, $S_{\theta}$. We can write again:

$$ (S50) \qquad \hat{\theta} = \bar{\theta} + \tilde{\theta}, $$

Noting that the error is independent from the true value, we have:

$$ (S51) \qquad \hat{\theta}^T \Sigma_X \hat{\theta} = \theta^T \Sigma_X \theta + \tilde{\theta}^T \Sigma_X \tilde{\theta} $$

To express $\tilde{\theta}^T \Sigma_X \hat{\theta}$, it is important to note that $S_{\theta,ij} = E(\tilde{\theta}_i \tilde{\theta}_j)$, since $E(\tilde{\theta}) = 0$. Then, we can note that, the error being unknown, we actually want to compute $E_r(\tilde{\theta}^T \Sigma_X \hat{\theta})$ taken across virtual repeats $r$ of the experiment:

$$ (S52) \qquad E_r(\tilde{\theta}^T \Sigma_X \hat{\theta}) = E_r(\sum_{ij} \tilde{\theta}_i \tilde{\theta}_j \Sigma_{X,i,j}) = \sum_{ij} E_r(\tilde{\theta}_i \tilde{\theta}_j) \Sigma_{X,i,j} = \sum_{ij} S_{\theta,ij} \Sigma_{X,i,j} = \text{Tr}(S_{\theta} \Sigma_X) $$

This is similar to the result of Brown and Rutemiller (1977). Finally, we have:

$$ (S53) \qquad V_{Plas} = \hat{\theta}^T \Sigma_X \hat{\theta} - \text{Tr}(S_{\theta} \Sigma_X). $$
Appendix F. Full results for the section “Perfect modelling of quadratic curves”

This section provides the full results corresponding to the section “Perfect modelling of qua-
dratic curves” in the main text. The results of all investigated values for the number of envi-
ronments (10 or 4) and number of genotypes (20 or 5 for the discrete case, 200 or 50 for the
continuous case) are provided for the discrete and continuous cases.
<figure><img src="image_12.png" /><figcaption>**Figure S1** – Distribution of the error (difference between the inferred and true value) for each of the inferred variance components for three discrete scenarios: $N_{Env}$: number of environments, $N_{Gen}$: number of different genotypes, $N_{Rep}$: number of replicates per genotype. Estimates are for $\hat{P}_{RN}^2$ (proportion of variance generated by plasticity after averaging across genotypes), $\hat{h}_{RN}^2$ (total heritability of the reaction norm), $\hat{h}^2$ (environment-blind heritability) and $\hat{h}_1^2$ (heritability from plasticity) for both the curve-parameter and character-state approaches. For the curve-parameter, the $\pi$-decomposition of $\hat{P}_{RN}^2$ into $\pi_{SI}$ (contribution of the slope) and $\pi_{CV}$ (contribution of the curvature); the $\gamma$-decomposition of $\hat{h}_{RN}^2$ into $\gamma_a$ (genetic contribution of the intercept), $\gamma_b$ (genetic contribution of the slope), $\gamma_c$ (genetic contribution of the curvature) and $\gamma_{ac}$ (genetic contribution of the covariance between the intercept and the curvature) and the $\ell$-decomposition of $\hat{h}_1^2$ into $\ell_b$ (slope) and $\ell_c$ (curvature) are also shown. The grey dots correspond to the average over the 1000 simulations.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

**Figure S2** – Distribution of the error (difference between the inferred and true value) for each of the inferred variance components for four continuous scenarios: $N_{\text{env}}$: number of environment tested per genotype, $N_{\text{Gen}}$: number of different genotypes. The character-state approach was impossible for the continuous environment scenario. Estimates are for $\hat{P}_{\text{RN}}^2$ (proportion of variance generated by plasticity after averaging across genotypes), $\hat{h}_{\text{RN}}^2$ (total heritability of the reaction norm), $\hat{h}^2$ (environment-blind heritability) and $\hat{h}_1^2$ (heritability from plasticity) for both the curve-parameter and character-state approaches. For the curve-parameter, the $\pi$-decomposition of $\hat{P}_{\text{RN}}^2$ into $\pi_{\text{SI}}$ (contribution of the slope) and $\pi_{\text{CV}}$ (contribution of the curvature); the $\gamma$-decomposition of $\hat{h}_{\text{RN}}^2$ into $\gamma_a$ (genetic contribution of the intercept), $\gamma_b$ (genetic contribution of the curvature) and the $\ell$-decomposition of $\hat{h}_1^2$ into $\ell_b$ (slope) and $\ell_c$ (curvature) are also shown. The grey dots correspond to the average over the 1000 simulations.

<figure><img src="image_14.png" /></figure>
## Appendix G. Comparison with the approach from Murren et al. (2014)

Murren et al. (2014) studied variation of the reaction norm shapes across different datasets, using their own metrics. We argue in the main text that our variance decomposition is more appropriate than the ones suggested by Murren et al. (2014), and we develop here why.

The first step in the approach of Murren et al. (2014) is to choose a reference reaction norm in each of the studies and compute contrasts (i.e. difference with) to that particular reaction norm. The contrasts are then analysed, rather than the reaction norms themselves. For the sake of simplicity, and because this does not (or marginally) impact our comments on this approach, we will overlook that step and consider reaction norms directly.

For each genotype k and from its given reaction norm (or contrast) $z_k = \{z_{k,1}, \dots, z_{k,n}\}$, Murren et al. (2014) compute four statistics (we removed the absolute values for the sake of simplicity here):

(1) The offset, $O_M$, measures the “location” of the reaction norm, i.e. its mean. Comparison of the offsets allows detecting whether reaction norms are “shifted” toward higher or lower values. It is computed, for each genotype k, as the absolute value of the average of the norm across environments:

$$ (S54) \qquad O_{M,k} = \frac{\sum_{i}^{n} |z_{k,i}|}{n}. $$

(2) The slope, $S_M$, measures the linear trend of the reaction norms. Formally, it is the absolute sum of the differences between two consecutive environments, divided by the number of intervals ($n-1$):

$$ (S55) \qquad S_{M,k} = \frac{\sum_{i}^{n-1} |z_{k,i+1} - z_{k,i}|}{n-1}. $$

(3) The curvature, $C_M$, is computed as the absolute value of the average change in phenotype between two consecutive pairs of environments:

$$ (S56) \qquad C_{M,k} = \frac{\sum_{i}^{n-2} |(z_{k,i+2} - z_{k,i+1}) - (z_{k,i+1} - z_{k,i})|}{n-2}. $$

(4) The wiggle, $W_M$, is, according to the authors the “the variability in shape not described by any of the previous three measures”:

$$ (S57) \qquad W_{M,k} = \frac{\sum_{i}^{n-2} |(z_{k,i+2} - z_{k,i+1}) - (z_{k,i+1} - z_{k,i})|}{n-2} - C_{M,k}. $$

Given the lower interest in this latter statistics, we will not comment on it any further. Most of the comments on the other statistics also apply to this one.

One strong assumption underlying the calculations above is that environmental values $\epsilon = \{\epsilon_1, \dots, \epsilon_n\}$ on which the reaction norms were evaluated are evenly spaced, e.g. that the differences $\epsilon_{i+1} - \epsilon_i$ are equal for all possible values of i. The assumption is actually that the space between two measures is equal to 1 (which, admittedly, is only a matter of rescaling when evenly-spaced values are already assumed). If this is the case, then there is indeed no loss in generality in using the number of components ($n, n-1$ and $n-2$) rather than actual values of $x$ in the denominator. Although it is common for studies on reaction norms to use evenly-spaced environmental values, it is an unnecessary assumption that shall not be satisfied by all studies.

Second, developing the sums in $S_M$ and $C_M$ above show that the intermediate values cancel each other out, leaving only the values at each extreme of the environmental range in the estimate:

$$ (S58) \qquad \begin{aligned} S_{M,k} &= \frac{z_{k,n} - z_{k,1}}{n-1}, \\ C_{M,k} &= \frac{(z_{k,n} - z_{k,n-1}) - (z_{k,2} - z_{k,1})}{n-2}. \end{aligned} $$

The issue here is double: (i) the estimation is highly sensitive to the random noise coming from a small number of values (two or three/four); and (ii) the intermediate values in the reaction norm are simply thrown out and not used for a more robust estimation. In other words, it would have
been exactly the same to not measure the reaction norm at these intermediate values, since they
are not accounted for in the calculation.

A final issue is that the approach uses the measured values of the reaction norms without ac-
counting for the uncertainty in their estimation (i.e. standard-deviation and sample size for each
genotype and environmental value) which poses the well-known issue of non-propagation of
the error when doing "statistics on statistics".

Although we also provide estimators of the impact of several aspects of reaction norms on the phenotypic variation, our approach differs from the one from Murren et al. (2014) by many aspects. First, our variance decomposition makes the explicit distinction between the average shape of the reaction norm and the genetic variance surrounding it. As such, to $O_M, S_M$ and $C_M$ corresponds not only the $\pi$-, but also the $\gamma$- and $\iota$-decomposition. We clearly delimit the domain of validity of each of these decomposition. We also account for possible correlation between those components. Second, we use the whole of the statistical inference to define our variance decomposition estimates. Third, we explicitly account for the uncertain estimation of reaction norms.

## Acknowledgements

We are grateful to Jarrod Hadfield, Thibaut Morel-Journel, Stéphane Robin and John Stinchcombe for useful discussions and/or comments that much improved the quality of the paper. Preprint version 5 of this article has been peer-reviewed and recommended by Peer Community In Evolutionary Biology (https://doi.org/10.24072/pci.evolbiol.100708 Jacob 2025).

## Fundings

We thank the European Research Council (Grant STG-678140-FluctEvol to LMC).

## Conflict of interest disclosure

The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article. The authors declare the following non-financial conflict of interest: both PdV and LMC are recommenders for PCI Evolutionary Biology.

## Data, script, code, and supplementary information availability

The code for the data simulation and analyses performed in this article is available at the following repository: github.com/devillemereuil/CodePartReacnorm (de Villemereuil 2025a). The R package Reacnorm can be found at the following repository: github.com/devillemereuil/Reacnorm (de Villemereuil 2025b).

## References

Albecker MA, Trussell GC, Lotterhos KE (2022). A Novel Analytical Framework to Quantify Co-Gradient and Countergradient Variation. Ecology Letters **25**, 1521-1533. https://doi.org/10.1111/ele.14020.

Angilletta MJ (2009). Thermal Adaptation: A Theoretical and Empirical Synthesis. OUP Oxford. 304 pp. https://doi.org/10.1093/acprof:oso/9780198570875.001.1.

Bates D, Mächler M, Bolker B, Walker S (2015). Fitting Linear Mixed-Effects Models Using Lme4. Journal of Statistical Software **67**, 48. https://doi.org/10.18637/jss.v067.i01.

Bonamour S, Chevin LM, Charmantier A, Teplitsky C (2019). Phenotypic Plasticity in Response to Climate Change: The Importance of Cue Variation. Philosophical Transactions of the Royal Society B: Biological Sciences **374**, 20180178. https://doi.org/10.1098/rstb.2018.0178.

Bradshaw AD (1965). Evolutionary Significance of Phenotypic Plasticity in Plants. In: Advances in Genetics. Ed. by EW Caspari and JM Thoday. Vol. 13. Cambridge (MA, USA): Academic Press, pp. 115-155. https://doi.org/10.1016/S0065-2660(08)60048-6.
Brown GG, Rutemiller HC (1977). Means and Variances of Stochastic Vector Products with Applications to Random Linear Models. Management Science **24**, 210-216. https://doi.org/10.1287/mnsc.24.2.210.

Bürkner PC (2018). Advanced Bayesian Multilevel Modeling with the R Package Brms. The R Journal **10**, 395-411. https://doi.org/10.32614/RJ-2018-017.

Charmantier A, McCleery RH, Cole LR, Perrins C, Kruuk LEB, Sheldon BC (2008). Adaptive Phenotypic Plasticity in Response to Climate Change in a Wild Bird Population. Science **320**, 800-803.
https://doi.org/10.1126/science.1157174.

Chevin LM, Collins S, Lefèvre F (2013). Phenotypic Plasticity and Evolutionary Demographic Responses to Climate Change: Taking Theory out to the Field. Functional Ecology **27**, 967-979.
https://doi.org/10.1111/j.1365-2435.2012.02043.x.

Chevin LM, Lande R, Mace GM (2010). Adaptation, Plasticity, and Extinction in a Changing Environment: Towards a Predictive Theory. PLOS Biology **8**, e1000357. https://doi.org/10.1371/journal.pbio.1000357.

de Jong G (1990). Quantitative Genetics of Reaction Norms. Journal of evolutionary biology **3**, 447-468. https://doi.org/10.1046/j.1420-9101.1990.3050447.x.

de Jong G (1995). Phenotypic Plasticity as a Product of Selection in a Variable Environment. The American Naturalist **145**, 493-512. https://doi.org/10.1086/285752.

de Jong G (1999). Unpredictable Selection in a Structured Population Leads to Local Genetic Differentiation in Evolved Reaction Norms. Journal of Evolutionary Biology **12**, 839-851. https://doi.org/10.1046/j.1420-9101.1999.00118.x.

Des Marais DL, Hernandez KM, Juenger TE (2013). Genotype-by-Environment Interaction and Plasticity: Exploring Genomic Responses of Plants to the Abiotic Environment. Annual Review of Ecology, Evolution, and Systematics **44**, 5-29. https://doi.org/10.1146/annurev-ecolsys-110512-135806.

Deutsch CA, Tewksbury JJ, Huey RB, Sheldon KS, Ghalambor CK, Haak DC, Martin PR (2008). Impacts of Climate Warming on Terrestrial Ectotherms across Latitude. Proceedings of the National Academy of Sciences **105**, 6668-6672. https://doi.org/10.1073/pnas.0709472105.

de Villemereuil P (2018). Quantitative Genetic Methods Depending on the Nature of the Phenotypic Trait. Annals of the New York Academy of Sciences. The Year in Evolutionary Biology **1422**, 29-47. https://doi.org/10.1111/nyas.13571.

de Villemereuil P (2025a). *Devillemereuil/CodePartReacNorm: Accepted Version*. Version accepted. Zenodo. https://doi.org/10.5281/zenodo.14674496.

de Villemereuil P (2025b). *Devillemereuil/Reacnorm: Initial Release on CRAN*. Version v0.2. Zenodo. https://doi.org/10.5281/zenodo.14704700.

de Villemereuil P, Charmantier A, Arlt D, Bize P, Brekke P, Brouwer L, Cockburn A, Côté SD, Dobson FS, Evans SR, Festa-Bianchet M, Gamelon M, Hamel S, Hegelbach J, Jerstad K, Kempenaers B, Kruuk LEB, Kumpula J, Kvalnes T, McAdam AG, et al. (2020). Fluctuating Optimum and Temporally Variable Selection on Breeding Date in Birds and Mammals. Proceedings of the National Academy of Sciences **117**, 31969-31978. https://doi.org/10.1073/pnas.2009003117.

de Villemereuil P, Morrissey MB, Nakagawa S, Schielzeth H (2018). Fixed-Effect Variance and the Estimation of Repeatabilities and Heritabilities: Issues and Solutions. Journal of Evolutionary Biology **31**, 621-632. https://doi.org/10.1111/jeb.13232.

de Villemereuil P, Schielzeth H, Nakagawa S, Morrissey MB (2016). General Methods for Evolutionary Quantitative Genetic Inference from Generalised Mixed Models. Genetics **204**, 1281-1294.
https://doi.org/10.1534/genetics.115.186536.

Falconer DS (1952). The Problem of Environment and Selection. The American Naturalist **86**, 293-298. https://doi.org/10.1086/281736.

Falconer DS, Mackay TF (1996). Introduction to Quantitative Genetics. 4th ed. Harlow, Essex (UK): Benjamin Cummings.

Gavrilets S, Scheiner SM (1993a). The Genetics of Phenotypic Plasticity. V. Evolution of Reaction Norm Shape. Journal of Evolutionary Biology **6**, 31-48. https://doi.org/10.1046/j.1420-9101.1993.6010031.x.
Gavrilets S, Scheiner SM (1993b). *The Genetics of Phenotypic Plasticity*. VI. Theoretical Predictions for Directional Selection. *Journal of Evolutionary Biology* **6**, 49–68. https://doi.org/10.1046/j.1420-9101.1993.6010049.x.

Gienapp P, Teplitsky C, Alho JS, Mills JA, Merilä J (2008). *Climate Change and Evolution: Disentangling Environmental and Genetic Responses*. *Molecular Ecology* **17**, 167–178. https://doi.org/10.1111/j.1365-294X.2007.03413.x.

Gomulkiewicz R, Kirkpatrick M (1992). *Quantitative Genetics and the Evolution of Reaction Norms*. *Evolution* **46**, 390–411. https://doi.org/10.1111/j.1558-5646.1992.tb02047.x.

Hammill E, Rogers A, Beckerman AP (2008). *Costs, Benefits and the Evolution of Inducible Defences: A Case Study with Daphnia Pulex*. *Journal of Evolutionary Biology* **21**, 705–715. https://doi.org/10.1111/j.1420-9101.2008.01520.x.

Jacob S (2025). *Phenotypic Plasticity across Traits, Organisms and Environments: Decomposing the Components of Phenotypic Variance*. *Peer Community in Evolutionary Biology* **1**, 100708. https://doi.org/10.24072/pci.evolbiol.100708.

Johnson JB, Omland KS (2004). *Model Selection in Ecology and Evolution*. *Trends in Ecology & Evolution* **19**, 101–108. https://doi.org/10.1016/j.tree.2003.10.013.

Johnson PC (2014). *Extension of Nakagawa & Schielzeth’s R2GLMM to Random Slopes Models*. *Methods in Ecology and Evolution* **5**, 944–946. https://doi.org/10.1111/2041-210X.12225.

King JG, Hadfield JD (2019). *The Evolution of Phenotypic Plasticity When Environments Fluctuate in Time and Space*. *Evolution Letters* **3**, 15–27. https://doi.org/10.1002/evl13.100.

Kirkpatrick M (2009). *Patterns of Quantitative Genetic Variation in Multiple Dimensions*. *Genetica* **136**, 271–284. https://doi.org/10.1007/s10709-008-9302-6.

Kirkpatrick M, Heckman N (1989). *A Quantitative Genetic Model for Growth, Shape, Reaction Norms, and Other Infinite-Dimensional Characters*. *Journal of Mathematical Biology* **27**, 429–450. https://doi.org/10.1007/BF00290638.

Lande R (1979). *Quantitative Genetic Analysis of Multivariate Evolution, Applied to Brain:Body Size Allometry*. *Evolution* **33**, 402–416. https://doi.org/10.2307/2407630.

Lande R (2009). *Adaptation to an Extraordinary Environment by Evolution of Phenotypic Plasticity and Genetic Assimilation*. *Journal of Evolutionary Biology* **22**, 1435–1446. https://doi.org/10.1111/j.1420-9101.2009.01754.x.

Lande R, Arnold SJ (1983). *The Measurement of Selection on Correlated Characters*. *Evolution* **37**, 1210–1226. https://doi.org/10.2307/2408842.

Landsman Z, Nešlehová J (2008). *Stein’s Lemma for Elliptical Random Vectors*. *Journal of Multivariate Analysis* **99**, 912–927. https://doi.org/10.1016/j.jmva.2007.05.006.

Landsman Z, Vanduffel S, Yao J (2013). *A Note on Stein’s Lemma for Multivariate Elliptical Distributions*. *Journal of Statistical Planning and Inference* **143**, 2016–2022. https://doi.org/10.1016/j.jspi.2013.06.003.

Lawley DN (1956). *Tests of Significance for the Latent Roots of Covariance and Correlation Matrices*. *Biometrika* **43**, 128–136. https://doi.org/10.2307/2333586.

Lynch M, Walsh B (1998). *Genetics and Analysis of Quantitative Traits*. *Sunderland, Massachussets (US): Sinauer Associates*.

Lynch M, Gabriel W (1987). *Environmental Tolerance*. *The American Naturalist* **129**, 283–303. https://doi.org/10.1086/284635.

Merilä J, Hendry AP (2014). *Climate Change, Adaptation, and Phenotypic Plasticity: The Problem and the Evidence*. *Evolutionary Applications* **7**, 1–14. https://doi.org/10.1111/eva.12137.

Mitchell DJ, Houslay TM (2021). *Context-Dependent Trait Covariances: How Plasticity Shapes Behavioral Syndromes*. *Behavioral Ecology* **32**, 25–29. https://doi.org/10.1093/beheco/araa115.

Moczek, Emlen (1999). *Proximate Determination of Male Horn Dimorphism in the Beetle Onthophagus Taurus (Coleoptera: Scarabaeidae)*. *Journal of Evolutionary Biology* **12**, 27–37. https://doi.org/10.1046/j.1420-9101.1999.00004.x.

Morrissey MB (2015). *Evolutionary Quantitative Genetics of Nonlinear Developmental Systems*. *Evolution* **69**, 2050–2066. https://doi.org/10.1111/evo.12728.
Morrissey MB, Liefting M (2016). *Variation in Reaction Norms: Statistical Considerations and Biological Interpretation*. Evolution **70**, 1944-1959. https://doi.org/10.1111/evo.13003.

Murren CJ, Maclean HJ, Diamond SE, Steiner UK, Heskel MA, Handelsman CA, Ghalambor CK, Auld JR, Callahan HS, Pfennig DW (2014). *Evolutionary Change in Continuous Reaction Norms*. The American Naturalist **183**, 453-467. https://doi.org/10.1086/675302.

Nakagawa S, Schielzeth H (2013). *A General and Simple Method for Obtaining R2 from Generalized Linear Mixed-Effects Models*. Methods in Ecology and Evolution **4**, 133-142. https://doi.org/10.1111/j.2041-210x.2012.00261.x.

Narasimhan B, Johnson SG, Hahn T, Bouvier A, Kiêu K (2023). *Cubature: Adaptive Multivariate Integration over Hypercubes*. https://doi.org/10.32614/cran/package.cubature.

Nussey DH, Wilson AJ, Brommer JE (2007). *The Evolutionary Ecology of Individual Phenotypic Plasticity in Wild Populations*. Journal of Evolutionary Biology **20**, 831-844. https://doi.org/10.1111/j.1420-9101.2007.01300.x.

Nussey DH, Postma E, Gienapp P, Visser ME (2005). *Selection on Heritable Phenotypic Plasticity in a Wild Bird Population*. Science **310**, 304-306. https://doi.org/10.1126/science.1117004.

Pélabon C, Hilde CH, Einum S, Gamelon M (2020). *On the Use of the Coefficient of Variation to Quantify and Compare Trait Variation*. Evolution Letters. https://doi.org/10.1002/evl3.171.

Pinheiro J, Bates D, DebRoy S, Sarkar D, the R Core team (2009). *Nlme: Linear and Nonlinear Mixed Effects Models*. https://doi.org/10.32614/CRAN/package.nlme.

Reed TE, Waples RS, Schindler DE, Hard JJ, Kinnison MT (2010). *Phenotypic Plasticity and Population Viability: The Importance of Environmental Predictability*. Proceedings of the Royal Society B: Biological Sciences **277**, 3391-3400. https://doi.org/10.1098/rspb.2010.0771.

Rice SH (2004). *Evolutionary Theory: Mathematical and Conceptual Foundations*. 1st edition. Sunderland (Mass.): Sinauer Associates is an imprint of Oxford University Press.

Robertson A (1966). *A Mathematical Model of the Culling Process in Dairy Cattle*. Animal Science **8**, 95-108. https://doi.org/10.1017/S0003356100037752.

Rovelli G, Ceccobelli S, Perini F, Demir E, Mastrangelo S, Conte G, Abeni F, Marletta D, Ciampolini R, Cassandro M, Bernabucci U, Lasagna E (2020). *The Genetics of Phenotypic Plasticity in Live-stock in the Era of Climate Change: A Review*. Italian Journal of Animal Science **19**, 997-1014. https://doi.org/10.1080/1828051X.2020.1809540.

Schaum CE, Collins S (2014). *Plasticity Predicts Evolution in a Marine Alga*. Proceedings of the Royal Society B: Biological Sciences **281**, 20141486. https://doi.org/10.1098/rspb.2014.1486.

Scheiner SM (1993a). *Genetics and Evolution of Phenotypic Plasticity*. Annual Review of Ecology and Systematics **24**, 35-68. https://doi.org/10.1146/annurev.es.24.110193.000343.

Scheiner SM (1993b). *Plasticity as a Selectable Trait: Reply to Via*. The American Naturalist **142**, 371-373. https://doi.org/10.1086/285544.

Scheiner SM, Lyman RF (1989). *The Genetics of Phenotypic Plasticity I. Heritability*. Journal of Evolutionary Biology **2**, 95-107. https://doi.org/10.1046/j.1420-9101.1989.2020095.x.

Schlichting CD, Pigliucci M (1998). *Phenotypic Evolution: A Reaction Norm Perspective*. 1st edition. Sunderland, Mass: Sinauer Associates Inc.

Stinchcombe JR, Function-valued Traits Working Group, Kirkpatrick M (2012). *Genetics and Evolution of Function-Valued Traits: Understanding Environmentally Responsive Phenotypes*. Trends in Ecology & Evolution **27**, 637-647. https://doi.org/10.1016/j.tree.2012.07.002.

Stirling G, Roff DA (2000). *Behaviour Plasticity without Learning: Phenotypic and Genetic Variation of Naïve Daphnia in an Ecological Trade-Off*. Animal Behaviour **59**, 929-941. https://doi.org/10.1006/anbe.1999.1386.

Suzuki Y, Nijhout HF (2006). *Evolution of a Polyphenism by Genetic Accommodation*. Science **311**, 650-652. https://doi.org/10.1126/science.1118888.

Teplitsky C, Mills JA, Alho JS, Yarrall JW, Merilä J (2008). *Bergmann's Rule and Climate Change Revisited: Disentangling Environmental and Genetic Responses in a Wild Bird Population*. Proceedings of the National Academy of Sciences **105**, 13492-13496. https://doi.org/10.1073/pnas.0800999105.
Tredennick AT, Hooker G, Ellner SP, Adler PB (2021). *A Practical Guide to Selecting Models for Exploration, Inference, and Prediction in Ecology*. Ecology **102**, e03336. https://doi.org/10.1002/ecy.3336

Tufto J (2000). *The Evolution of Plasticity and Nonplastic Spatial and Temporal Adaptations in the Presence of Imperfect Environmental Cues*. The American Naturalist **156**, 121–130. https://doi.org/10.1086/303381.

Tufto J (2015). *Genetic Evolution, Plasticity, and Bet-Hedging as Adaptive Responses to Temporally Autocorrelated Fluctuating Selection: A Quantitative Genetic Model*. Evolution **69**, 2034–2049. https://doi.org/10.1111/evo.12716.

Vedder O, Bouwhuis S, Sheldon BC (2013). *Quantitative Assessment of the Importance of Phenotypic Plasticity in Adaptation to Climate Change in Wild Bird Populations*. PLOS Biology **11**, e1001605. https://doi.org/10.1371/journal.pbio.1001605.

Via S, Lande R (1985). *Genotype-Environment Interaction and the Evolution of Phenotypic Plasticity*. Evolution **39**, 505–522. https://doi.org/10.1111/j.1558-5646.1985.tb00391.x.

Wilson AJ, Réale D, Clements MN, Morrissey MM, Postma E, Walling CA, Kruuk LEB, Nussey DH (2010). *An Ecologist's Guide to the Animal Model*. Journal of Animal Ecology **79**, 13–26. https://doi.org/10.1111/j.1365-2656.2009.01639.x.

Wolterock R (1909). *Weitere Experimentelle Unter-suchungen Uber Artveranderung, Speziell Uber Das Wesen Quantitativer Artunterschiede Bei Daphniden*. Verhandlungen der Deutschen Zoolоgischen Gesellschaft **1909**, 110–172.