Hybrid larch heterosis: for which traits and under which genetic control?

Despite the interest foresters have for inter-specific hybrid trees, still little is known about their quantitative genetics. This is especially true for the hybrid (HL) between Larix decidua (EL) and L. kaempferi (JL). Long-term, well-designed, multi-site experiments are necessary to estimate the parameters required for HL breeding programs. This paper presents the results from a diallel mating trial between nine EL and nine JL, set up in three contrasted sites. Growth traits (height, circumference), quality traits (wood density, stem form, heartwood proportion), and bud flush were measured from plantation to up to 18 years after plantation. Wood density and heartwood proportion were assessed using increment cores. We did a spatial analysis to take into account environmental heterogeneity at the tree level, and we fitted a multi-trait, Bayesian MCMC (Markov chain Monte Carlo) genetic model. Our study confirmed, in most situations, that HL expressed heterosis over its best parent for growth traits taking advantage of an early faster growth, with no loss in wood density. However, growth traits showed low levels of heritability. On the other hand, bud flush and stem flexuosity had high heritabilities, and wood density was clearly under JL control. Site-dependent heritabilities were expressed by EL. Additive genetic correlations were presented. The traits with high heritabilities showed high correlation between their performances in pure species and in hybridization, as well as high across-site correlations. The discussion focused on the interest of these genetic parameters for the hybrid larch breeding programs.


