Multi‐trait genomic selection for weevil resistance, growth, and wood quality in Norway spruce

Abstract Plantation‐grown trees have to cope with an increasing pressure of pest and disease in the context of climate change, and breeding approaches using genomics may offer efficient and flexible tools to face this pressure. In the present study, we targeted genetic improvement of resistance of an introduced conifer species in Canada, Norway spruce (Picea abies (L.) Karst.), to the native white pine weevil (Pissodes strobi Peck). We developed single‐ and multi‐trait genomic selection (GS) models and selection indices considering the relationships between weevil resistance, intrinsic wood quality, and growth traits. Weevil resistance, acoustic velocity as a proxy for mechanical wood stiffness, and average wood density showed moderate‐to‐high heritability and low genotype‐by‐environment interactions. Weevil resistance was genetically positively correlated with tree height, height‐to‐diameter at breast height (DBH) ratio, and acoustic velocity. The accuracy of the different GS models tested (GBLUP, threshold GBLUP, Bayesian ridge regression, BayesCπ) was high and did not differ among each other. Multi‐trait models performed similarly as single‐trait models when all trees were phenotyped. However, when weevil attack data were not available for all trees, weevil resistance was more accurately predicted by integrating genetically correlated growth traits into multi‐trait GS models. A GS index that corresponded to the breeders’ priorities achieved near maximum gains for weevil resistance, acoustic velocity, and height growth, but a small decrease for DBH. The results of this study indicate that it is possible to breed for high‐quality, weevil‐resistant Norway spruce reforestation stock with high accuracy achieved from single‐trait or multi‐trait GS.

Most of the genomic selection studies in tree species focused on modeling growth and wood traits, which are quantitative traits likely to be controlled by a large number of genes of small effects (Namkoong, Kang, & Brouard, 1988). Despite the fact that many tree breeding programs screen for biotic resistance against pests or disease (Mullin et al., 2011), to our knowledge, there is no study on the accuracy of GS for insect resistance in trees. One study tested different GS models to predict resistance to a pathogen in loblolly pine (Pinus taeda L.) (Resende, Munoz, et al., 2012;. The authors reported that fusiform rust resistance was likely controlled by a few genes of large effects since it was best predicted with Bayesian regression models that allowed for different variance of marker effects. Not only is the genetic architecture of the trait an important consideration in the choice of GS model, but one must also deal with the nature of the phenotypic data. Screening data for pest resistance is most often qualitative or semiquantitative, which needs appropriate statistical approaches that can handle non-normality of modeling errors. Thus, there is an opportunity for testing different GS approaches with weevil resistance data in Norway spruce, for example the Bayesian generalized linear regression models that support binary or ordinal data (Perez & de los Campos, 2014) or the threshold GBLUP model developed for ordinal data (Montesinos-López et al., 2015).
In practice, breeders generally need to consider multiple traits simultaneously in their selections for improved genetic stock. Besides reducing insect attack, improving growth and wood volume are major goals for Norway spruce, as it is for several plantation-grown conifers (Mullin et al., 2011). However, reduction of wood quality was observed in the past under selection for accelerated growth (Chen F I G U R E 1 [Box 1] . Genomic selection modeling and integration in tree breeding (a) Build GS models  Lenz, Cloutier, MacKay, & Beaulieu, 2010). Therefore, important wood traits for mechanical applications such as wood stiffness should be considered when performing selections (Lenz, Auty, Achim, Beaulieu, & Mackay, 2013). Mottet et al. (2015) concluded that selection for weevil resistance would not reduce height growth, but could affect diameter and therefore volume. However, the genetic relationships between weevil resistance and intrinsic wood properties such as wood density and stiffness have never been investigated. Thus, there is a need to understand the genetic correlations between weevil resistance, wood quality, and growth traits and to look at the possibility to combine them in a multi-trait genomic selection framework.
Multivariate genomic selection models can improve the accuracy of predictions by taking advantage of the genetic correlations between traits (Calus & Veerkamp, 2011). This ability is especially advantageous for prediction of traits that are costly or difficult to measure by conventional means on a large number of candidate trees, such as weevil resistance or wood quality, by using available correlated indicator traits. Simulation studies showed that multi-trait models can increase the accuracy of a target trait of low heritability when it is modeled together with genetically correlated indicator traits harboring high heritability (Guo et al., 2014;Jia & Jannink, 2012). In addition, the benefits of multi-trait GS are increased for target traits with scarce phenotypic records when they are coupled with intensively phenotyped indicator traits (Guo et al., 2014;Jia & Jannink, 2012;Schulthess et al., 2016). Few studies have applied multi-trait genomic selection methods to real plant breeding datasets (e.g., Bao, Kurle, Anderson, & Young, 2015;Fernandes, Dias, Ferreira, & Brown, 2018;Schulthess et al., 2016). In tree species, the accuracy of multi-trait versus single-trait GS models was only tested on a few traits using bivariate models in loblolly pine (Cheng, Kizilkaya, Zeng, Garrick, & Fernando, 2018;Jia & Jannink, 2012) and in Eucalyptus , and these studies did not investigate the effects of missing phenotypic data. Given the observed positive genetic correlations between weevil resistance and height growth (Mottet et al., 2015), and the possible correlations with wood quality traits, multi-trait GS models could improve the accuracy of predictions for traits that are difficult and expensive to assess on a large number of trees.
Once accurate genomic-estimated breeding values for each trait have been obtained from either single-trait or multi-trait models, it is possible to combine them into a selection index (SI; Hazel, 1943) and to rank individuals based on their overall performance across all traits. This strategy is especially useful in the presence of negative genetic correlations in order to select material that achieve a good balance in their performance for all traits of interest. In addition, GS is especially well suited to allow identifying correlation breakers in sufficient number, given the higher selection intensities that can be achieved by screening larger numbers of candidates than with conventional methods (Park et al., 2016). Hence, multi-trait genomic selection models and index selection are two different tools that can be combined to improve accuracies of breeding values and optimize genetic gains, respectively, in a multi-trait breeding program.
Here, we present a comprehensive genetic study of weevil resistance in Norway spruce and its relationships with growth and wood quality traits in the context of establishing a multi-trait genomic selection breeding program. Our objectives were to (a) better understand the genetic relationships between weevil attack and other growth and wood traits, in particular intrinsic wood quality traits; (b) evaluate the performance of different single-trait genomic selection models, especially for weevil resistance; (c) test the performance of multi-trait genomic selection models for predicting a target trait (weevil resistance or wood quality) when coupled with genetically correlated indicator traits (e.g., height growth); and (d) develop multitrait genomic selection indices for the production of high-quality and weevil-resistant seedling stock in Norway spruce.

| Genetic material and phenotyping
The data analyzed in this study are a subset of a larger breeding population derived from a partial diallel mating design (see Mottet et al., 2015 for more details). We focused our phenotyping and genotyping efforts on the trees planted on two sites affected by white pine weevil Each parent was crossed on average 2.3 times ( Figure S1). The two trials were established according to a randomized complete block design, each with five blocks and three-tree row plots for each family (tree spacing: 2.5 m × 2 m). Because of tree mortality, not all families were represented in every block (i.e., incomplete block design).
Tree height (Height 15 ) and diameter at breast height (DBH 15 ) were measured at age 15. The height-to-diameter ratio (Height 15 / DBH 15 ) was calculated as a proxy of stem taper. Three wood quality traits related to mechanical properties were assessed: average wood density and cellulose microfibril angle were determined from wood increment cores using X-ray densitometry and diffractometry at age 15 (Density 15 and MFA 15 ) as described in Lenz et al. (2017); acoustic velocity was measured at age 16 (Velocity 16 ) with the ST300 Hitman tool. Acoustic velocity is a proxy for wood stiffness or modulus of elasticity measured at standing trees Desponts, Perron, & DeBlois, 2017;Lenz et al., 2013).
For each of two surveys at ages 10 and 15, the presence/absence of weevil damage in the current year (coded 0/1) and in previous years (coded 0/1) were recorded. Current year damage was detected by inspecting the terminal shoot for emergence holes or death of the leader shoot. Weevil damage in previous years was visible by forks, curves, bayonets, or multiple stems. We calculated the cumulative number of attacks (CWA) as: where WA current10 and WA current15 are the presence/absence of attack at age 10 and 15, respectively; WA previous10 is the presence/absence of attack prior to age 10; and WA 11-14 is the presence/absence of attack between ages 11 and 14. The variable WA 11-14 , which was calculated to avoid double-counting previous attacks, was equal to 1 (presence) if there was an attack prior to age 15, but no attack at age 10 or prior to age 10. The resulting CWA variable is ordinal (ordered categories of 0, 1, 2, and 3 weevil attacks). Table 1 and Figure S2 provide the summary statistics and violin plots, respectively, for the traits assessed in this study. Figure S3 shows the spatial distribution of CWA values in the test sites.

