Tree breeding and silviculture_ Douglas-fir volume gains with minimal wood quality loss under variable planting densities

Validating performance of genetically-selected trees under realistic planting scenarios is essential for confidence in tree breeding programs. Quantifying the relative impact of genetic selection and initial planting density on tree size and quality can further guide operational forest practices. We evaluate volume gains, survival, stem quality and wood quality traits on 20-year old trees representing three levels of genetic selection that were grown under four initial planting densities. Working in a realized gain trial for coastal Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco) on five replicated sites, we ask: (1) Do predicted stand productivity levels materialize as expected under different planting densities? (2) Are single-tree plot designs, simulating progeny trials, capable of producing reliable results relative to large-block designs, simulating realistic planting scenarios? (3) If trees selected for volume gain show declines in wood and stem quality relative to wild-stand controls, can this be effectively managed by altering stocking densities? Because young progeny trials are used to estimate genetic gain in tree volume at a rotation age of 60, we use a growth and yield model calibrated to wildstand controls to assess whether genetically-selected families meet projections at age 20. On average, observed stand volumes exceeded projections on four out of five sites, and in three out of four initial planting densities. At the level of site by planting density, the moderate genetic-gain population (mid-gain) exceeded projections 13 out of 20 times while the top genetic-gain population (top-crosses) exceeded projections 11 out of 20 times. This fits expectations for breeding values, which are designed to reflect general performance averaged across all environments. Using a different validation approach, large-block designs showed better performance relative to simulated progeny trial designs. Very high planting densities (1890+ stems/ha) may minimize wood quality losses of genetically-selected planting stock but effects are relatively minor, while good performance among all traits was observed at operational planting densities (~1189 stems/ha). Wood density and microfibril angle proxy measures in top-crosses showed relatively minor and non-significant losses (−1.1 to −4.0%) compared to major and significant gains in volume per hectare at age 20 (29.0%) when averaging values at 1189 and 1890 stems/ha. Altogether, the genetic selection systems produce reliable results across a range of site qualities.


