Eucalyptus population and growth data

The Eucalyptus population employed in this study belongs to the breeding program of International Paper do Brasil. It involved a progeny trial encompassing almost exclusively E. grandis × E. urophylla F₁, backcrosses, and F₂ full-sib families. The trial was planted in randomized complete blocks, with five trees per plot comprising 58 full-sib families generated by an incomplete half-diallel mating design, totaling 2,784 trees planted in 2006 in Brotas, São Paulo State, Brazil (22°17′05″S 48°07′38″W). A subset of 1,000 trees out of the 2,784 in the trial was sampled for measuring growth and wood properties and DNA marker genotyping (Fig 1). Sampling was carried out in a stratified manner to keep equivalent representation of the majority of families, with approximately 20 to 25 trees per family and avoiding off-type trees (e.g. branched, diseased, deformed, etc.). The sampled set of 1,000 trees included 45 full-sib families derived from 46 parents considered genetically unrelated based on the fact that they were originally selected in unrelated families. This sample included 610 trees (61%) of 25 F₁ hybrid families of E. grandis and E. urophylla, 366 trees (37%) of 19 backcross and F₂ hybrid families of E. grandis and E. urophylla, and 24 (2%) trees of one E. grandis × E. camaldulensis hybrid family (Table 1). For brevity, the population will be called an urograndis hybrid population. Diameter at breast height (DBH) and tree height (H), were measured in 2011 at age five, and wood volume (VOL) in cubic meters was estimated using a taper factor of 0.45 [49]. Mean annual increment (MAI) in m³.ha^-1.year^-1 was calculated by multiplying the total tree volume by 1,200 trees per hectare and dividing the result by 5 years of growth.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Summary flowchart of the experimental steps employed to estimate genetic parameters for growth and wood traits in Eucalyptus hybrid urograndis.

Growth data and wood sample NIRS spectra were collected from the sampled subset of 1,000 trees across full-sib families. A subset of between 200 and 350 trees selected based on maximizing NIRS spectra distance was phenotyped (wood chemical and physical traits) and data used to develop acceptable NIRS calibration models used to predicted lignin and wood density for the remaining 650 trees. Pedigree, genotypes (SNPs and DArT-seq), growth and wood trait data, either directly measured for 200 trees (cellulose, hemicellulose, microfibril angle, fibers, coarseness) or predicted for 1000 trees (lignin, wood density) were employed for genetic parameter estimation. Block arrows indicate step processes, thin dashed arrows indicate data or sample use.

https://doi.org/10.1371/journal.pone.0218747.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Number of families and trees sampled for the different mating types and their respective species and hybrids involved.

https://doi.org/10.1371/journal.pone.0218747.t001

Wood sampling and sample selection for wet lab analyses

Sampling of wood shavings for all 1,000 trees was performed when trees were 5-year-old. Bark was removed from both sides of the tree, and wood shavings were sampled through the entire stem using a driller (12 mm diameter) at breast high (1.3 m), always in a north-south direction (S1 Fig). Wood shavings were stored in paper envelopes dried at room temperature and ground using a Willey mill. The 40–60 mesh (0.297–0.420 mm) woodmeal portion of the sample was used in all subsequent wet chemistry lab and NIRS analyses (S2 Fig). To select a robust set of samples for wet lab chemical and physical trait measurements and for subsequent NIRS model calibration, spectra of woodmeal for all 1,000 samples were obtained using a NIRSystems 5000 equipment (FOSS, Hillerød, Denmark), reading every second wavelength, from 1,100 to 2,500 nm. Each sample was read 16 times, using the average of each one of 700 wavelengths. In order to reduce the cost and time needed for the wet lab chemical and physical wood analysis, a representative subset of 350 samples was selected constituting a NIRS calibration/validation set (Fig 1). Samples were selected by the Kennard and Stone sampling algorithm [50], based on Euclidean distances of samples spectra aiming at maximizing sample variation and representativeness of the range variation for predicted chemical variation.

Wood chemical analyses

