Accounting for competition in multi-environment tree genetic evaluations: a case study with hybrid pines

Abstract

Key message

A novel multi-environment competition individual-tree mixed model resulted in better fit, and greater individual narrow- and broad-sense heritabilities than the model without competition, notably for traits showing competition. In multi-environment tests, the proposed model would allow tree breeders to select genotypes with the best performance in both additive direct and competition breeding values, increasing forest productivity.

Context

Genetic merit of trees is known to be affected and interact with local competition effects as well as changes across environmental conditions. Recent studies showed that competition genetic effects can affect the genetic variance and bias the tree breeding values, and its covariance with direct breeding values has been variable across traits.

Aims

The present paper extends a mixed-model methodology to the problem of accounting for competition in a multi-environment set of forest genetic trials and exploring its impact on genetic variances as well as the multi-environment genetic correlation.

Methods

The proposed model is illustrated using data from two full-sib trials of Pinus elliottii var. elliottii × Pinus caribaea var. hondurensis F₁. Dispersion parameters and (co)variance of total breeding values were estimated for diameter at breast height, total tree height, and stem straightness at age 10.

Results

For traits showing competition effects (diameter at breast height and total tree height), the proposed multi-environment competition model gave better fit than the simpler model. Accounting for competition increased the direct additive variance, reduced the residual variances, and did not change significantly the across-site additive genetic correlation. However, for diameter at breast height, top 5% best genetic rankings showed differences.

Conclusion

When traits are strongly affected by inter-tree competition, the use of the proposed model in multi-environment analyses can efficiently identify the phenomenon with general benefits in the fitting of genetic components and open the door to select on the basis of competitiveness.

1 Introduction

The productivity of plantations depends on the genetic potential of the species and the interaction of the species with edaphic and climatic conditions, as well as other site management decisions (McKeand et al. 2006). Therefore, genotypes superior in one environment may not be correspondingly superior elsewhere, or the degree of difference among genotypes may vary according to the environment (Wu and Matheson 2005). Hence, tree improvement programs usually involve evaluating similar sets of genotypes across a range of sites (Cappa et al. 2012), with implications in testing, deployment, and subsequent selection strategies (Zas et al. 2004). Developing efficient plant breeding strategies requires a knowledge of the magnitude and structure of the genotype by environment interaction (G × E; e.g., Cullis et al. 2014). Moreover, the relevance of G × E might increase in a context of continuous climatic changes that can generate pronounced differences between environments (Morais Júnior et al. 2017). The magnitude of G × E can be examined by genetic correlations among environments as obtained from a mixed-model analysis of multi-environment forest genetic trials (MET thereof; e.g., Dungey et al. 2000; Costa e Silva et al. 2005; Cappa et al. 2012).

Competition is defined as the negative impact on plant growth or fitness caused by the presence and interplay with its neighbors, usually by reducing the availability of resources. In trees, it may be mediated by the genotypic composition and spatial arrangement of a tree’s neighbors (Hinson and Hanson 1962). Furthermore, it can be divided into genetic and environmental sources (Magnussen 1989). Under competition, the superiority of best genotypes tend to be overestimated due to greater use of available resources at a cost of the worst genotypes which exhibit greater sensitivity to competition (Resende et al. 2005). Considering the long rotation periods of forest trees, competition is likely to be an important factor to consider in genetics analysis (Costa e Silva and Kerr 2013). Recent studies in forest tree breeding have used a univariate single-environment competition mixed model to investigate direct and indirect genetic relationships among trees (Cappa and Cantet 2008; Costa e Silva et al. 2013; Costa e Silva and Kerr 2013; Cappa et al. 2015, 2016). Also using a single-environment approach, Costa e Silva et al. (2017) and Hernández et al. (2019) fitted a multivariate competition mixed model to address the direct and indirect genetic correlations among traits.

These studies have shown that ignoring genetic and/or environmental competition effects tended to underestimate the additive genetic effects and overestimate the residual estimated variances (e.g., Cappa and Cantet 2008; Costa e Silva and Kerr 2013), hence biasing predictions. The same approach could be extended to multi-environment analyses where across site estimates could also be affected by inter-tree competition. Moreover, as tree improvement involves often selections across many environments, it is important to consider also the impact on covariances among trials for direct and competition genetic effects to better understand to what extent the inter-tree competition affects the within- as well as the across-sites selection accuracies. Consequently, tree breeders could be interested in determining to what extent multi-environment ranking might be affected by ignoring inter-tree competition. Therefore, the multi-environment standard mixed model must be adapted to be used with genetic and environmental competition effects, which can be referred to as “competition MET analysis.” The use of a multi-environment competition individual-tree mixed model has not been reported in forest breeding programs.

The current study extends the standard MET individual-tree mixed model to a competition setting. Developments are illustrated for growth and stem straightness traits from two full-sib progeny tests of Pinus elliottii var. elliottii Engelm × Pinus caribaea var. hondurensis (Senecl) Barrett and Golfari F₁ hybrid pine assessed at age 10. The resulting parameter estimates (including across-site additive and dominance genetic correlations) for the proposed multi-environment competition individual-tree mixed model are compared with corresponding estimates from the multi-environment standard individual-tree mixed model without competition effects. We also investigated for each trait the variance and across-site correlation of total breeding values (i.e., additive direct plus additive competition), and the influence of inter-tree competition on the additive and dominance genetic correlations among environments (G × E). The rankings of selection candidates calculated for the multi-environment standard mixed model and multi-environment competition mixed model were also compared to determine the importance of accounting for competition in multi-environment analyses.

2 Materials and methods

2.1 Genetic material, description of trials, and quantitative traits

The Pinus elliottii var. elliottii Engelm (PEE) × Pinus caribaea var. hondurensis (Senecl) Barrett and Golfari (PCH) F₁ hybrid data set used in the present study was obtained from the tree improvement program initiated by the National Institute of Agricultural Technology (INTA, Argentina) in 2004. The program currently comprises 17 trials of F₁ hybrid progenies, where two of the oldest were chosen here for the competition MET analysis. Both trials were established in 2007 in Puerto Laharrague and Colonia Delicia (Misiones, Argentina), and for the present study they were identified as sites 1 and 2, respectively (INTA 2020).

The full-sib families used in these trials were originated by crossing 13 PEE mothers from INTA Clonal Seed Orchard (CSO) with 8 PCH fathers to form a factorial array. However, given that not all crosses produced enough seeds, only 66 full-sib families out of possible 104 were included in these trials. The number of common parents between both trials was 20 (i.e., only one mother was not represented at trial 1) and the number of common full-sib hybrid families was 37 (out of a total of 66). All PCH parents were plus-trees situated in Misiones Province, four of which were selected from two provenance trials planted by INTA, while the remaining four parents corresponded to selections from commercial plantations belonging to PINDO S.A. of unknown provenances of origin. In addition, the two trials used here included five check lots of unknown parentage which represented about 10% of the total of phenotypic data: (1) PEE from the CSO of INTA (San Antonio, Misiones), (2) PEE from CSO of INTA (Cerro Azul, Misiones), (3) PCH from Brazil, (4) PEE × PCH F₂ hybrid from the Queensland Forestry Research Institute (QFRI, Australia), and (5) Pinus taeda Marion from CSO of Arauco Argentina S.A. (Puerto Libertad, Misiones).

The F₁ hybrid trees were planted in a randomized complete block design with five and six replicates, respectively, for trials 1 and 2, involving five-tree row plots per full-sib family, and with an initial spacing of 3.0 × 3.0 m between trees. All surviving trees at both trails were measured at age 10 for diameter at breast height in centimeters (DBH, 1.3 m above ground), total tree height (TH) in meters, and stem straightness (STR) in a 6-point subjective scale with 1 indicating the most crooked trees and 6 indicating the straightest trees. The categorical traits (STR) were transformed into normal scores (NS; Gianola and Norton 1981) to meet the requirements for normal distributions and renamed as NSTR. Table 1 provides general information about the two trials. A detailed description of the genetic materials and trials used in the present study can be found in Belaber et al. (2018).