Introduction
Selecting the most appropriate seed sources for planting after harvesting is one of the most important decisions a forester can make to maintain long-term forest health and productivity. In addition to species selection and climate suitability for a planting site, foresters must consider the genetic quality of planting material. Planting seedlings derived from breeding and seed orchard programs (collectively, tree improvement) can help return the landscape to a forested state more efficiently and reliably (Stoehr et al., 2004), reduce rotation times (Rehfeldt et al., 1991;Stoehr et al., 2004), increase potential for economic gain when planted extensively across the landscape (Schreiber and Thomas, 2017;Chang et al., 2019), and increase operational flexibility. In addition to genetic quality, initial planting density is another key silvicultural consideration: It can influence site establishment by affecting early growth and survival (e.g., Scott et al., 1998;Woodruff et al., 2002), stocking at free-growing, and long-term growth and yield (Smith and Reukema, 1986;Reukema and Smith, 1987). It also impacts timber value by altering branch thickness, crown recession and juvenile to mature wood ratios (Jozsa and Middleton, 1994), which directly affect the proportion of log volume that can be recovered for lumber as well as lumber grade or quality (Fahey et al., 1991).
In Canada's western-most province of British Columbia (BC), tree breeding programs aim to increase the potential economic value of seedlings used in reforestation programs. This is done using recurrent selection that begins with seed and typically scions collected from phenotypes with desirable traits in natural stands. Many tree breeding programs emphasize volume gain and do so by selecting for height at a relatively young age in progeny trials and applying age-age correlations and multipliers to estimate timber gain at rotation age (Lambeth, 1980;Lambeth et al., 1983;Magnussen, 1988;Rehfeldt et al., 1991;Vargas-Hernández and Adams, 1992;Xie and Yanchuk, 2003;Ye and Jayawickrama, 2012). The longest-running tree breeding program in BC is for coastal Douglas-fir (Pseudotsuga menziesii var. menziesii (Mirb.) Franco), which is now entering its fourth generation and is currently expected to produce gains in volume of 20 to 30% by a rotation age of 60 years. Operationally, trees have been planted at initial planting densities ranging from 600 to 1300 stems/ha (Arnott, 1986). Planting densities have varied over time, however, which likely affects lumber recovery and value given documented impacts on Douglas-fir growth, crown and branch attributes (Briggs et al., 2007;Hein et al., 2008;Newton et al., 2012;Lowell et al., 2014;Grace et al., 2015).
Experiments at the intersection of silviculture and genetics are rare yet exceptionally informative. Realized-gain trials offer an ideal experimental situation to examine how genetic gains interplay with initial planting density. They are designed to test the performance of genetically-selected planting stock relative to wild-stand seed sources and several examples are well-documented in the literature (see St Clair, 1993;Carson et al., 1999;Lambeth and Springs, 2000;St. Clair et al., 2004;Vergara et al., 2007;Weng et al., 2008;Stoehr et al., 2010;Ye et al., 2010). Although realized-gain trials evaluate past selections, and breeding programs must necessarily continue making selections after trial establishment, they provide an important means to assess the accuracy of estimates on past selections. Unlike small plots in progeny trials, realized-gain trials test trees of similar gain levels in large block designs. Large blocks reduce the potential inflation of growth differences among gain levels, thereby simulating more realistic operational planting scenarios. Because these trials are spatially explicit, often incorporating different initial planting densities into the design, they enable assessments of the gains estimated from genetic selection in single-tree-plot progeny trials on an area basis. With increasing age, they also facilitate an evaluation of wood quality attributes under different stand density regimes.
Realized gain trials are therefore well-suited for answering two main questions about tree improvement programs. The first asks if gains estimated from progeny testing materialize as expected on an area-yield basis. Because gains are estimated from small plots in progeny trials, where inter-tree competition among plots or over-sampling of microenvironments may amplify differences among trees and families, genetic gains could be over-estimated (Magnussen, 1989(Magnussen, , 1993St. Clair and Adams, 1991;Carson et al., 1999). Inter-genotypic competition may, however, be less of a concern if final selections occur at young ages before significant competition occurs (St. Clair and Adams, 1991). A second question is whether genetic selection for volume gain reduces log value due to negative genetic correlations between growth and wood quality attributes (Zobel and Jett, 1995;Beaulieu et al., 2006;Baltunis et al., 2007;Jayawickrama et al., 2009Jayawickrama et al., , 2011El-Kassaby et al., 2011;Lenz et al., 2013). Wood density and microfibril angle properties are less desirable in juvenile wood in trees (Jozsa and Middleton, 1994;Barnett and Bonham, 2004;Baltunis et al., 2007;Donaldson, 2008). Faster-growing trees and shorter rotation times may therefore be associated with a higher proportion of juvenile wood (Jozsa and Middleton, 1994;Kennedy, 1995;Beaulieu et al., 2006;Ukrainetz et al., 2008;Jayawickrama et al., 2009). Thus, there can be compromises between genetically-based volume gains and wood quality.
Since realized-gain trials are valuable long-term experiments, wood quality measurements must rely on non-destructive methods to avoid damaging trees. A suitable proxy measure for evaluating wood density is provided by the Resistograph (Rinn et al., 1996;Gao et al., 2017). The Resistograph is a specialized drill that measures the resistance to torque, every 0.1 mm of drill depth, as a relative, unitless measure. Because it is non-destructive, fast and inexpensive, it is becoming more common in forest genetic field trials (Isik and Li, 2003;Gao et al., 2017;Fundova et al., 2018). A surrogate metric for assessing stiffness is the velocity of a sound wave measured longitudinally on standing trees. This method involves inserting two probes into a tree, hitting one probe with a hammer, and detecting the resulting sound wave with a sensor on the second probe. The longer the transit time (time of flight), the higher the inferred microfibril angle. Acoustic velocity methodologies assess standing tree wood quality and can predict wood quality attributes in final lumber (Wang et al., 2007;Jones and Emms, 2010). Because they are non-destructive, these well-tested methods have been useful in tree breeding programs and long-term genetic studies (Isik and Li, 2003;Knowles et al., 2004;Jayawickrama et al., 2011;Lenz et al., 2013;Moore et al., 2015). They are therefore also well-suited to making relative comparisons among genetic classes in realized gain trials.
We quantify the influence of genetic selection and initial planting density on tree and stand volume as well as survival, crown ratio, branch thickness and wood quality proxies. We work within a 20-year old realized-gain trial for coastal Douglas-fir, based on selections from first-generation tests. This trial represents trees from three levels of genetic selection grown under four different initial planting densities, replicated twice on each of five sites. Each site additionally contains two replicates of single-tree plots of the same genetic populations, simulating progeny trials. This design therefore enables us to assess the reliability of selection for volume via progeny trials. We aim to answer the following specific questions: (1) Do predicted stand productivity levels materialize as expected under different planting densities? (2) Are single-tree plot designs, simulating progeny trials, capable of producing reliable results relative to large-block designs, simulating realworld planting scenarios? (3) If trees selected for volume gain show declines in wood and stem quality relative to wild-stand controls, can this be effectively managed by altering stocking densities? Because approximately 15 million coastal Douglas-fir seedlings are planted annually in BC, of which almost all are sourced from seed orchards containing selected parents, there is potential for foresters to increase productivity without sacrificing wood quality.

Experimental design
The realized-gain trials include trees representing three genetic classes for volume that were grown under four initial planting densities on five replicated sites (Fig. 1). The three genetic gain levels include: (1) Wind-pollinated wild-stand (WS) seed sources represent 0% gain in tree volume. Seven local and operationally-collected seed sources were mixed at the seed level and serve as a control representing baseline planting stock used for reforestation in the absence of tree improvement.
(2) Fully-pedigreed mid-gain (MG) crosses from 24 parent trees of intermediate genetic quality. This class represents an average breeding value of 10% gain in tree volume at rotation age of 60 years, relative to wild-stand controls. (3) Fully-pedigreed top crosses (TC) of elite genetic quality represent controlled pollinations among 10 parents. This class is expected to attain 18% additional volume gain by rotation age of 60 years. According to Woods (1993), these classes were developed by selecting from 560 parents in first-generation tests. The base breeding populations represented selections from natural stands in the breeding zone as well as from western Washington, chosen based on form and size relative to other trees in even-aged stands. Other considerations included broad coverage of the breeding zone and minimum distances between trees to avoid co-ancestry (Woods 1993). These gains represent average parental breeding values based on 12-year height and diameter measurements from previous progeny trial series established using a non-overlapping half-diallel mating design, where parents were selected for broad adaptation, as they were tested across 11 environments (Woods, 1993;Stoehr et al., 2010). The mid-gain and top-cross gain levels are average mid-parent values of the created full-sib families and express the anticipated volume gain at a rotation age of 60, above wild-stand controls.
The five sites were chosen to represent a range of site qualities for productivity. These range from low to high site index or SI (Supplemental Table S1). Site index is a measure of potential site productivity, or the capacity to grow trees on a given area of land. Here, we use the accepted definition for site index of average height of the largest diameter at breast height (DBH, 1.3 m) undamaged tree in a 0.01 ha plot after 50 years of growth above breast height (British Columbia Ministry of Forests Lands Natural Resource Operations and Rural Development, 2018). The sites range in SI from 43.65 (Norrish Creek) to 36.14 (Spirit Lake). All sites are characterized by mild temperatures, high annual precipitation but relatively low summer rainfall (Supplemental Table S1). All sites were regularly brushed to reduce competition from non-test trees, and some sites were fenced. One-yearold container-grown seedlings were planted on the sites in multiple-tree plots under four initial planting densities. These planting densities include: 1.6 m × 1.6 m (equivalent to 3906 stems/ha), 2.3 × 2.3 m (equivalent to 1890 stems/ha), 2.9 m × 2.9 m (equivalent to 1189 stems/ha) and 4.0 m × 4.0 m (equivalent to 625 stems/ha). The 1189 stems/ha approximate an operational planting density and the 1890 stems/ha would be a very high planting density. The 3906 stems/ha and 625 stems/ha test extreme scenarios that anchor the results for modelling purposes. Each site contains two replicates per treatment combination of spacing by genetic level.
The multiple-tree plots are represented by 144 trees, planted in 12 rows by 12 columns. These large plots were separated by two rows of buffer trees where initial planting density changes from plot to plot. To provide additional buffer, only the inner 100 trees (10 rows by 10 columns) were included in the analysis (Fig. 1). Thus, the main experimental design included three genetic entry levels, four spacing treatments, five sites, two replications/site, and 100 trees/plot, for approximately 12,000 trees. The number of living trees with measured values varies by plot and are reported in Supplemental Table S2. In addition to the larger, multiple-tree plots, all sites contain two replicates of single-tree plots to allow for the evaluation of potential inflation due to design effects and to provide a direct comparison between the area-based plots and the single-tree plot designs used for progeny testing. Single-tree plots were established at 2.9 m × 2.9 m spacing in two 90-tree blocks that included 30 randomly sub-selected trees representing each genetic level (Fig. 1). The map on the right shows the distribution of Douglas-fir (Pseudotsuga menziesii) in the Pacific Northwest and site locations in coastal British Columbia. A hypothetical site layout for one planting site is shown on left. The experimental design relies on 3 genetic entry levels × 4 spacing treatments × 5 sites × 2 replications/site × 144 trees/plot. The green colour of the plots illustrates the genetic level of the treatment. Analyses were conducted on the inner 100 trees in each plot, providing a buffer of one row of trees, while buffer trees were also installed between spacing densities. Two replicates of single-tree plots (STPs) were installed per site and are highlighted in blue. The map was produced using ArcMap (ESRI ArcGIS version 10.6). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Tree size and stem quality measurements
In 2015, when trees were 20 years old, tree size and stem quality measurements were collected. Tree DBH was measured with a steel diameter tape. Tree height as well as height to the base of the live crown were measured with a Vertex IV hypsometer (Haglöf, Sweden), where distances were shot a minimum of two times. Where possible, shorter trees were measured using a pole graduated in centimeters for greater speed and accuracy. Tree height and DBH were used to calculate volume per tree according to equations published in Omule et al. (1987). Volume per hectare was calculated by summing the live tree volume in a plot, dividing it by plot area and multiplying it by 10,000 (more detail in Supplementary Information file).
Height to the base of the live crown is a measure of crown recession, which is an important characteristic affecting branch development that in turn affects log quality through knots (Briggs et al., 2007). For this study, the live crown base was defined as the lowest live branch that is part of the uninterrupted live crown. Crown ratio was calculated for each tree as a percentage of live crown to total tree height. Crown ratio is an indicator of tree relative photosynthetic capacity, stem quality and potential juvenile-to-mature wood ratios (Jozsa and Middleton, 1994). The average diameter of the largest branch in a whorl at DBH (or directly above) was recorded in mm using calipers. This is a proxy for overall branch thickness, which is a stem quality parameter affecting knot size and value of lumber (Fahey et al., 1991;Briggs et al., 2007). A binary response of mortality of this largest branch was further recorded since branch mortality can also affect stem quality through loose knots (Jozsa and Middleton, 1994).

Wood quality proxy field measurements
A subset of plots was measured with two acoustic hammer methods: The IML Hammer (IML GmbH, D-69168 Wiesloch, Germany) and the Director ST300 (fibre-gen, 7140 SW Fir Loop, Tisgard, Oregon, USA). ST300 data were collected on one treatment combination replicate (spacing by genetic level) on three sites; readings were averaged on a minimum of two measurements per tree. IML Hammer data were also collected on one block per treatment combination on two sites.
We used the Resistograph manufactured by IML (IML GmbH, D-69168 Wiesloch, Germany) on a subset of plots on three sites on the 100 inner trees per plot (10 columns × 10 rows) in 2017. Bark was not removed prior to resistance drilling, as it is not considered necessary and may damage trees (Gao et al., 2017). Drilling resistance readings for each tree include local maxima and minima that often reflect annual growth rings and other density variations within the tree. To obtain one value per tree, we obtained an average across the diameter of the tree. This differs from the recommendations to use bark-to-pith measurements (Gao et al., 2017). This may be because resistance values can sometimes increase slightly over drilling distance due to increasing friction from accumulation of sawn material (Gao et al., 2017). In the field, the tip of the drill bit was cleaned after every three to four trees, to avoid accumulation of resin or sawdust, and to ensure zero values until contact with the wood. In the measurements, a mid-point cut-off was not always obvious in our measurements. We therefore chose to use both radii and did no post-measurement manipulation of the data other than standard data quality control measures. The high sample size in this study may help overcome anomalies at the individual-tree level.
A consistent direction for drilling (e.g. North-South) was not considered necessary here. A recent multi-national publication further suggests that there are no systematic biases based on sampling direction (Gut et al., 2019). Measurements were obtained when air temperatures were above freezing and within several weeks rather than across an entire season, helping reduce variations in air temperature or wood moisture content. This may help minimize measurement error (Ukrainetz and O'Neill, 2010). We also did not measure volumetric wood density because the literature suggests correlations with drilling resistance are moderate to high (Isik and Li, 2003;Bouffier et al., 2008). For coastal Douglas-fir from the same tree breeding program in British Columbia, El-Kassaby et al. (2011) further found genetic correlations of 0.85 (±0.11) of Resistograph values with wood density values derived from x-ray densitometry. Here, we rely on this well-established proxy methodology to make relative comparisons among genetic classes.

Evaluating genetic gain estimates: Experimental design, stand age and growth and yield modeling
We perform three analyses to evaluate the reliability of genetic estimates. First, we compare single-tree plot values to multiple-tree plot values to gauge the potential for inflated results in tree-size and survival estimates. We compare the gain differences for genetically-selected seed sources relative to controls in multiple-tree plots and in single-tree plots. Second, we compare age 12 and age 20 data to evaluate if genetic gains diminish over time (Rehfeldt et al., 1991;Carson et al., 1999;Xie and Yanchuk, 2003;Stoehr et al., 2004). Third, we gauge whether genetically-selected sources, with total stem-volume gains estimated for age 60, are meeting expectations at age 20. We compare the observed values at age 20 to a growth and yield model calibrated according to the performance of the wild-stand controls on the five sites in this study. While realized-gain trials can be used to update genetic gain modules in growth and yield models, the purpose here is to retrospectively assess tree genetic selection systems. Since theory suggests that age 20 gains are expected to be larger than at a rotation age of 60, the growth and yield model provides a general indication as to whether gains are ontrack. It is recognized that such models are themselves considered to be not fully accurate as they operate with their own assumptions and simplifications, which are described in greater detail in the Supplementary Information file.
To calibrate the projected growth curves, we ran simulations with no genetic gain (genetic worth entry = 0%) as a base using the Tree and Stand Simulator (TASS 2.07.76). We compared these values at age 20 to the performance of the wild-stand control. The site index for the model simulations was adjusted until the projections matched the control based on plot means for total volume. This was done for all wild-stand control plots at each spacing treatment across all five sites. Volume-per-hectare projections for each plot were then produced incorporating 0%, 10% and 18% genetic worth into the simulation. The total stem volume per hectare in TASS was calculated based on all living trees, including volumes of stumps and tops to ensure consistency with the plot volumes calculated using the local volume equations (Omule et al., 1987). The runs were done on one-hectare size plots. Predictions at the plot level were then averaged at the site-by-spacing level and at the spacing level for reporting and visualization. More information about the growth and yield modeling is provided in the Supplementary Information file.

Statistical analyses
Several statistical analyses were performed to study the influence of genetics and silviculture on tree size, survival, wood quality, and stem quality. Analyses were primarily conducted using the ASReml-R version 4 package (Butler et al., 2018) in the R statistical programming environment (R Core Team, 2019). To quantify the relative influence of genetics and initial planting density as well as the influence of the experimental design, p-values were derived with an F-test on a fixed effects model (model 1) incorporating the pedigree. Residuals were carefully checked and outliers were removed. Variance components of design and interaction effects were subjected to a random effects model (model 2), averaging output at the site level. Acoustic velocity variables were fit with modified models, where replicate-within-site acts as a surrogate for site and the interaction term is removed. (1) where y represents the observed values for a particular phenotypic trait for a series of individuals; µ is the overall mean; g is a vector of fixed genetic entry level effects; d is a fixed vector of planting density treatment effects; g × d is a vector of fixed effects for the interaction of genetic entry level with planting density; s is a random vector of site effects, with s~MVN(0, σ 2 s I); b s is a vector of random replicate-withinsite effects, with b s~M VN(0, σ 2 bs I); s × d is a random vector representing the interaction between site and planting density, with s × d~MVN(0, σ 2 sd I 20 ); p s is a vector of random plot effects at each site, with p s~M VN(0, D ps ); r s is a random vector of row-within-plot effects (post-hoc blocking), with r s~M VN(0, D rs ); c s is a random vector of column-within-plot effects (post-hoc blocking), with c s~M VN(0, D cs ); a × s × d is a vector of random additive effects by site and planting density, with a × s × d~MVN(0, A ⊗ D s ⊗ I d ); and e represents a vector of errors, with e~MVN(0, D e ). The element 1 is a vector of ones, and the letters X and Z represent incidence matrices for their associated effects; A is the additive genetic relationship matrix derived from the pedigree; G is a 20 × 20 variance-covariance matrix between genotypes across sites and initial planting densities; I x is an identity matrix of its corresponding order; and ⊗ is the Kronecker product. ( 2) where y represents the observed values for a particular phenotypic trait for a series of individuals; µ is the overall mean; s is a random vector of site effects, with s~MVN(0, σ 2 s I s ); b s is a vector of random replicatewithin-site effects, with b s~M VN(0, σ 2b bs I bs ); d is a random vector of planting density treatment effects, with d~MVN(0, σ 2 d I d ); s × d is a random vector representing the interaction between site and planting density, with s × d~MVN(0, σ 2 sd I sd ); a is a vector of random additive effects, with a~MVN(0, A); a × s is a vector of random additive by site interaction effects, with a × s~MVN(0, A ⊗ D as ); a × d is a vector of random additive by planting density interaction effects, with a × d~MVN(0, A ⊗ D ad ); a × s × d is a vector of random additive by site and planting density interaction effects, with a × s × d~MVN(0, A ⊗ D asd ⊗ I d ); and e represents a vector of errors, with e~MVN(0, D e ). All matrix terms are as previously described.
Stem quality variables like branch thickness and crown ratio are influenced by stand development, e.g., after canopy closure, crowns begin to recede, which also alters branch properties (Briggs et al., 2007). Therefore, planting stock that has been genetically-selected for growth rate may progress through stand development stages earlier than wild-stand sources: Any potential difference in mid-gain or topcross classes relative to the wild-stand control may not be due to genetics for these variables, but due to indirect effects from faster growth. To determine the relative influence of genetics and stand development on these traits, we used a modified model 1 including height as a fixed covariate effect as well as its interaction with genetic level. Should the interaction between height and genetic level for each stem quality variable not be significant, it would suggest that stand development is playing a larger role than genetics for that given response variable.
In addition, best linear unbiased estimates (BLUEs) for fixed effects of genetic level and initial planting density were extracted from the simplified model: where y represents the observed values for a particular phenotypic trait for a series of individuals; µ is the overall mean; g is a vector of fixed genetic entry level effects; d is a fixed vector of planting density treatment effects; s is a fixed vector representing the site; s × g is a random vector representing the interaction between site and genetic entry, with s × g~MVN(0, σ 2 sg I sg ); s × d is a random vector representing the interaction between site and planting density, with s × d~MVN(0, σ 2 sd I sd ); p s is a vector of random plot effects at each site, with p s~M VN(0, σ 2 ps I ps ); and e represents a vector of errors, with e~MVN(0, D e ). All matrices were previously defined.
Models for variables calculated at the plot level (e.g. percent mortality) dropped the plot term, while acoustic velocity traits used a replicate-within-site term as a surrogate for site. Estimated BLUEs were later used for plotting, calculation of realized gains and comparisons of multiple and single-tree plots. Differences among genetic levels for all variables were tested with the least significant difference (LSD) at the 5% significance level and multiple comparisons were adjusted according to a false discovery rate (FDR).
Realized gains, or percent deviations of top-cross and mid-gain genetic levels relative to wild-stand controls, were calculated for all response variables according to: where x̅ is the mean of the genetic-level and spacing treatment combination; MG is the mid-gain genetic level; TC is the top-cross; and WS is the wild-stand control. To determine whether gains are inflated in single-tree plots (STPs) relative to multiple-tree plots (MTPs), we subtracted the STP gains from MTP gains.

Size and survival of selected trees higher relative to controls
Planting stock selected for faster volume gain outperformed wildstand seed sources in all tree size characteristics observed (Figs. 2 and  3). At age 20, volume per tree was the highest in top crosses planted under lower densities (Fig. 2). Average volume gain per tree was 29.3% (±4.3%) in the top crosses and 21.7% (±4.2%) for mid-gain trees (Fig. 3). Both genetic class and planting density were significant influences, but the interaction term was non-significant (Table 1); genetic gains were maintained across planting densities (Fig. 2) and sites (Fig. 4).
Genetic gains in tree volume are comprised of gains in height and diameter. Genetic effects were pronounced for height, with top-cross and mid-gain trees outperforming wild-stand controls by 13.1% (±1.2%) and 8.4% (±1.2%), respectively. Moderate increases in height were observed at the lowest planting densities (Fig. 2), but the effect of planting density on height was not significant (Table 1). The performance of genetically-selected seed sources was maintained across planting densities (Fig. 2) and the genetic group by planting density interaction term was not significant in the mixed model (Table 1). Average diameter gains for the top-cross and mid-gain genetic entries were 9.9% (±1.2%) and 7.8% (±1.2%), respectively. Initial planting density had a stronger effect on tree diameter than on height (Fig. 2). The planting density term was significant for diameter (Table 1). The genetic gains in diameter were maintained across the four planting densities (Fig. 2) and there was no significant interaction between genetic level and planting density (Table 1).
Survival was 92% on average, across the experiment. Relative to wild-stand controls, average survival of mid-gain and top-cross trees differed by approximately −2.2% (±5%) and 10.7% (±5.1%), respectively (Fig. 3), but these differences were not significant (Table 2). Survival, however, varied by initial planting density (Fig. 2). Survival was generally highest at 1189 stems/ha, where the top-cross class represented the highest survival of all treatment combinations at 96%. The lowest survival was observed under the highest planting density tested (3906 stems/ha). At this initial planting density, the wild-stand had the poorest survival rate at 85% while the mid-gain and top-cross genetic levels maintained higher survival rates, at 87% and 92%, respectively (Fig. 3).
Although volume/ha generally increased with higher planting densities (Fig. 2), the genetic gains relative to wild-stand controls varied (Fig. 3)

Selected trees show losses in wood quality traits under some initial planting densities
The Resistograph proxy shows a trend of lower wood density in genetically-selected sources relative to wild-stand controls (Fig. 3). The average difference in resistance was −3.5% (±2.1%) for top-cross trees and −2.7% (±2.1%) for mid-gain trees, although differences were not significant ( Table 2). The ST300 acoustic velocity values, representing inferred microfibril angle, indicated little to no loss (Fig. 3). The midgain and top-cross levels showed non-significant deviations of 0.4% (±4.9%) and −1.3% (±4.9%), respectively, relative to the wild-stand control (Table 2). However, the IML acoustic velocity method, another surrogate for inferred microfibril angle, suggested losses that were significant. The top-cross trees showed deviations from the control of −3.9% (±4.1%) and the mid-gain class showed deviations of −1.9% (±4.1%) ( Table 2).
While planting density was not significant for average drilling resistance, there was a significant interaction of planting density and genetic class (Table 1). Relative to the wild-stand control trees, genetically-selected trees showed variable losses in this wood density proxy under different planting densities. For example, this wood density proxy showed the largest deviations at the lowest planting density (625 stems/ha), where the top-cross trees deviated from the wild-stand control by −9.5% (±2.0%) and mid-gain trees deviated by −8.9% (±2.0%) (Fig. 3). The wood density proxy showed less pronounced deviations relative to wild-stand controls at 1189 and 1890 stems/ha. Specifically, mid-gain trees showed deviations in average resistance of −0.6% (±2.1%) at 1189 stems/ha and 0.5% (±2.1%) at 1890 stems/ ha relative to the wild-stand control. Relative to wild-stand controls, average resistance for the top-cross trees was 0.2% (±2.1%) at 1189 stems/ha and −2.4% (±2.1) at 1890 stems/ha (Fig. 3). Averaged across the 1189 and 1890 stems/ha planting density treatments, Resistograph deviations were −1.1% for the top-cross trees and −0.02% for the mid-gain trees.
The planting density term was significant for both acoustic velocity parameters, while the ST300 acoustic velocity parameter additionally showed a significant genetic class by planting density interaction ) and stem quality (panel C) as influenced by genetic selection and initial planting density in coastal Douglas-fir (Pseudotsuga menziesii var. menziesii) in a 20-year-old realized-gain trial in British Columbia. Colours represent the three genetic levels tested in this study, where the light green represents the wild-stand control (0% gain), a slight darker shade of green represents the mid-gain level (average expected 10% gain in volume at rotation age relative to the control) and the darkest green represents the top-cross level (average expected 18% gain in volume at rotation age relative to the control). Averages across all sites were derived from model 3. Error bars are standard errors of the predicted means. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (Table 1). For IML acoustic velocity readings, mid-gain and top-cross trees showed the lowest values relative to wild-stand controls at the 1189 stems/ha density, with deviations of −5.3% (±0.6%) and −5.9% (±0.6%), respectively. At the 1890 stems/ha density, however, the mid-gain showed deviations of −1.5% (±0.6%) and the top-cross showed deviations of −2.1% (±0.6%) (Fig. 3). Averaging values between 1189 and 1890 stems/ha suggested IML readings of −4.0% for the top-cross level and −3.4% for the mid-gain level. The 1890 stems/ ha planting density also showed increased values in the ST300 acoustic velocity readings, with the top-cross 0.3% (±0.7%) and the mid-gain 2.2% (±0.7%) higher than the wild-stand, respectively (Fig. 3). Averaging values between 1189 and 1890 stems/ha showed that ST300 deviated from wild-stand controls by 0.7% in the mid-gain level and −1.4% in the top-cross level.

Stem quality characteristics of selected trees affected by stand development
The top-cross and mid-gain genetic levels had higher heights to the base of the live crown and lower crown ratios than wild-stand controls (Figs. 2 and 3). On average, height to live crown for the top-cross was 32.7% (±4.1%) greater than the wild-stand control, and the mid-gain was 15.8% (±4.1%) greater. Similarly, crown ratio differed by −20.0% (±7.3%) for top cross and −7.8% (±7.3%) for mid gain when compared to wild-stand controls. The top-cross and mid-gain trees differed significantly from wild-stand trees for both variables (Table 2). Although both the recession of the crown and crown ratio are strongly influenced by initial planting density, the interaction terms were not significant and no apparent genotype-by-environment interaction was observed in the data (Table 1, Fig. 3). Branch thickness is primarily controlled by initial planting density (Table 1, Fig. 3) while differences among genetic levels were not significant (Table 2). Branch mortality was greater in the selected trees, significantly differing from the control by 54.6% (±20.7%) for the top-cross level and 15.3% (±20.7%) for the mid-gain level (Fig. 3). Differences between the topcross and mid-gain levels were not significant (Table 2). To determine whether the differences in stem quality traits in genetically-selected trees are due to genetic mechanisms or to faster stand development, a model incorporating height and an interaction term with height was applied. The interaction term of genetic level by height was significant for branch thickness but was not significant for crown ratio (Supplemental Table S3).

Performance on large-blocks relative to single-tree plots
The potential for inflated estimates due to the use of single-tree-plot experimental designs for progeny testing is illustrated in Fig. 5. Greater gains in the large-block designs are observed for height, diameter and volume per hectare in both mid-gain and top-cross levels at both age 12 and 20 (Fig. 5). At age 12, height gains in multiple-tree plots are 1.2% and 1.1% greater than in single-tree plots for mid-gain and top-cross levels, respectively. At age 20, gains in multiple-tree plots are 1.7% and 1.4% greater than in single-tree plots for the mid-gain level and topcross levels, respectively (Fig. 5). For diameter, the difference in gains between multiple-tree plots and single-tree plots depended on the genetic level. Gains in diameter for the mid-gain level are slightly greater in the multiple-tree plots relative to the single-tree plots, at 1.1% at age 12 and 1.0% at age 20 (Fig. 5). The gains in the top-cross diameters are 1.0 and 0.8% age 12 and 20, respectively, in multiple tree plots relative to single tree plots (Fig. 5). For survival, the mid-gain levels showed lower survival than wild-stand controls in multiple-tree plots relative to single-tree plots (Fig. 5).

Gains diminish with age as expected from the Lambeth model
The gains we found closely match those predicted by the Lambeth Fig. 3. Percent change of the mid-gain and top-cross genetic levels relative to the wild-stand control for tree size and survival, wood quality and stem quality traits in a 20-year realized-gain trial for coastal Douglas-fir (Pseudotsuga menziesii var. menziesii) in British Columbia. Colours represent the genetic levels in this study, where the darkest shade of green represents the top-cross level (average expected 18% gain in volume at rotation relative to the control) while the lighter shade of green represents the mid-gain level (average expected 10% gain in volume at rotation relative to the control). More information about response variables and proxy measurements are provided in the Methods section. Realized genetic gains were calculated according to models 3 and 4. Error bars are standard error of plot differences relative to wildstand controls. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) age-age correlation equation (Lambeth, 1980). We observed 29.3% gain in volume for the top-cross level relative to the wild-stand control at age 20, while the Lambeth equation predicted 29.6%. For the mid-gain level, the Lambeth equation prediction was 16.5% while 21.7% were observed. Since these gains are expected to diminish with time, however, we illustrate the age trend in Fig. 5. The advantage of the selected trees over the wild-stand control diminishes at age 20 relative to age 12 across all planting densities for height, diameter, tree volume and volume/ha (Fig. 5). The higher survival increased relative to the wildstand control due to mortality by age 20 (Fig. 5).