Cell wall chemical composition of wood samples was measured using a subsample of the same woodmeal samples previously scanned with NIRS. Cell wall carbohydrate and total lignin were measured as follows: woodmeal sample (40–60 mesh) was extracted with acetone in Soxhlet apparatus for 12 hours. Cell wall carbohydrates, namely cellulose (CEL) and hemicellulose (HC), and lignin traits, acid-soluble lignin (SL) and acid-insoluble lignin (IL), in combination forming the total lignin (TL) were determined as described [51] using the Klason method, with small modifications. Cell wall carbohydrates were quantified with a high-performance liquid chromatography system (HPLC) using a Dionex (DX-600, Sunnyvalle, CA, USA) equipped with a PA1 (Dionex) column, detector with a gold electrode and SpectraAS3500 auto injector (Spectra-Physiscs, Santa Clara, CA, USA). Carbohydrates amounts were quantified relative to monomeric cell wall-associated carbohydrates (glucose, xylose, mannose, galactose, rhamnose and arabinose) as standards [35]. The amounts of Klason lignin and cell wall sugars were obtained in percentages, relative to the initial weight of dry wood sample analyzed. The ratio of lignin monomer subunits (syringil/guayacil) (S/G), was determined from acetone-extracted ground wood and analyzed via gas chromatography (GC) on a Hewlett Packard 5890 series II equipment (Agilent Technologies, Santa Clara, CA, USA), equipped with an autosampler, splitless injector, flame ionizing detector and a 30-m 5% diphenyl-95% dimethyl polysiloxane-coated RTX-5MS capillary column (inner diameter, 0.25mm) [52]. Lignin chemistry traits were measured on the 350 samples, while cellulose and hemicelluloses were ultimately measured on a smaller set of 200 samples due to resource limitations (Fig 1).

Wood physical analyses

Increment wood cores were collected for the 350 selected trees based on NIRS spectra variation. Increment cores (12 mm) were collected at breast height (1.3 m), in a north-south direction (S1 Fig). The northern half of the core was used for wood density (WD) and microfibril angle (MFA) analysis and the southern half for fiber length (FL), fiber width (FW) and coarseness (COA). To measure wood density the northern half of the wood cores were precision cut in 1.67 mm-thick sections, using a custom-built twin-blade pneumatic saw [35], and acetone extracted in Soxhlet apparatus for 12 hours. The wood sections were acclimated to 7% moisture content and then scanned by X-ray densitometry (QTRS-01X; Quintek Measurement Systems Inc., Knoxville, TN, USA), from pith to bark. The measurements across the section were averaged to determine the sample density. With the purpose of establishing a regression for X-ray densitometry, ten samples were randomly selected and had precisely recorded their weight over volume and then scanned in the equipment, to estimate the remaining samples of the phenotyping calibration set. Unlike subtropical and temperate tree species, tropical eucalypts have a continuous growth pattern, and annual growth rings are unclear. Thus, MFA was measured 10 mm adjacent to the bark, instead of choosing a certain year, following procedures described earlier [35]. Briefly, precision-cut samples, used for wood density determination, were also used for MFA in a Bruker D8 Discover (Bruker AXS Inc., Madison, WI, USA) X-ray diffraction instrument equipped with an area detector (GADDS) to collect diffraction patterns, which contain reflection information of the microfibril orientation in the wood sample. The southern part of the increment cores were used for fiber traits analysis (FL, FW, COA). Samples were incubated in Franklin solution (30% hydrogen peroxide and glacial acetic acid; 1:1 ratio) at 70°C for 48 hours and macerated. Afterward, samples were washed in deionized water until the samples had been neutralized. The samples were filtered and oven-dried at 105°C. A part of the sample had weight recorded and resuspended in distilled, deionized water and analyzed on a Fiber Quality Analyzer (FQA; Optest Equipment Inc., Hawkesbury, Ontario, Canada). Fiber length and width were taken as the average of all fibers measured and coarseness measured by the dry fiber mass per unit length (mg.100m^-1). Wood physical traits were ultimately measured on between 337 (MFA) and 350 (WD) samples as some samples were lost in the procedure.

NIRS model calibration

Calibrations of the NIRS models were attempted for all chemical and physical traits based on spectra outputs of the 40–60 mesh (0.297–0.420 mm) woodmeal isolates of wood shavings and their corresponding wet chemical and physical measurements. The calibration population, composed by trees with both spectra and wet-lab data, was randomly divided into two subsets, one for model estimation and one for model validation. The estimation/validation sample sizes were either N = 250/N = 100 for those traits for which 350 samples were measured, or N = 150/N = 50 for those traits for which only 200 samples could be measured (Fig 1). The Unscrambler software (v.9.0; CAMO A/S, Oslo, Norway) was used to estimate the model parameters for each wood trait using Partial Least Squares (PLS) analysis for phenotype prediction, based on each sample spectra. For each trait, different spectral transformations were tested in order to obtain the highest possible accuracies. Root mean squared error of prediction (RMSEP; the difference between the true and estimated compositional value expressed in units of the phenotype), bias (the average value of the difference between predicted and measured values), and the coefficient of determination of prediction (R²_p) following transformations, were calculated for the external validation set and used to compare model estimates. The phenotypes of the remaining 650 samples of the population were predicted only for those traits for which NIRS models showed satisfactory performance. In the case of traits for which NIRS models were deemed poor (see below), only data for the directly measured trees were used in quantitative parameter estimation.