Introduction
Inter-specic hybridization is a common breeding strategy used in forest tree improvement. Commercial hybrids are produced for several species including poplars and aspens (Li et al, 1993), eucalypts, sub-tropical pines (reviewed by Nikles and Grin, 1991), and larches (reviewed by Pâques, 2013). For foresters, hybrids are often interesting because of their superiority over their parental species for several traits of ecological and economic importance (Zobel and Talbert, 1984).
Although there is an abundant literature on the topic of heterosis, this term remains somewhat ambiguous. Some authors prefer to associate heterosis to non-additive gene eects expressed in hybrids, and sometimes add the complementarity between parental traits under additive control to the denition of hybrid vigor (Nikles and Grin, 1991). Nevertheless, 'heterosis' was originally dened as a free-of-hypothesis concept, synonymous of 'hybrid vigor' (Shull, 1948). For a given quantitative trait, heterosis is considered as the dierence between the hybrid and the mid-parent value (Falconer and Mackay, 1996). The best-parent heterosis, i.e. the dierence between the hybrid and its best parent, is often preferred because it highlights perspectives in terms of breeding (Kearsey and Pooni, 1996;Dungey, 2001). In the present paper, heterosis will also be assumed a free-of-hypothesis concept, and will refer by default to best-parent heterosis. An extra ambiguity comes from the denition of 'hybrid'. Indeed, a 'hybrid' may refer to the cross between inbred lines, between populations, or between dierent species. Mechanisms underlying the heterosis of hybrids between inbred lines (e.g. maize lines) start to be well understood (Gallais, 2009;Baranwal et al, 2012;Chen, 2013); whereas the mechanisms behind the heterosis of inter-population and inter-species hybrids, i.e. hybrids between heterozygous parents, remain poorly documented (Perron, 2008)(see Li and Wu, 1996, though).
Despite the sparse information about the mechanisms, forest tree geneticists have well documented the superiority of some inter-populations and interspecic forest trees hybrids (Nikles and Grin, 1991;Li et al, 1993;Pâques, 2013). However, this superiority over some commercial references, which is the nal aim of a hybrid breeding program, does not provide information on the genetic architecture of the performances of the hybrid. At best, the study of heterosis by means of a proper mating design, and therefore proper parental references, allows the estimation of genetic parameters such as variance components (Hinkelmann, 1974) which can in turn be used to choose and develop a breeding strategy (Falconer and Mackay, 1996). In addition, the expression of heterosis may change with time and environmental condi-tions and its study, therefore, would require long-term multi-site experimentation. As an illustration, it can be demonstrated that genotype-by-environment interactions (G×E) alone can theoretically lead to heterosis including best-parent heterosis (Knight, 1973). Practical examples of the link between heterosis and G×E are also available. Thus, heterosis may be conditional to the site as shown by Li and Wu (1997) for aspen and by Weng et al (2014) for Eucalyptus urophylla × E. tereticornis. The G×E may not only manifest on the means but also on the components of variance. For instance, the parental contribution of Pinus caribaea var. hondurensis and P. tecunumanii to their hybrid heritability depends on the testing site (Mutete et al, 2015). Unfortunately, the creation of proper mating designs and the establishment of multi-site eld trials over contrasted environments are often complicated by the species biology (crossing of 2 dierent species), expensive and time-consuming. Therefore such assays of true heterosis are rare in the forestry literature.
Hybridization between several species of larches (Larix decidua Mill., L. kaempferi (Lamb.) Carr. and L. laricina (Du Roi) K. Koch) is possible and the hybrid between European (L. decidua ) and Japanese (L. kaempferi ) larches (respectively EL and JL) particularly shows promising prospects for forestry in Europe and North-America (Pâques, 1992b;Baltunis et al, 1998;Greenwood et al, 2015). Hereafter, 'hybrid larch' (HL) will exclusively refer to L. decidua × L. kaempferi or its reciprocal L. kaempferi × L. decidua. The superiority of HL for height over its parents has been described in the literature, as reviewed by Pâques (1989Pâques ( , 2013. Depending on the site, HL can outperform either EL, JL, or both. This superiority justies the existence of several hybrid larch breeding programs, in Europe and abroad (Perron, 2008). In addition to the heterosis on growth traits, HL inherits canker resistance from JL (Sylvestre-Guinot et al, 1999), which is of critical importance for lowland cultivation of EL in Western Europe.
However, the comparison of hybrid to its true parental references, i.e. heterosis sensu stricto, is poorly documented for larches (Pâques, 1989(Pâques, , 2013. Heterosis on growth traits was conrmed for HL at a juvenile stage (Baltunis et al, 1998;Pâques, 2002) and, more recently, on 22-years-old trees (Greenwood et al, 2015). In a breeding perspective, there is also a need for identifying the role of each parent species in the genetic control of the hybrid traits. Using 2 factorial mating designs, Pâques (2004) could identify the parental contribution to the phenotype of 16-years-old hybrid larches. He distinguished traits whose variability was more under Japanese parent control, such as stem volume and exuosity.
In summary, while hybrid superiority has been shown for hybrid larch, true heterosis is still poorly docu-

Manuscrit d'auteur / Author
Manuscript mented because of the absence of proper experiments, as discussed before. Multisite, long-term, well designed experiments are of primary importance in order to gather reliable estimations of the genetic parameters needed to plan a hybrid larch breeding strategy (Pâques, 1989). To our knowledge, hybrid larch heterosis has never been studied with such an experiment yet. The rst part of this study attempted to show for which trait and which amplitude one can expect heterosis from hybrid larch. True heterosis was assessed by the mean of a diallel mating design. The question of whether heterosis would be a juvenile character was addressed together with that of the impact of the environment on its expression. A better knowledge of the quantitative genetics of the traits should help us to better dene a breeding strategy, and to dene where and how eorts should be emphasized. For this reason, we investigated the importance of additive and non-additive genetic variances, the magnitude of genetic correlations, as well as the relative contribution of parental species in the genetic control of hybrid traits. We made a choice of solving such a complex genetic model in a Bayesian framework using MCMC techniques. Indeed, this framework brings the desirable feature of integrating estimation, prediction, and decision into a single analytical process (Gianola et al, 1989). Finally, the hypothesis of a better buering of environmental constraints by hybrids vs. parental species was tested through a preliminary study of stability of species and hybrids across environments (G×E interaction).

Experimental population
Nine (9) EL parents from Sudetan Mountains (Czech Republic) and 9 JL parents from various Hondo Island origins (Japan) were sampled from INRA (Institut National de la Recherche Agronomique) EL and JL breeding populations and crossed following a full diallel mating design producing pure species and hybrid full-sib progenies. Parents involved in the experiment were selected at random among trees able to ower and to produce cones. Out of the 171 possible families (that is 324 crosses, a family being a cross plus its reciprocal), self-crosses mostly failed and none of them were eld tested, and for 8 families (7 EL families and 1 HL family) no ospring were available. The remaining 145 families were represented by one cross and/or its reciprocal. In each site where the diallel was set up there were between 18 and 29 EL families, between 25 and 28 JL families, and between 56 and 71 HL families; the families were represented by 1 to 63 trees per site, with a mean of 23.2 trees. Thirty-six (36) extra full-sib and half-sib families, each having one of the parents of the diallel, were added to the main design. Details of available crosses are summarized in Appendix 1, Table 3-5. All seeds were sown at INRA nursery in Orléans (France) in spring 1995. Then 2-years-old bare-root seedlings were planted in 3 ecologically contrasted sites in spring 1997: Saint-Appolinaire (abbreviated SA), Saint-Saud (SS) and Cumières (not considered in this study because of a too high mortality), following a single-tree randomized incomplete blocks (RIB) design. To replace the Cumières site, 1311 trees were randomly sampled in SS, and 3 scions were collected from each of these trees and grafted on 2-year-old seedling. Because of the scarcity of available rootstocks, 2 scions were grafted on HL rootstocks (half-sib EL×JL hybrids) and the last one was grafted on a JL rootstock (from a commercial seedlot). The grafted trees were grown for 2 years in PNRGF (PÃ´le National des Ressources Génétiques Forestières) nursery at Peyrat-le-Château (abbreviated PC). Resulting material was planted in spring 2004 in PC vicinity, following a RIB design with 3-tree row plots (the 3 scions from the same ortet together). On all 3 sites, progenies from the 3 taxa were randomly allocated to each block. We favored this option in order to better take into account micro-environmental eects.
The site of SA is established at the eastern edge of the Massif Central range in a mountainous area on a steep southern-aspect slope. On the contrary, the sites of SS and PC are located at lower elevations on the western edge of Massif Central on mostly at areas, under oceanic inuence. The site of SS was a former grassland, frequently fertilized, whereas the other sites were already forest lands in the past. Description of the sites is provided in Appendix 1, Table 6. Two thinning treatments were applied to SS: in early 2004 and in 2010. During the rst thinning, 30% of the poorest performing trees (based on a selection index combining volume and exuosity) within each family were removed. The second thinning favored an even spatial distribution of remaining trees, like in forest production. No thinning was applied to the 2 other sites. Counts of trees per family and per site, as well as the evolution of the number of trees due to mortality and thinning, is provided in Appendix 1, Fig. 4.

Measurements
Trees were measured several times from plantation to winter 2010 in SA, and to winter 2014 in SS and in PC. Measurements included breast-height circumference (BHC), total height (HT), stem exuosity (FLX) and bud ush (BUD), and concerned all available trees at each assessment (except for SA in 2005). Stem exuosity and terminal bud ushing were visually assessed using subjective scoring scales, following Keiding and Olsen (1965)(see Pâques, 1992b)  Manuscript straight stem to 1: stems with several severe crooks) and Gauchat and Pâques (2011) (Fig. 1) for bud ush (from 0: dormant bud to 5: elongated needles, 1-2 cm in open rosette; plus an extra class 6: stem elongation). Each observation of bud ush at any given site was done in a single day. Diameter increment cores were harvested at breast height in the 3 sites over two or three collecting periods depending on site (Appendix 1, Table 7). All individuals in SA and SS trials were sampled, with a special eort in SS for harvesting increment cores from trees before being thinned. In PC, only 1 ramet per ortet was randomly sampled. Increment cores were dried o, sawn in 2mm-thick axial boards and then X-rayed on negative lms. Out of the two, only the pith to bark radius with the least defects (resin pockets, knots, compression wood) was kept for subsequent analysis. Microdensitometric proles were produced from high resolution scans of the X-rayed lms using WinDENDRO TM (Regent Instruments Canada Inc., 2008). From microdensitometric proles we visually delimited annual ring increments, which were used to estimate yearly BHC values. Inferred BHC from annual ring increments (∆) for year n + 1 was: with δ being an individual correction term, calculated as the ratio between dierence in radius calculated from BHC measures and dierence in radius calculated from ring increments, assuming a circular stem section and a constant bark thickness. Overall wood density (DEN n ) over years was calculated from the ring mean densities (d) weighted by the corresponding ring surfaces (S): we used as the rst ring the one for which enough observations were available (i 0 = 2000 in SA and SS and i 0 = 2007 in PC, that is constantly 4 years after plantation, for which at least 50% of the ring observations were available in each site).
Delimitation between sapwood and heartwood was visually assessed on the boards based on color dierences, and heartwood surface proportion (HWP) was calculated for 2009 and 2011 in SS, for 2013 in SA, and for 2014 in PC. All the traits (BHC, HT, DEN, HWP, FLX and BUD) were analyzed without any transformation.

Models
The data were analysed using a two-step approach. Firstly, we tted a model M1 that aimed to estimate the taxa means, from which we assessed best-parent heterosis at the species level. Model M1 was tted for each combination of site, age and trait; it took into account genetic and spatial environmental sources of variation. All 4 taxa (EL, JL, EL×JL and JL×EL) were analyzed simultaneously, including all available genetic material, at the tree level, so that (i) the contrasts between taxa were directly estimated and (ii) the spatial variance estimation used all available information within each trial. Estimating simultaneously genetic and spatial eects was computationally costly. For this reason M1 was a uni-trait model, meaning it ignored information across traits. However, traits are usually correlated both at genetic and environmental levels. This can be leveraged with a multi-trait analysis, gaining more insight on the genetic merits for each individual trait (e.g. Marchal et al, 2016), but also allowing for inference on the genetic cross-correlations. For these reasons, as a second step, we tted a genetic model M2 that was multi-trait. The spatial eects estimated in M1 were subtracted from the phenotype before processing to the second step.
Model M2 was tted for each combination of site and taxon; with the goal to estimate the genetic variance parameters such as heritabilities and genetic correlations between traits. Model M2 was also used to predict the genetic performances of the 18 parents of the diallel. The stability of these performances across sites, as well as their stability in pure species vs. in hybridization, were assessed by the mean of Pearson correlations. At this second step of M2, we used a Bayesian framework in order to infer directly the genetic parameters using posterior distributions. Attention must be drawn here on the terminology and concepts attached to the use of Bayesian approaches that dier from those in M1 that followed a frequentist approach. While the former makes decisions about the magnitude and relevance of estimates based on posterior distributions, the latter uses hypothesis testing and condence intervals. With the choice of a Bayesian framework in M2, we focused on the magnitude of the genetic parameters and the uncertainty around their estimations, because these are the pieces of information that ultimately can guide breeders in their selection decisions.

Model M1 -heterosis estimation
Each trait at each age in each site (y) was analyzed with the model M1 that took the form: with f ixed the xed eects at the taxon level, genet the individual random genetic eects, and envt the individual random environmental eects. The xed eects took the form: (M1a) The xed eects of M1 were β 0 , the EL×JL taxon mean used as a reference; and β E , β J , β J×E , the contrasts between each taxon t and the reference EL×JL. The indicative variable ι P C was dened such as the hybrid rootstock eect γ t:r were only accounted in the site PC. The hybrid rootstock eects γ t:r were calculated as contrasts between the trees whose rootstock r was HL compared to the trees whose rootstock r was JL; this eect included an interaction between rootstock species and scion species. The rootstock eects estimates in PC were presented in Appendix 2, Fig. 5. The genetic eects took the form: (M1b) For pure species, the additive eect or breeding value (BV) of the genotype i (that is, an individual tree in SA and SS and a clonal genotype in PC) was a E ∼ N (0, σ 2 AE A E ) or a J ∼ N (0, σ 2 AJ A J ) for EL or JL larches respectively, with A E and A J the additive relationship matrices. Dominance eects were d E ∼ N (0, σ 2 DE D E ) or d J ∼ N (0, σ 2 DJ D J ) for EL or JL larches respectively, with D E and D J the dominance relationship matrices. For hybrids, we used a Stuber and Cockerham (1966) model. The general hybridization abilities (GHA) from the mother j and the father k were h E ∼ N (0, σ 2 HE I HE ) and h J ∼ N (0, σ 2 HJ I HJ ), assuming unrelated, non-inbred parents (coancestry matrices for parents were half-identity matrices). This model also assumed that L. kaempferi (2n = 24) and L. decidua (2n = 24) contributed to the same extent to the hybrid nuclear genome, as Nkongolo and Klimaszewska (1995) showed in in vitro embryogenic lines of EL×JL karyotypes ranging between 2n = 24 and 2n = 25. The dominance variance was captured using a hybrid family eect, or specic hybridization ability (SCA): f ∼ N (0, σ 2 DH F) (Lo et al, 1997;Cros, 2014). As self-crosses were absent and diallel parents were assumed to be unrelated and non-inbred, F = 1 4 I F . A random clonal eect c was added (i) in PC for all taxa, in order to capture the epistasis and the nongenetic eects associated with each ortet (except for traits measured from the increment cores in PC, as no clonal repetition was available for them) and (ii) for all sites with HL only, in order to capture the within-family Mendelian segregation genetic variance and thus to separate it from the residual. The clonal eects had taxon specic variances.
Environmental heterogeneity was assumed to be spatially structured, so that the spatial eect at location l was s ∼ N (0, σ 2 S S) with S = {ρ distance(a,b) } where ρ had to be estimated for each combination of trait, age and site. Spatial correlation matrices S were built using the R package elds (Nychka et al, 2015;R Core Team, 2015). Finally e ∼ N (0, σ 2 E I R ) were the unstructured residuals, common to the 4 taxa; its variance was σ 2 E . All single-trait models were tted using restricted ML (REML) and the R package breedR .
For computational reasons ρ was estimated by Maximum Likelihood (ML) before the genetic eects were included in the model. Once ρ estimated, the spatial variance σ 2 S was obtained at the same estimation step as the other variance components. Estimates of ρ and spatial variance proportion (SVP) that we dened as Table 9. For 2 traits in PC (FLX at 11-years-old and BUD at 3-years-old), ρ was estimated to be low (0.69 and 0.63 respectively), correspondingly with very high SVP over 99% in both cases (results not shown). Such a parameterization of ρ may lead to a confusion between σ 2 S and the residual variance σ 2 E . Therefore, for these traits in PC we did not use the ML estimation for ρ but we xed ρ = 0.99 instead.

Taxa means and heterosis
Using model M1 (including all taxa simultaneously), we estimated directly the contrast between each taxon mean and the reference EL×JL. Doing this way, we measured the best-parent heterosis where appropriate, or the parental superiority otherwise; and we obtained a standard error for each contrast estimation. We used the standard errors of these contrasts to estimate condence intervals and we computed p-values (null hypothesis: âthe contrast is nullâ) under Gaussian distribution assumption.

Model M2 -variance parameters estimation
Model M2 was independently tted for each combination of site and taxon: EL, JL, and HL (EL×JL and JL×EL pooled together, assuming no reciprocal eects), on data previously adjusted for spatial heterogeneity. The response variables y of M2 were the observations of all p traits from which were subtracted the spatial eects predicted with model M1. In PC, we also subtracted the rootstock species eect γ estimated with M1. Multi-trait analysis is also known to account

Manuscrit d'auteur / Author
Manuscript properly for individuals that are censored (the case in SS due to thinning; and non-sampled individuals for increment cores traits in PC) in the estimation of genetic parameters (Wei and Borralho, 1998). Therefore in SS some traits (BHC, HT, DEN and FLX) were also included at a young stage in order to take into ac-count the thinning censure; they were then considered as new traits. All included BHC were eld measured, not inferred from ring increments. Exact ages for all traits involved in M2 are given in Appendix 1, Table  8. For pure species EL and JL, model M2 took the form: The overall means were m. A complete variance-covariance matrix was used to structure the additive eects distribution, that is a ∼ N (0, Σ A ⊗ A) with: The additive genetic correlation between traits r A were thus estimated as the model was tted. Dominance eects were independent between traits, that is d ∼ N (0, Σ D ⊗ D) with: Properly conditioned inverse matrices A −1 and D −1 were built directly using R packages MCMCglmm (Hadeld, 2010) and nadiv (Wolak, 2012). A complete variance-covariance matrix was also used to structure the residual eects distribution e ∼ N (0, Σ R ⊗ I R ) so that residual correlations between traits r R were estimated. In PC only, we added an unstructured, independent between traits, clonal eect c ∼ N (0, Σ C ⊗ I C ). Again, Stuber and Cockerham (1966) model was applied for hybrids, following model M2b: Additive variances in hybridization were structured using complete variance-covariance matrices h E ∼ N (0, Σ HE ⊗ I HE ) and h J ∼ N (0, Σ HJ ⊗ I HJ ) so that r HE and r HJ , the components of additive genetic correlation in the hybrids due to each of the parental species, were estimated. The family eects were independent f ∼ N (0, Σ D ⊗ F) Similarly to M2a, the residual effects distribution was r ∼ N (0, Σ R ⊗ I R ) Once again, we added an unstructured, independent between traits, clonal eect in PC.
The models were tted using Markov Chain Monte Carlo (MCMC) with the R package MCMCglmm (Hadeld, 2010). Means priors were vague with m ∼ N (0, 10 10 I p ). Variance-covariance matrices priors were inverse-Wishart (W −1 ) distributions. We used parameter expansion (PX) (Liu et al, 1998) for the additive components: . For hybrids, we used PX for Σ HE and Σ HJ with the same hyperparameters. On one hand, PX allows a better mixing of the chains. On the other hand, it avoids the modal shape of W −1 and it allows close-to-zero estimation for marginal variances (Gelman, 2006;Hadeld, 2010). The hyperparameters V = I p and ν = p + 1 were set following Gelman and Hill (2007) Manuscript priors remains undocumented. For this reason and to check any eventual impact we sampled 10 5 values in the parameter expanded inverse-Wishart distribution, with the same parameterization, and we used Monte-Carlo to draw the marginal correlation prior distribution. We visually assessed that it remained uniform on [−1, 1] (result not shown). Parameter expansion was also set for dominance variances so that for each trait σ 2 D α ∼ W −1 (V = 1, ν = 2) with working parameters α D ∼ N (0, 10 3 I p ). We used Sorensen and Gianola (2007)(p. 579) improper uniform prior for the residuals Σ R ∼ W −1 (V = 10 −10 I p , ν = −(p + 1)). Marginal residual correlations estimations are provided in Appendix 2, Fig. 7. In PC and for the hybrids, the clonal eect variance prior was also set improper uniform so that for each trait σ 2 C ∼ W −1 (V = 10 −10 , ν = −2). All chains were at least 1.5 × 10 6 iterations long (up to 10.5 × 10 6 long, depending on the taxon and site combination) with 5 × 10 5 iteration burn-in, and a regular thinning so that 10 3 uncorrelated samples from each chain were nally kept for inference.

Variance parameters
The genetic variance parameters were estimated from M2 only. For pure species, the phenotypic variance was σ 2 P = σ 2 A + σ 2 D + σ 2 R as in this case σ 2 R the residual variance from the model was equal to σ 2 E the environmental variance, i.e. the variance of non-genetic effects (epistasis was neglected). The narrow-sense heritability was h 2 = σ 2 A /σ 2 P . For hybrids we tted a model at the family level, therefore due to Mendelian segregation part of the genetic variance expressed at the individual level was captured by the residual variance σ 2 R . The individual level phenotypic variance was showing an analogy with pure species variance decomposition. In PC, we included the clonal variance Stuber and Cockerham, 1966;Cros, 2014), these heritabilities being analogous to half their pure species counterparts. Degree of dominance d 2 = σ 2 D /σ 2 P , residual correlations (r R ) and coecients of variations CV = σ P /m along with taxa means (m) were provided for all taxa in Appendix 2 (Figure 6, Figure 7 and Table 9 respectively). We used the variance components Markov chains to produce Markov chains of h 2 , d 2 and the genetic and residual correlations and we computed 95% credibility intervals (CIs) of these parameters.

Parents performances
For each trait and site, we predicted BV and GHA for the 18 parents of the diallel using model M2. Individual performances of both intra-and inter-specic crosses of the 9 EL and the 9 JL parents are presented with 95% CIs in Appendix 3, Fig. 8-10 for the most heritable traits (bud ush, exuosity and density). We assessed the across-sites stability of BV and GHA by the mean of Pearson correlations between the parents' performances in each site; moreover, a visual representation of the G×E stability of the 18 parents is provided for the most heritable traits in Appendix 3, Fig. 11. In the same way, we calculated the Pearson correlations between BV and GHA (r BV,GHA ) in order to assess the stability of performances in pure species vs. in hybridization. For both the across-site stability and the stability in hybridization, we used the parents' performances Markov chains to produce Markov chains of the Pearson correlations and we checked if their 95% CIs exceeded some given thresholds (0, 0.3, 0.6 and 0.9).

Importance of reciprocal eects and rootstock eects in the model M1
Although some reciprocal eects (β J×E , designed as a contrast between JL×EL and EL×JL, Cf equation M1a) were signicant (density in SS, 2 and 3 year height in SA) (p < 0.05 to p < 0.01), their magnitude was low (up to â8.29 kg.m-3 for density in SS and up to +6.23 cm for height in SA) (Fig 1, d, g, h). For all other traits (BHC (Fig 1, a-c), height and density in the other sites (Fig 1, e-f, g, i), exuosity (Fig 1, j-l), bud ush (results not shown) and HWP (results not shown)) reciprocal eects were very small and not signicant. Therefore, we took EL×JL as the reference (HL) for heterosis assessment, with negative contrasts (β J and β E in equation M1a) indicating heterosis, and positive contrasts parental superiority.
In the particular case of PC, we found strong, positive, signicant rootstock species eects (γ in equation M1a) for growth traits (Appendix 2, Fig. 5 a-b) (p < 10 −5 at the older age, for both BHC and height). Depending on the scion species, the eect of hybrid rootstock compared to Japanese rootstock varied between 7.3 cm and 5.5 cm for BHC at age 11. For height, there might be an interaction between the scion species and the rootstock species, JL scions seemed to undergo a lower HL rootstock eect. Thus for height the rootstock eect varied between 90.4 cm for EL scions and 48.3 cm for JL scions at age 11. Rootstock eects for the other traits were low though sparsely signicant (Appendix 2, Fig. 5

Heterosis
Growth traits (height and BHC) showed a clear bestparent heterosis (Fig 1, a-f) with the exception of BHC in PC. In SA, the best parental reference for growth traits was EL. In this site, at age 14, HL exceeded the best parental reference EL by 112.6 cm for height and 7.2 cm for BHC. Circumference contrasts evolved very slowly after age 14. In SS, at age 18, EL and JL were not signicantly dierent for growth traits; the contrast between HL and JL (respectively HL and EL) was 12.2 cm (respectively 17.3 cm) for BHC and 185.5 cm (respectively 197.7 cm) for height. Still in SS, growth traits contrasts showed a plateau after age 13. In both SA and SS, heterosis for growth traits was highly signicant (p < 10 −3 ) at the oldest age. In PC, the best parental reference for growth traits was JL. For trait BHC at age 11, JL signicantly exceeded HL of 2.8 cm (p < 0.05) and HL signicantly exceeded EL of 13.0 cm (p < 10 −5 ). Hybrid larch height exceeded both JL height of 55.6 cm (p < 0.05) and EL height of 254.8 cm (p < 10 −5 ).
Japanese larch was straighter (higher exuosity score) than the hybrid in all sites and ages (Fig 1, j-l). This dierence was signicant at age 14 in SA (dierence of 0.41, p < 0.05), at any age from age 4 in SS (maximum dierence of 0.58, p < 0.01), and at age 3 in PC (dierence of 0.43, p < 0.01). We also found a significantly higher straightness for EL than for HL in SS at age 4 (dierence of 0.31, p < 0.05). Density pattern was consistent across sites: EL was the densest, JL the least dense, and HL was intermediate (Fig 1,  g-i). These results were signicant in SA (p < 0.01), in SS for JL vs. HL (p < 10 −3 ) and EL vs. HL at age 5 then from age 8 onwards (p < 0.01), and in PC only for JL vs. HL until age 8 (p < 0.05). There was no clear trend with age for density contrasts, except in PC where dierences between taxa diminished with time.
Bud ush occurred slightly but signicantly earlier for HL than for EL at age 12 in SS (dierence of 0.39, p < 0.001) and later for HL than for JL at age 3 in SA (dierence of 0.32, p < 0.05). The other dierences between taxa for this trait were small and non-signicant (results not shown). In the same way, we found no clear pattern for HWP but a few contrasts were signicant: HL had a higher HWP than EL in PC (+7.16%, p < 0.05) and than JL in SA (+3.55%, p < 0.05); these were the highest contrasts for this trait.
In summary, best-parent heterosis was found in most sites for growth traits (HT and BHC). In PC, for BHC, we only showed average-parent heterosis. Most of these heterotic advantages were attained at early ages (usually before 15), reaching a plateau in subsequent later ages. Hybrids looked intermediate for density, and poorer than the parental references for exuosity.

Heritabilities
As previously detailed, genetic parameters were estimated in a Bayesian framework (no p-values and test of signicance). We considered heritabilities lower than 0.2 as low, between 0.2 and 0.4 as moderate, and higher than 0.4 as high. European larch pure species heritabilities pattern strongly depended on the site (Fig 2, a-c). In SA, all EL heritabilities were moderate to high (ĥ 2 E ∈ [0.232; 0.440]) (ĥ 2 E : maximum a posteriori estimate for h 2 E ); in SS only bud ush heritability was high (ĥ 2 E = 0.660), all the other heritabilities were low (ĥ 2 E ≤ 0.189); in PC heritability was high (ĥ 2 E = 0.471) for exuosity and moderate for height and bud ush (ĥ 2 E = 0.288 and 0.310 respectively). On the contrary, JL pure species heritabilities pattern was better preserved across sites compared to EL (Fig 2, d- There was also a weak contribution from JL to HWP heritability in SS (ĥ 2 HJ = 0.160). All other HL heritabilities (excluding young stage heritabilities in SS) were low (Fig 2, a-f).
Therefore, bud ush and exuosity resulted in the highest heritabilities across sites and taxa. Density heritability was high for JL and was under JL control in hybridization. Growth traits and HWP attained generally the lowest levels of heritability. The ranking in heritabilities among traits was better preserved in Manuscript JL across sites compared to EL. In particular, EL heritabilities both in pure species and in hybridization were generally higher in SA than in the other sites, even for growth traits. In general, EL had higher contributions to hybrid heritabilities for exuosity and bud ush than JL.

Degrees of dominance
Degrees of dominance were very low for all traits for all taxa at all sites (Appendix 2 Fig. 6). The highest degree of dominance we estimated was 0.136 for wood density of HL in PC; any other degree of dominance was under 0.089.

Additive correlations
Only a few additive genetic correlations CIs excluded 0 (Fig. 3), whereas a lot of residual correlations CIs were above or below 0 (Appendix 2, Fig. 7). Additive correlation between bud ush and exuosity was negative across sites and taxa (ranging from -0.581 to -0.073), meaning the later the bud ush, the straighter was the stem. This correlation was markedly negative in PC for pure EL species, with most of the posterior distribution below 0 (r AE = −0.396, 95% CI below 0 meaning P (r AE > 0) < 2.5%) and similarly for pure JL species (r AJ = −0.380, 90% CI below 0 meaning P (r AJ > 0) < 5%). On the other hand, in SA, bud ush and BHC additive correlation was positive for EL pure species (r AE = 0.501, P (r AE < 0) < 5%), meaning the earlier the bud ush, the larger was the stem.
In PC, additive correlation between height and density was positive for the pure EL and the pure JL species (r AE = 0.862, P (r AE < 0) < 2.5% andr AJ = 0.390, P (r AJ < 0) < 2.5% respectively). In the same way, in SS, we found a positive additive correlation between height and density in the pure JL species (r AJ = 0.514, P (r AJ < 0) < 5%). We also estimated positive additive correlations, though with some more uncertainty, between height and density for the pure EL species (r AE = 0.404) and both species contributions to hybridization (r HJ = 0.356 andr HE = 0.292). For JL in pure species as well as in hybridization, we always found negative additive correlation between BHC and density (ranging between -0.066 and -0.671). On the contrary, a positive correlation between BHC and density was found on the pure species EL side in PC (r AE = 0.769, P (r AE < 0) < 5%). Finally, additive correlation between BHC and height was positive for both pure species in SA (r AE = 0.712, P (r AE < 0) < 2.5% andr AJ = 0.484, P (r AJ < 0) < 5%); as well as in PC (r AE = 0.888, P (r AE < 0) < 2.5% and r AJ = 0.488, P (r AJ < 0) < 5%).

Stability of performances in hybridization
Stability (assessed as a Pearson correlation) between performances in pure species (BV) and in hybridization (GHA) showed a similar pattern to that of the heritability. For European larch, r BV,GHA were strongly site-dependent (Table 1). Indeed, for EL, r BV,GHA were moderate to high for all traits in SA, from 0.592 for height to 0.809 for HWP, 0.819 for exuosity and 0.909 for bud ush. In PC, r BV,GHA was high for exuosity (0.825) and for bud ush (0.722). Only bud ush had a high correlation in SS (0.801). For JL, still following the same trend as heritabilities, r BV,GHA were high for bud ush, exuosity and density in any site (between 0.807 and 0.918), except for exuosity in PC where r BV,GHA was only moderate (0.596). A graphical representation of performances correlations between BV and GHA, including performances CI, for bud ush, density and exuosity is provided in Appendix 3, Fig 8-10. In short, BV appeared to be a good proxy of GHA for the highly heritable traits: bud ush and exuosity, and density for JL. For the other less heritable traits, correlations were weak or, for EL, site-dependent.

Stability of performances across sites
Like the stability in hybridization, the stability of parents performances across sites (also assessed as Pearson correlations) were linked to the heritability of the traits; the most heritable traits were often the most stable across sites (Table 2). In pure species, EL showed high stability across sites for bud ush (between 0.763 and 0.918) only. However, the across-site stability of EL performance in hybridization was high for bud ush (between 0.894 and 0.911) and also for exuosity (between 0.726 and 0.854), for HWP and for density in SA vs. SS (respectively 0.785 and 0.772), and for height in SS vs. PC (0.834). Following heritabilities trend, JL pure species performances were highly stable across sites for bud ush, exuosity and density (between 0.687 and 0.955); across-site stability for height was also high in SS vs. PC (0.738). Across-site stabilities of JL performances in hybridization were high for density (between 0.817 and 0.899), but they showed a surprising pattern for exuosity and bud ush: although they were high in SA vs. SS (respectively 0.927 and 0.923), they were only moderate in the comparisons involving PC (i.e. SA vs. PC and SS vs. PC) (between 0.643 and 0.681). Finally, JL across-site stability for height was moderate to high in hybridization in SA vs. SS (0.805) and SS vs. PC (0.672). A graphical representation of the across-site stability of the performances in hybridization is provided in Appendix 3, Fig. 11. Briey, JL showed more stable performances across sites than the EL counterpart for most of the assessed traits. Across-site correlations were particularly high in JL for highly heritable traits: bud ush, exuosity and density. This trend of JL was paralleled in its contribution to hybrid performance across sites, with similar correlations and for the same set of stable traits. European larch, with less stable performances except for bud ush, showed a higher across-site stability in its contribution to the hybrid performances.

Discussion
Interspecic hybrids between European and Japanese larches present many desirable features in terms of growth and stability that are not observed simultaneously at parental level. Many hybrid varieties ex-ist, but hybrid larch as many other hybrids lacks to a large extent recurrent breeding eorts, notably by the absence of a formal breeding program. The rst steps towards fullling this absence are the production of quantitative genetic knowledge on hybrid larch performance. Our study was based on one of the very few intra-/inter-specic diallel mating designs existing for forest trees.
Our study showed clear heterosis for growth traits (height and BHC). This advantage occurred already at a juvenile stage, with apparently little or no extra heterosis gained after age 15-16. Those traits that beneted the most from heterosis tended to show low to medium heritabilities, and were among those in the study being less stable across sites and with relatively low correlation between performances in pure and hybrid crosses. Heterosis over the best parent was absent from the traits showing the highest heritabilities  in the study, namely bud ush, stem exuosity and wood density. These traits were stable across sites and showed generally high correlation between performances in pure and hybrid crosses. These results, although based on a relatively narrow base population, are of key importance when planning hybrid breeding for larch.

Hybrid larch heterosis
Best-parent heterosis in growth traits (BHC and total height) was validated in all sites, except in PC for BHC where HL almost reached its best parent (JL), though it overpassed it for height. As a result, regardless if the site was more favorable to EL or to JL, the hybrid yielded more than or as much as the best parent species for growth traits. For BHC, the relative best-parent global heterosis was 15.7% in SA (best parent: EL), the least 'fertile' site and 12.0% in SS (best parent: JL), the most 'fertile' site. In PC (similar 'fertility' as SS), the hybrid BHC was 6.3% lower than that of JL and 25.8% higher than that of EL ( Fig. 1 and Appendix 2, Table 9). Using the multi-site diallel mating design trial we described here (though, including Cumière site instead of PC), heterosis was already shown for growth traits at a juve-  (Pâques, 2002). Near-rotation age (at year 22) heterosis between EL and JL has also been shown using a diallel mating design trial, in a single site in Central Maine, USA (Greenwood et al. 2015). Like us, Greenwood et al (2015) found no reciprocal eect between JL×EL and EL×JL, and they demonstrated best-parent heterosis at the taxon level for height, BHC and volume over proper parental references. The novelty of our study was to show heterosis at near-rotation age in several contrasted sites, using a constant genetic material. In SA and in SS, after age 15-16, hybrid superiority over its parental species reached a plateau for growth traits and more clearly for BHC. In SS only, the last thinning campaign could play a role in this phenomenon. Nevertheless, this result comforts the hypothesis that heterosis would be expressed at a juvenile stage, and that the acquired superiority would be conserved along the tree's life. Indeed, it is now well-known by foresters that hybrid larch in pure plantations reaches its maximum mean annual increment

Manuscrit d'auteur / Author
Manuscript for stem volume 5 to 20 years earlier than its parental species (Rondeux, 2003). Also, in trials where taxa are closely inter-mixed as in our study, one can expect that hybrid with its rapid growth could have overgrown its parental counterparts in the neighborhood and prevented them from catching up at later ages. As an illustration, a link between genetic performances and competitiveness has been shown by Cappa et al (2015) in loblolly pine, with a high correlation (CI: 0.77 to 0.92) between BHC breeding value and competition ability. We failed to t such a competition model (Muir, 2005;Muñoz and Sánchez, 2015) to our data. The model led to unreliable results, likely due to the mix of species and the over-parametrization, and it was thereby disregarded (result not shown). Presumably, competition could have interfered with the phenomenon of heterosis for growth traits, notably for BHC which is known to be aected by competition.
We did not nd a clear pattern of best-parent heterosis for the other investigated traits. Hybrid wood density was intermediate between both parental values; this trend was clear and signicant in SA and SS but more confused in PC. Similar results for hybrid wood density were also found in Pinus hybrids (Dungey, 2001). In larch, Charron et al (2003) already reported a lower hybrid wood density when compared to EL counterparts. Nevertheless, the fact that faster growing HL still produced higher density than its worst parent for this trait (JL) is promising for further breeding and for the sustainability of larch wood products quality. Concerning the other traits, hybrid seemed to be less straight than its parents, notably when compared against JL. For bud ush and HWP, the ranking of taxa was not conclusive (result not shown).
Some singularities aected PC site, namely the absence of best-parent heterosis for BHC, and the confusion of density ranking that was otherwise clear in the other sites. Given that PC was the only site with grafted trees, one can suspect complex interactions between rootstocks and scions to be one of the underlying causes of this singularity. Rootstock species eects were found for growth traits (Appendix 2, Fig.  5). This nding indicates that the use of grafted trees may introduce an additional unwanted source of variation that should be taken into account in any case. On the other hand, we showed for growth traits that any species grafted on a hybrid rootstock may express some heterosis that led to a gain in height and in BHC after 11 years. This nding suggests thus that heterosis results not only from favorable functional properties of the aerial part but also from the below-ground part of the tree.

Heritabilities
Some traits showed high heritabilities: bud ush, stem exuosity and, for JL only, wood density. The gain in selection relies on high heritabilities as well as on phenotypic variances (Falconer and Mackay, 1996). Across sites and taxa, the estimated phenotypic standard deviations ranged from 0.79 to 1.22 for exuosity on the 1-5 scale; from 0.66 to 1.18 for bud ush on the 0-6 scale (Appendix 2, Table 9). Therefore, it is possible to improve these traits through breeding. However, wood density standard deviation was quite low (ranging from 28.6 to 42.0 kg.m −3 ) so the perspectives of gain for this trait remain low.
For all taxa, height heritability was often intermediate, and for BHC and HWP, heritabilities were mostly low. The general trend was not perfectly followed by EL however, showing a site-conditional expression of its heritabilities. In particular, BHC heritability was moderate for EL in SA (ĥ 2 E = 0.356) and accompanied by a phenotypic standard deviation of 101.9 mm. We also found a very low EL heritability for stem exuosity in SS. This may be explained by the thinning that occurred in this site. Indeed, thinning censure could not be recovered from young stage observations because of the quite low additive correlation between young and old stage exuosity for EL in SS (r AE = 0.455). On the other hand, JL heritabilities were generally better preserved across sites, with high magnitudes for density, bud ush and exuosity. Interestingly, considering all taxa, the most heterotic traits (BHC and height) were associated with low to moderate heritabilities, whereas non-heterotic traits such as bud ush, stem exuosity and wood density were highly heritable. This suggests that biomass related traits (BHC and height) could mostly be improved at the hybridization stage, whereas bud ush, density and exuosity could be improved by recurrent selection within the pure species. Indeed, these later traits showed generally high stability of the performances in hybridization. Hybrid wood density was clearly under JL additive control. Only HWP showed a dicult path towards genetic improvement, with no heterosis and a very low level of heritability. However, as absolute heartwood size follows stem size (Pâques and Charpentier, 2015), total heartwood size may be improved along with BHC through hybridization.
We compared our narrow-sense heritabilities estimates (Fig. 2) to the ones from the literature which were reviewed by Pâques (2013) (Appendix 4, Fig.  12). Only EL and HL heritabilities were available for a few traits in this review. The EL heritabilities (h2E) that we estimated were well within the range of variation of the other authors' estimates for bud ush, height and exuosity, except for exuosity in SS for which our estimation was noticeably low. On

Manuscrit d'auteur / Author
Manuscript the other hand, BHC heritabilities were lower than the ones found in the literature, except in SA where it was consistent with published values. Hybrid heritabilities (h 2 HE + h 2 HJ ) were close to the literature references for exuosity and density, but our study obtained lower estimates for height and BHC than those already published. Only Hering (1990, in Pâques, 2013 found also low heritabilities for these growth traits. More generally, moderate heritabilities for height and high ones for wood density could be expected for tree species (e.g. see Cornelius, 1994;Kain, 2003), but according to these studies our estimated exuosity heritabilities were unusually high and the BHC ones were unusually low.
Heritability underestimation, notably for growth traits, could be due to the inter-specic competition in our experimental design. Indeed, a potential limitation of our study came from our single-(or in PC three-) tree plot design, which may have enhanced competition between taxa. Some authors found a link between competition and depletion of heritability. Lin et al (2013) showed a depletion on breast-height diameter heritability at age 5 in Pinus radiata under competition induced by a high plantation density. Also at high plantation density, Bouvet et al (2003) observed a depletion of both BHC and height heritabilities in Eucalyptus hybrids. Although this latter depletion was not signicant, this eect remained stable from age 4 onwards and it was validated in 3 independent trials. The question remains nevertheless on how far interspecies vs. intra-species competition aects heritability estimates. Some of our results had no correspondence to those in Pâques (2004). The wood density, which we found mostly under Japanese larch control, was estimated to be under European larch control by Pâques (2004). Unlike this latter publication, we found little heritability for growth traits (height and BHC), and stem exuosity being under European or both parents control. Several reasons may contribute to explain this. On one hand, our ndings highlight the need to assess these parameters in several environments, while Pâques (2004) results come from a single site, which was dierent from ours. On the second hand, the mating design was more unbalanced in Pâques (2004) and included dierent genetic basis with EL from the Alps. Pâques and Charpentier (2015) found HL HWP heritability (full-sib family mean h 2 ) to be lower (0.49) than heartwood diameter heritability (0.68), but still higher than our narrow-sense heritability estimates (ĥ 2 E , h 2 J andĥ 2 HE +ĥ 2 HJ ranging between near 0 and 0.232). This could suggest that for this trait there was either little genetic variance in our sampling, or lack of accuracy in our measurement method.

Dominance
The dominance variances and the subsequent degrees of dominance we estimated were very low (d 2 ≤ 0.136), even for traits showing heterosis. Kain (2003) also found very low d 2 (between 0 and 0.1) for 17 traits in the 2 Pinus pure species he investigated, however he obtained high dominance to additive variance ratio in Pinus elliottii var. elliottii × Pinus caribaea var. hondurensis hybrid for some traits including breast height diameter and volume. Dungey (2001) pinpointed that in several species the relative importance of dominance variance tends to diminish with age. This trend was also found in hybrid larch (Pâques, 1992a). In the present article, we presented the analysis of data collected on old trees, generally older than in many of the reviewed articles in Pâques (2013), which overall leads to particularly reliable estimations of the additive genetic performances.

Additive correlations
Bud ush was genetically correlated with exuosity. We found a negative additive correlation between these two traits, meaning the later the bud ush, the straighter the stem was. We also found a positive additive correlation between bud ush and BHC in SA, pinpointing to the need for a trade-o between BHC and exuosity if selecting on bud ush. We also showed positive additive genetic correlation between height and density in SS and in PC. This could be an eect of the competition within these 2 sites. Indeed, though regularly thinned, SS is from far the most fertile site, and PC is also more fertile than SA (see Appendix 1, Table 6, and taxa means ranking in Appendix 2, Table 9). Therefore we could expect competition to be stronger in SS and PC than in the somehow poorer SA site. However, the hypothesis that competition leads to positive correlations between height and density would deserve further investigation. Wood density variability was under JL control, and for this taxon we found a negative correlation between density and BHC, in pure species as well as in hybridization. This trend was only supported by moderate statistical evidence in hybridization in SA; but it might pinpoint at the possible trade-o between wood volume and wood quality, an issue raised by Kain (2003) with hybrid pines. However, on one hand, gain in wood productivity using best-parent heterosis at the taxon scale was not accompanied by a loss of density, as hybrid density was intermediate between its parental species. On the other hand, based on heritability dierences alone, a genetic improvement of wood density would be more eective than that of wood volume.

Stability of performances
Correlations between BV and GHA and across-sites correlations followed the same pattern as heritabilities. This could be expected as low heritability traits lack the genetic information that is needed to properly assess these correlations. Conversely, we found high correlations between BV and GHA for high heritability traits. This result was in line with Dieters and Dungey (2000)'s suggestion that the correlation between BV and GHA may be positively related to additive to dominance variance ratio, as we found very low d 2 and h 2 d 2 = Correlations between BV and GHA and across-site correlations are both valuable parameters for guiding the choice of breeding strategies. On one hand, the high correlations between pure species and hybrid performances indicate that breeding eorts at the parental species level can be easily translated into comparable progress at the hybrid stage. On the other hand, we expect the genetic gains to be stable across sites, especially for exuosity and bud ush in hybridization. Moreover, the high proportionality between JL density performances in pure species and in hybridization together with the stability of JL density performances across sites are very valuable; they allow a eld screening of JL trees for their genetic value for this trait in hybridization.

About the method
The MCMC algorithms used in this study are rarely used in quantitative genetic studies of forest tree species. However, these algorithms present many computational advantages. Overall, they tend to be robust solvers even when low level of genetic information is available. Using PX prior, we could handle close-to-zero variances, an issue sometimes met with average-information REML algorithm (e.g. see Mutete et al, 2015). The prior parametrization we used was exible enough for variances as well as for covariances components. Heritabilities estimates ranged from almost 0 to 0.660, showing that phenotypic variance could be properly partitioned according to data information. In the same way, correlations estimates ranged between close-tozero and close-to-one (plus or minus one, max(|r A |) = 0.888 and max(|r R |) = 0.895); for the PX (genetic) posteriors (Fig. 3) as well as for the improper (residual) posteriors (Appendix 2, Fig. 7). The possibility to make direct inference on derived parameters such as heritabilities and genetic correlations was very powerful. We could also measure the reliability of our correlations between BV and GHA, an issue raised by Dungey (2001).
In order to get rid from the inter-taxa competition interference, it could be interesting to re-think our experimental design for further study of such contrasted taxa. Pure taxa micro-plot design might be an alternative solution. However, border trees would still be exposed to inter-taxa competition, so they would have to be excluded -raising an extra-cost. Moreover, the spatial eect would be in confusion with the taxon eect. Spatial analysis methods that are less susceptible to high frequency eects, such as splines , could be a solution for the analysis of such an experimental design. Another solution would be the use of large pure taxa macroplot design. However, taxa and environment eects would be confounded; using clonal repetitions as control in each macro-plot would produce an extra-cost because of border trees; and an extra-diculty because of the potential interactions between the genetic material and the environment.