Deviations of gains from growth and yield model projections
Predicted average gains at rotation age met expectations most of the time when averaged by site or initial planting density (Table 3). By site, selected sources met or exceeded projections on four out of five sites for the top-cross level and on all sites for the mid-gain level (Table 3, Fig. 4). When averaged by initial planting density, gains underperform on the lowest planting density (625 stems/ha) for both genetic classes but meet expectations at higher planting densities (Table 3, Fig. 4). At the level of site by planting density, the mid-gain level exceeded predictions 13 out of 20 combinations (five sites by four planting densities) while the top-cross level exceeded predictions 11 out of 20 times (Table 3). At the operational planting density of 1189 stems/ha, the top-cross exceeded projections on two of five sites. Across the five sites, the top-cross deviated from growth and yield model projections by 0.3% on average but ranged from −12.6% to +20.0% (Table 3). The mid-gain exceeded projections on four out of five sites, averaging 8.5% and ranging from −0.2% to +23.4% (Table 3).

Implications of genetic gains in volume under different initial planting densities
The theoretical predictions of growth gains from genetically-selected seed derived from tree breeding and seed orchard programs are validated to age 20. The mid-gain and top-cross genetic classes consistently out-performed non-selected seed sources in terms of productivity and there was no genetic interaction with spacing, as was also found in other experiments for Douglas-fir (St. Clair and Adams, 1991;Johnson, 1997;St. Clair et al., 2004). At the tree level, slower individual tree growth to age 20 is observed at higher planting densities, which is consistent with the effect of increasing competition (Woodruff et al., 2002). On an area basis, however, greater yield was achieved with higher planting densities. This is congruent with realized gain trials for other conifers (Rehfeldt et al., 1991;Weng et al., 2008). Although 1890 and 3906 stems/ha result in greater volume gains per area, low individual tree diameter could reduce lumber recovery (Fahey et al., 1991). At very low planting densities (625 stems/ha), however, low volume per hectare is also combined with other undesirable characteristics in stem and wood quality, discussed further below.