Genotypic data

Genomic DNA was extracted from xylem scrapings isolated from the sapwood at breast height (1.3m) using an optimized Sorbitol+CTAB method for high quality DNA from wood samples [53], quantified with a Nanodrop 2000 spectrometer and adjusted to concentrations between 20 and 40ng.uL^-1. DNA samples were genotyped at Geneseek (Lincoln, NE) using an Infinium (Illumina) custom made chip for Eucalyptus (EucHIP60k.br) [54]. DNA samples were also genotyped using DArT-seq, a sequence-based genotyping method developed by Diversity Arrays Technology Pty Ltd (DArT P/L, Canberra, Australia) [46]. From the sampled population of 1,000 trees, 970 were ultimately genotyped with both marker types, co-dominant SNPs and dominant (presence/absence) DArT-seq markers, and used for the subsequent quantitative genetics analyses.

Quantitative genetics analyses

The single-trait analysis was based on the following univariate individual-tree mixed model with additive and dominance genetic effects: [1] Where, y is the vector of phenotypic data, β is the vector of fixed effects (block design effects); p is the vector of random plot effects following where I is the identity matrix and is the plot variance; a is the vector of random additive genetic effects following , where A is the average numerator relationship matrix and is the additive genetic variance; d is the vector of random dominance effects following where D is the average dominance relationship matrix and is the dominance genetic variance; and e is the vector of the random residual effect following where is the residual error variance. X, Z_p Z_a, and Z_d are incidence matrices relating fixed and random effects to measurements in vector y.