| SNP genotyping
The 726 trees were genotyped using an Infinium iSelect SNP array (Illumina), which was assembled from a catalog of high-confidence gene SNPs obtained from exome capture and sequencing (Azaiez et al., 2018). The genotyping reproducibility rate was high (99.94%), as estimated from two positive controls replicated on each genotyping plate. From 5,660 successfully manufactured SNPs representing as many distinct gene loci well distributed over the 12 spruce linkage groups (Pavy et al., 2017), we retained a total of 3,914 SNPs with call rate ≥90% (average call rate of 99.6%), minor allele frequency (MAF) ≥0.005, and a fixation index |F IS | < 0.50. SNPs were well distributed across MAF classes with 86% of the SNPs with MAF ≥0.05 ( Figure S4). Missing genotypes (only 0.4% of genotypes) were imputed using a k-nearest neighbor method based on linkage disequilibrium (LD-kNNi) with the software LinkImpute (Money et al., 2015). The software estimated an accuracy of 0.83 for imputed genotypes by randomly masking genotypes.

| Relationship matrices and pedigree verification
All analyses were performed in the R v.3.3.1 environment (R Core Team, 2016). A pedigree-based relationship matrix (A) and its inverse were first computed based on the registered pedigree information using the function "asreml.Ainverse" of the R package ASReml-R v.3.0 (Butler, Cullis, Gilmour, & Gogel, 2007). For use in genomic selection (GS) models, the realized genomic relationship matrix (G, Figure S5) was computed from the marker data with the "A.mat" (1) CWA = WA previous10 + WA current10 + WA 11−14 + WA current15 TA B L E 1 Phenotypic means, standard deviations (SD), and coefficients of variation (CV) for each site and across sites for the 714 trees retained for analyses  function of the R package rrBLUP (Endelman & Jannink, 2012) with the default options, which was equivalent to the formula described by VanRaden (2008). A comparison of the A and G matrices revealed 11 "misclassified" trees that presumably did not belong to their expected cross. In addition, one tree was deemed an outlier because it had abnormally high average wood density (Density 15 ), after checking model residuals (see Equation (2) below). After removing misclassified and outlier trees, 714 trees genotyped on 3,914 SNPs were used in subsequent analyses. The resulting A and G matrices were highly correlated with a Pearson r = 0.94 ( Figure S6). Large values of G within each class of A are mostly due to two inbred families ( Figure S7).