Genetic gains congruent with modelled expectations at the landscapelevel
Our study found that relative genetic gains in tree size and areabased yield diminish with age. This builds on earlier work by Stoehr et al. (2010), who found realized height gains for the mid-gain and topcross classes to be 10.4% and 16.1%, respectively. These authors found volume gains per tree to be 28.4% for the mid-gain and 48.5% for the top-cross, while volume gains per hectare were 27.8% and 52.8%, respectively. However, the gains diminish as expected according to Lambeth model predictions. This aligns with research by Ye and Jayawickrama (2012), who studied 68 progeny trials for coastal Douglas-fir up to 41 years of age, and found good fit with the Lambeth equation. It is also consistent with findings by St. Clair et al. (2004), who found that gains materialized as expected in a coastal Douglas-fir realized-gain trial in the United States Pacific Northwest.
When comparing volume gains at age 20 to a growth and yield model calibrated by wild-stand sources, however, results depended on the scale of analysis. Counting occurrences of positive deviations at the level of site by spacing (20 combinations), the mid-gain and top cross met or exceeded growth and yield model projections 65% and 55% of the time, respectively. Averaged by site, however, the selected sources met or exceeded projections on four out of five sites. Averaged by initial planting density, only the lowest density (625 stems/ha) showed performance lower than projected by the growth and yield model, suggesting genetic gains may not materialize as predicted at these low planting densities or that the model used over-estimates gains at low planting densities. Low planting densities may be sub-optimal regardless, due to undesirable crown ratios and thicker branches. While genetic gains in volume per hectare materialize across a range of higher planting densities, they appear most congruent with projections at operational planting densities (1189 stems/ha) and very high operational planting densities (1890 stems/ha). Tree size and quality traits were measured at age 20 (in 2015) and wood quality traits were measured at age 22 (in 2017). Significance for fixed effects is based on model 1. Variance component proportion is expressed as a percent (VC %) and is based on model 2. Wood quality measurements use a concatenated replicate-withinsite variable as a surrogate for site.
The generally positive deviations relative to a growth and yield model fit expectations for breeding values, which are designed to reflect general performance averaged across all environments. We note that these are based on a calibrated but imperfect growth and yield model. Furthermore, it is uncertain how gains will materialize over the next 40 years to a rotation age of 60. The projections illustrate possible final outcomes resulting from different levels of genetic gain, planting density and site productivity based on the information to date. For this reason, realized gain trials are very valuable experiments worthy of long-term monitoring. Fig. 4. Observed values for volume per hectare of all living plot trees relative to calibrated growth and yield model projections for coastal Douglas-fir (Pseudotsuga menziesii var. menziesii) for three genetic classes: wild-stand (light green), mid-gain (medium green; representing 10% volume gain at rotation age) and top-cross (representing 18% volume gain at rotation). Panel A shows the observed volume per hectare values, representing straight means derived from age 12 and 20 measurements in a realized-gain trial in British Columbia. Panel B represents growth projection curves to a rotation age of 60 under the four planting densities tested here. Growth projections are estimates derived from the Tree and Stand Simulator (TASS), British Columbia's regenerated stand growth and yield model. More information about growth and yield modeling is provided in the Methods and Supplementary Information sections. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Genetic selection system: Simulated progeny trials provide reliable estimates
Gains in tree size and volume per hectare for genetically-selected material were greater in multiple-tree plots representing more operational planting scenarios, relative to single-tree plots, which test mixtures of families representing different gains levels (simulated progeny designs). However, the survival of trees from genetically-selected levels was generally higher in single-tree plots than in block plots. This may be because trees with poorer performance are more likely to be outcompeted when grown alongside top performers. On multiple-tree plots, where genetic classes are grown in pure stands, differences in survival among the genetic classes are not as great because the relative competitive advantage is lower. Regarding tree size variables, though, it was interesting to find that gains in multiple tree plots were generally comparable to single-tree plot values. We had expected that gains in the large-block designs would be smaller relative to single tree plot designs due to expectations of inflated estimates from single-tree plots. The data therefore suggest that offspring from parents selected in single-tree-plot progeny tests perform well as a group when planted together. Altogether, this suggests that the existing selection systems for coastal Douglas-fir, based on single-tree plots at age 12 on a 2.9 × 2.9 m spacing, produces reliable results, and that inter-tree competition up to at least age 20 does not introduce bias in breeding value estimates for stem volume.