Table 1 Location, design information, survival, and means (standard errors in brackets) for diameter at breast height (DBH), total tree height (TH), and stem straightness (STR), measured at age 10 across the two Pinus elliottii var. elliottii × Pinus caribaea var. hondurensis F₁ hybrid progeny trials

2.2 Statistical analysis

Preliminary analyses of the two F₁ hybrid pine trials indicated that the most competitive trait (or the one most sensitive to competition) was DBH, with competition effects being found at both non-genetic (i.e., residual) and genetic levels, especially at site 2. By contrast, TH showed moderate to strong competition at the residual and genetic levels only at site 2, and NSTR revealed competition at the residual level at both sites and strong competition at the genetic level only at site 2. A detailed diagnosis of competition effects for each trait–site combination is given in the Annex (see Diagnosis of competition effects and Fig. 1)

Fig. 1

Isotropic empirical semivariograms of residuals after fitting a multi-environment individual-tree mixed model with fixed effect of overall mean, random check lot effects, random additive (breeding values) and dominance genetic effects, and random errors, for diameter at breast height (DBH), total tree height (TH), and normal score of stem straightness (NSTR) at age 10 for sites 1 and 2

.

The three traits assessed (DBH, TH, and NSTR) were analyzed independently using the following multi-environment standard and competition individual-tree mixed models:

• 1. Multi-environment standard individual-tree mixed model (MSM)

