Abstract
The link between phenotypic plasticity and heterosis is a broad fundamental question, with stakes in breeding. We report a case-study evaluating temporal series of wood ring traits of hybrid larch (Larix decidua × L. kaempferi and reciprocal) in relation to soil water availability. Growth rings record the tree plastic responses to past environmental conditions, and we used random regressions to estimate the reaction norms of ring width and wood density with respect to water availability. We investigated the role of phenotypic plasticity on the construction of hybrid larch heterosis and on the expression of its quantitative genetic parameters. The data came from an intra-/interspecific diallel mating design between both parental species. Progenies were grown in two environmentally contrasted sites, in France. Ring width plasticity with respect to water availability was confirmed, as all three taxa produced narrower rings under the lowest water availability. Hybrid larch appeared to be the most plastic taxon as its superiority over its parental species increased with increasing water availability. Despite the low heritabilities of the investigated traits, we found that the expression of a reliable negative correlation between them was conditional to the water availability environment. Finally, by means of a complementary simulation, we demonstrated that random regression can be applied to model the reaction norms of non-repeated records of phenotypic plasticity bound by a family structure. Random regression is a powerful tool for the modeling of reaction norms in various contexts, especially perennial species.
The link between phenotypic plasticity and heterosis has been well documented for annual plants, notably maize (Janick 1999; Gallais 2009). Usually, hybrids are able to stay productive in environments where one or both parents’ performances drop. This enhanced stability is referred to as ’hybrid homeostasis’, and can be responsible for the environmental conditionality of the expression of heterosis (Knight 1973). Comparing larch growth in several contrasted sites, we recently highlighted homeostasis for hybrid larch (HL, Larix decidua Mill. × L. kaempferi (Lamb.) Carr., and the reciprocal cross), a highly productive conifer cultivated for wood in Western Europe and North America (Marchal et al. 2017). Because the plasticity of some traits can be a tool for the homeostasis of the whole organism, phenotypic plasticity is suspected to play a key role in the construction of hybrid larch heterosis.
Phenotypic plasticity, in its narrow-sense definition, is the ability of a genotype to produce several phenotypes depending on the environmental conditions. It can be studied by means of reaction norms (Schlichting and Pigliucci 1998). A reaction norm is an equation, or simply a graphical representation, of the value taken by a phenotypic trait along an environmental gradient. Estimating reaction norms is not a trivial task. In particular, exposing the same genotype to different environments may be experimentally challenging, depending on the biological model. In some cases, the only solution is to expose related individuals to the different environments, relying on their genetic relationship to draw the common, additively inherited, component of the reaction norms (e.g., Gibert et al. 2004; Valladares et al. 2006). The random regression model is an extension to the classical quantitative genetics model (Kirkpatrick and Heckman 1989), and as such it can predict the additive component of reaction norms, and give access to causal components of population variation. In addition, Murren et al. (2014) demonstrated with a large meta-analysis that when comparing reaction norms of close species or populations, changes in shapes (i.e., slope, curvature) were generally higher than the changes in the intercepts (i.e., the taxon means), suggesting the need of high-order modeling. In that view, random regression modeling has also the valuable capability to fit complex curves for reaction norms.
Random regression is a special case of covariance function (Meyer 1998). Covariance functions present a particular interest in quantitative genetics, as they allow the representation of quantitative genetics parameters as functions. For instance, the cattle breeding literature is rich in illustrations of heritabilities and genetic correlations estimated as functions of time in the milk production context (e.g., Miglior et al. 2007; Muir et al. 2007; Jamrozik et al. 2010). Although covariance functions and random regression have been often suggested for the modeling of reaction norms (Kirkpatrick and Heckman 1989; De Jong and Bijma 2002; Schaeffer 2004), their application is still rare for most taxa (Morrissey and Liefting 2016). For instance, in the forestry context, random regression has been used to model growth over time (Apiolaza and Garrick 2001; Wang et al. 2009), but not tree growth reaction norms over environmental gradients until very recently (Carnwath and Nelson 2016; Marcatti et al. 2017). Li et al. (2017) advocated for the use of random regression in forest trees, but the growth reaction norms they reviewed were only performed at the population scale, with fixed effects, without covariance functions. Nevertheless, because they are sessile, long-living, and because they record radial growth increments in the form of annual rings, trees are remarkable biological models for the study of phenotypic plasticity. The use of random regression has been suggested for the study of tree rings plasticity (Sánchez et al. 2013).
For a given tree, the succession of wood annual rings constitutes an archive of its plastic responses to the succession of past climatic environments. Indeed, trees growth is conditioned by solar radiation, temperature and precipitation. They are particularly sensitive to drought, as this can eventually lead to lethal embolism of their hydraulic system. As a response to water stress, trees have a range of plastic responses including modification of the cambial activity to control the sap flow rate and xylem resistance to embolism (Bréda et al. 2006; Rennenberg et al. 2006). Under a temperate climate, water availability is usually high in spring with cool temperatures and lower in summer with higher temperatures. Trees respond to this seasonal succession by producing an early wood, made of large cells with thin walls, progressively followed by a late wood, made of narrow cells with thick walls, allowing a rapid change in the trunk of the water conductance along the season (Sánchez-Vargas et al. 2007; Martinez-Meier et al. 2009; Bryukhanova and Fonti 2013). This succession of early and late wood is recorded in the wood as an annual growth ring.
In summary, trees present the remarkable feature of recording long repeated series of phenotypes that are known to respond plastically to the climatic environment, and differentially among individuals. This provides a gold mine of information for the study of phenotypic plasticity. First, based on a random regression approach, we leveraged this information to address the question of the role of phenotypic plasticity in the construction of HL heterosis for radial growth. Indeed, integrative heterosis as observed for HL stem circumference necessarily arises from the cumulation of heterosis occurring at the annual ring scale: this might be expressed differentially with respect to parental genotypes depending on water availability. Second, we described how the water availability affected the quantitative genetic parameters of HL and of its parental species. Finally, since most biological models do not archive repeated series of phenotypic plasticity records, we addressed the question of the generalizability of the method. Using simulations, we evaluated the robustness of the random regression approach for the modeling of reaction norms in the case where no repeated series of phenotypes per individual are available.
Materials and Methods
The data are part of a multi-site progeny trial established in early 1997 on two environmentally contrasted sites at Saint-Appolinaire (SA, 45∘58’N 4∘26’E, 784 m a.s.l.) and Saint-Saud (SS, 45∘31’N 0∘48’E, 307 m a.s.l.) in central France. The site of SA is a relatively high-elevation site, on a steep slope with a southern aspect formerly planted with Douglas-fir; whereas the site of SS is a low elevation site on a former meadow, under oceanic influences.
Progenies were produced by control-crossing in the frame of a diallel mating design between 9 European larch (Larix decidua, EL) parents and 9 Japanese larch (L. kaempferi, JL) parents, producing pure species and HL full-sib progenies. A total of 327 EL, 472 JL and 1199 HL genetically distinct trees were ultimately available for analysis in SA; and 311 EL, 706 JL and 1261 HL in SS. The mating design is detailed in Supplementary 1, Tables S1 and S2, with the number of genotypes available for each combination of parents. The set-up is also described further in Marchal et al. (2017).
Phenotypic data: wood formation records
One breast-height diameter increment core was collected from each tree from each site. For each diameter increment core, only one radius (the one exhibiting the fewest defects) was kept for further analysis. These radial increment cores were sawed in 2 mm thick boards and X-rayed to obtain microdensitometric profiles (Figure 1). From the alternating of early wood and late wood, the year of formation for each ring was identified (Regent Instruments Canada Inc. 2008). Ring width (RW) and ring mean density (RMD) were measured from the microdensitometric profiles. A total of 1998 increment cores in SA and 2278 increment cores in SS were collected. In SS, the increment cores were collected at three different periods: before the first thinning in 2003 from trees to be felled, and two later ones before and after the second thinning, that is in 2009 and in 2011. In SA, the collection was done in 2013. The average number of rings available per core was 13.2 in SA and 9.7 in SS (Supplementary 1, Fig. S1).
Microdensitometric profile from one wood core sample in Saint-Appolinaire. Bottom: wood core sample (pith on the left). Middle: X-ray radiography of the sample. Top: wood density variation (kg/m3) along the core (mm). The peaks correspond to late wood and delimit the annual growth rings.
Environmental data: soil water availability
The soil daily relative extractable water (REW) content was estimated using a water balance model. The REW varies between 1 (field capacity) and 0 (permanent wilting point, meaning the water is no longer available for the plants). We used an adapted, simplified implementation of Granier et al. (1999)’s daily water balance model. This implementation allowed two important features. First, it accounted for the temporal evolution of the stand. Indeed, during the growth of the trees, the canopy leaf area increases and subsequently increases the transpiration, whereas understorey shrubs and grass vegetation decrease. Second, the model could be run at the individual tree level. Indeed, the intra-site spatial heterogeneity of the soil storage capacity causes variation in the water availability at the individual tree level. Thus, the daily REW was estimated for each combination of individual tree within each site and year for which ring data were available. The water balance model is described in Supplementary 2.
To summarize a year of fluctuation of REW, we selected the lowest decile of the whole series of the year’s daily REW values (D1rew). This index was considered because of its easy interpretation (the REW below which lay the ∼35 driest days of the year), because of its likelihood to carry the information of drought events and therefore to affect the growth, and finally because it ensures a proper coverage of the environmental gradient (Supplementary 1, Fig. S2), unlike e.g., indexes based on a drought threshold that may not be reached some humid years. The environmental gradient coverage guarantees the stability of the regression parameters and the quality of the subsequent analyses.
The index D1rew, derived from a complex water balance model, was compared to a simpler index MJJA defined as the sum of rain from May to August (inspired from Fallour-Rubio et al. (2009)), from which was subtracted the daily potential evapotranspiration (Supplementary 1, Fig. S2). Moreover, Bryukhanova and Fonti (2013) showed that for European larch, several traits including RW were more strongly correlated to the soil water deficit of year 
than to the deficit of the current year t. Therefore, we modeled the reaction norms along D1rew for year t, for year 
, and MJJA for year t as three separate environmental gradients, and we compared the coefficient of determination (
) (Nakagawa and Schielzeth 2013; Johnson 2014) of the resulting reaction norms.
Modeling of reaction norms by random regression model
Reaction norms were modeled using orthogonal Legendre polynomials (Kirkpatrick et al. 1990; Schaeffer 2004). Let 
be the 
order Legendre polynomial of 
with 
the standardization of x on [-1, 1], that is 
We fitted the following model for each taxon (EL, HL or JL):
where 
was the 
observation of individual k, of genotype j, from site i, with environment x. The site’s effect s was fixed. The additive effects a were random, and depended on the species. For EL and JL pure species, respectively, 
and 
, where ⊗ indicates a Kronecker product, A were the additive relationship matrices computed from the pedigree, and 
were the estimated variance-covariance matrices for additive effects. For the hybrid, 
with g the additive contributions from each parental species. Thus, 
on the EL side and 
on the JL side, with 
the estimated variance-covariance matrices for additive effects in hybridization (Stuber and Cockerham 1966). Given that the parents were supposed outbred and unrelated, the resulting relationship matrices 
and 
reduced to identity matrices. The permanent environment, i.e., the similarity from non-additive genetic origin between the repeated measures of the same individual tree (Mrode and Thompson 2005), was 
. The residual was unstructured 
We chose to fix 
so that an ’order M’ applied to the whole model. Models of different orders are characterized by how much information on plastic responses is included: order 0 does not estimate any plasticity, order 1 calculates slopes, order 2 additionally fits parabolas, and so on.
This model was fitted for both traits RW and RMD separately (univariate models), and also for RW and RMD simultaneously (multivariate model). In this latter case, the effects s, a and p were estimated for each trait; the variance-covariance matrices Σ gathered the variances of each combination of trait and order and the covariances between them all; the residual variance 
was independent for each trait.
The model was fitted by Markov chain Monte Carlo (MCMC) with the same priors as in Marchal et al. (2017) (Hadfield 2010; R Core Team 2017): parameter expansion was used on the genetic variance-covariance matrices, and flat improper priors were set on the permanent environment variance-covariance matrices as well as on the residual variances. All chains were 
iterations long, with 
iterations burn-in and a thinning of 
. Point estimations from chains were maximum a posteriori, and 95% credible intervals (CIs) were computed where appropriate. The quality of fitting depending on the order was assessed using coefficients of determination 
.
Estimation of genetic parameters
The additive variances (
and 
) and covariance between RW and RMD (
) were calculated from the multivariate random regression model as functions of the environmental gradient x. The additive variance-covariance matrix decomposes into 
where 
and 
are the variance-covariance sub-matrices for each trait, 
is the across-traits covariance sub-matrix, and the operator: t indicates a transposition. The additive variance of a trait T is the variance of the additive performances within the population for this trait, that is the variance of a linear combination:
Where 
The permanent environment variance 
was computed in the same way. Similarly, the additive covariance between the 2 traits was computed as:
Finally, the permanent environment covariance 
was computed in a similar way. Narrow-sense heritabilities and additive correlation were then computed respectively:
and:
For hybrids, narrow-sense heritabilities and additive correlations were computed for each of the parental contributions 
and 
These heritabilities 
and 
quantify the proportion of hybrid phenotypic variance that is due to the additive inheritance from each parental side (Marchal et al. 2017). In the same way, the additive correlations 
and 
are the correlations between the additive contributions for the 2 traits from each parental side. Therefore, on the EL side and given that 
, for a trait T:
from which:
and for the covariance parameters between the 2 traits:
from which:
and idem on the JL side.
Simulation: evaluation of the random regression model with single record per individual
We used the software Metagene to simulate data. The functioning of this simulator to derive in silico populations with phenotypic plasticity records is detailed in Supplementary 3. Basically, the genotypic effect at a given locus was set as a function of the environment 
where parameters 
and δ defined the parabola that was associated to each genotype in a set of X diallelic loci constituting the genome. Then, for each independent simulation, the simulator randomly sampled a genome for each founder, produced the mating between founders, the new offspring genomes, and returned their phenotypic plasticity in the form of longitudinal records over Y environments. We parameterized the allelic effects in such a way that the additive reaction norms were very interactive, that is, the ranking of the parents varied along the environmental gradient due to slopes and parabolas.
The simulated mating design consisted of a full diallel between 10 monoecious founders, excluding selfs. Unlike the real mating design previously described, the simulated mating design involved a single species so there was a single additive variance to estimate. Each combination of parents (A×B) produced 
sibs, so that the size of a full-sib family (A×B + B×A) was n. The environmental gradient was divided in n random positions, each position being defined by an environmental value x. Each of these environments hosted one sib per family, and as many individuals as families. The progenies’ phenotypes, their pedigree, and the environmental values x were included in the analysis. As there were no repeated observations of the individuals, the permanent environment was unidentifiable, it was neither simulated neither included in the analysis.
We tested 4 different scenarios, resulting from the combination of low (0.1) and high (0.6) heritabilities with small (20) and large (120) family sizes. Thus, these scenarios measured the importance of the quantity (i.e., the number of progenies) and quality (i.e., the heritability) of information for genetic inference. For each scenario, 100 independent simulations were run and analyzed. The analysis was performed using pure species univariate random regressions:
where 
was the only observation for individual i, from environment x, μ were the parameters for the mean reaction norm, a were the additive effects such as 
and r was the unstructured residual.
The chains were 
iterations long, with 
iterations burn-in and a thinning of 
We measured the ability of the random regression to infer additive components of the parental reaction norms. To do so, we defined the accuracy of the model as the correlation between the parents’ predicted additive performances and their true simulated additive performances (Mrode and Thompson 2005), and we computed it at each point of the environmental gradient. The reaction norms of the parents were only inferred for the ranges on which the progenies were tested.
Prediction of the larch reaction norms with single records per individuals
Finally, we also investigated the possibility to infer reaction norms from single records per individuals using the empirical larch data. We randomly sampled one observation per tree in the dataset, and used it to fit the random regression model. The model we used was the same as for the complete empirical larch dataset, except without the permanent environment component that became unidentifiable.
Data availability
Our datasets are available at the INRA repository GnpIS https://data.inra.fr/dataset.xhtml?persistentId=doi:10.15454/REMJWZ&version=DRAFT (https://doi.org/10.15454/REMJWZ).
The Metagene simulator is available on the NOVELTREE project page https://data.inra.fr/dataset.xhtml?persistentId=doi:10.15454/ELCSZI (https://doi.org/10.15454/ELCSZI). Supplemental material available at Figshare: https://doi.org/10.25387/g3.7211597.
Results
Quality of model fitting to real data
The order 0 model (fixed and random intercepts) explained between 8.9% and 12.1% of RW variance depending on taxa and site combination (Figure 2). The addition of fixed and random slopes (1st order) along D1rew of the current year (t) greatly improved the model, allowing it to account for 39.9% (EL in SA) to 55.8% (HL in SS) of the variance. Analyzing higher orders (order 2 and order 3 along D1rew of the current year) slightly increased the 
(by 5.1% on average compared to order 1). Such high 
were not reached using D1rew of the previous year, nor using the simpler index MJJA. For RMD, on the contrary, order 0 or higher orders led to very close 
(Figure 2). Using D1rew of the previous year instead of that of the current year led to a marginal improvement of 
that happened only in SA, reaching up to 39.4% for HL with order 3.
Coefficients of determination (
) for ring width (RW) (a) and ring mean density (RMD) (b), depending on the model specification. The polynomial order (M) varied from 0 to 3. Three covariates were compared: the first decile of the soil daily relative extractable water (D1rew) for the current year (year t, implicit) and for the previous year (year 
, specified), and the sum of daily differences between precipitation and potential evapotranspiration from May to July (MJJA) for the current year. The last model specification (’Multi.’) corresponds to 1st order regression with the covariate D1rew of the current year t, but in this case the model was multivariate and RW and RMD were analyzed simultaneously.
Heterosis expressed mostly in radial growth, that was better explained by D1rew of the current year. For this reason, we focused on this environmental gradient only, and all the subsequent results refer to random regression models along D1rew at year t. The gain in 
with orders beyond 1 was overall due to more variance explained by the fixed component of the models, as shown in Figure 3. Therefore, we present in Figure 4 and Fig. S10 (Supplementary 4) the taxon scale reaction norms (i.e., the fixed component) estimated from the order 3 univariate random regression models. Then, for parsimony consideration and because increasing the order did not increase the portion of 
captured by the random effects (Figure 3), the multivariate random regression models (RW and RMD analyzed simultaneously), from which the genetic variance and covariance parameters were estimated, were fitted with order 1. The 
s of the multivariate models were also presented in Figure 2 and were very similar to their univariate counterparts.
Decomposition of the coefficients of determination (
) obtained with the order 1 model (a) and with the order 3 model (b), for each trait: ring width (RW) and ring mean density (RMD), and for each combination of site (SA or SS) and taxon (EL, HL or JL). The environmental gradient was the first decile of the daily relative extractable water for the current year. The proportion of variance explained by each component is presented: genetic additive effects (in black), permanent environment (in white), and fixed terms (in gray).
Reaction norms of ring width (a, b) and ring mean density (c, d) along the first decile of the daily relative extractable water (D1rew), in the sites SA (a, c) and SS (b, d), for each taxon: European larch (in blue, EL), Japanese larch (in red, JL) and their hybrid (in green, HL). These average reaction norms at the taxon level represent the fixed components of the order 3 random regressions. Dashed lines: 95% credible intervals.
Stability of the heterosis and its dependence on the water availability
The ranking in RW performance between EL and JL varied depending on the site: EL performed better in SA while JL did better in SS (Figure 4). Meanwhile, the superiority of the hybrid over both its parental references, hence best-parent heterosis, occurred in the two sites and over the whole range of D1rew.
For any taxon and in any site, RW was plastic as it increased with increasing water availability (Figure 4). The three taxa showed however different curves along the D1rew, being close to each other when water availability was minimal (D1rew close to 0), and splitting apart with increasing water availability, where the superiority of HL over its parental references was the highest. Thus, HL had the steepest reaction norm in any site, as shown in Table 1. The gain in superiority for RW of HL over its parental references due to increasing D1rew ranged between +3.26 mm (HL vs. EL in SA) and +1.04 mm (HL vs. JL in SS). For high D1rew, the CIs of the HL reaction norm were not overlapping with those of the parental species.
On the opposite, the trait RMD showed neither heterosis nor plasticity. The hybrid ranged between both its parents, and all the norms of reaction were almost flat for this trait, showing no conspicuous pattern of variation along the water availability gradient D1rew.
Heritabilities and genetic performances along the water availability gradient
All the narrow-sense heritabilities we estimated were very low (Figure 5). Low heritabilities were overall due to high residual variances in comparison to the lower additive and permanent environment variances (Figure 3). Despite this residual noise, some extreme parental performances of contrasting genotypes were different for both traits with statistical credibility (Supplementary 4, Fig. S6 - S7).
Narrow-sense heritabilities along the first decile of the daily relative extractable water (D1rew), for European larch (a, c) and Japanese larch (b, d), for the traits: ring width (a, b) and ring mean density (c, d). These heritabilities were computed from the variance parameters estimated with the multivariate, order 1 random regressions. Black line: narrow-sense heritability in pure species crosses; and gray line: narrow-sense heritability in hybridization (i.e., proportion of hybrid phenotypic variance that is due to the additive inheritance from each parental side). Dashed lines: 95% credible intervals.
Ring width heritabilities were close to 0 (Figure 5). The signal for performance contrasts for RW in pure species was also very weak, but both species showed contrasted performances in hybridization as the water availability increased; some of these contrasts were supported by non-overlapping 95% CIs when D1rew was high (Supplementary 4, Fig. S7).
Heritabilities for RMD were higher than those of RW, especially on the JL side for which they varied between 0.08 and 0.12 both in pure species and in hybridization (Figure 5). Pure species heritability and heritability in hybridization were very close for JL, and they both increased with water availability, though the increase was not supported by statistical confidence (the CIs were large and encompassed the whole variation). Increasing contrasts between individual performances, supported by some non-overlapping 95% CIs for high D1rew (Supplementary 4, Fig S7), may reflect this possible gain in heritability along the D1rew gradient for RMD on the JL side. Noteworthily, the ranking of the 9 JL parents’ performances for RMD was consistent in pure species and in hybridization (Supplementary 4, Fig S6).
Correlations along the water availability gradient
The additive genetic correlation between RMD and RW showed similar patterns between pure species and their respective contributions in hybridization (Figure 6). The additive correlation decreased slightly along the gradient from around 0 to -0.41 in EL pure species (-0.35 in hybridization). On the JL side, it started from a positive correlation of 0.48 in pure species (0.14 in hybridization) and it steeply switched to a negative correlation of -0.60 (-0.41 in hybridization). Uncertainty was high around this parameter, however, for D1rew high enough, the 95% CIs excluded 0 (completely or mostly) on the JL side (in pure species or in hybridization, respectively) (Figure 6). Due to the higher genetic variance for RW for JL in hybridization than in pure species, the negative correlation pattern was especially visible in the JL performances in hybridization: for the highest D1rew, the ranking between RW and RMD was almost inverted (Supplementary 4, Fig S6).
Additive correlations between ring width and ring mean density along the first decile of the daily relative extractable water (D1rew), for European larch (a) and Japanese larch (b). These correlations were computed from the variance and covariance parameters estimated with the multivariate, order 1 random regressions. Black line: correlation in pure species crosses; and gray line: correlation in hybridization (i.e., correlation between the additive contributions for the 2 traits from the corresponding parental side). Dashed lines: 95% credible intervals.
The permanent environment correlation between RW and RMD was negative for HL and JL, and did not vary along the environmental gradient. It was also null to negative for EL (Supplementary 4, Fig S8). This means that the sum of effects that the model did not explicitly account for (i.e., micro-environment, competition between trees, non-additive genetic effects, etc.) tended to induce a negative correlation between radial growth and wood density.
Simulation: accuracy of the random regression model With single record per individual
Using simulated data, we evaluated the ability of the random regression model to predict the additive component of reaction norms from family series of single observations per environment. The accuracy of the model depended on the scenario and on the order of the random regression (Figure 7). The 
scenario (
progenies, 
) showed very poor predictive abilities, independently of the order. When a larger number of progenies (
) or a higher heritability (
) was available, the accuracy was still low for order 0 (fixed and random intercept model) but it greatly increased with order 1 (addition of fixed and random slopes). However, in each case (
or 
) no further accuracy was gained from order 1 to order 2. Only the accuracy for the last scenario (both 
and 
) increased with order 2 (addition of fixed and random parabolas). The accuracy for the last scenario analyzed with order 2 model ranged between 0.905 and 1 depending on the simulation and on the position along the environmental gradient.
Accuracy of the predictions of parents’ additive reaction norms in each simulated scenario: (1) 
and 
(a, e, i), (2) 
and 
(b, f, j), (3) 
and 
(c, g, k), and (4) 
and 
(d, h, l); for 100 simulations in each scenario, and for each order of random regression: order 0 (a-d), order 1 (e-h) and order 2 (i-l). Each gray curve is the accuracy of 1 simulation.
The order 0 models were not able to estimate the additive variance properly (Supplementary 4, Fig. S9). From order 1 and over, the ability of the model to estimate the additive variance appeared overall dependent on the heritability. Indeed, the estimated additive variance was close to the true one in scenarios 3 and 4 (both 
) with order 1 or 2. With lower heritability (
) and order 1 or 2, the additive variance was overestimated.
Finally, we also fitted random regression using only a single observation per tree from the empirical dataset. As show in Supplementary 4, Fig. S10, the average taxon-scale reaction norms thus estimated were close to the ones estimated using the whole dataset (Figure 4). The uncertainty was much higher using a single observation per tree, however, both the slope of the reaction norms and the ranking of the taxa were well preserved.
Discussion
In this study, we constructed the reaction norms of annual wood-formation traits along a water availability gradient in a larch multi-site diallel mating experiment. The reaction norms were fitted using random regression modeling, allowing to investigate the changes in heritability and genetic correlations along the gradient. Our study was complemented by using simulations involving the same analytical approach, where we evaluated the ability of the random regression model to estimate the additive component in a frequent phenotypic plasticity experimental setting: that of parental reaction norms from non-repeated observations of their progenies as data.
The annual ring width was plastic and as expected increasing water availability allowed a higher radial growth (Figure 4). At the taxon level, the hybrids performed better than their parents in each of the two sites of the study. This alternative modeling confirmed our previous findings (Marchal et al. 2017), in that hybrid larch had a stable superiority across sites. However, within each site, HL demonstrated more plasticity than its parental references: indeed, under water stress all taxa produced a similarly narrow ring, whereas in favorable conditions of higher water availability the HL expressed superiority over its parental references. In other words, the HL reaction norms were steeper than those of the parental counterparts over the same environmental gradient.
On the contrary, the second trait under study, RMD, was not plastic for this water availability gradient and expressed no heterosis in any site. Although RMD displayed higher narrow-sense heritabilities than RW, it did not reach high values, with a global maximum of only 0.12. Despite the weak genetic signal, we represented the genetic parameters and the performances of the 18 parents of the diallel as functions of the water availability gradient. Some performance functions were distinct with statistical credibility, only when the water availability index was high enough. Interestingly, on the Japanese larch side, the additive correlation switched from being uncertainly positive to reliably negative as water availability increased (Figure 6). The seemingly emergence of this negative correlation may be explained by a genetically variable increase in the early wood / late wood ratio with increasing water availability, with early (spring) wood being generally less dense than that of late (summer) wood (Figure 1). However, we need to look more carefully to other ring traits (as did, e.g., Bryukhanova and Fonti 2013) and their respective correlations before proposing any causal explanation, with the goal to better understand the structure of the genetic variability of larch wood plasticity.
We obtained fairly high 
for the reaction norms models, suggesting that water plays an important role in the tree ring phenotypic plasticity. We evaluated a simpler environmental factor (sum of daily rain minus potential evapotranspiration from May to August) but it showed a lower 
highlighting the relevance of our water availability index ’D1rew’. It should be noted, however, that our water balance model has not been field-calibrated, and it should then be considered with care if generalizations are to be made. Moreover, the index D1rew gives no indication on the distribution of the driest days along the year. The timing of a water deficit, in spring or in summer, could have more or less effect on different ring traits; for instance, those relative to early or late wood, or to the transition between the two. It also has to be said that the relation between water balance and radial growth that we showed in the present study does not necessarily imply a direct causality. Indeed, other factors may play a role in the observed plastic response. For instance, heat affects directly the photosynthetic efficiency and the resources that may be allocated to growth (Rennenberg et al. 2006). Heat and drought being highly correlated, their effects could well be confounded to some extent. More broadly, the environment has a multivariate nature. Soil, climate, but also competition with neighboring trees are known to affect the tree’s growth. We isolated what we expected to be one of the most important environmental factor for radial growth, yet the existence of an important site effect pinpoints the fact that some other environmental factors might be involved in the tree ring phenotypic plasticity. Identifying relevant environmental factors of plant plastic reactions is an open area of research, notably in the context of global warming. The present study did not aim explicitly at the identification of relevant environmental triggers, rather it presented an approach that could help in such identification.
We studied phenotypic plasticity at two levels. The first level was spatial, at the across-site scale, and the second level was longitudinal, at the individual scale. At the across-site scale, heterosis was shown to be stable, supporting the common statement that hybrids are more stable across macro-environmental sites than their parental counterparts (Gallais 2009). Specifically, the ranking of the parents species varied across sites whereas hybrid was invariably the highest performing taxon, in what could be qualified as hybrid homeostasis according to the theory developed by Knight (1973). The spatial plasticity is generally, and historically, the one that interests breeders the most because of the operational implications for the deployment of varieties. Using non-linear random regressions, Marcatti et al. (2017) fitted eucalyptus growth reaction norms along a gradient of spatial environments. The spatially distributed climatic environment was described with principal component analysis to account for its multivariate nature. This method is very appealing in order to deal with spatial phenotypic plasticity in tree breeding. However, Marcatti et al. (2017)’s approach could be further improved by accounting for genetic information, such as the one from a pedigree and from which quantitative genetic parameters could be estimated (Hinkelmann 1974; Falconer and Mackay 1996). In this context, the use of random regressions as defined by Kirkpatrick and Heckman (1989) and presented in this paper, or other covariance functions, would be relevant.
The second level of phenotypic plasticity that we addressed was longitudinal, along the year-to-year water availability gradient. The resulting reaction norms were strongly influenced by the sites. Indeed, our study steps in the direction that tree ring longitudinal phenotypic plasticity should be seen as a plastic trait in itself, varying spatially (De Luis et al. 2013; Natalini et al. 2016), varying with long-term trends such as global warming (Natalini et al. 2016), and varying with the level of competition between neighboring trees in wet years (Carnwath and Nelson 2016). Growth recovery, the ability for trees to produce large rings the years following a drought event, is also site-dependent (Gazol et al. 2017). Besides the site effect, a substantial part of the individual variation was shown to occur between taxa. Unlike what we initially expected, hybrid larch reaction norms for radial growth along the water availability gradient were not flatter than the ones from its parental species: on the contrary, they displayed more plasticity and were steeper than the parental ones. Interestingly, this increase in longitudinal plasticity for HL incrementally contributes to the construction of heterosis that was observed on an integrative scale across years, namely in the total circumference (Marchal et al. 2017).
The molecular and physiological mechanisms behind the increased plasticity of hybrid larch remain open questions. Longitudinal plasticity as shown here is still a novel approach, with less immediate application to current plant breeding, unlike spatial plasticity. It certainly opens up new possibilities whenever long-time series are available, where extreme events, such as the 2003 drought in France (Bréda et al. 2006; Rennenberg et al. 2006), are recorded. Some initiatives, for instance, compared dead trees vs. alive neighbors immediately after extreme climatic events for their past wood records (Martinez-Meier et al. 2008), finding that both classes had long-term distinctive patterns of reaction. This kind of study could well be undertaken with a random regression approach to gain insight in the quantitative genetics of such longitudinal patterns and their environmental drivers, and be of potential use ultimately for breeders.
One eventual problem with longitudinal data are autocorrelation. In order to minimize its effects, several authors propose an extra step that consists in fitting an autocorrelation model on the tree ring series. The resulting residuals, that are more independent than the raw data, are then used as the response variable in the subsequent phenotypic plasticity models (Bryukhanova and Fonti 2013; De Luis et al. 2013). In this study, we did not do so because the chronology was somehow already involved in the environmental variable. Indeed, as the trees grew older, the stand LAI increased, and so did the transpiration, making water generally less available. Though this trend was not so strong (rain and potential evapotranspiration during the growing season were the main drivers of D1rew, as seen on Supplementary 1, Fig. S2), we did not want to account twice for the same chronological effect, and therefore we chose to work with raw data instead. The much weaker explicative power of D1rew of the previous year compared to that of D1rew of the current year supports our decision. We acknowledge though that there remains a risk of non-controlled autocorrelation in the data, in particular due to the trees’ ontogeny and to the onset of competition between trees (Sánchez et al. 2013).
The simulations showed that it was possible to estimate the additive component of reaction norms using only singular observations of related individuals. However, the quantity and the quality of information (i.e., respectively, the number of related individuals and the heritability) were key factors to estimate properly the additive components of reaction norms. Although the simulation was not meant to mimic the real case in the genomic layout of effects, it pinpointed the eventuality of potential biases in the estimation of the genetic variances (Supplementary 4, Fig. S9). As emphasized by Misztal et al. (2000), a limitation of the model used in the present study is the lack of covariance function for the residuals. This limitation could be a source of bias in the estimation of variance components and of heritabilities. Unfortunately, this feature was not available yet with the software we used. Indeed, we fitted linear mixed models in which the covariance functions were implicit and computed from the covariance between the regression coefficients. On the other hand, despite this issue with the variance, the simulated additive reaction norms could be properly estimated in most scenarios. Moreover, we also showed that the average reaction norm at the taxon level could be estimated from a single observation per genotype using the empirical larch data.
Random regression is already used for the modeling of reaction norms, especially in dairy cattle for which industry produces a large flow of longitudinal data (e.g., Kolmodin et al. 2002; Windig et al. 2006; Santana et al. 2017). However, plasticity has been suggested (Bradshaw 1965) and demonstrated (Murren et al. 2014) to be of special importance in plants. Indeed, because they are sessile, plants have to face their environment in a different way than animals that are capable of behavioral responses and locomotion. Manifestations of phenotypic plasticity have been reported in several perennial crops. For instance, grape vine manifests phenotypic plasticity in terms of fruit weight and chemical composition (Dai et al. 2011). Even in equatorial regions, oil palm is able to react to subtle variations in photoperiod and drought events by changing its bunch productivity (Legros et al. 2009). Cherry tree phenology reacts promptly to climate, notably heat, with global warming expected to bring flowering a month forward (Allen et al. 2014). Ismaili et al. (2016) showed a significant genotype-by-year interaction in apricot tree, and recommend the use of mixed models for the analysis of perennial plants’ longitudinal data. All these examples could be good candidates for analysis based on covariance functions. Indeed, the possibility to define quantitative genetic parameters as functions of the environment and to model the additive contributions to phenotypic plasticity opens wide perspectives in terms of selection (De Jong and Bijma 2002), especially in a global warming context (Koski 1996). In that sense, Chevin and Hoffmann (2017) pinpoint the importance of inferring the complex functions underlying plastic responses, notably those modeled under extreme environments, as essential steps to address evolutionary questions related to the amount of genetic differences in plasticity to extreme events, and to what extent these reflect heritability. The same authors also argue that comprehending the evolutionary dynamic of phenotypic plasticity would require a fine knowledge on the genetic constraints that could operate across environments. Our study contributed to illustrate methodologically several of the perspectives highlighted by Chevin and Hoffmann (2017).
Like trees, other organisms naturally accumulate growth records that reflect their reaction to past environmental conditions, etched in hard organs that grow incrementally: for instance, fish otoliths, mollusc shells, corals, whale ear plugs, ibex horns, etc. (reviewed by Morrongiello and Thresher 2015). Morrongiello and Thresher (2015) advocate for the use of random regression for the analysis of such natural records of longitudinal data, especially with regards to these species’ phenotypic plasticity. On a prospective review concerning new integrative ways to assess phenotypic plasticity, Stinchcombe et al. (2012) already highlighted the potential advantages of using random regression in the context of non-linear continuous reaction norms, notably to study growth curves matched by chronological age.
The random regression framework as developed by Kirkpatrick and Heckman (1989) exhibits two particular strengths that deserve to be emphasized once more. First, the use of orthogonal base functions, such as Legendre polynomials, allows the fit of virtually any shape of growth curves or reaction norms. Although the example we presented here did not illustrate this need, neglecting curvature when studying evolutionary divergence in reaction norms leads to a risk of missing critically important information (Murren et al. 2014). The traditional use of x and 
as covariates should be avoided in any case, given the high correlation that binds the identity and the square (and any power) functions.
A second point of interest is the fact that genetic information can be taken into consideration in the model, which can be of relevance not only for breeding but also in ecology studies looking for drivers and patterns of natural selection (Brommer et al. 2005). As we illustrated with the results of our simulations, this possibility opens up the use of a random regression framework to species whose individuals do not cumulate in any known form longitudinal records of their plastic responses. In this sense, our simulation provided an example of random regression being an alternative to traditional methods: related individuals can indeed give access to the additive component of the reaction norms, and this can likely be extrapolated to isofemale lines (Gibert et al. 2004) or half-sib families (Valladares et al. 2006). Finally, the pedigree information can be conveniently replaced by molecular information (e.g., Ly et al. 2018), extending the potential of the random regression framework beyond the limitation of our capability to realize time-consuming, sometimes impossible, artificial mating. This also opens up the possibility to realize genome-wide association studies in order to detect the genetic variants involved in phenotypic plasticity.
Acknowledgments
The authors sincerely acknowledge the technical staff of INRA experimental units (UE GBFOR and UE BIOGECO) who have established, maintained and assessed the field trials as well as collected the increments cores. The authors thank Genobois platform technical staff who have prepared and X-rayed wood samples and provided microdensitometry profiles. We thank Philippe Rozenberg for his implication all along the paper. We thank the think-tank ’PlasPhen’ for providing opportunities to move forward in our thinking on this paper, and in particular Patricia Gibert and Vincent Debat who organized it. We thank Isabelle Cousin and Ghislain Girot for their help with the interpretation of the soil data. We thank Vincent Ducrocq for his advices with the random regression model. We thank Météo-France and the platform INRA CLIMATIK for the availability of the climatic data. F. Muñoz was partially funded by research grant MTM2016-77501-P from the Spanish Ministry of Economy and Competitiveness. F. Muñoz and L. Sánchez received funding from the European Union’s Seventh Framework Program for research, technological development, and demonstration under grant agreement no. 284181 (’Trees4Future’, coord. L.E.Pâques).
Footnotes
Supplemental material available at Figshare: https://doi.org/10.25387/g3.7211597.
Communicating editor: E. Josephs
- Received August 24, 2018.
- Accepted October 10, 2018.
- Copyright © 2019 Marchal et al.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