| Heritability and genotype-by-environment interactions
First, variance components, heritability, and breeding values were estimated using the conventional pedigree-based (ABLUP) or the genomic-based (GBLUP) individual-tree mixed models (the so-called "animal model") in ASReml-R v.3.0: where y is the phenotype; μ is the overall mean; s is the fixed site effect; b(s) is the random effect of block within site, with b (s) ∼ N(0, 2 b I b ); a is the random additive genetic effect, with a ∼ N 0, 2 a A ; sa is the random interaction of site with additive genetic effects, with sa ∼ N(0, 2 sa I s ⊗ A); and e is the residual term, assuming homogeneity across sites with e ∼ N(0, 2 e I e ). For the ABLUP method, the matrix A is the pedigreebased relationship matrix, which was replaced by the realized genomic relationship matrix G for the GBLUP method (a ∼ N 0, 2 a G and sa ∼ N(0, 2 sa I s ⊗ G)). The X and Z matrices are incidence matrices of their corresponding effects, and I x is an identity matrix of its proper dimension. The symbol ⊗ refers to the Kronecker product. To test the hypothesis of greater than zero variance for each effect (H 0 : σ 2 = 0; H 1 : σ 2 > 0), we performed a likelihood-ratio test with one degree of freedom between the full model in Equation (2) and a reduced model without the effect to be tested. The dominance effect was not included in models because it was not significant for all traits under study (Table S1), and, according to BIC, the fit of models including dominance was similar or worse than the models including the additive effect only (Table S2). We compared model (1) with a model fitting a different residual variance for each site and found that the latter model was slightly better for three out of the seven traits studied (Table S3).
However, the breeding values between both approaches were highly correlated for all traits, with r > 0.99. Thus, we opted for the simpler model in Equation (2) to keep ABLUP and GBLUP models comparable with marker-based genomic selection models. Finally, inspection of the residuals from Equation (2) versus the distance in X and Y in the trials showed no detectable spatial patterns (not shown).
Narrow-sense individual heritability was estimated as: The size of genotype-by-environment interaction (GxE), or type-B correlation (r B ), was estimated as: Standard errors of heritability and type-B correlation estimates were obtained using the delta method (pin function from the R package nadiv; Wolak, 2012).

| Correlations between weevil resistance, wood quality, and growth traits
To estimate single-site phenotypic and genetic correlations between traits, bivariate models were run for all pairs of traits in ASReml.
Bivariate models were run for each site separately because of convergence problems of the multi-site model due to large GxE for some traits (e.g., DBH 15 , see Results). The following model was fitted: where y i and y j are the stacked vectors of phenotypic observations for trait i and trait j, respectively; t is the vector of fixed effects of traits (i.e., the grand mean for each trait); b(t) is the random effect of block The matrix A (ABLUP) was replaced by the realized genomic relationship matrix (G) for the GBLUP method. The matrices V B , V A , and V R are 2 x 2 variance-covariance matrices defined by the correlation of effects between traits (r b , r a , and r e , respectively) and unique variances for each trait (i.e., CORGH in ASReml). To facilitate convergence, we provided starting values for the variance components in V B , V A , and V R matrices that were taken from the results of the single-trait models (Equation (2), Table S4). For r b , r a , and r e , the starting value was set to 0 (no correlation). The genetic correlation between traits was directly provided by the estimated parameter r a and the phenotypic correlation was calculated as: where COV(i,j) p is the phenotypic covariance between traits, and ̂2 pi , ̂2 bi , ̂2 ai , and ̂2 ei are the estimated phenotypic, block, additive, and residual variance of trait i (same for trait j), respectively. The significance of the genetic correlation (H 0 : r a = 0; H 1 : r a ≠ 0) was tested by performing a likelihood-ratio test with one degree of freedom between the full model in Equation (5) and a reduced model assuming r a = 0 (i.e., a diagonal V A matrix). The significance of the phenotypic correlation (H 0 : r p = 0; H 1 : r p ≠ 0) was tested by performing a likelihood-ratio test with three degrees of freedom between the full model in Equation (5) and a reduced model assuming no correlation between traits (r b = 0, r a = 0, and r e = 0). The single-site heritability of trait i was given by:

| Single-trait genomic selection models
We evaluated four single-trait GS methods: GBLUP, Bayesian ridge regression (BRR), BayesCπ, and, in addition for CWA, threshold GBLUP (TGBLUP) developed for ordinal traits. GBLUP was implemented as described in Equation (2) in ASReml. The G matrix describes the realized genomic relationships between trees, which better account for within-family Mendelian sampling, as well as deeper (unknown) pedigree relationships. GBLUP relies on the "infinitesimal" model of quantitative genetics, assuming that the genetic control of complex traits is equally distributed across many (infinite) loci with small effects (Falconer & Mackay, 1996 (2015). To implement TGBLUP, the same model as in Equation (2) was fitted, but the response type was set to "ordinal" where y is the phenotype; μ is the overall mean; s is the fixed site effect (i.e., modeled with a flat prior); b(s) is the random effect ; a m is the random additive effect of markers; and e is the residual term, with e ∼ N(0, 2 e I e ). Note that compared with ABLUP and GBLUP models, the genotype-by-environment interaction (GxE) component was not fitted in BRR and BayesCπ models, thus assuming that marker effects were stable across environments. To account for different distributions of marker effects, the prior for a m changes depending on the method (see Appendix S1 for a full description). Briefly, the method BRR is a Bayesian version of ridge regression, in which marker effects are normally distributed (i.e., Gaussian prior) and have identical variance (a m ∼ N(0, 2 m I m )). In BRR, all markers have a nonzero effect, and so this method is appropriate for traits controlled by a large number of genes with small effects. In contrast, the method BayesCπ takes into account that only a proportion π of markers have an effect, while a proportion (1 − π) of marker effects are shrunk toward zero (Habier, Fernando, Kizilkaya, & Garrick, 2011). For the BRR and BayesCπ methods, the response type was set to "ordinal" (probit link function) for the trait CWA , and to "gaussian" for all other traits. BGLR was run for 50,000 iterations and a thinning interval of 20, with the first 15,000 iterations discarded as a burn-in. Genomic-estimated breeding values (GEBVs) were obtained by summing over the effects of all markers, with GEBV i = ∑ m j=1 Z' ijâj , where â j is the estimated effect of the j th marker, and Z ′ ij is an indicator of the genotype of individual i at the j th marker.

| Multi-trait genomic selection models
Genomic selection models that incorporated the information of multiple correlated traits into a single analysis were evaluated using GBLUP multivariate models in ASReml. To facilitate the convergence of multivariate models, the modeling was done in two steps. First, phenotypes were adjusted for block and site effects (y*) by taking the residuals (e) of a model that included a fixed site effect (s) and a random block within site effect (b(s)): y = + Xs + Zb (s) + e. After adjusting phenotypes, the portion of GxE due to rank changes in different sites remains, but the "level-of-expression GxE" (i.e., spread of breeding values across environments) is controlled for (Li, Suontama, Burdon, & Dungey, 2017). Second, a multivariate model with p traits was fitted: where y * i are the stacked vectors of adjusted phenotypes from trait i to trait p; t is the vector of fixed effects of traits (i.e., the grand mean for each trait); a(t) is the random additive effect within trait, with a (t) ∼ N 0,G ⊗ V A ; and e is the residual error, with e ∼ N 0,I e ⊗ V R . The matrices V A , and V R are p × p variance-covariance matrices, defined by correlations between all pairs of traits and unique variances for each trait. GxE (as rank-changes interaction) was not fitted to simplify the model and facilitate convergence. We provided starting values for the variance components in V A and V R matrices that were taken from the results of the single-trait models. We obtained GEBVs of individual trees for each trait separately from the BLUPs of the random additive effect (a) within trait.
We evaluated the performance of multi-trait GS models for predicting a target trait, when coupled with genetically correlated indicator traits. We chose three target traits, the cumu- a two-trait model with one of the target traits and Velocity 16 ; and (d) a three-trait model with one of the target traits, Height 15 / DBH 15 ratio, and Velocity 16 . To simulate a situation where the target trait was only measured for a smaller subset of the trees, the percentage of missing data in the training set (see below) was varied from 0% to 90% by randomly adding missing values to the phenotypes of the target trait. We repeated this process 100 times for each trait and level of random missing values in crossvalidation (see below).
The full set of individual trees was randomly split into tenfold, each containing ~10% of the trees from each family. For each round of CV, ninefold (~642 trees or 90%) was used in model training, which was used to predict the breeding values for the remaining fold (~71 validation trees or 10%). This tenfold CV was repeated ten times, for a total of 100 models for each trait.
The predictive ability (PA) of the models was evaluated as the Pearson correlation coefficient between the predicted breeding values of the validation trees and the adjusted phenotypes (y*) for block and site effects. The predictive accuracy (PACC) of models was estimated from PA as PACC = PA∕ √ĥ 2 ind (Dekkers, 2007;Legarra, Robert-Granie, Manfredi, & Elsen, 2008). For all the methods tested (ABLUP, GBLUP, TGBLUP, BRR, BayesCπ, and multi-trait GBLUP), we used the ĥ2 ind estimated from the GBLUP model (Equations 2 and 3) as our best estimate of heritability for the calculation of predictive accuracy (PACC). PA and PACC were calculated within each fold to avoid including a fold effect, then averaged across folds and repetitions to obtain standard errors (Isik, Holland, & Maltecca, 2017).
For multi-trait models, the prediction ability (PA) and prediction accuracy (PACC) of multi-trait GBLUP models were compared with the equivalent single-trait GBLUP models using adjusted phenotypes as a response variable (see Appendix S2).

| Genetic gains and multi-trait selection indices
Expected genetic gains from single-trait selection were calculated as the mean estimated breeding value of the top 5% trees for each trait separately, as estimated from the single-trait ABLUP (EBVs) or GBLUP (GEBVs) analysis (Equation 2). These estimated gains represented the maximum possible gain for each trait. To estimate genetic gains in a multi-trait selection context, we combined four traits of economic interest into a SI as follows: where Height 15EBV , CWA 6.15EBV , Velocity 16EBV , and Density 15EBV are the BLUP estimated breeding values from the single-trait ABLUP or GBLUP analysis (Equation 2) for the corresponding trait; and w i are the relative weight given to each trait, with the restrictions: 0 ≤ w i ≤ 1 for all traits, and We assigned a negative sign to CWA since a decrease in the value of this trait (fewer weevil attacks) represents an improvement. Breeding values for each trait were scaled to a variance of one (already centered) prior to SI calculations. DBH 15 is an economically important trait, but was excluded from the SI scenarios because its heritability was not significantly different from zero (see results). The four weight coefficients were varied between 0 and 1 by intervals of 0.05, resulting in 1,771 different selection indices. For each SI, trees were ranked according to decreasing values of the index (I) and the top 5% trees were selected to calculate the expected genetic gain. For each trait, the relative genetic gain (%) was calculated as the ratio of the expected gain to the maximum possible gain from single-trait selection. We present, a first SI scenario (SI-1) that corresponded to the priorities for the Norway spruce breeding program in Québec, which put more emphasis on weevil resistance (w 2 = 0.6), followed by growth, represented here by height growth (w 1 = 0.3), and acoustic velocity (w 3 = 0.1). We chose to present two further SIs that maximized the total relative gain of the (10) The model fitted is described in Equation (2). b Level of statistical significance: *p < 0.05; **p < 0.01; ***p < 0.001. c See Table 1 for full description of traits.
TA B L E 2 Individual narrow-sense heritability (ĥ 2 ind ) and type-B genetic correlation (r B ) estimates (standard errors in parentheses) using the single-trait ABLUP and GBLUP methods for the across-site analyses a . For heritability estimates, the significance of the additive variance component is shown (see Table  S4). For type-B genetic correlations, the significance of the site × additive variance component is shown b . A significant site × additive variance indicates significant genotype-by-environment interaction (i.e., smaller values of r B ) following traits: Height 15 , CWA, and Velocity 16 (SI-2) and in the other case Height 15 , CWA, Velocity 16 , and Density 15 (SI-3).

| Heritability and genotype-by-environment interaction
In the across-site analyses, heritability ranged from 0 to 0.47 using the ABLUP method and from 0 to 0.29 using GBLUP (

| Correlations between weevil resistance, growth, and wood quality traits
Phenotypic and genetic correlations between traits using the GBLUP method are presented for each site separately, in Table 3 for GPI, the most severely affected site by weevil attacks, and in Table 4 for STM, the site that was moderately affected. Results using ABLUP were similar (GPI : Table S5; STM: Table S6), and so we present below only the results using GBLUP. Weevil resistance was significantly correlated with growth and wood quality traits.
On both sites, we found a moderate negative genetic correlation TA B L E 3 Site GPI: phenotypic (r p , above diagonal) and genetic correlations (r a , below diagonal) between traits calculated with the GBLUP method a . Diagonal elements indicate the single-site narrow-sense heritability (ĥ 2 The model fitted is described in Equation (5)  r a = 0.55). Overall, these results indicated that weevil-resistant genotypes were taller, with larger height-to-DBH ratio and higher wood stiffness as measured by acoustic velocity.

| Accuracy of single-trait genomic selection models
The four tested single-trait genomic selection (GS) methods (GBLUP, TGBLUP, BRR, BayesCπ) and the conventional pedigree-based method (ABLUP) resulted in similar predictive abilities and predictive accuracies. The PA, which was defined as the correlation between the predicted breeding values for the validation trees and the phenotypic values, ranged from 0.10 for DBH 15 to 0.46 for Velocity 16 ( Figure 4) For the cumulative number of weevil attacks (CWA), the three methods that considered ordinal data, namely TGBLUP, BRR, and BayesCπ, had PA and PACC similar to those of GBLUP, which assumed normality of residuals ( Figure 4). In addition, genomic-estimated breeding values obtained from TGBLUP and GBLUP were highly correlated (Pearson r = 0.997), and the heritability estimates were within the same range (GBLUP: ĥ2 ind = 0.27 (0.07); TGBLUP: ĥ2 ind = 0.30 (0.08)). Thus, the assumption of the normality of residuals in GBLUP ( Figure S8) did not appear to affect heritability estimates and the performance of the model in our dataset.

| Accuracy of multi-trait genomic selection models
Multi-trait GBLUP models were used to predict the breeding values of a target trait (CWA, Density 15 , or MFA 15 ), when combined with genetically correlated indicator traits (Height 15 , Velocity 16 , Height 15 / DBH 15 ratio). The cumulative number of weevil attacks (CWA), average wood density at age 15 (Density 15 ), and microfibril angle at age 15 (MFA 15 ) were chosen as target traits because they are rather cumbersome and expensive to assess on a large number of trees, whereas indicator traits are easier to track for the majority of trees in a breeding population. The multi-trait models were compared with the single-trait GBLUP models. As expected, when the percentage of missing phenotypic data increased from 0% to 90% for the target traits CWA, Dens 15 , and MFA 15 , the predictive accuracy (PACC) of the single-trait models (dashed gray line in Figure 5) sharply dropped from 0.83 to 0.59, from 0.71 to 0.43, and from 0.91 to 0.55, respectively. Similar trends were found for the PA ( Figure S9).
The PACC of multi-trait models was not improved over singletrait models when all of the phenotypic data for the target trait were included for model training (0% missing data) ( Figure 5; standard F I G U R E 4 (a) Predictive ability (PA) and (b) predictive accuracy (PACC) of the single-trait genomic selection models (GBLUP, BRR, BayesCπ) and the conventional pedigree-based model (ABLUP) tested in this study. For the cumulative number of weevil attacks (CWA), three models accounted for ordinal data type, namely the threshold GBLUP model (TGBLUP), BRR, and BayesCπ, while ABLUP and GBLUP assumed that errors were normally distributed. Error bars indicate the standard errors of the estimates. The PACC of models for the trait DBH 15 was not calculated because the estimated heritability was null. See Table 1 for full description of traits (a) (b) errors are given in Table S7). However, for the target trait CWA with 40% or more missing data, the accuracy of multi-trait models was improved as compared with that of the single-trait model (Figure 5a). F I G U R E 5 Predictive accuracy (PACC) of GBLUP multi-trait genomic selection models for predicting the target traits: (a) the cumulative number of weevil attacks (CWA); (b) Density 15 ; and (c) MFA 15 . The different colored lines represent different multi-trait models with different indicator traits. The dashed gray line is the single-trait GBLUP model. The percentage of missing phenotypic data for the target trait in the training sets was varied from 0% to 90% (x-axis), while 100% of the training data was retained for the indicator traits. See Table 1   For the target traits Density 15 and MFA 15 (Figure 5b,c), multitrait models did not improve PACC over single-trait GBLUP, even when the target trait had large amount of missing data. For MFA 15 , some of the multi-trait models showed even a lower PACC than the single-trait models.

| Genetic gains and multi-trait selection indices
The expected genetic gains, expressed as a percentage of the phenotypic mean when selecting the top 5% trees for each trait separately (single-trait selection), are given in Table S8. For the cumulative number of weevil attacks (CWA), a positive gain represents a reduction of the number of attacks (i.e., higher resistance). CWA could be improved by as much as 78% (ABLUP) or 66% (GBLUP), while other traits showed much smaller gains in the 4%-20% range.
The large observed gains for CWA are likely due to a moderate-tohigh heritability of CWA, a large coefficient of phenotypic variation (Table 1), and the non-normal distribution of the trait.
The expected gains from the multi-trait selection indices (SIs) are shown in Table 5. Results from ABLUP and GBLUP were similar, and so, we present only the GBLUP results below. The first selection index scenario (SI-1) presents the planned breeding focus and put most emphasis on weevil resistance (w 2 = 0.6), followed by Height 15 (w 1 = 0.3), and by Velocity 16 (w 3 = 0.1). This scenario achieved the largest gains in terms of reduction of cumulative weevil attacks (CWAs) (83% of maximum possible gains from single-trait selection in

| Genetic control of weevil resistance and its relationships to growth and wood quality traits
Tree breeding uses the existing natural intraspecific variation to identify superior genotypes with desirable attributes. In this context, the heritability of natural resistance to pests may be seen as an indicator for the evolutionary potential of the species (Charmantier & Garant, 2005;Geber & Griffen, 2003). In a recently introduced species such as Norway spruce to North America, natural selection for resistance to native pests may have acted only for one or two generations at best. Nevertheless, we showed that resistance to white pine weevil attack was under moderate-to-high genetic control. Individual heritability estimates (ABLUP) were slightly higher than in an earlier study by Mottet et al. (2015), who combined data from more families and tests than this study, but values were in the same range than those observed in the native interior spruce (King et al., 1997) and slightly higher than those in the native Sitka spruce (King, 2004). In its native range, Norway spruce suffers damages from another weevil species, the large pine weevil (Hylobius abietis), mostly at the seedling stage after clearfelling operations (Day & Leather, 1997). Although both weevil species do not attack trees at the same developmental stage, both feed on the bark and phloem, so it is likely that resistance mechanisms to different weevils are partly related. Zas et al. (2017) found moderate family heritability and low GxE for resistance to the pine weevil in Norway spruce. Norway spruce has also faced outbreaks of bark beetles, such as Ips imitinus and Ips typographus, after major storms such as those that occurred in Central Europe in the 1990s (Wermelinger, 2004). Resin canal traits relevant for constitutive resistance against bark beetles were found to be under strong genetic control (Rosner & Hannrup, 2004).
Thus, it is likely that genetic variation at resistance genes for the North American white pine weevil was already available at the time of introduction of Norway spruce in North America, as the result of natural selection against insect pests in Europe, with opportunities for rapid change in allele frequencies at resistance loci.
In our study, both CWA and tree height at age 15 were under significant genetic control and the genetic correlation between these traits was negative, but significant only for the site that suffered the most frequent weevil attacks (GPI). These results are not surprising given that height growth is mechanistically stunted in consequence of attack, while the tree continues to grow in DBH and volume. Indeed, with increasing age, a higher height/DBH ratio and thus a lower stem taper was observed in more weevilresistant trees (Holst, 1955;Mottet et al., 2015; this study). King et al. (1997) also reported negative genetic correlations between CWA and tree height for interior spruce in British Columbia, both before and after the occurrence of weevil attacks. They concluded that inherently faster growing families have higher level of genetic resistance. In Norway spruce, Mottet et al. (2015) found negligible genetic correlations between weevil attacks and tree height measured before the majority of weevil attacks (age 5) and concluded that genetic improvement for resistance to white pine weevil would not adversely affect growth. On the other hand, the strong positive genetic correlation between resistance to weevil attacks and height at age 15 found on site GPI in this study may indicate that they are controlled by common genes, which was also suggested by an earlier QTL study in interior spruce (Porth et al., 2012). In addition, resistance mechanisms to European bark beetles such as resin canal traits were found to be positively genetically correlated with both height growth and DBH (Rosner & Hannrup, 2004). Overall, our results suggest that accelerated breeding of resistant seed stock through genomic selection tools will result in taller trees, either because the leaders of resistant trees will be less affected by attacks, or because alleles underlying growth genes will be simultaneously favored.
The genetic control of variation in wood quality traits, such as average wood density and acoustic velocity as a proxy for wood stiffness, was in the same range as that for weevil resistance. Compared to recent wood quality studies of Norway spruce from Scandinavia (Chen et al., 2014, our heritability estimates were lower for average wood density, but higher for acoustic velocity. Heritability estimates for MFA were low compared with those reported in previous studies (Chen et al., 2014;Lenz et al., 2010), which is most likely related to the measuring approach used in the present study where only the last ring was assessed. The moderate heritability and the resulting sizeable genetic gain observed here for acoustic velocity make it a promising quick-assessment trait for improving wood stiffness and product quality (Lenz et al., 2013). In addition, wood stiffness can be improved together with weevil resistance since acoustic velocity was negatively genetically correlated with the cumulative number of weevil attacks. Moreover, weevil resistance and wood traits showed low genotype-by-environment interactions, which should facilitate reforesting selected planting stock across larger breeding zones. Average wood density did not appear to be the best trait for overall improvement of Norway spruce in the tested conditions given its negative correlations with growth traits and slightly higher genotype-by-environment interaction. Overall, given their positive genetic correlations and moderate-to-high heritability, weevil resistance, acoustic velocity, and height growth appear as excellent candidate traits for simultaneous genetic improvement.

| Genomic selection models for accurate and hastened selection of best genotypes
Genomic selection can significantly shorten breeding cycles through the prediction of breeding values of nonphenotyped material using their genomic profiles and allows screening more candidates for increase selection intensity or multi-trait selection. Compared with the conventional pedigree-based approach (ABLUP), genomic selection models using GBLUP, BRR, or BayesCπ had comparable PA and accuracy, confirming early proof-of-concept studies in other conifers and spruces (Beaulieu, Doerksen, Clément, et al., 2014;Gamal El-Dien et al., 2015;Lenz et al., 2017;Ratcliffe et al., 2015). In BRR and BayesCπ, we did not fit genotype-by-environment interactions as opposed to ABLUP and GBLUP methods, but this did not affect PA nor predictive accuracy.
Predictive ability was related to heritability and was lowest for the low heritability traits MFA 15 and DBH 15 . This is because the proportion of phenotypic variation that can be explained by additive genetic effects is smaller for these traits. Comparisons of predictive accuracy of breeding values across studies are difficult because of the absence of a standard way of determination. For similar growth and wood quality traits, predictive accuracy was slightly higher than previously observed for native Norway spruce (Chen et al., 2018) or for white spruce (Beaulieu, Doerksen, Clément, et al., 2014;Gamal El-Dien et al., 2015;Ratcliffe et al., 2015), and it was lower and more variable than that observed for black spruce (Picea mariana [Mill.] B.S.P.) (Lenz et al., 2017). In contrast to those previous GS studies, we assumed herein that true breeding values were unknown and calculated the predictive accuracy by dividing the PA by the square root of heritability (Dekkers, 2007;Legarra et al., 2008). Hence, predictive accuracy should not be correlated with heritability, but depends on other characteristics of the dataset, such as the accuracy of phenotypic measurements, the effective population size, and the genetic architecture of traits (Grattapaglia & Resende, 2011). However, the precision of the present estimates of predictive accuracy may be affected by the precision of heritability estimates.
In theory, GS should perform better than pedigree-based models because it allows correcting pedigree errors and capturing the within-family variation resulting from Mendelian segregation (Grattapaglia et al., 2018). Thus, in a forward selection scenario with nonphenotyped young material, the selection of the best individuals within families becomes possible with GS. In the present context of small size of the breeding population, we conclude that the estimated genomic predictions for resistance to weevil attack, tree height, height/DBH ratio, and wood quality traits allow for accurate and hastened selection of best candidates based on their genomic profiles.

| Evidence for polygenic control of weevil resistance in spruces
To our knowledge, the present study is the first one applying genomic selection to breed for insect resistance in conifers. We therefore had particular interest in testing different algorithms that considered the ordinal distribution of this trait and that reflected different distributions of marker effects. First, we found no advantages of using methods that accounted for ordinal data (threshold GBLUP, BRR, and BayesCπ), as compared with approaches that assumed normality of residuals (ABLUPs and GBLUPs). Second, the BayesCπ algorithm, which assumed that only a portion of genes had an effect, did not lead to improvement of PA or accuracy compared with methods assuming that all genes had small effects (GBLUPs, BRR). In BayesCπ, the estimated proportion of markers having an effect ( ≅ 0.50) for weevil resistance was in the same range as that for other traits (Table S9). In contrast, Resende, Munoz, et al. (2012) and  obtained a higher predictive accuracy using BayesCπ for fusiform rust resistance in loblolly pine, which suggested the presence of large effect genes. Our findings can be explained by two possible phenomena: (a) weevil resistance is effectively controlled by many genes of small effects, with no detectable major gene effects; or (b) none of the sampled SNPs was in close linkage with genes having effects for resistance. It is difficult to discriminate between both hypotheses, especially since GS models with the current genome coverage and small size of the breeding population should mostly retrace relatedness between trees from long-range linkage disequilibrium (e.g., Lenz et al., 2017) and could thus provide high predictive accuracies following both sets of conditions. However, polygenic control of weevil resistance seems plausible since earlier reports associated resistance to weevil attack with different constitutive and induced mechanisms. Physical defense barriers to weevil were described through sclereids and stone cells in the bark of sitka and hybrid spruces (King, Alfaro, Lopez, & Van Akker, 2011;Whitehill et al., 2016). In addition, chemical defense strategies are also at play through the presence of abundant constitutive resin ducts (King et al., 2011;Rosner & Hannrup, 2004) and the additional formation of traumatic resin ducts (Poulin, Lavallee, Mauffette, & Rioux, 2006), together with increase in terpenoid metabolite production (Robert et al., 2010) following insect feeding. Porth et al. (2012) identified many candidate genes for weevil resistance in interior spruce, including some master regulatory genes. Moreover, transcriptome studies in Sitka spruce revealed several thousand differentially expressed genes between resistant and sensible genotypes, as well as following weevil wounding and feeding (Ralph et al., 2006;Whitehill et al., 2019). Hence, weevil resistance in Norway spruce is most likely polygenic given the complex nature of physical and chemical defense mechanisms, but whether some larger gene effects exist remains to be elucidated.

| Multi-trait GS as a tool to improve accuracy of scarcely phenotyped traits
The joint modeling of multiple traits can benefit from genetic correlations between traits and increase predictive accuracy (Calus & Veerkamp, 2011;Guo et al., 2014;Jia & Jannink, 2012). In Eucalyptus, Cappa et al. (2018) reported modest improvements in the predictive accuracy of breeding values from multi-trait over single-trait GS models (~2%-4%) when a low heritability target trait (tree height) was coupled with a highly genetically correlated trait (DBH, r a = 0.92). Other empirical plant or tree breeding studies found little benefits of using multi-trait GS to predict nonphenotyped selection candidates when 100% of the individuals in the training set were phenotyped (Bao et al., 2015;Cheng et al., 2018;Fernandes et al., 2018;Jia & Jannink, 2012;Schulthess et al., 2016). Similarly, we found that multi-trait GBLUP did not outperform single-trait models when all individuals in the training set were phenotyped for the target trait. Thus, in cases when phenotypic information is balanced across traits, single-trait modeling is the recommended method given that they harbor reduced model complexity (Schulthess et al., 2016 (Schulthess et al., 2016). In this study, convergence problems precluded us to fit simultaneously more than three traits. Multi-trait Bayesian models are expected to be more efficient than multi-trait GBLUP when the number of traits increases (Calus & Veerkamp, 2011). Overall, the benefits of adding more traits seems limited, but this needs to be further tested with different GS models and with traits combinations of different levels of heritability and genetic correlations.

| Index selection provided positive gains for most focus traits
Index selection is a commonly used tool in tree improvement to se- Norway spruce in the Canadian breeding and forestry contexts, and the relative weights that we have chosen for the SI-1 scenario are currently the most likely to be implemented in the short term but are still of an indicative nature.
In this study, favorable correlations between weevil resistance, height growth, and acoustic velocity resulted in only minor tradeoffs among those traits in the tested selection indices. However, the improvement for DBH or the closely related volume (not estimated in this study) appears more challenging. The low genetic control and the important genotype-by-environment interaction observed for DBH (see also Mottet et al., 2015) make it difficult to predict the effect of selecting for weevil resistance on this trait. However, the possible loss in lumber volume would likely be compensated by less log defects due to increased resistance of plantations to weevil attacks (Daoust & Mottet, 2006). Thus, improving for weevil resistance would undoubtedly benefit Norway spruce as an exotic plantation species in eastern Canada.

| CON CLUS IONS
In this study, we investigated the genetic control of weevil resistance and its relationship with wood and growth traits in Norway spruce as an exotic plantation species. We further integrated these traits into a multi-trait genomic selection (GS) framework. We found that in such a realistic context, it was possible to improve significantly for weevil resistance using GS, given that weevil resistance was moderately to highly heritable and that it was positively genetically correlated with the height/DBH ratio and wood stiffness (acoustic velocity). By taking advantage of these existing genetic relationships, we showed that multi-trait genomic selection models could improve the accuracy of the prediction for a scarcely phenotyped target trait (weevil resistance) by using the information from a readily available indicator trait. Finally, by combining multiple correlated traits into a selection index, we obtained the best compromise for all traits of interest that corresponded to the priorities of the breeders. Thus, this integrated approach showed how genomic selection can be used to breed simultaneously for taller, stiffer, and more weevil-resistant Norway spruces. We conclude that single/multi-trait GS models and index selection are efficient selection tools that can be integrated into operational breeding programs to accelerate the realization of genetic gains for most traits of interest.
Another advantage of integrating genomic selection to a breeding program is that, once the breeding population has been genotyped, models can easily be recalibrated for additional traits, such as pest resistance or resilience to episodic climate extremes such as droughts. Such traits may be costly and difficult to measure in forest trees (e.g., Housset et al., 2018) and the use of correlated indicator traits in multi-trait models would reduce phenotyping costs. Mass phenotyping using remote sensing technologies could also help identify indicator traits that are correlated with economically important ones (Dungey et al., 2018). The increase of information from multiple correlated traits will provide opportunities to test novel multi-trait models and to improve genomic selection accuracies. Finally, given that large numbers of candidate trees can be genotyped at the very juvenile stage, multi-trait GS will result in highly significant reductions in testing time for such long-lived species, while allowing to increase selection intensity compared to more conventional selection approaches.