$$ {\displaystyle \begin{array}{l}\left[\begin{array}{c}{\boldsymbol{y}}_1\\ {}{\boldsymbol{y}}_2\end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{X}}_1& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{X}}_2\end{array}\right]\left[\begin{array}{c}{\boldsymbol{\beta}}_1\\ {}{\boldsymbol{\beta}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{b_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{b_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{b}}_1\\ {}{\boldsymbol{b}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{p_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{p_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{p}}_1\\ {}{\boldsymbol{p}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{s_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{s_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{s}}_1\\ {}{\boldsymbol{s}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{a_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{a_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{a}}_1\\ {}{\boldsymbol{a}}_2\end{array}\right]\\ {}\kern.75em +\left[\begin{array}{cc}{\boldsymbol{Z}}_{d_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{d_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{d}}_1\\ {}{\boldsymbol{d}}_2\end{array}\right]+\left[\begin{array}{c}{\boldsymbol{e}}_1\\ {}{\boldsymbol{e}}_2\end{array}\right]\end{array}} $$

(1)

where y₁ and y₂ are the vectors of individual tree observation for the same trait at sites 1 and 2, respectively. The matrices X₁ ⊕ X₂, ${\boldsymbol{Z}}_{b_1}\oplus {\boldsymbol{Z}}_{b_2}$, ${\boldsymbol{Z}}_{p_1}\oplus {\boldsymbol{Z}}_{p_2}$, ${\boldsymbol{Z}}_{s_1}\oplus {\boldsymbol{Z}}_{s_2}$, ${\boldsymbol{Z}}_{a_1}\oplus {\boldsymbol{Z}}_{a_2}$, and ${\boldsymbol{Z}}_{d_1}\oplus {\boldsymbol{Z}}_{d_2}$ relate the observations to the fixed effects of overall site means in $\left[{\boldsymbol{\beta}}_1^T|{\boldsymbol{\beta}}_2^T\right]$, the random effects of block in $\left[{\boldsymbol{b}}_1^T|{\boldsymbol{b}}_2^T\right]$, of plot in $\left[{\boldsymbol{p}}_1^T|{\boldsymbol{p}}_2^T\right]$, of check lot in $\left[{\boldsymbol{s}}_1^T|{\boldsymbol{s}}_2^T\right]$, of additive genetic (breeding value) effects in $\left[{\boldsymbol{a}}_1^T|{\boldsymbol{a}}_2^T\right]$, and of the dominance genetic effects in $\left[{\boldsymbol{d}}_1^T|{\boldsymbol{d}}_2^T\right]$. The vector $\left[{\boldsymbol{e}}_1^T|{\boldsymbol{e}}_2^T\right]$ is the residual vector. The symbols ⊕ and ^T indicate the direct sum of matrices and the transpose operation, respectively. The vector of block effects was assumed distributed as

$$ \left[\begin{array}{c}{\boldsymbol{b}}_1\\ {}{\boldsymbol{b}}_2\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\end{array}\right],\left[\begin{array}{cc}{\sigma}_{b_1}^2& \mathbf{0}\\ {}\mathbf{0}& {\sigma}_{b_2}^2\end{array}\right]\otimes \boldsymbol{I}\right) $$

where ${\sigma}_{b_1}^2$ and ${\sigma}_{b_2}^2$ are the block effect variances of sites 1 and 2, respectively, and I is the identity matrix. The vector of plot effects was assumed distributed as

$$ \left[\begin{array}{c}{\boldsymbol{p}}_1\\ {}{\boldsymbol{p}}_2\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\end{array}\right],\left[\begin{array}{cc}{\sigma}_{p_1}^2& \mathbf{0}\\ {}\mathbf{0}& {\sigma}_{p_2}^2\end{array}\right]\otimes \boldsymbol{I}\right) $$

where ${\sigma}_{p_1}^2$ and ${\sigma}_{p_2}^2$ are the plot effect variances of sites 1 and 2, respectively. The vector of check lot effects was assumed distributed as

$$ \left[\begin{array}{c}{\boldsymbol{s}}_1\\ {}{\boldsymbol{s}}_2\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\end{array}\right],\left[\begin{array}{cc}{\sigma}_{s_1}^2& {\sigma}_{s_{1,2}}\\ {}{\sigma}_{s_{2,1}}& {\sigma}_{s_2}^2\end{array}\right]\otimes \boldsymbol{I}\right) $$

where ${\sigma}_{s_1}^2$ and ${\sigma}_{s_2}^2$ are the check lot effect variances of sites 1 and 2, respectively, and ${\sigma}_{s_{1,2}}$ is the check lot covariance between sites 1 and 2. The vector of individual breeding values was assumed distributed as

$$ \left[\begin{array}{c}{\boldsymbol{a}}_1\\ {}{\boldsymbol{a}}_2\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\end{array}\right],\left[\begin{array}{cc}{\sigma}_{a_1}^2& {\sigma}_{a_{1,2}}\\ {}{\sigma}_{a_{2,1}}& {\sigma}_{a_2}^2\end{array}\right]\otimes \boldsymbol{A}\right) $$

where ${\sigma}_{a_1}^2$ and ${\sigma}_{a_2}^2$ are the additive genetic variances of sites 1 and 2, respectively, ${\sigma}_{a_{1,2}}$ is the additive genetic covariance between sites 1 and 2, and A is the average numerator relationship matrix (Henderson 1984). The vector of dominance genetic effects was assumed distributed as

$$ \left[\begin{array}{c}{\boldsymbol{d}}_1\\ {}{\boldsymbol{d}}_2\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\end{array}\right],\left[\begin{array}{cc}{\sigma}_{d_1}^2& {\sigma}_{d_{1,2}}\\ {}{\sigma}_{d_{2,1}}& {\sigma}_{d_2}^2\end{array}\right]\otimes \boldsymbol{D}\right) $$

where ${\sigma}_{d_1}^2$ and ${\sigma}_{d_2}^2$ are the dominance genetic variances of sites 1 and 2, respectively, ${\sigma}_{d_{1,2}}$ is the dominance genetic covariance between sites 1 and 2, and D is the dominance relationship matrix (Mrode 2014). Finally, the residual vector is distributed as

$$ \left[\begin{array}{c}{\boldsymbol{e}}_1\\ {}{\boldsymbol{e}}_2\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\end{array}\right],\left[\begin{array}{cc}{\sigma}_{e_1}^2& \mathbf{0}\\ {}\mathbf{0}& {\sigma}_{e_2}^2\end{array}\right]\otimes \boldsymbol{I}\right) $$

where ${\sigma}_{e_1}^2$ and ${\sigma}_{e_2}^2$ are the residual variances of sites 1 and 2, respectively.

• 2. Multi-environment competition individual-tree mixed model (MCM)

The MSM (Eq. (1) was extended to include the direct and competition additive genetic effects as well as environmental competition effects, as follows:

$$ \left[\begin{array}{c}{\boldsymbol{y}}_1\\ {}{\boldsymbol{y}}_2\end{array}\right]=\left[\begin{array}{cc}{\boldsymbol{X}}_1& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{X}}_2\end{array}\right]\left[\begin{array}{c}{\boldsymbol{\beta}}_1\\ {}{\boldsymbol{\beta}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{b_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{b_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{b}}_1\\ {}{\boldsymbol{b}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{p_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{p_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{p}}_1\\ {}{\boldsymbol{p}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{s_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{s_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{s}}_1\\ {}{\boldsymbol{s}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{a_{d_1}}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{a_{d_2}}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{a}}_{d_1}\\ {}{\boldsymbol{a}}_{d_2}\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{a_{c_1}}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{a_{c_2}}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{a}}_{c_1}\\ {}{\boldsymbol{a}}_{c_2}\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{d_1}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{d_2}\end{array}\right]\left[\begin{array}{c}{\boldsymbol{d}}_1\\ {}{\boldsymbol{d}}_2\end{array}\right]+\left[\begin{array}{cc}{\boldsymbol{Z}}_{p_{c_1}}& \mathbf{0}\\ {}\mathbf{0}& {\boldsymbol{Z}}_{p_{c_2}}\end{array}\right]\left[\begin{array}{c}{{\boldsymbol{p}}_{\boldsymbol{c}}}_1\\ {}{{\boldsymbol{p}}_{\boldsymbol{c}}}_2\end{array}\right]+\left[\begin{array}{c}{\boldsymbol{e}}_1\\ {}{\boldsymbol{e}}_2\end{array}\right] $$

(2)

where the matrices X₁ ⊕ X₂, ${\boldsymbol{Z}}_{b_1}\oplus {\boldsymbol{Z}}_{b_2}$, ${\boldsymbol{Z}}_{p_1}\oplus {\boldsymbol{Z}}_{p_2}$, ${\boldsymbol{Z}}_{s_1}\oplus {\boldsymbol{Z}}_{s_2}$, ${\boldsymbol{Z}}_{d_1}\oplus {\boldsymbol{Z}}_{d_2}$, and the vectors $\left[{\boldsymbol{y}}_1^T|{\boldsymbol{y}}_2^T\right]$, $\left[{\boldsymbol{\beta}}_1^T|{\boldsymbol{\beta}}_2^T\right]$, $\left[{\boldsymbol{b}}_1^T|{\boldsymbol{b}}_2^T\right]$, $\left[{\boldsymbol{p}}_1^T|{\boldsymbol{p}}_2^T\right]$, $\left[{\boldsymbol{s}}_1^T|{\boldsymbol{s}}_2^T\right]$, $\left[{\boldsymbol{d}}_1^T|{\boldsymbol{d}}_2^T\right]$, and $\left[{\boldsymbol{e}}_1^T|{\boldsymbol{e}}_2^T\right]$ are specified as described above. In Eq. (2), the matrices ${\boldsymbol{Z}}_{a_{d_1}}\oplus {\boldsymbol{Z}}_{a_{d_2}}$, and ${\boldsymbol{Z}}_{a_{c_1}}\oplus {\boldsymbol{Z}}_{a_{c_2}}$ related the observations to the individual additive direct genetic breeding values in $\left[{\boldsymbol{a}}_{d_1}^T|{\boldsymbol{a}}_{d_2}^T\right]$, and the additive competition genetic breeding values of the neighboring trees in $\left[{\boldsymbol{a}}_{c_1}^T|{\boldsymbol{a}}_{c_2}^T\right]$. The i^th row of matrices ${\boldsymbol{Z}}_{a_{c_1}}$ and ${\boldsymbol{Z}}_{a_{c_2}}$ has all elements equal to zero except in the position corresponding to the m_i neighbors of the tree i, with positive values interpreted as the intensity of competition (IC, f_ij, j = 1, …, m_i) that each neighbor exerts over the phenotype of the i^th tree. A detailed explanation of these IC values can be found in Cappa and Cantet (2008). The stacked vector of breeding values for the MCM (Eq. 2) is distributed as

$$ \left[\begin{array}{c}\begin{array}{c}{\boldsymbol{a}}_{d_1}\\ {}{\boldsymbol{a}}_{d_2}\\ {}{\boldsymbol{a}}_{c_1}\end{array}\\ {}{\boldsymbol{a}}_{c_2}\end{array}\right]\sim N\left(\left[\begin{array}{c}\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\\ {}\mathbf{0}\end{array}\\ {}\mathbf{0}\end{array}\right],\kern0.5em \left[\begin{array}{cc}\begin{array}{cc}{\sigma}_{a_{d_1}}^2& {\sigma}_{a_{d_{1,2}}}\\ {}{\sigma}_{a_{d_{2,1}}}& {\sigma}_{a_{d_2}}^2\end{array}& \begin{array}{cc}{\sigma}_{a_{d_1{c}_1}}& {\sigma}_{a_{d_1{c}_2}}\\ {}{\sigma}_{a_{d_2{c}_1}}& {\sigma}_{a_{d_2{c}_2}}\end{array}\\ {}\begin{array}{cc}{\sigma}_{a_{c_1{d}_1}}& {\sigma}_{a_{c_1{d}_2}}\\ {}{\sigma}_{a_{c_2{d}_1}}& {\sigma}_{a_{c_2{d}_2}}\end{array}& \begin{array}{cc}{\sigma}_{a_{c_1}}^2& {\sigma}_{a_{c_{1,2}}}\\ {}{\sigma}_{a_{c_{2,1}}}& {\sigma}_{a_{c_2}}^2\end{array}\end{array}\right]\otimes \boldsymbol{A}\right) $$

where ${\sigma}_{a_{d_1}}^2$ and ${\sigma}_{a_{d_2}}^2$ are the additive direct genetic variance for sites 1 and 2, respectively; ${\sigma}_{a_{c_1}}^2$ and ${\sigma}_{a_{c_2}}^2$ are the additive competition genetic variance for sites 1 and 2, respectively; ${\sigma}_{a_{d_{1,2}}}$ and ${\sigma}_{a_{d_{2,1}}}$ are the covariance between direct breeding values of sites 1 and 2; ${\sigma}_{a_{c_{1,2}}}$ and ${\sigma}_{a_{c_{2,1}}}$ are the covariance between the competition breeding values of sites 1 and 2; ${\sigma}_{a_{d_1{c}_1}}$ and ${\sigma}_{a_{d_2{c}_2}}$ are the covariance between the additive direct and competition breeding values within sites 1 and 2, respectively; and ${\sigma}_{a_{d_1{c}_2}}$ and ${\sigma}_{a_{d_2{c}_1}}$ are the covariance between the direct breeding values of site 1 and competition breeding values of site 2 and between direct breeding values of site 2 and competition breeding values of site 1, respectively. The vector of dominance genetic effects was specified as described above. Finally, Eq. (2) included the matrices ${\boldsymbol{Z}}_{p_{c_1}}\oplus {\boldsymbol{Z}}_{p_{c_2}}$ related to the environmental competition effects (or permanent environmental competition effects; Cappa and Cantet 2008) in $\left[{\boldsymbol{p}}_{c_1}^T|{\boldsymbol{p}}_{c_2}^T\right]$ and is distributed as

$$ \left[\begin{array}{c}{\boldsymbol{p}}_{c_1}\\ {}{\boldsymbol{p}}_{c_2}\end{array}\right]\sim N\left(\left[\begin{array}{c}\mathbf{0}\\ {}\mathbf{0}\end{array}\right],\left[\begin{array}{cc}{\sigma}_{p_{c_1}}^2& \mathbf{0}\\ {}\mathbf{0}& {\sigma}_{p_{c_2}}^2\end{array}\right]\kern0.5em \otimes \boldsymbol{I}\right) $$

where ${\sigma}_{p_{c_1}}^2$ and ${\sigma}_{p_{c_2}}^2$ are environmental competition variances of sites 1 and 2, respectively.

2.3 Parameter estimation and model comparison

Restricted maximum likelihood (REML, Patterson and Thompson 1971) was used to estimate variance components for models (1) and (2), which were obtained using the Expectation Maximization algorithm in R with the R (R Core Team 2015) with the function remlf90 from the package breedR (Muñoz and Sanchez 2019), followed by one round with Average Information (AI) algorithm to compute the standard errors of the variance components and heritabilities (as in Chateigner et al. 2020).

The total heritable (direct and competition) effect of an individual’s genes on trait values in a population, or total tree breeding value of tree i (TBV_i), from the MCM model (2) was calculated following Costa e Silva et al. (2013; Eq. 5), i.e., ${\mathrm{TBV}}_i={\boldsymbol{a}}_{d_i}+\left({\overline{n}}_R{\overline{f}}_{ijR}+{\overline{n}}_C{\overline{f}}_{ijC}+{\overline{n}}_D{\overline{f}}_{ijD}\right)\ {\boldsymbol{a}}_{c_i}$, where ${\overline{n}}_R$, ${\overline{n}}_C$, and ${\overline{n}}_D$ denote the average numbers of their neighbors, and ${\overline{f}}_{ijR}$, ${\overline{f}}_{ijC}$, and ${\overline{f}}_{ijD}$ are the average IC values in the row (R), column (C), and diagonal (D) directions, respectively. The values of the expression $\left({\overline{n}}_R{\overline{f}}_{ijR}+{\overline{n}}_C{\overline{f}}_{ijC}+{\overline{n}}_D{\overline{f}}_{ijD}\right)$ were 1.99 for site 1 and 1.79 for site 2. The total heritable variances of traits DBH, TH, and NSTR were also estimated for each site following accordingly Costa e Silva et al. (2013); Eq. 6, i.e., ${\hat{\sigma}}_{TBV}^2={\hat{\sigma}}_{a_d}^2+2\left({\overline{n}}_R{\overline{f}}_{ijR}+{\overline{n}}_C{\overline{f}}_{ijC}+{\overline{n}}_D{\overline{f}}_{ijD}\right){\hat{\sigma}}_{a_{dc}}+{\left({\overline{n}}_R{\overline{f}}_{ijR}+{\overline{n}}_C{\overline{f}}_{ijC}+{\overline{n}}_D{\overline{f}}_{ijD}\right)}^2{\hat{\sigma}}_{a_c}^2$. Finally, the across-site additive genetic correlation of the TBV (i.e., additive direct plus competition breeding values) was estimated using the eight immediate neighbors of a focal tree following Costa e Silva et al. (2017; Supplementary Note S7, Eq. S_9 and S_10).

An important aspect for breeding selection is to examine the impact on the ranking of candidates across a multi-environmental setting in the presence of competition effects. In that sense, Spearman rank correlations were calculated for the top 5% (111 trees) and for each trait to assess the extent by which the ranking of predicted breeding values for individuals (offspring) from the MSM differed from the additive direct breeding values predicted from the MCM model. In addition, we plotted the ranking of 111 trees from the two models and the three studied traits.

The single-site individual narrow- and broad-sense heritability from the MSM (${\hat{h}}_s^2$ and ${\hat{H}}_s^2$, respectively) and from the MCM (${\hat{h}}_c^2$ and ${\hat{H}}_c^2$, respectively) were estimated as

$$ {\displaystyle \begin{array}{c}{\hat{h}}_s^2=\frac{{\hat{\sigma}}_a^2}{{\hat{\sigma}}_a^2+{\hat{\sigma}}_d^2+{\hat{\sigma}}_e^2};\kern0.5em {\hat{H}}_s^2=\frac{{\hat{\sigma}}_a^2+{\hat{\sigma}}_d^2}{{\hat{\sigma}}_a^2+{\hat{\sigma}}_d^2+{\hat{\sigma}}_e^2}\\ {}{\hat{h}}_c^2=\frac{{\hat{\sigma}}_{a_d}^2+{\hat{\sigma}}_{a_c}^2}{{\hat{\sigma}}_{a_d}^2+{\hat{\sigma}}_{a_c}^2\kern0.5em +{\hat{\sigma}}_d^2+{\hat{\sigma}}_{p_c}^2+{\hat{\sigma}}_e^2};\kern0.5em {\hat{H}}_c^2=\frac{{\hat{\sigma}}_{a_d}^2+{\hat{\sigma}}_{a_c}^2+{\hat{\sigma}}_d^2}{{\hat{\sigma}}_{a_d}^2+{\hat{\sigma}}_{a_c}^2\kern0.5em +{\hat{\sigma}}_d^2+{\hat{\sigma}}_{p_c}^2+{\hat{\sigma}}_e^2}\end{array}} $$

where ${\hat{\sigma}}_a^2$ is the estimated additive genetic variance, ${\hat{\sigma}}_{a_d}^2$ is the estimated additive direct genetic variance, ${\hat{\sigma}}_{a_c}^2$ is the estimated additive competition genetic variance, ${\hat{\sigma}}_d^2$ is the estimated dominance genetic variance, ${\hat{\sigma}}_{p_c}^2$ is the estimated environmental competition variance, and ${\hat{\sigma}}_e^2$ is the estimated residual genetic variance. Note that for the MCM model, the additive variance was replaced with the estimated total additive genetic variance, ${\hat{\sigma}}_{a_d}^2+{\hat{\sigma}}_{a_c}^2$ (Cappa et al. 2017). In addition, for MCM the denominator of the above expression also included the estimated variance of the environmental competition effects, ${\hat{\sigma}}_{p_c}^2$.

The across-site additive genetic correlation (${\hat{r}}_{a_{1,2}}$) and across-site dominance genetic correlations (${\hat{r}}_{d_{1,2}}$) from models (1) and (2) were calculated as

$$ {\hat{r}}_{a_{1,2}}=\frac{{\hat{\sigma}}_{a_{1,2}}}{\sqrt{{\hat{\sigma}}_{a_1}^2\times {\hat{\sigma}}_{a_2}^2}};\kern0.75em {\hat{r}}_{d_{1,2}}=\frac{{\hat{\sigma}}_{d_{1,2}}}{\sqrt{{\hat{\sigma}}_{d_1}^2\times {\hat{\sigma}}_{d_2}^2}} $$

In addition, the across-site additive direct–direct (${\hat{r}}_{a_{d_{1,2}}}$) and competition–competition (${\hat{r}}_{a_{c_{1,2}}}$) genetic correlations, and the across-site additive direct–competition (${\hat{r}}_{a_{d_1{c}_2}}$ and ${\hat{r}}_{a_{d_2{c}_1}}$) genetic correlations from model (2) were calculated as

$$ {\displaystyle \begin{array}{c}{\hat{r}}_{a_{d_{1,2}}}=\frac{{\hat{\sigma}}_{a_{d_{1,2}}}}{\sqrt{{\hat{\sigma}}_{a_{d_1}}^2\times {\hat{\sigma}}_{a_{d_2}}^2}};{\hat{r}}_{a_{c_{1,2}}}=\frac{{\hat{\sigma}}_{a_{c_{1,2}}}}{\sqrt{{\hat{\sigma}}_{a_{c_1}}^2\times {\hat{\sigma}}_{a_{c_2}}^2}}\\ {}{\hat{r}}_{a_{d_1{c}_2}}=\frac{{\hat{\sigma}}_{a_{d_1{c}_2}}\ }{\sqrt{{\hat{\sigma}}_{a_{d_1}}^2\times {\hat{\sigma}}_{a_{c_2}}^2}};{\hat{r}}_{a_{d_2{c}_1}}=\frac{{\hat{\sigma}}_{a_{d_2{c}_1}}}{\sqrt{{\hat{\sigma}}_{a_{d_2}}^2\times {\hat{\sigma}}_{a_{c_1}}^2}}\end{array}} $$

The Akaike Information Criterion (AIC, Akaike 1974) was used to compare models (1) and (2) and assess which one yielded the best fit to the data. A smaller AIC value indicates a better trade-off between goodness-of-fit and parsimony. In addition, to determine the statistical significance, we applied the one-tailed likelihood-ratio tests with nine degrees of freedom that represented the number of additional parameters in model MCM versus model MSM (LRT; Stram and Lee 1994). The significances from zero of the across-site genetic correlations for additive and dominance effects, and genetic additive correlations for direct–direct, competition–competition, and direct–competition effects, were evaluated by the two-tailed likelihood LRT with one degree of freedom.

3 Results

3.1 Model comparison

When comparing the AIC values from the multi-environment standard individual-tree mixed model without (MSM) and with competition being included (MCM), we found an improvement (better fit) for the MCM analysis for DBH and TH (Table 2). By contrast, for NSTR the simpler MSM model had a better fit. MCM showed best fit for traits with strong to moderate genetic inter-tree competition effects (Table 3), DBH and TH, where the average of ${\hat{r}}_{a_{d.c.}}$ across both sites was − 0.43, and − 0.19, respectively. Contrarily, the trait for which MSM had the best fit, straightness (NSTR), showed to be less sensitive to genetic competition effects with average of ${\hat{r}}_{a_{d.c.}}$ across both sites being equal to 0.08. Since the AIC statistic is composed of a measure of total goodness of fit and a penalty on the complexity of the model, the MSM had a smaller AIC than MCM for NSTR. When considering only the measure of the goodness of fit (i.e., LRT), the MCM improved the fit significantly for traits DBH (p = 0.0089) and TH (p = 0.0000). However, NSTR showed that there were no differences between the two models (logL = − 3992), suggesting that the MCM did not differ significantly from the MSM.

Table 2 Akaike information criterion (AIC), estimated additive genetic variance (${\hat{\sigma}}_a^2$), additive direct genetic variance (${\hat{\sigma}}_{a_d}^2$), additive competition genetic variance (${\hat{\sigma}}_{a_c}^2$), dominance genetic variance (${\hat{\sigma}}_d^2$), check lot effect variance (${\hat{\sigma}}_s^2$), block effect variance (${\hat{\sigma}}_b^2$), plot effect variance (${\hat{\sigma}}_p^2$), environmental competition effects variance (${\hat{\sigma}}_{p_c}^2$), residual variance (${\hat{\sigma}}_e^2$), and single-site individual narrow- (${\hat{h}}_s^2$ and${\hat{h}}_c^2$) and broad-sense (${\hat{H}}_s^2$ and ${\hat{H}}_c^2$) heritability for diameter at breast height (DBH), total tree height (TH), and normal score of stem straightness (NSTR) from the multi-environment standard (MSM) and competition (MCM) individual-tree mixed models of sites 1 and 2 (the lowest AIC are highlighted in bold)

Table 3 Across-site genetic correlations for additive (${\hat{r}}_{a_{12}}$), dominance (${\hat{r}}_{d_{12}}$), additive direct–direct (${\hat{r}}_{a_{d_{12}}}$), additive competition–competition (${\hat{r}}_{a_{c_{12}}}$), and additive direct–competition (${\hat{r}}_{a_{d1c2}}$and ${\hat{r}}_{a_{d2c1}}$) effects for diameter at breast height (DBH), total tree height (TH), and normal score of stem straightness (NSTR) from the multi-environment standard (MSM) and competition (MCM) individual-tree mixed models

3.2 Variance components and heritabilities from the MSM and MCM

In general, ignoring competition’s genetic and environmental effects by fitting a MSM resulted in lower additive genetic variances, compared with the results of MCM. These differences were more noticeable for DBH, especially at site 2, where the additive variance was 45.3% larger in MCM than in MSM (Table 2). The exceptions were TH at site 2, where the additive genetic variance was 1.9% smaller with MCM than MSM, and for NSTR at site 2, where additive genetic variances were the same with MSM and MCM. The standard errors of additive genetic variances were in general lower than half of the magnitude of parameters estimated from both models, having the MCM in general lower standard errors than MSM. Regarding the residual variance estimates (${\hat{\sigma}}_e^2\Big)$, the MCM model resulted in lower ${\hat{\sigma}}_e^2$ at both sites for DBH (an average drop of 9.9%) and for TH (average drop of 18.1%) compared to MSM. For NSTR, the ${\hat{\sigma}}_e^2$ decreased slightly with MCM compared with MSM at both sites (average of 1.6%). The standard errors were always less than a quarter of the residual variance estimates for the two models.

Regarding the dominance genetic variance (${\hat{\sigma}}_d^2$), the MCM estimated lower values than MSM (from 6.7 to 14.3%), except for DBH and TH at site 1, where the ${\hat{\sigma}}_d^2$ increased slightly (< 4.5%) from MSM to MCM (Table 2), but differences were minor in all cases. In general, the standard errors were about half of the estimated dominance genetic variance, except for TH at site 1, where they were very high for MSM and MCM. This result, along with the lower values of additive genetic variances from MSM than from MCM mentioned above, caused the ratio of ${\hat{\sigma}}_a^2/{\hat{\sigma}}_d^2$to be always smaller for MSM (average = 1.61) than (${\hat{\sigma}}_{a_d}^2/{\hat{\sigma}}_d^2$) for MCM (average = 1.95). On the other hand, for both models, the ${\hat{\sigma}}_d^2$ was lower than ${\hat{\sigma}}_{a_d}^2$ for all traits and sites (0.27 ≥ ${\hat{\sigma}}_d^2/{\hat{\sigma}}_{a_d}^2$ ≤ 0.88), except for DBH at site 1, where the ratio of ${\hat{\sigma}}_d^2/{\hat{\sigma}}_{a_d}^2$ was equal to 1.52 and 1.42 for MSM and MCM, respectively.

In addition to the competition’s genetic effects, the proposed multi-environment competition mixed model also considered a permanent competition’s environmental effect (${\boldsymbol{p}}_{c_1}$ and ${\boldsymbol{p}}_{c_2}\Big)$, fitted as an additional random term to model (2). For growth traits, site 1 showed a greater variance of p_c than site 2, while for NSTR the variance of p_c was similar at both sites (Table 2). However, in all cases, the standard errors of permanent competition environmental effects were very high.

In general, the MCM model resulted in greater individual narrow- and broad-sense heritability values than the MSM model for the three traits evaluated (Table 2). Narrow-sense heritability values ranged from 0.15 to 0.28 for the MSM model, and from 0.14 to 0.40 for the MCM model. On the other hand, broad-sense heritability values ranged from 0.21 to 0.47 for the MSM model and from 0.21 to 0.51 for MCM model. The greater differences between the heritabilities from the two models were observed at site 2 for DBH trait, where narrow-sense heritability was 0.28 for the MSM model and 0.40 for the MCM model, and broad-sense heritability was 0.40 for the MSM model and 0.51 for the MCM model. The standard errors of heritability were in general lower than half of the heritability estimates for both models, except for trait TH at site 1, where they were about half.

In our study, the competition additive genetic variance (${\hat{\sigma}}_{a_c}^2$) represented on average less than 10.6% of the corresponding direct additive genetic variance (${\hat{\sigma}}_{a_d}^2$) for growth traits, being 0% for NSTR. In all cases, however, the standard errors for ${\hat{\sigma}}_{a_c}^2$ were very high. The effect of combining together ${\hat{\sigma}}_{a_c}^2$ and ${\hat{\sigma}}_{a_d}^2$was a reduction in the total heritable variance (${\hat{\sigma}}_{TBV}^2$) for DBH, with reductions at site 1 of 13.0% (i.e., ${\hat{\sigma}}_{a_d}^2$ = 4.91 vs. ${\hat{\sigma}}_{TBV}^2$ = 4.28) and 48.0% at site 2 (i.e., ${\hat{\sigma}}_{a_d}^2$ = 10.52 vs. ${\hat{\sigma}}_{TBV}^2$ = 5.48). The greater reduction in ${\hat{\sigma}}_{TBV}^2$ at site 2 can be explained by the negative sign and high absolute value of covariance between direct additive and competition genetic effects (${\hat{\sigma}}_{a_c}^2$ = − 1.67), relative to the small value of additive competition genetic variance (${\hat{\sigma}}_{a_c}^2$ = 0.29), which gives more weight to the second than to the third term in the equation of total heritable variance (${\hat{\sigma}}_{TBV}^2$, see Eq. 6 in “Materials and methods” section). By contrast, for TH our result revealed that ${\hat{\sigma}}_{TBV}^2$ increased at site 1 by 50% (i.e., ${\hat{\sigma}}_{a_d}^2$ = 1.15 vs. ${\hat{\sigma}}_{TBV}^2$ = 1.72), due to a greater contribution of ${\hat{\sigma}}_{a_c}^2$ (0.20) to the equation of total heritable variance than the negative covariance (${\hat{\sigma}}_{a_{dc}}$ = − 0.06). Nevertheless, at site 2, ${\hat{\sigma}}_{TBV}^2$ showed a slight decrease (1%; i.e., ${\hat{\sigma}}_{a_d}^2$ = 1.05 vs. ${\hat{\sigma}}_{TBV}^2$ = 1.06) due to a comparable contribution of the ${\hat{\sigma}}_{a_c}^2$ (0.04) and the negative ${\hat{\sigma}}_{a_{dc}}$ (− 0.03) to the ${\hat{\sigma}}_{TBV}^2$. Finally, the competition effects showed no consistent trend on ${\hat{\sigma}}_{TBV}^2$ of NSTR, increasing by 20% at site 1 (i.e., ${\hat{\sigma}}_{a_d}^2$ = 0.18 vs. ${\hat{\sigma}}_{TBV}^2$ = 0.22) and decreasing the ${\hat{\sigma}}_{TBV}^2$ by 19% at site 2 (i.e., ${\hat{\sigma}}_{a_d}^2$ = 0.12 vs. ${\hat{\sigma}}_{TBV}^2$ = 0.10).

3.3 Across-site genetic correlations from the MSM and MCM

Estimates of across-site additive genetic correlations from MSM $\left({\hat{r}}_{a_{12}}\right)$ and additive direct genetic correlations from MCM $\left({\hat{r}}_{a_{d12}}\right)$ were strong for TH and NSTR (> 0.89 for both models) and moderate for DBH (0.77 for MSM and 0.79 for MCM), and all significantly greater than zero (Table 3). When we compared these estimates across models, for growth traits, slightly higher ${\hat{r}}_{a_{d12}}$ values were obtained from the MCM model compared with the MSM model (a 2.6% drop), while for NSTR ${\hat{r}}_{a_{d12}}$ the MCM yielded lower correlations than from MSM (a drop of 5.1%). In contrast to ${\hat{r}}_{a_{d12}}$, estimates of the across-site dominance genetic correlation $\left({\hat{r}}_{d_{12}}\right)$ for growth traits were lower (< 0.55 for both models) than for the NSTR trait (> 0.88).

The estimates across-site correlations between additive competition–competition genetic effects were small and, in general, not significantly different from zero for the three traits studied (${\hat{r}}_{a_{c12}}$ < 0.24, Table 3). The across-site correlations between additive direct–competition genetic effects $\Big({\hat{r}}_{a_{d1c2}}$ and ${\hat{r}}_{a_{d2c1}}\Big)$ were small and, in general, negative for all the studied traits, except for DBH where the ${\hat{r}}_{a_{d1c2}}$ was moderate (− 0.66) and statistically significant (Table 3).

We also estimated a single across-site additive direct and competition genetic correlation from MCM by computing the correlation between the TBV for sites 1 and 2 (${\hat{r}}_{TBV_{12}}\Big)$. The ${\hat{r}}_{TBV_{12}}$values were positive and moderate to high for the three studied traits (i.e., ${\hat{r}}_{TBV_{12}}$ = 0.83 for DBH, ${\hat{r}}_{TBV_{12}}$ = 0.57 for TH, and ${\hat{r}}_{TBV_{12}}$ = 1.00 for NSTR). For DBH and NSTR, ${\hat{r}}_{TBV_{12}}$ values were close to the ${\hat{r}}_{a_{d12}}$ equivalents (Table 3), while for TH the differences between the two correlations were much larger. The TBV correlation depends on how strong the across-site covariance between direct–direct, direct–competition, and competition–competition genetic effects are, plus on the magnitude of the intensity competition factor and total heritable variance $\left({\hat{\sigma}}_{TBV}^2\right)$ between sites 1 and 2. In this sense, the low values of ${\hat{r}}_{TBV_{12}}$ observed for TH was mainly associated to the negative direct–competition covariance (${\sigma}_{a_{d_2{c}_1}}$ = − 0.11) which reduced the covariance of TBV, plus the increment of ${\hat{\sigma}}_{TBV}^2$ (1.72) with respect to ${\hat{\sigma}}_{a_d}^2$ (1.15) at site 1.

Finally, the Spearman rank correlations between offspring predicted additive breeding values (BVs) from the MSM model and additive direct breeding values predicted from the MCM model were positive and very high (≥ 0.93) for TH and NSTR (Fig. 2),

Fig. 2

Ranking for the top 111 offspring trees (5 %) based on additive breeding values from the multi-environment standard individual-tree mixed model (MSM) and additive direct breeding values from the multi-environment competition individual-tree mixed models (MCM), for diameter at breast height (DBH), total tree height (TH), and normal scores of stem straightness (NSTR). The 10 most competitive trees (i.e., those with highest—negative—breeding values) for DBH are highlighted in blue color

indicating that similar selection for both traits will result from both models. However, for the trait most affected by competition, DBH, the Spearman correlation was only moderate (0.71), suggesting that different selection decisions would be made following either model.

4 Discussion

This research presents an extension of the multi-environment standard individual-tree mixed model to account for competition in multi-environment forest genetic trials. This model provides a formal genetic and environmental structure to study the genetic additive direct and competition effects by environment interactions (G × E) by estimating across-site direct and competition genetic correlations. In the present research, the (co)variance components and their functions (including across-site additive and dominance genetic correlations), and the variances (and across-site correlations) of total breeding values were estimated using the proposed model for growth and stem straightness traits in two full-sib genetic trials of Pinus elliottii var. elliottii × Pinus caribaea var. hondurensis F₁ assessed at age 10.

4.1 Variance components and heritabilities from the MSM and MCM

In general, the analysis reported here showed that the MCM increased the estimated additive variance and reduced the residual variance compared with results from the MSM. These results are in agreement with those reported by Cappa and Cantet (2008) and Cappa et al. (2015) for DBH in Pinus taeda. For instance, Cappa et al. (2015) observed that the inclusion of competition effects in the mixed model increased the additive variance by 40% and reduced the residual variance by 34%. The same tendency was reported by Costa e Silva and Kerr (2013), who indicated that with a medium relatedness level (i.e., full-sib mating design and multi-tree row plots), the single-trait models that ignored competition either at genetic or residual levels tended to bias downwardly the additive variance and upwardly the residual variance. Cappa et al. (2016) observed the same trend on additive and residual variances for DBH and TH in Douglas fir when comparing a single-environment competition-spatial mixed model to a standard mixed model. As pointed out by Cappa and Cantet (2008), the differences in the size of the additive and residual variances estimated based on the standard and competition models depend on the sign and the magnitude of the within-site additive direct–competition correlation (${\hat{r}}_{a_{d.c.}}\Big),$compared with the magnitude of the within-site additive competition genetic variance $\left({\hat{\sigma}}_{a_c}^2\right)$. In our study, for DBH, the values of within-site correlation from MCM were ${\hat{r}}_{a_{dc}}$ = − 0.36 (and ${\hat{\sigma}}_{a_c}^2$ = 0.13) and − 0.95 (0.29) at sites 1 and 2, respectively, where the ratio of ${\hat{r}}_{a_{dc}}$/${\hat{\sigma}}_{a_c}^2$ was − 2.8 at site 1 and − 3.3 at site 2, generating a greater increase in additive variance and reduction in residual variance at site 2 than at site 1.

On the other hand, the greater proportion of additive variance compared with dominance variance observed in the present study is in agreement with the results of univariate standard analyses using the same hybrid reported in the present study (PEE × PCH; Belaber et al. 2018) and other Pinus hybrids (Dieters et al. 1997; Gwaze et al. 2000; Dungey et al. 2000; Brawner et al. 2005; Mutete et al. 2015). For instance, Dieters et al. (1997) reported a ratio between additive variance and dominance variance for DBH between 1.0 and 1.9 in hybrids between Pinus oocarpa Schiede and Pinus tecunumanii (Schw.) Eguiluz and Perry at age 5. The high proportion of ${\hat{\sigma}}_{a_d}^2$ compared to ${\hat{\sigma}}_d^2$ in this population of F₁ pine hybrid suggests that breeding strategies that maximize the use of additive genetic variance could be efficient. However, the ${\hat{\sigma}}_d^2$ was not negligible, representing more than 50% of the ${\hat{\sigma}}_{a_d}^2$ in the 66.7% of studied trait-site combinations. This would open up opportunities for exploiting the best specific combining abilities to deploy a family forestry via vegetative multiplication.

We have not found reports of competition models that include dominance effects. However, in a preliminary analysis to the present study, we evaluated the interplay between dominance and competition by comparing the variance components from the proposed MCM model with those obtained from the same model without dominance effects. Leaving dominance out of the MCM model mainly increased plot and residual variances but did not change the additive competition genetic variance (results not shown), indicating that there is no interaction between dominance and competition genetics effects.

Regarding the permanent competition’s environmental effect, Cappa and Cantet (2008) reported that leaving these environmental effects out of the competition model may bias upwardly the ${\hat{\sigma}}_{a_d}^2$ and downwardly the ${\hat{\sigma}}_e^2$. However, a preliminary analysis with present data that used a single-environment competition model showed that the omission of competition environmental term caused only a slight increment in the ${\hat{\sigma}}_{a_d}^2$ (mean = 5.2%) and a reduction in the ${\hat{\sigma}}_e^2$ (mean = 4.5%) across all the trait–site combinations, except for TH at site 1, where the ${\hat{\sigma}}_{a_d}^2$ decreased (14.9%). These differences could be related to lower survival of the two trails analyzed in the present study (mean = 82.0%) than that analyzed by Cappa and Cantet (2008) (95.7%). Therefore, as concluded by Cappa and Cantet (2008), the presence of dead trees promotes a less stressful environment for the surrounding plants, which may bias the prediction of direct breeding values, if environmental competition effects are not accounted for in the model.

Our results showed that ignoring a strong inter-tree competition may downwardly the estimation of narrow- and broad-sense heritabilities. A similar tendency was observed by Cappa et al. (2017) in Eucalyptus grandis (Hill ex Maiden) for the trait DBH, and by Hernández et al. (2019) in Corymbia citriodora subsp. variegata (F.Muell.) for DBH, TH, and wood density traits, when they compared narrow-sense individual heritability estimated from single-site single-trait individual-tree standard and competition mixed models. In short, we observed that the MCM model displayed a consistent increase in the estimated additive genetic variance and a reduction in the residual variance, with a consequent increase of the heritability with respect of MSM model. On the other hand, for the same traits and trials evaluated in the present study, but at age 7, Belaber et al. (2018) observed, in general, lower narrow- and broad-sense heritability values from a standard individual-tree mixed model than the values reported here from the MSM model, with the exception of the trait NSTR. We could suggest that the MCM model reveals the highest heritability values whenever the trait shows strong to moderate genetic inter-tree competition effects, which is often the case in older trials.

The heritable component is not only the result of the focal individual’s own genes but also from the effects of genes carried by its neighboring conspecifics (Bijma 2011). In other words, each individual’s breeding value include the gene’s direct additive effects but also the genetic component of the susceptibility to neighboring conspecifics (viz. a heritable “social” effect). Even though the indirect genetic variance is often considerably smaller than the direct genetic variance, such tree genetic interactions may substantially alter a population’s ability to respond to selection (Costa e Silva et al. 2013). In our study, the inter-tree genetic competition decreased the ${\hat{\sigma}}_{TBV}^2$ in respect to the ${\hat{\sigma}}_{a_d}^2$ of DBH at site 1 and 2. The drop in heritable variance due to competition effects is in agreement with single-site analyses in Eucalyptus globulus by Costa e Silva et al. (2013), who found that inter-tree genetic competition decreased the ${\hat{\sigma}}_{TBV}^2$ for DBH by 75% and 85% at ages 2 and 4, respectively. However, unlike in Costa e Silva et al. (2013), our results found that the ${\hat{\sigma}}_{a_d}^2$ contribute more to the ${\hat{\sigma}}_{TBV}^2$ of DBH than ${\hat{\sigma}}_{a_c}^2$. Nevertheless, this reduction in ${\hat{\sigma}}_{TBV}^2$ observed in the current study for DBH, due to inter-tree genetic competition, limits the potential for response to selection. In other words, the genotype growing rapidly in DBH will have (on average) a substantial negative heritable effect on the DBH growth from its neighbors’ competition. On the other hand, for TH, the increase and slight decrease in ${\hat{\sigma}}_{TBV}^2$ at sites 1 and 2, respectively, means that faster TH genotypes will have (on average) favorable (site 1) or null effects (site 2) on the TH of its neighbors. Finally, the increase in ${\hat{\sigma}}_{TBV}^2$ at site 1 was associated with the small positive contribution of ${\hat{\sigma}}_{a_c}^2$ (0.003) plus the positive additive-competition covariance (${\hat{\sigma}}_{a_{dc}}$ = 0.01). By contrast, the negative contribution of competition genetic effect to the ${\hat{\sigma}}_{TBV}^2$ at site 2 was due to a slightly negative additive-competition covariance (${\hat{\sigma}}_{a_{dc}}$ = − 0.01) compared with the small and positive value of ${\hat{\sigma}}_{a_c}^2$ (0.002) (Table 2). We have not found reports on the effects of inter-tree genetic competition for stem straightness, which added to our uncertain results, suggesting the need for new studies to investigate this issue.

4.2 Across-site genetic correlations from the MSM and MCM

An evaluation of the importance and consequences of genotype by environment interaction (G × E) in a breeding program is essential for decision-making about testing, deployment, and selection strategies (Zas et al. 2004). Previous studies on forest trees highlighted the potential of multi-environment individual-tree mixed models to investigate the genetic basis of G × E (e.g., Dungey et al. 2000; Costa e Silva et al. 2005; Cappa et al. 2012). However, despite the relevance of this issue, we have not found any reports on forest tree improvement that account for genetic competition in the estimation of across-site genetic correlation. Therefore, this is the first study to our knowledge reporting estimates of across-site genetic correlations from a MCM and comparing across-site additive and dominance genetic correlations from the MCM with those obtained from the simple MSM without genetic and environmental competition effects.

When we compared the across-site genetic correlations estimates from the two models, our results showed that slightly higher (growth traits) and lower (NSTR trait) ${\hat{r}}_{a_{d12}}$ were obtained from the MCM model compared with the MSM model. These results suggest that ignoring competition effects is not likely to change the magnitude of the across-site additive genetic correlation estimates for growth and stem straightness. In contrast to ${\hat{r}}_{a_{d12}}$, estimates of the across-site dominance genetic correlation $\left({\hat{r}}_{d_{12}}\right)$ for growth traits were lower than for the NSTR trait. This means that for growth traits, trees at site 1 do not have the same performance at the dominance genetic level as trees at site 2. On the contrary, the dominance effects do not change the stem straightness performance of the genotypes across the two sites evaluated. When we compared the ${\hat{r}}_{d_{12}}$ across models, the DBH and NSTR traits showed slightly higher values for the MCM model with respect to the MSM model (1.9% for DBH and 5.7% for NSTR), while TH trait showed lower value for the MCM model than the MSM model (36.4%). Nevertheless, the conclusion of ${\hat{r}}_{d_{12}}$ must be taken with caution since accurate estimation of non-additive dominance parameters requires large data files with a high proportion of full-sibs (Misztal 1997; Toro and Varona 2010). In that sense, the requirement of a large number of families to estimate dominance effects may not have been fulfilled in our study since only 37 of the 66 full-sib families evaluated are common to both sites, and 90% of these contain a maximum of 20 full-siblings at each site.

The multi-environment competition analysis based on MCM also allows the study of the across-site correlations between additive competition–competition genetic effects $\left({\hat{r}}_{a_{c12}}\right)$ and the across-site correlations between additive direct–competition genetic effects $\Big({\hat{r}}_{a_{d1c2}}$ and ${\hat{r}}_{a_{d2c1}}\Big)$. The values of ${\hat{r}}_{a_{c12}}$ were small and, in general, not significantly different from zero, indicating that competition additive genetic effects at one site would not be a good predictor of these effects at the other site (i.e., trees that are competitive at site 1 may not have the same performance at site 2). In general, the ${\hat{r}}_{a_{d1c2}}$ and ${\hat{r}}_{a_{d2c1}}$ were low and negative for all studied traits, excepting ${\hat{r}}_{a_{d1c2}}$ for DBH (Table 3). This suggests that there is a high level of uncertainty when predicting the competition breeding values in a non-tested environment based on the additive direct breeding values from a tested environment. Therefore, these results show that it is not feasible to predict the additive (or competition) genetic effects in a site based on the additive (or competition) genetic effects evaluated elsewhere.

Inter-tree competition may affect the efficiency of multi-environmental genetic selection. Therefore, a relevant question for breeding purposes is whether or not the predicted tree BVs from MSM and the predicted BVs from MCM lead to comparable rankings. In our study, positive and very high Spearman correlations between the BVs of both models were obtained for TH and NSTR traits. However, the moderate Spearman correlation observed for DBH could be indicating that for more competitive traits, the selected trees from any of these models could be very different. It would be difficult to identify which of the two predictors (MSM versus MCM) is closer to the real (unknown) values without a proper validation. However, we could assume that predicted additive direct breeding values from MCM are supposedly less affected by bias due to the fact of a proper account of competition, which is ignored under MSM. Such prediction could, therefore, lead to more accurate rankings of tree genetic merit. In that sense, ignoring competition in multi-environment analyses may substantially compromise selection decisions and genetic progress for growth in DBH. As this is the first study to use a multi-environment competition individual-tree mixed model in forest genetics, no comparison with other studies is possible. However, from a single-site univariate standard and competition-spatial mixed models, Cappa et al. (2016) showed a similar result for the top 10% ranked individuals for DBH. In addition, trees with the most competitive advantage over its neighbors for DBH (i.e., those with the most negative breeding values) were, in general, the best ranked and stable trees across the two models (Fig. 2), implying that the selection based only on the additive (or direct) breeding value could increase the level of competitiveness in DBH growth for the next population of improvement. As suggested by Hernández et al. (2019), a breeding strategy based on selecting less competitive genotypes might reduce the number of suppressed individuals in forests, increasing yields per unit area of future commercial plantations.

5 Summary and conclusions

Given the likely presence of genetic competition when testing trees across environments, we presented an extension of the single-environment competition individual-tree mixed model (MSM). To our knowledge, this is the first study that uses a multi-environment competition individual-tree mixed model (MCM) in forest tree breeding. We illustrated the use of the proposed model by fitting growth (DBH and TH) and stem straightness (NSTR) data from two full-sib progeny trials of Pinus elliottii var. elliottii × Pinus caribaea var. hondurensis assessed at age 10. The proposed MCM gave a better fit (lower AIC value) than the simpler model (i.e., one without competition effects or MSM). This was noticeable especially for traits showing the highest competition genetic effects in the study such as DBH and also to a lesser extent in TH. In general, accounting for competition increased the estimated direct additive variance and reduced the dominance genetic and residual variances compared with results from the MSM. This tendency was stronger for DBH due to its greater susceptibility to inter-tree competition. For TH and NSTR traits, the trend was also detectable, but weaker. Although ignoring competition genetic effects did not result in a significant change in the magnitude of the across-site additive genetic correlation, it led to different predictions of breeding values for DBH, compared to those of the model with competition MCM. Finally, across-site additive competition–competition genetic correlations were low for all traits (${\hat{r}}_{a_{c12}}$ < 0.24), suggesting that genetic component of competition among trees may depend on environmental circumstances. In summary, when traits are strongly affected by inter-tree competition, the use of the multi-environment competition individual-tree mixed model in MET analyses of forest genetic trials can improve the fit of the model, which could further benefit accuracy and bias in the prediction of breeding values. Moreover, this kind of analyses open the door to the use of competition in selection decisions with the aim of improving global at the level of the plantation.

Annex

1.1 Diagnosis of competition effects

The presence of competition effects for each trait–site combination was examined by using the isotropic empirical semivariogram of residual (see Fig. 1) from a multi-environment individual-tree mixed model with fixed effect of overall site means, and random plot, check lot, additive (breeding values) and dominance genetic effects, and random residuals (see details in Eq. (1). In the presence of inter-tree competition, the semivariogram depicts a decreasing variation at very short range (i.e., a sudden initial drop). For the first trial analyzed (site 1), while DBH and NSTR exhibited a moderate competition at the residual level, with a tendency toward a decreasing variation at very short range and a non-general trend at medium and large distances, TH displayed a slight tendency toward increasing variation at close distances. Site 2 showed a strong competition at the residual level with a sudden initial drop of the variation for the three traits analyzed. To identify competition effects at the genetic level, we then studied the correlation between additive direct and competition genetic effects (${\hat{r}}_{a_{dc}}$) for each trait within each site using a multi-environment individual-tree competition mixed model (see details in Eq. 2). A high and negative correlation between both types of genetic effects (${\hat{r}}_{a_{dc}}$ > − 0.3) suggests strong genetic competition (e.g., Resende et al. 2005). In our study, this correlation was negative for the three traits studied at both sites, except for NSTR (0.41) at site 1, while it was strong for DBH at sites 1 and 2 (− 0.36 and − 0.95, respectively) and NSTR at site 2 (− 0.42), and moderate for TH at sites 1 and 2 (− 0.12 and − 0.17, respectively).