Low genetic effect on crown and branch thickness: Best managed through silviculture
Stem quality traits like crown ratio and branch thickness were more influenced by initial planting density than by genetics. Higher planting densities were associated with smaller crown ratios, which may improve lumber recovery and quality (Fahey et al., 1991;Jozsa and Middleton, 1994;Schimleck et al., 2018). Although the mid-gain and top-cross trees appear to have faster crown recession relative to wildstand controls, this is likely due to advanced stand development resulting from faster-growing and larger trees: The lack of a significant interaction between genetic level and height suggests that differences in the genetic levels in crown ratio are due to different rates of stand development rather than differences among the genetic groups. However, the significant interaction of genetic level and height for branch thickness indicates a genetic influence. Although the mid-gain and top-cross trees do not show poorer stem quality characteristics in general, our results indicate that these traits remain more effectively managed through initial planting density.

Productivity-wood quality trade-offs: Avoid low initial planting densities
The large gains in productivity found in this study were associated with relatively minor reductions in wood-quality. At age 20, the topcross genetic group produced an additional 29.0% in volume/ha relative to the baseline scenario when averaged at 1189 and 1890 stems/ ha planting density levels. In contrast, wood-quality deviations ranged from 0.7 and −4.0% at the 1189 and 1890 stems/ha planting density levels, varying by trait and genetic-gain level. The IML acoustic velocity technique generally showed larger losses (−3.7% average but down to −6.0%) than the ST300 approach (−0.5% average but down to −3.2%). These traits provide a proxy for wood stiffness and reductions in stiffness down to −5% may have minimal impacts in construction (Jayawickrama et al., 2009). While reductions were observed here, it is worth noting that the parent trees chosen for the top-cross and mid-gain populations used in this study were selected solely based on expected stem-volume breeding values, with no consideration for wood quality. This contrasts with breeding and seed orchard populations that are selected on the basis of both growth and wood quality traits. Although previous selections may not have incorporated wood quality as a selection criterion, the evaluation of earlier selections here seems to indicate that the minor wood quality losses are manageable.
The lowest initial planting density tested here (625 stems/ha) produces the lowest wood quality, large and undesirable branches and crowns, and the lowest volume per hectare at age 20. These findings are congruent with the literature: Adverse crown, branch and wood quality attributes for coastal Douglas-fir are observed under low initial planting densities (Briggs et al., 2007;Hein et al., 2008;Newton et al., 2012;Lowell et al., 2014;Grace et al., 2015). The acoustic velocity values were generally lowest under lower planting densities, indicating higher microfibril angle and lower wood quality. This aligns with results found in other conifer species, where lower stocking densities or thinning treatments can sometimes lead to lower stiffness inferred from acoustic velocity measurements (Lasserre et al., 2005;Lowell et al., 2014;Moore et al., 2015). In the present study, we find these responses are generally amplified in the genetically-selected planting stock, as indicated by realized gains across different planting densities.
While higher planting densities could increase volume per hectare and simultaneously help avoid minor wood quality losses, it is important to evaluate economic feasibility. Higher stocking may incur higher planting costs but may be associated with reduced competition control and stand tending costs. The higher upfront costs are a reason why planting densities greater than 1200 stems/ha are currently less common in coastal BC. Anecdotally, however, planting densities were higher in the 1980s (1200 + stems/ha) when survival was lower due to poorer planting stock. When stock quality improved and survival increased, initial planting densities were as low as 800 stems/ha to reduce planting costs. The present study supports concerns over quality at lower planting densities. Future economic analyses could explore value and return on investment, but in the meantime, higher operational planting densities (1189 + stems/ha) seem suitable for all traits studied here for both mid-gain and top-cross genetic levels.