In the marker-based approach, the pedigree-based relationship matrices for additive A and dominance D effects of the previous mixed model [1] were substituted by the corresponding marker-based additive (G_A; based either on SNPs or DArT-seq data) and dominance (G_D; based on SNPs data only) relationship matrices. The additive genomic relationship matrix (G_A) from the co-dominant SNPs was calculated using the function ‘A.mat’, in R (http://www.R-project.org/) package rrBLUP, that uses the formula proposed earlier [55]. where, W is the n × m (n = number of individuals, m = number of SNPs) rescaled genotype matrix following M—P, where M is the genotype matrix containing genotypes coded as 0, 1, and 2 according to the number of alternative alleles, and P is a vector of twice the allelic frequency, p_i. The additive genomic relationship matrix (G_A) from the dominant DArT-seq markers was calculated following the formula proposed by [56] where, S is a n × m matrix (n = 970, m = number of DArTseq markers) rescaled genotype matrix following Z—P, where Z is the genotype matrix that specifies the genotypes expressed as 1/0 denoting the presence/absence of the DArT-seq marker, and P is a matrix containing the allelic frequency of the code 1 at locus i, p_i.

We examined the effect of progressively reducing the number of SNP or DArT-seq markers on the estimate of relatedness coefficients and narrow-sense heritability using an additive-only model, i.e. model in Eq [1] without the dominance component. An additive-only model was used for comparing the two types of markers because DArT-seq are dominant and, as such, do not allow building a bona-fide dominance (G_D) relationship matrix. To that end, subsets of 500 (05K), 1,000 (1K), 3,000 (3K), 5,000 (5K), 10,000 (10K), 20,000 (20K), and 30,000 (30K) randomly selected SNPs, and subsets of 500 (05K), 1,000 (1K), 3,000 (3K), 5,000 (5K), 10,000 (10K), and 20,000 (20K) randomly selected DArT-seq markers were used to build the corresponding genomic additive relationship matrices. Then, a product-moment correlation coefficient was used to evaluate the connection between pairs of marker-based additive relationship matrices using each combination of marker (SNPs and DArT-seq) and the pedigree-based relationship matrix.

The dominance pedigree-based (D) and SNPs-based (G_D) relationship matrices were calculated in R (http://www.R-project.org/) using the package “AGHmatrix [57]. Two different parameterization approaches [58,59] were employed to build the dominance relationship matrices G_D, herein named G_DVitezica and G_DSu, respectively.

Additive genetic correlations between two different traits measured from the same individual were estimated based on the following individual-tree mixed model: [2] where, y_i and y_j are, respectively, the vectors of individual-tree observations on trait i and j. Matrices X_i⨁X_j, and relate observations to elements of , plot effects in , and additive genetic effects in , respectively, and is the error vector. The symbols ⨁ and ' indicate the direct sum of matrices and transpose operation, respectively. Finally, the expectation and variance-covariance matrix for plot effects in model [2] are respectively equal to: where, , and are the plot variances for the traits i and j, respectively. The expected value and variance-covariance matrix of the additive genetic effects in model [2] are respectively equal to: where, and are the additive genetic variances for the traits i and j, respectively, whereas is the additive covariance between traits i and j. The symbol ⊗ indicates the Kronecker products of matrices. The expected value and variance-covariance matrix of e are respectively equal to:

The residual variances for the traits i and j are , and , respectively, and is the residual covariance between the two traits.

Estimation of genetic parameters

Restricted maximum likelihood (REML) [60] was used to estimate variances and covariances for the random effects in the mixed models [1] and [2], and were obtained with the ASREML program [61]. Estimates of pedigree- and marker-based variances for the plot, additive, and dominance effects, and residual errors, i.e. , respectively, were re-parameterized to additive genetic correlations (r), and individual trait narrow- and broad-sense heritability ( and , respectively) as follows:

Additionally, pairwise trait Pearson’s correlations (and its significance expressed in probability levels) were calculated using phenotypes, and breeding values estimated using pedigree (pedigree-based numerator relationship matrix), SNP-based (33,398 SNPs) and DArT-seq-based (24,001 markers) genomic relationship matrices, using the cor() function in R and the function corr.test() in R (http://www.R-project.org/) package “physch”. Finally, Spearman correlations were also calculated to compare whether the ranking of predicted breeding values differed whether estimated using pedigrees and the two marker types used.

Phenotypic data and NIRS models

Extensive phenotypic variation was observed for all measured traits in the breeding population sample (Table 2). Range variation for mean annual increment in growth ranged from 29 to 122 m³.ha^-1.year^-1, total lignin from 24 to 32%, cellulose content from 42 to 55%, and S/G ratio from 1.8 to 4.2. Predictions with NIRS data once calibrated using the subset of 350 trees were satisfactory for part of traits. Overall, considerably better predictions were obtained for wood chemical than physical traits (S1 Table). For S/G ratio, for example, the coefficient of determination of prediction following external validation was R_p² = 0.86, while for hemicellulose it was only R_p² = 0.36. For WD density, the NIRS model had a poorer prediction (R_p² = 0.60), but was still considered effective for predicting values and therefore increase sample size for the subsequent quantitative analyses. For all other traits (CEL, HC, MFA, FL, FW and COA) only between 200 and 339 lab-measured samples were used in subsequent quantitative genetic parameters analyses (Fig 1). Trees for which direct wood trait measurements were obtained and used in quantitative parameter estimation, involved between 40 and 46 parents, 39 and 45 families with an average of 7.7±3.7 trees per family, a fully representative sample from the genetics standpoint.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Range variation for the 15 phenotypic traits assessed in the Eucalyptus urograndis hybrid population.

Number of trees for which trait values were ultimately used in the quantitative parameters analyses (n), and statistics: mean, median, standard deviation (SD), phenotypic coefficient of variation (CV), minimum (Min.), and maximum (Max.) values observed.

https://doi.org/10.1371/journal.pone.0218747.t002

Genomic relationships

A total of 33,398 SNPs with call rate ≥ 0.90 and MAF (minor allele frequency) ≥ 0.01, and 24,001 DArT-seq dominant (presence/absence) variants at a call rate ≥ 0.80 and estimated MAF ≥ 0.02 were used for the genetic analyses providing good genome-wide coverage of the 11 Eucalyptus chromosomes (S2 Table). To capture the realized relationship structure in the population sample, individual pairwise relatedness were estimated using either genome-wide marker data or pedigrees of unrelated (0.00 relatedness), half-sibs (0.25) and full-sibs (0.50). For the 970 trees ultimately genotyped with both marker types, we gathered 469,965 pairwise estimates of pairwise relationships. The value distribution involved 87.2% estimates of unrelated individuals, while half-sibs represented 10.3% and full-sibs 2.5%. Genomic realized relationships estimated using SNPs or DArT-seq markers showed bell shaped distributions with a better approximation centered at zero for SNPs and a left skewed distribution toward negative values for DArTseq markers (Fig 2). The average genomic relationships for unrelated, half-sibs and full-sibs were always lower than the expected values based on pedigree information. They were estimated at -0.026, 0.125 and 0.34 using SNPs, and -0.029, 0.145 and 0.376 using DArTseq markers, and essentially the same average relationship was apparent once >1,000 markers were used (S3 Table). A comparison of the pedigree expected and genomic realized relationship matrices is also depicted using heatmaps showing similar patterns, indicating good pedigree control in the production of the full-sib families, although marker data, by capturing the realized genetic relationships, provides considerably more refined estimates of the continuous distribution of true relatedness in the population sample (Fig 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Pedigree and genomic relationships.

Distribution of the number of pairwise additive relationships (excluding the diagonal elements) (left) and heat maps of the pairwise relationship matrices (right) among the 970 trees of the Eucalyptus population, estimated using the expected pedigrees, 33,398 SNPs and 24,001 DArT-seq markers (top to bottom). The heat map scales show the continuum of the realized genetic relationships between pairs of individuals, from no relationship (dark blue areas corresponding to values below and up to zero), increasing to half-sib relationships (light blue shades around 0.25) up to full-sib relationships (yellow areas corresponding to values around 0.5).

https://doi.org/10.1371/journal.pone.0218747.g002

The median and interquartile ranges of genomic relationships in the box plots spanned consistently lower values than the pedigree expected relationships for half-sibs and full-sibs, irrespective of the type and number of DNA marker used to estimate relationships, while the whiskers revealed a considerable number of relationships outside the expected range (S3 Fig). Pearson correlations between the A matrix and the different GRM (i.e., G) matrices built with variable numbers and types of DNA markers varied from 0.640 to 0.774, showing a very small variation with decreasing numbers of markers (S4 Table). Correlations between GRMs built with variable numbers of SNPs or DArTseq were generally high, above 0.9. A random set of 1,000 SNPs or DArT-seq markers recovered essentially the same GRM built with 33,000 SNPs or 24,000 DArT-seq (r = 0.95; S4 Table). Finally, the correlations of GRM built with SNPs and DArT-seq markers, varied from 0.742 (between the smaller marker subsets) to 0.907 (between GDART10.5K and GSNP30K). Slightly higher correlations were observed for the GRMs built with the genome positioned 10,500 DArTseq when compared to using all 24,000 DArTseq markers.

Heritabilities and genetic correlations

Heritability estimates for lignin traits and wood density were higher than for growth traits with all three approaches (i.e., A, G_DVitezica and G_DSu), while for the remaining wood traits heritabilities showed variation according to the estimation approach. Growth traits showed higher broad sense () than narrow sense () heritabilities with all three estimation approaches (A matrix and the two G_D approaches, G_DVitezica and G_DSu), revealing a substantial effect of non-additive genetic control on these traits (Table 3). For all wood traits, on the other hand, narrow and broad-sense heritabilities were essentially the same, except for FL where dominance also seemed relevant. Additive variance is therefore the main driver of genetic variation in chemical and physical traits, in contrast to what was observed for growth. While for lignin traits heritabilities using pedigree or genomic data were relatively similar, for CEL and HC approximately 2-fold larger heritabilities were estimated using genomic data with both approaches. For FW narrow-sense heritability estimates of zero to 0.08 were obtained, while broad-sense varied between 0.08 and 0.33 suggesting an implausible picture of only dominance variance controlling this trait, a result that will require further scrutiny. Estimates of narrow-sense heritabilities obtained from the realized genomic relationships using varying numbers and types of DNA markers were compared to the estimates obtained with pedigree data under an additive-only model (S5 Table). For growth traits and wood density the estimates obtained using genomic data were again smaller than those obtained from pedigree-based estimates, while for wood chemical and fiber traits they varied. Heritabilities estimates using the two types of markers were essentially the same, and increased with increasing numbers of markers, only stabilizing when more than 10,000 markers were employed, indicating that sparser genome coverage below this density does impact the accuracy of parameter estimation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Narrow-sense () and broad-sense () heritabilities and their approximate standard error (SE), and the ratio of dominance to additive variance for each growth, chemical and physical wood trait.

Heritabilities were estimated using the additive relationship matrix based on the pedigree (A) and genomic relationship matrix using for the dominance component the parametrization of Vitezica et al. (2013, G_DVitezica) and Su et al. (2012, G_DSu) constructed from all available SNPs (~ 33K).

https://doi.org/10.1371/journal.pone.0218747.t003

Pairwise Pearson genotypic (from breeding values) and phenotypic correlations were estimated amongst the fifteen traits using a univariate analysis (Table 4). Genotypic correlations estimated with the two types of markers were equivalent, and generally slightly higher when compared to estimates based on pedigrees relationships and higher than phenotypic correlations. However, this varied depending on which pairwise traits were compared. Genotypic correlations, as expected, were high amongst the traits within each category of correlated traits, i.e., growth, lignin, cellulose and fiber. Overall, significant but low correlations, positive or negative, were apparent between growth and wood traits. A slightly higher negative correlation was estimated between growth traits and WD. Higher negative or positive correlations were also seen between the different wood traits, including CEL × TL (-0.65), WD × S/G (-0.45), WD × CEL (0.32), WD × HC (-0.28), WD × FL (0.44), WD × COA (0.55), MFA × TL (-0.39), FL × CEL (0.35). Additive genetic and phenotypic correlations among the 15 traits were also estimated (S6 Table). Worth mentioning are the consistently high negative additive correlations of CEL with TL and IL, WD with MAI, DBH and S/G, and the high positive correlations of WD with CEL, FL and COA, of MFA with HC, and of FW and DBH, irrespective of whether a SNP, DArT-seq or pedigree matrix was used. Finally, high Spearman correlations (0.71 to 0.99) were observed among all three pairwise comparisons of breeding values estimated using pedigrees, SNPs and DArT-seq markers (S7 Table). However, correlations between breeding values derived from observed phenotypes with those obtained from pedigrees or marker data were slightly lower for growth and wood physical traits, consistent with the observed differences between narrow and broad-sense heritabilities.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Pearson correlations between different growth, chemical and physical wood traits from the phenotype and breeding values of the univariate analysis of the Eucalyptus grandis × E. urophylla breeding population.

In each cell from top to bottom: genotypic correlation based on SNP-based realized relationship matrix (~33K); genotypic correlation based on DArTseq-based realized relationship matrix (~24K); genotypic correlation based on pedigree-based relationship matrix; phenotypic correlations.

https://doi.org/10.1371/journal.pone.0218747.t004

NIRS phenotyping and range variation of wood traits

In addition to confirming the outstanding volume growth of urograndis hybrids, a key finding of our report is the significant phenotypic variation observed for wood chemical and physical traits. Trees with excellent growth rate, wood density in the range of 550–600 kg. m^-3 and lignin content and S/G ratio, as low as 24% and above 4.0 respectively (Table 2), were observed in a relatively small set of 1,000 individuals sampled in 45 full-sib families from 46 unrelated parents. The combination of high growth rate, with average density and high S/G ratio wood is what pulpwood tree breeders usually target. These properties ultimately result in high mean annual cellulose increment per hectare, because wood with higher S/G tends to be easier to delignify during chemical pulping [62], consuming less chemical and energy and producing higher pulp yield.

Our study emphasizes that a better use of the existing phenotypic variation in breeding populations depends on efficient ways to collect data for wood property traits for large numbers of samples. Direct wet chemistry methods for thousands of tested trees in a progeny trial are usually not an option, in spite of faster analytical methods for some traits. In most Eucalyptus breeding programs, only wood density and sometimes pulp yield are measured in a subsample of trees of a trial that passed the initial cutoff of volume growth assessment. Additional wood traits are then measured on a very limited set of trees at the final clonal testing stage, precluding a wider sampling of the available variation for wood properties during progeny trials. Near-infrared reflectance spectroscopy (NIRS) applied to wood properties in Eucalyptus is changing this scenario [22,24,30,63,64], allowing the prediction of wood traits by non-destructive sampling of very large numbers of samples. In our study, a good calibration model was obtained for S/G ratio (R_p² = 0.86) and reasonable models for the other lignin traits with R_p² ≥ 0.70 (S1 Table). Very good prediction models for S/G ratio have recurrently been reported for Eucalyptus with R_p² up to 0.97, a particularly valuable result to breeding for cellulose, as S/G ratio and Kraft pulp yield have been shown to be strongly positively correlated both in temperate and tropical Eucalyptus species [21,30,64–66]. A NIRS calibration below the threshold for a ‘‘good” model [17] was obtained for WD (R_p² = 0.60), possibly impacted by the fact that NIRS spectra were collected on woodmeal samples and not on solid samples. Model for WD was, however, deemed adequate for the purpose of this study including the estimation of quantitative parameters and tree ranking, as both these data applications involve relative traits values. Better calibrations for WD in eucalypts were developed when NIRS readings were obtained on more elaborate samples of radial wood surfaces [67], although good models could also be developed with woodmeal material when larger sample sizes were used [66]. For cellulose and physical fiber traits, our NIRS models were not as strong as those reported in previous studies [24,30]. As a result, only the actual 200 to 339 directly wet-lab measured trees were used for subsequent quantitative genetics analyses. These sample sizes and their source in terms of the number of families and parents involved nevertheless constituted a representative sample of the population for quantitative parameter estimation. Recently, global NIRS calibrations for wood chemistry in Eucalyptus were reported using the same NIRS instrument as the one we employed in our study [30]. Such a development opens opportunities for a follow-up study by simply inputting the spectra data collected in our work into those models, to potentially improve our estimates and include individual carbohydrate traits not yet contemplated here.

Genetic control of growth and wood traits

Wood properties, particularly wood density and chemical traits, showed a moderate to stronger genetic control than growth traits (Table 3), in agreement with previous reports [14,15] and a number of studies in different Eucalyptus species [22,28,63]. On the other hand, FL, MFA and COA showed a moderate to low genetic control, consistent with studies in E. globulus and E. urophylla [27,63], while FW showed a poor heritability. The use of additive-dominant models and the inclusion of genome-based data in our work further highlights the significant role of non-additive sources of variance in the control of growth in this Eucalyptus hybrid, confirming earlier studies that showed either a balance or a higher weight of dominant relative to additive effects [7,15,38]. Furthermore, with the exception of height growth, the ratio of non-additive to additive variance became higher when a genomic relationship matrix was used (Table 3). Interestingly, the key role of dominance variation in growth seems to be a peculiarity of this hybrid, as dominance variance was largely insignificant for growth in all other eucalypts where this was investigated, including E. urohylla × E. tereticornis hybrids [29], E. globulus [20], E. nitens [41] and E. globulus × E. nitens hybrids [68]. These results help elucidate the long-lasting success in achieving significant volume gains by capturing the non-additive portion of the genetic variation into clonally propagated selected eucalypt urograndis trees [69]. For all wood properties, the dominance variance was either unimportant for chemical traits and WD ( ≈ 0), or only slightly relevant for FL and COA with a ratio around 0.2 to 0.6 when a genomic matrix was used. Equivalent results were recently reported for growth, wood density and pulp yield when using a genomic relationship matrix in an additive-dominant model in a similar urograndis hybrid population [15].

It is important to point out that in estimating genetic parameters, we used an individual-tree mixed model (i.e., progeny model) with additive and dominance genetic effects instead of a parent mixed model with male and female effects and their interaction. The individual-tree mixed model used here assumes that the alleles controlling the traits would be common in the parental lines studied, that the variance genetic components in the segregating eucalypt hybrid population used are similar, and that epistasis is negligible. Identifying which genetic model provides the best description of a forest tree hybrid is a challenge [70], as the assumptions of the infinitesimal model may not be fully appropriate [71]. Studies in eucalypt hybrids have estimated genetic variance components using either an individual-tree mixed model [72,73], a parent mixed model (e.g., [11,29,74], or both [38]. When comparing parental and progeny mixed models involving pedigree- and marker-based information in a urograndis hybrid population with 13 female and 9 male parents, Bouvet et al. (2016) showed that progeny models with genome-wide information improve the variance estimates and the prediction of genetic values, and that a larger numbers of parents could improve model performance even more. Our study included a relatively large number parents (N = 46) and the progeny additive-dominance model using marker-based data reduced the overestimation of the additive variance observed when using pedigrees and consequently the narrow-sense heritabilities, mainly for DBH, volume and MAI (S3 Table). Particularly for such traits where the non-additive component is important, the pedigree-based analysis cannot capture the full genetic relationship among the individuals, failing in disentangling the non-additive genetic component from the additive one [37]. Genome-based relationships therefore not only provide considerably more realistic estimates of narrow-sense heritability, but also more accurate genetic variance decomposition into additive and non-additive factors. Of relevance is the fact that the use of a GRM allows such variance decomposition even if breeding is carried out by incomplete pedigrees [41], an easy, economical alternative that now becomes an appealing option for advancing eucalypt populations in light of the availability of low-cost genotyping.

Genetic correlations

Reported additive correlations between growth and wood traits have varied significantly across different Eucalyptus species, germplasm sources, and sites such that it is difficult to find general patterns for the relationships among these traits. Correlations from SNP-based relationships indicate that selecting for increased growth in urograndis will result in very modest negative genetic responses in wood density (-0.26), cellulose content (-0.12) and fiber length (-0.18), and positive for soluble and insoluble lignin and S/G ratio (0.15) (Table 4). These findings are consistent with earlier reports showing low negative correlations of growth with WD and pulp yield [14]. Although statistically significant, likely due to the large sample size used, these low correlations indicate that simultaneous selection in any desired direction for these wood properties can be applied with minor effects on the genetic gain for higher growth rate. Reports on correlations between growth and wood density have been somewhat conflicting, most of them showing negative correlations in E. globulus [20,63] and E. urophylla [27] consistent with our study, while others detected either positive correlation in E. urophylla × E. tereticornis hybrids [29], or no correlation in E. globulus [22,75]. Taken together, these results suggest that growth and wood density in eucalypts might be in fact largely independent, with their correlation impacted by species, germplasm source and environment. We also found low genotypic correlations between volume growth and chemical and fiber traits (positive with lignin, S/G and hemicellulose; negative with cellulose and fiber dimensions), again indicating opportunities for flexible concurrent selection for high growth and wood properties. Results from other studies either agree or contrast ours, confirming that no general patterns seem to exist when trying to establish any standard correlation between growth and wood traits in eucalypts [23,28]. On the other hand, correlations between wood chemical traits and density have been more consistent across eucalypt studies, likely due to the evident physical content interdependency of cellulose and lignin in cell wall structure and the impact of cellulose content on wood density. Our results agree with a general pattern of moderate to strong negative correlations between cellulose and lignin contents, density and S/G ratio, and positive between density with cellulose content and fiber traits [22,27,29,63].

Implications for Eucalyptus urograndis breeding

Urograndis hybrids have been and will continue to be the mainstay of the vast majority of eucalypt forest based industrial operations in the tropics, with an increasingly outstanding role in the world supply of sustainable forest-based products. By studying a typical urograndis breeding population, our report confirms some known patterns of variation for growth traits and brings novel quantitative genetics data for an array of wood properties. While increased volume growth is a common objective to all breeding programs, the relative importance and direction in which wood traits are selected varies. For example, a specific wood density window of approximately 500–550 kg.m³ is the target for trees for maximal efficiency in cellulose production, but higher density is aimed for energy (specifically, calorific value). While low lignin is usually desired to improve pulp yield in cellulose production, high lignin content to enhance the calorific power of wood is pursued by industries that use eucalypt for charcoal or biomass. Our analysis showed that ample genetic variation is available for all traits in a typical urograndis hybrid breeding population. However, while additive variance represents essentially the exclusive source of genetic variation for all wood property traits, dominance variance plays at least an equally important role as the additive variance when growth traits are considered. Additionally, our results indicate a general pattern of low genetic correlations between growth and wood traits, indicating good potential for simultaneous genetic improvement of individual trees for higher volume growth and wood properties in any desired direction. Wood chemical traits, however, do show significant genetic correlations, and given their important influence on final production traits, they call for unique breeding tactics to manage the connected responses to selection.

Our results have useful consequences for the current and future prospects of breeding urograndis hybrids. Reciprocal recurrent selection (RRS) between E. grandis and E. urophylla was tentatively adopted by a number of industrial breeding programs in the 1980’s and 1990’s, based on models derived from hybrid breeding of annual crops. It proved, however, to be time-consuming, costly and inefficient. Recently, simple recurrent selection (SRS) in synthetic hybrid population has been adopted in urograndis programs due to its simplicity, higher speed and lower cost, with genetic gain per year projected to far exceed that of RRS [69,76,77], although little empirical data is yet available to back its expected performance. The mostly additive control of wood properties described in our study would warrant using SRS, but the relevance of dominance would question its full efficiency for growth. However, as pointed out earlier [69], while non-additive effects would apparently be disregarded when choosing SRS, they would be efficiently captured nonetheless during selection of elite individual trees followed by their deployment by large-scale clonal propagation.

Finally, our study adds to the recent stream of promising reports showing the positive impact of using genomic realized relationships to estimate heritability, genetic correlations and breeding values in essentially all mainstream forest trees [45]. By effectively tracking the Mendelian sampling term, the use of genomic data gets closer to the true relatedness among individuals within and across families, allowing better genetic variance decomposition when compared to using the expected documented pedigrees. Both SNP and dominantly inherited DArT-seq markers efficiently recovered relationships; however, co-dominant SNPs were superior in providing more accurate estimates of heritability and allowing the estimation of dominance variance. Of note is the fact that the DArT-seq genotyping method currently also provides co-dominant SNP data. Although this study was carried out on fully pedigreed families, genomic relationships from genome-wide SNPs would allow advancing breeding populations by partial pedigrees followed by sib-ship reconstruction with the same or better precision of fully pedigreed populations with the added advantage of close control of inbreeding and improved management of genetic diversity. Costs of SNP genotyping have fallen drastically in recent years making this a viable option for eucalypts where controlled crosses are still relatively cumbersome, expensive, and prone to errors. Such flexibility in operational breeding practice and the prospects of accelerating breeding cycles by genomic selection for multiple traits at ultra-early ages, convincingly point to a new era of genomically-informed eucalypt tree improvement.