Conclusions
The genetic gains in volume that are realized under simulated operational planting scenarios met expectations based on estimates from prior progeny trials as well as from single-tree progeny trials simulated P-values adjusted for multiple comparisons according to a false discovery rate (Benjamini and Hochberg, 1995).
in this study. On average, they also met expectations according to a calibrated growth and yield model. These results support confidence in protocols used for genetic-selection programs and the growth gains predicted from small or single-tree-plot progeny trials. In the context of large volume gains, the genetically-selected trees used in this study were associated with relatively minor losses of wood quality even though the parent trees used in this study were selected for growth traits only, with no selection for wood-quality. This is in contrast to seed orchard populations that are selected for both growth and wood quality traits. Results from this study support the use of operational planting densities between about 1100 and 1800 trees per ha. Within this range, genetic gains in volume per ha are optimized, with minimal losses in wood quality.

Materials statement
Genetic stock is potentially available upon request yet subject to availability of controlled crosses among specific parents. Institutional Material Transfer Agreements that specify their use is for non-profit research work only is required.

Data statement
Datasets are available from the British Columbia Ministry of Forests, Lands, Natural Resource Operations and Rural Development upon request.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 5.
Deviations of genetic gains in tree size and survival of genetically-selected mid-gain (10% predicted gain in volume at rotation age) and top-cross (18% predicted gain in volume at rotation age) relative to wild-stand controls (0% gain) at age 12 and 20. Panel A, in blue, represents the comparison of multiple tree plots (MTP) to single tree plots (STP) in terms of deviations of genetic gains of the two genetically-selected groups relative to the wild-stand control. Single-tree plots were established at an initial spacing of 2.9 m × 2.9 m and were compared to multiple-tree plots of the same planting density. Panel B, in green, highlights deviations from wild-stand controls across four initial planting densities tested on large block plots only (multiple-tree plots or MTPs). These panels show the expected reduction in gains over time, across the four initial planting densities tested here. For both single tree and multiple-tree plots, averages are based on 2 replicates on 5 sites. Measurements were derived in 2006 and 2015, when trees were 12 and 20 years old. Trees were grown in British Columbia in a realized-gain trial for coastal Douglas-fir (Pseudotsuga menziesii var. menziesii). The realized gains were calculated according to models 3 and 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 3 Percent deviations of observed volume per hectare for top-cross (TC) and mid-gain (MG) genetic levels relative to projections from a growth and yield model that was calibrated to wild-stand controls. Values are reported for five sites representing different levels of stand productivity and four initial planting densities from a 20-year old realized-gain trial of coastal Douglas-fir (Pseudotsuga menziesii var. menziesii). An operational planting density is simulated by the 1189 stems/ha scenario, the average for which is highlighted in bold.