Meta-Modelling to Quantify Yields of White Spruce and Hybrid Spruce Provenances in the Canadian Boreal Forest

: Tree improvement programs can improve forest management by increasing timber yields in some areas, thereby facilitating conservation of other forest lands. In this study, we used a meta-analytic approach to quantify yields of alternative white ( Picea glauca (Moench) Voss) and hybrid spruce ( Picea engelmannii Parry ex Engelmann x Picea glauca (Moench) Voss) stocks across planting sites in the boreal and hemiboreal forests of Canada. We extracted meta-data from published tree improvement program results for five Canadian provinces covering 38 planting sites and 330 white and hybrid spruce provenances. Using these meta-data and a random-coefficients nonlinear mixed-effects modelling approach, we modelled average height over time trajectories for varying planting site characteristics, as well as climate transfer distances between planting sites and provenances. Climatic transfer distances had strong effects on the height trajectory parameters. In particular, the asymptote parameter had a nonlinear increasing trend with planting site versus provenance mean annual temperature differences. We incorporated the height trajectory meta-analysis model into an existing growth and yield model to predict volume yields. Overall, this research provides a mechanism to quantify yields of alternative provenances at a particular planting site, as a component of decision support models for evaluating evaluate forest management investment into improved planting stocks alternatives under current and possible future climates.


Introduction
Forest industries in Canada and elsewhere are facing increased competition for forest land from urban expansion and other resource sectors including agriculture, oil, and gas, resulting in a diminishing fiber supply [1,2]. At the same time, climate changes are expected to impact forests via changes in forest growth rates and in natural disturbances caused by fire, diseases, and insects [3]. Tree improvement programs can provide improved seedling stocks with faster growth rates and/or better disease and insect resistances [4][5][6][7][8]. Using improved stocks can ameliorate timber supply shortages and facilitate a more competitive forest sector, while continuing to provide a broad range of ecosystem services by supporting sustainable forest management over the entire forest land area [9][10][11]. However, a mechanism for quantifying the impacts of using improved stocks at different sites is needed to support this initiative.
The boreal forest covers an extensive expanse of the world's forests, circling the northern hemisphere of the earth across Scandinavian and Nordic countries, Russia, Canada, and USA [12]. The Canadian boreal forest represents 28% of this biome, covering 552 million hectares [13]. The hemiboreal forest is south of the boreal forest and is transitional with temperate forests, except in Alberta and British Columbia, where it is generally to the south and west of the boreal forest and is transitional to mountainous forests ( Figure 1). Coniferous species are generally dominant in the boreal forest, including several Picea species [14]. In Canada, white spruce (Picea glauca (Moench) Voss) and hybrid spruce (Picea engelmannii Parry ex Engelmann x Picea glauca (Moench) Voss) are important commercial species that occupy much of this spatial extent [15] (Figure 1).

Figure 1.
Boreal and hemiboreal forests of Canada and USA [15].
Canadian tree improvement programs began in the 1950s with the overall goal of providing improved stocks with superior traits [16]. Initially, seeds were collected from phenotypically superior individuals called "plus trees" from natural stands with unknown male parents; performances were then assessed using provenance or progeny trials representing multiple seed sources to establish seed planning zones. Subsequent selections of seeds through both wind (i.e., natural-stand mating) or controlled-pollinated crosses of these first-generation trees were used to establish progeny trials to assess the extent of genetic control over selection attributes, to rank parents and offspring, to establish production populations (seed orchards), and to quantify the expected genetic gain from improved stocks [10,16,17].
Quantifying gain from improved tree stocks is a challenge for a number of reasons. First, trees are long-lived and, as a result, the time to maturity (and subsequent generations) is long, unlike most agricultural crops. Even the earliest provenance trials in Canada represent only short periods of time relative to biological rotation ages reported for natural stands. Second, growth rates can increase with subsequent selections of best performing trees [18], and, therefore, it has been considered to be a "moving target" [19]. Third, there often is an environmental effect and there may also be a provenance by environment effect because provenance performances vary with planting site characteristics [20,21]. Finally, in quantifying improvements, a baseline for comparison must be used. This baseline could be wild stocks [22,23], seeds from the planting location (i.e., local provenance) [19,24,25], or the average of all provenances in the trials [26,27].
In spite of these difficulties, a number of researchers have quantified yields of improved stocks for some Canadian commercial tree species at a provincial or more localized spatial extent. Yields of white spruce provenances and their breeding performances were compared based on trials in Québec [28]. The 20-year growth increase of black spruce (Picea mariana (Mill.) B.S.P.), white spruce, tamarack (Larix laricina (Du Roil) K. Koch), and jack pine (Pinus banksiana Lamb.) was reported using provenance trials in New Brunswick [17]. A number of ways to quantify genetic gain was discussed for commercial species of British Columbia [23]. Growth and yield trends of first-generation jack pine and black spruce in New Brunswick were modelled using 20-year measurements [18]. Yields of improved white spruce were quantified as a means of determining the economic benefits of the tree improvement program in Québec [10].
Meta-modelling (also called meta-regression) can be used as a mechanism for extending the temporal and spatial extents of these studies by combining results of many provenance trials. Metamodelling is related to meta-analysis of experimental trials [29] and to using combined estimators to localize existing models [30], each of which provide more precise estimates by synthesizing data from different sources [31]. Meta-modelling refers to the combined use of aggregate data (termed macrodata in some papers) in a model previously reported in separate studies [32], but may also incorporate individual data (termed micro-data in some papers) [33].
The use of meta-modelling to combine information across provenance trials has been reported in a limited number of studies. Mean heights of different provenances of Larix sukaczewii Dylis, L. sibirica Ledeb. and L. gmelinii (Rupr.) Rupr. were predicted at four planting sites in Alberta using simple models and climate transfer variables (i.e., differences in climate between planting site and provenance locations) [34]. Yield responses for a number of conifer species in the Canadian boreal forest were modelled using reported results for several provenance trials from Canada and the United States [27]. Tree-level information from several lodgepole pine (Pinus contorta var. latifolia Douglas ex Louden) common garden experiments were combined in Washington State, USA and British Columbia, Canada [35] using a random-coefficients mixed-effects modelling approach to estimate the basal area growth of different provenances (i.e., seed source location) under current and forecasted future climates. A mixed-effects model and plot-level data from several plantations in the north-west USA were used to predict three-year height response of Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) [36]. Influence from climatic and competition variables on tree growth responses were examined by developing a nonlinear mixed effects model of black spruce in three forest types in the western boreal forest of Canada [37].
In this study, we used meta-modelling to quantify changes in yields by using stocks from white and hybrid spruce tree improvement programs across the boreal and hemiboreal forests of Canada. White spruce is considered a "plastic" species because of its ability to adapt to extreme climates [38]. The hybrid of white spruce and Engelmann spruce (Picea engelmannii Parry ex Engelm.) occurs naturally at mid-elevations [39] over a large extent of the hemiboreal forest in British Columbia [40] and in the foothills in Alberta [21]. While heights of white spruce can reach 30 m in natural stands, this hybrid can reach up to 60 m in height. To predict these yield changes, an average height trajectory meta-analysis model was developed using published results for white and hybrid spruce provenance trials across a large range of planting sites. A random-coefficients approach was used [35,36], where climatic transfer distances between provenance and planting site locations were used to modify the meta-analysis model parameters. The new average height trajectory meta-analysis model was then incorporated into an existing growth and yield model to forecast volume per ha, basal area per ha, and average diameter trajectories for different provenance and planting site climates. The developed model provides a mechanism to quantify yields of alternative provenances at a particular planting site. As a component of a decision support model, this research facilitates evaluating alternative forest management investment in using improved planting stocks alternatives under current climates, as well as under possible future climates [9].

Meta-Dataset
The meta-dataset used to develop the average height trajectory meta-analysis model included white spruce provenance trials from the Canadian boreal forest and hybrid spruce progeny trials from the hemiboreal forest of British Columbia. This meta-dataset was collected from published articles, internal reports, government reports, and other sources principally using the Web of Science ® (Science Citation Index expanded) search engine following an approach similar to that used in other research [27,41]. Only articles from well-known journals and research agencies were included in the final meta-dataset. Since a longer temporal extent would improve the accuracy and usefulness of the meta-model, provenance trials generally included open-pollinated and first-generation stocks. From this search, 24 studies were selected representing 38 planting sites in Canada with provenances from Canada and the United States (Table 1). For each selected study, the following variables were obtained: (1) species, (2) planting site label, (3) latitude and longitude of the planting site, (4) planting site elevation (m), (5) plantation start date and measurement dates, (6) initial planting density (stems per hectare), (7) plantation age (years), (8) average height of each provenance at each measurement time (m), and (9) latitude, longitude, and elevation of each provenance ( Table 2). The plantation start date, measurement dates, and densities varied across planting sites. In some cases, authors were contacted for more information, such as additional re-measurement data and unreported variables. For many of these studies, the location information for each provenance was identified from other sources. In cases where the planting density was not reported, this was calculated from the reported inter-tree spacing between all regularly-spaced trees in the plantation. Although repeated measures were preferable for our study, some provenances had only a single time measurement. The average height trajectories by provenance showed a relatively regular pattern ( Figure 2).   a MAT (°C) is the mean daily temperature; MAP (mm) is the mean annual precipitation; DD (days) is the mean degree days greater than 5 °C; and DMAT (°C) is mean daily temperature distance, DDdif (days) is the degree days distance, and DMAP (mm) is the mean annual precipitation distance between the site and the provenance.

Climate Data
For each provenance at every planting site, climate normals information was added to the metadataset. To obtain these data, the process followed was: 1. Weather station data for the 1981-2010 period were extracted from the Environment Canada climate normals database (www.climate.weatheroffice.ec.gc.ca, accessed 10th September, 2013) and from the National Oceanic and Atmospheric Administration's (NOAA) National Climatic Data Center (NCDC) (ftp://ftp.ncdc.noaa.gov/pub/data/normals/1981-2010, accessed 4th April, 2020). The USA data were then converted to equivalent SI units. 2. For each weather station, the following climate variables were calculated for each year and then averaged over the 1981-2010 period: (i) mean daily temperature (MAT); (ii) mean annual precipitation (MAP); and (iii) mean annual number of days with a temperature greater than 5 °C (i.e., "degree days", DD). 3. These weather station climate variables were then used to interpolate climate normals for each provenance location and planting site using inverse distance weighting (IDW).
Using the climate normals for each provenance and planting site, climate transfer distances were calculated for each provenance at a planting site as: 1. Mean daily temperature difference (DMAT) = site MAT-Provenance MAT 2. Mean annual precipitation difference (DMAP) = site MAP-Provenance MAP 3. Degree days difference (DDdif) = site DD-Provenance DD Summaries of the climate transfer distance variables for the meta-dataset are presented in Table  2. Of note, the DD for the provenances had a wide range. Since this is a threshold metric, the provenances from the USA sites had much higher values for DD as more days reached the threshold of 5 °C. This also resulted in a wide range for the DDdif.

Base Model and Parameter Prediction
For the base model representing the average height over plantation age trajectory for all provenances and sites, five models were considered based on recommendations in the literature, in particular, recommendations given by Burkhart and Tomé (pp. 114-123) [48]. The models considered were: 1. Chapman-Richards model (introduced into forestry by Pienaar and Turnbull [49]): where H is the average height of planting site i and provenance j at measurement time t; age is the plantation age; and is the error term. This model has three parameters: the asymptote ( 1 ) specifies the maximum height and the other two are shape parameters ( 2 and 3 ).

Schumacher model [50]
: This model has two parameters: the asymptote ( 1 ) specifies the maximum height and the other is a shape parameter ( 2 ). [51,52]:

Lundqvist-Korf function
This model has three parameters: the asymptote ( 1 ) specifies the maximum height and the other two are shape parameters ( 2 and 3 ). 4. Exponential model [48]: This model has three parameters: the asymptote ( 1 ) specifies the maximum height and the other two are shape parameters ( 2 and 3 ). 5. Monomolecular model [53]: This model has three parameters: the asymptote ( 1 ) specifies the maximum height and the other two are shape parameters ( 2 and 3 ).
We fitted each candidate model using the meta-data and used the commonly used Akaike's information criterion (AIC) [54], Pseudo R 2 , and root mean squared error (RMSE) to select a base model.

Pseudo 2 :
where = 1, … , planting sites; = 1, … , provenances within planting sites; = 1, … , measurement times for a provenance with a planting site; is the sum of squares of residuals; ( ) is the corrected total sum of squares; H ̅ is the average of all provenances average height and Ĥ is the predicted height.
2. Root mean squared error (RMSE, m): is the total number of observations and is the number of fixed-effects parameters in the model.
3. Akaike's information criterion (AIC): where L is the maximum likelihood reported for the fitted model.
We also superimposed the average trajectory over each of the trajectories in the meta-data. Since this base model was then modified by the parameter prediction methods described next, the results and justification for the selection of a base model are presented here.
The Chapman-Richards and Exponential models were the best of these five models, with similar Pseudo 2 (0.906 and 0.907, respectively), RMSE (1.091 m and 1.087 m, respectively), and AIC (6308 and 6294) values. The remaining three models also fit the meta-data well, with Pseudo 2 values from 0.85 to 0.884, RMSE values from 1.217 to 1.381, and AIC values from 6764 to 7294, but were not further considered. The Exponential model resulted in a predicted average height of more than 50 m at age 80 years, which is not biologically possible at this age for the forests of this study. Given the nature of the Chapman-Richards model, predicted average heights were constrained to a maximum of 30 m matching previously published silvics for white spruce. Finally, the fitted Chapman-Richards model function better reflected the measured height yield trajectories (Figure 2) than the Exponential model.
Given our results, we selected the Chapman-Richards model as the base model. This model is biologically based; it was developed by integrating Von Bertalanffy's model, which mimics organism growth rates as the result of metabolism rate differences [55]. Further, although model accuracy will change among datasets [53], the Chapman-Richards model has been commonly selected as the preferred model for yields of trees and stands (p. 116, [7]). For example, the Chapman-Richards model fitted the data better than other alternative models for predicting the dominant height of slash pine (Pinus elliottii Engelm. var. elliottii) [56]. Finally, the parameters of the Chapman-Richards model are easy to interpret, an important consideration for parameter prediction using a base model [57,58]. Wang et al. [59] used parameter prediction and the Chapman-Richards model for estimating dominant heights of Eucalyptus (Eucalyptus globulus Labill) in southeastern Australia. However, although the asymptotic nature of the Chapman-Richards model is considered a model strength, this may present a problem when forecasting yields of trees and stands that may decline over longer time periods. We return to this later in the discussion section of this paper.
Random-coefficients modelling (aka parameter prediction modelling) (pp. 54-56 of [60], pp. 403-530 of [61], pp. 317-342 of [62]) was then used to modify the base model for the meta-dataset using planting site and provenance climates and planting site characteristics. The parameters of the base model were replaced with functions of source (i.e., genetic population) and location (i.e., planting site) climate variables, along with other variables, following prior research on basal area increment [35] and on height growth response of Douglas-fir (Pseudotsuga menziesii (Mirb.) Franco) populations [36] that used climate transfer distances variables.
In particular, for this study, the Chapman-Richards model (Equation (1)) was modified to predict average height of each provenance at each site by replacing the parameters with functions of climate transfer distances and other variables: where H , age , and were previously defined in Equation (1); and 1 , 2 and 3 are error terms. Several random-coefficients models were fitted using combinations of predictor variables in submodels 2a to 2c. Possible predictor variables were climate transfer distance variables (DMAT, DMAP, DDdif), planting density, site elevation, and provenance elevation. Transformations of some variables were included to represent expected nonlinear trends with climatic transfer distances. In particular, the climate transfer model by Rehfeldt et al. [34] was used as a guide for the asymptote parameter sub-model Equation (2a). PROC MODEL of SAS software version 9.3 was used to fit the models [63]. To increase assurance that a global optimum was found for each model solution, several starting parameter sets were included and predicted average height-age trajectories were carefully compared to measured values.

Model Selection
As with the base model, the accuracy of alternative parameter prediction models was compared using the Pseudo R 2 , RMSE, and AIC; however, the numbers of estimated parameters used in the AIC calculations were much larger than the base model given each parameter of the base model was predicted using a linear model. As well as these fit statistics, graphs of predictions superimposed on the observed data for each provenance and planting site combination were used to check for any model lack-of-fit.

Model Validation
Model validation was applied only to the final model selected for the meta-dataset. First, the selected meta-analysis model was evaluated using the "leave-one-out" model validation approach. For this purpose, one planting site along with all provenances at that site was excluded from model fitting, and then the model was applied to the excluded planting site. This was repeated for each planting site and the results were summarized into the root mean square predicted error (RMSPE) and mean prediction error (MPE): Further validation of the selected meta-analysis model was done using the "0.632 bootstrap" validation approach [64,65]. For this, the model was fitted using 30 planting sites selected using simple random sampling with replacement and the fitted model was applied to the non-sampled sites (i.e., test set). Results were then summarized to obtain RMSPEc and error variance ( 2 ): where k is the number of observations in the test set and is the number of non-sampled planting sites. After repeating this process 1000 times, the average RMSPEs (RMSPE ̅̅̅̅̅̅̅̅̅ ) and average error variances were calculated ( ̅ 2 ). Next, the model was fitted using all meta-data and applied to the same meta-data; these results were also summarized into a root mean squared error (RMSPEM) and error variance ( 2 ). Then, the weighted averages of these (i.e., using the test datasets and using the full meta-dataset) were used to calculate the 0.632 bootstrap root mean square predicted error (RMSPE0.632) and error variance (Err0.632) as follows: Err0.632 = 0.368 2 + 0.632 ̅ 2 RMSPE0.632 = 0.368 RMSPE + 0.632 RMSPE ̅̅̅̅̅̅̅̅̅

Using the Height Trajectory Meta-Analysis Model to Forecast Yields
To examine the sensitivity of the three meta-analysis model parameters (θ 1 , θ 2 , and θ 3 , Equation (2a)-(2c)) to climatic transfer distances, estimated average height trajectories to plantation age 80 years were obtained using the developed meta-analysis model for eight simulated provenances, representing a wide spatial range (i.e., from northern and southern locations across the boreal and hemiboreal forests of Canada) planted at the Petewawa National Forestry Institute site in Ontario. This planting site was chosen from the meta-data since this was the oldest plantation (last measured at 44 years since planting) with the largest number of provenances (25 provenances). These were compared to those for the Chalk River provenance, the local provenance at this planting site. Two planting densities were also simulated (500 and 2500 stems per ha).
Since the motivation for developing this meta-analysis model was to forecast yields using improved planting stocks (i.e., best provenances or improved seedlots from progeny testing), the selected model was then incorporated into an existing white spruce plantation growth and yield model [66]. This model was chosen since it is a stand-level model, was developed for spruce in the boreal forest, and is a "height driven" model. The developed average height trajectory meta-analysis model replaced the average height model within this growth and yield model, but no other model components were changed (Figure 3, implemented in R, version 3.3.1 [67]). Because this growth and yield model is height-driven, replacing the average height model is expected to provide accurate yield predictions for volume per ha and other variables. However, no accuracy assessment was possible since information on these other variables was not commonly published for provenance trials.

Figure 3.
A previously developed white spruce stand-level growth and yield model [66] was modified for yields of different provenances and planting sites using the developed height trajectory model.

Average Height Trajectory Meta-Models
Using the meta-dataset, several combinations of climate transfer distance variables (DMAT, DMAP, DDdif), planting density, site elevation, and provenance elevation were used to model the parameters of the Chapman-Richards model (Equation (2) and (2a)-(2c)). Of these, the six top models are listed in Table 3. Model 6 had the highest Pseudo 2 (0.9234), smallest RMSE (0.9866 m), and smallest AIC (5894), and no lack-of-fit was noted on the graph of predicted average height trajectories superimposed on the observed data (see Figure A1 in Appendix A for example graphs, and location information in Table A1 and Figure A2). In addition, all parameters associated with the predictor variables were statistically significantly different from 0 ( = 0.05) ( Table 4). The RMSPE was 1.1 m and the MPE was 0.10 m using the leave-one-out model validation procedure. Using the 0.632 bootstrap model validation approach, RMSPE0.632 was 1.2 m and Err0.632 was 1.3 m 2 .   Table 3 replaces Equation (2a)-(2c)). For the selected meta-model, DMAT and DMAT 2 were included to model the nonlinear trend of the asymptote with DMAT. For provenance and site elevations of 500 m, the asymptote peaks at a DMAT of around −5 °C. However, within the range of DMATs within the meta-dataset (−6.83 to 8.66 °C; Table 2), the asymptote generally declined with an increase in DMAT, meaning that the asymptote is lower when a provenance from a cold location is planted at a warmer location. In terms of the two shape parameters, planting density was included in the parameter prediction model for 2 indicating that average height is affected by planting density. Although dominant height is commonly not affected across a broad range by planting density, average height may be affected, as later discussed.

Using the Height Trajectory Meta-Analysis Model to Forecast Yields
The predicted height trajectories for the eight simulated provenances over a wide spatial range showed a large variation, with a range of about 12 to 27 m at age 80 for both 500 and 2500 stems/ha planting densities (Figure 4). This was largely due to differences in DMAT corresponding with latitudinal differences. Provenances from northern areas of Alberta, BC, Ontario, and Yukon had lower predicted height trajectories (shown in grey on Figure 4) (DMATs from 5 to 8 °C). Southern provenances (shown in black on Figure 4) had the same or higher predicted height trajectories relative to the Chalk River provenance, the local provenance at the Petawawa National Forestry Institute. Similar results were obtained for both densities indicating that density has only a small effect on 2 .
Using the modified growth and yield model (Figure 3), the quadratic mean DBH (diameter at breast height) and volume per ha trajectories were predicted for the simulated provenances and two planting densities. Although the planting density had little impact on the average height trajectories, the other yield variables were responsive to planting density as would be expected (Figure 4). Most importantly, the predicted volume yields for the simulated provenances ranged from around 200 to 600 m 3 /ha at 80 years for the 2500 stems/ha planting density compared to the Chalk River local provenance with about 425 m 3 /ha, indicating the wide range of possible productivity.

Discussion
Tree improvement programs can provide options for ameliorating fiber shortages; however, forecasting possible changes in yields is challenging. In this study, a meta-modelling approach was used as a mechanism for providing yield forecasts by combining meta-data from white and hybrid spruce provenance trials across the wide expanse of the Canadian boreal forest. Height trajectories were modelled using the widely-used Chapman-Richards model (e.g., [7,68,69]) as the base model.
The derivation of the model was based on growth of organisms [55], translating the complex biological yield of organisms into a simplified mathematical form with only three parameters representing a wide range of curve shapes [49,70].
The model also reaches an asymptote. This feature has received some criticism in that stand yields may decline. For example, a weak relationship between dependent and predictor variables can cause the asymptote to occur at an unrealistically low age using the Chapman-Richards function [69,71]. However, the asymptotic nature likely represents average heights of high-value plantations. Average heights would only decline if: (1) regeneration/ingrowth trees were retained (i.e., no stand tending to remove non-crop trees) and included in the average height; (2) taller trees died at a higher rate than shorter trees; or (3) some trees were damaged and damaged heights were included in the average height. The first case is unlikely for higher value plantations; even if regeneration/ingrowth trees were present, often only crop-trees are included in reported average heights. The second case is also unlikely since little structural diversity is expected in plantations with controlled spacing and limited genetic stocks. In the final case, the plantation would be compromised due to catastrophic damage and would not likely be included in the meta-data. Therefore, the leveling off of average height to an asymptote is a reasonable conjecture for these plantations.
We used random-coefficients modelling and replaced the three parameters of the Chapman-Richards model with functions of climatic transfer distances between the provenance and planting site locations, since these variables can represent growth potential [20,34]. In particular, we found that DMAT was a very important variable, showing a nonlinear trend for the asymptote function (Equation (2a)). Other studies also found that differences in mean temperature between provenances and planting sites were important. For example, a simple quadratic surface was used to model mean height growth of different provenances and planting sites of Larix sukaczewii Dylis, Larix sibirica Ledeb., and Larix gmelinii (Rupr.) Rupr [34]. Of the large number of climate variables, mean annual temperature was chosen as the best single predictor variable, but other effective single predictor variables included mean temperature in the coldest month, the degree days less than 0 °C, annual moisture index, and the difference between warmest and coldest months mean temperatures. The authors initially concluded that the largest height growth occurred when the planting site and provenances had similar climates (i.e., climate transfer distances close to 0); however, when they excluded one planting site (Brooks) that had been irrigated, they found that the optimal height growth was obtained when provenances were moved from warmer to colder temperatures. In a study of lodgepole pine provenances, summer and winter temperatures, along with precipitation, were used in a random-coefficients model for basal area growth [35]. Basal area growth was larger when provenances came from a warm temperature and planted at a warmer and wetter site. For Douglasfir provenances in USA, the mean temperature of the coldest month transfer distance was important for modelling three-year height (i.e., seedling height at age three after planting) [36]. In that study, the final model indicated that provenances coming from colder climates had improved growth if planted at locations with warmer winter temperatures. Other studies have shown inconsistent results with climatic transfer distances. For example, the height growth of several species common in eastern North America decreased when provenances from colder areas were transferred to warmer sites for some species [72]. Conversely, other researchers [73] found that height growth decreased substantially when a provenance from a warmer temperature (4 °C more) was planted at a colder site. For our study, we found that the average height asymptote for white spruce and spruce hybrids was generally higher for provenances that were planted in a colder location (i.e., a negative DMAT), consistent with the Larix spp. study [33], but unlike the Douglas-fir study [35]. However, this asymptote declined for very large negative DMATs.
Planting density and elevations of the planting site and provenance were also included in the sub-model for the asymptote, but these variables had less impact. Planting density was also included in one of the two shape parameters sub-models of our study but, again, it had only a minimal effect on predicted average heights. These results were not unexpected since many authors have found that planting density has little effect on dominant height [56,60,74], and average and dominant heights are quite similar in plantations with similar genetic stocks. As further support, no statistically significant effects of planting density were found for white spruce average heights using Nelder plots [72]. However, even though planting density did not have a large impact on average height trajectories, there would be impacts on volume per ha and on quadratic mean DBH as reflected in the other sub-models of a growth and yield model.
Finally, to provide the desired yield estimates for differing provenances, the height model of an existing white spruce growth and yield simulator [64] was modified using our average height trajectory meta-model. Since models for basal area per ha, volume per ha, and stems per ha were linked with the height model, the trajectories of these other variables were altered indirectly. An alternative approach that is often used is to use a single percentage (or multiplier) at all ages to modify the growth and yield of regular stocks for changes in genetic stock. For example, Petrinovic et al. [10] modified previously published yield tables of another study [75] by increasing the dominant average height trajectory by 20%. Other researchers projected height, volume, and diameter using the STAMAN stand growth model [76] and obtained the percent gain at age 40 by comparing improved and unimproved stocks [16]. However, using a simple multiplier assumes that any changes in yield are anamorphic, whereas we allowed the shape parameters to also change for the height trajectory meta-analysis model, thus allowing for polymorphic changes over time. We directly modified the average height trajectory model in the growth and yield model only, which affects other yields (i.e., basal area, volume, and quadratic mean DBH) in turn. Genetic gain of radiata pine (Pinus radiata D. Don.) was incorporated into an existing growth and yield model by adjusting site productivity in New Zealand [77]. Another study modified both the mean top height and basal area per ha model components to predict yields of improved radiata pine [19]. In our research, we used a metamodelling approach using primarily published data from many provenance trials; unfortunately, information on variables other than average height was not commonly available. However, since models of stand-level variables are linked in growth and yield simulators to reflect biological reality [7], our change in the average height trajectory might be expected to produce realistic forecasts of the other yield variables.
The main contribution of this study is the use of a meta-modelling approach to predict yields using improved stocks. Our approach expanded the spatial and temporal extents of the individual studies, resulting in a larger range of climate transfer distances and plantation ages. However, as is often the case in other studies, the climate transfer distance variables were spatially interpolated using the closest weather stations. Further, the meta-data extended to only 44 years of age, whereas the harvest ages of white spruce plantations are more commonly about 80 to 100 years [23]. We used the Chapman-Richards model as the base model to constrain predicted heights to published biologically based limits beyond the meta-data age range. However, these limits should be validated as data become available for these plantations. In addition, we did not conduct uncertainty analyses to evaluate the risks associated with using the developed meta-model in any decision-making process in this paper. However, the results of this research were incorporated into a decision support model and used to evaluate alternative forest management investment strategies under possible future climates [9]. A final limitation of this study is that the average height trajectory meta-analysis model would be expected to underestimate successive generation yields [6].

Conclusions
Using meta-analysis, an average height trajectory meta-analysis model for white and hybrid spruce was developed to predict yields using improved stocks for the boreal and hemiboreal forests of Canada. The fitted model closely followed the observed patterns of average height yield over time. Further, the trajectories are polymorphic, allowing for the possibility that a provenance that performs well at an early stage may not, ultimately, be the best provenance at rotation. The average height yield trajectory model was incorporated in a growth and yield simulator to obtain the yield trajectories of other variables, notably volume per ha, as a means of estimating the potential of tree improvement programs for mitigating fiber shortages. However, caution should be used since rotations commonly used for white spruce are considerably longer than the 44 years represented in the meta-data.
Author Contributions: S.A. developed the meta-data, analyzed the data and prepared the original draft; S.A. and V.L. conceptualize the methodology; V.L reviewed and edited the manuscript; A.Y. contributed in providing more repeated measures in the meta-data; S.A., V.L. and A.R. selected methodologies for the model validation; S.A, V.L, A.Y., A.R., P.M, G.B. contributed to the manuscript editing. All authors have read and agreed to the published version of the manuscript.
Funding: Funds were received by G.B. through Genome Canada, Genome Québec, Genome British Columbia and Genome Alberta for the SMarTForests Project (www.smartforests.ca).

Acknowledgments:
We thank Barry Jaquish for sharing the BC spruce provenance trials. Also, thanks go to Jean Beaulieu for providing data for a further measurement of Québec spruce provenance trials. Finally, special thanks go to Yousry El-Kassaby for providing very helpful suggestions on an earlier draft of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.  Figure A1 for provenance locations; b MAT (°C) is the mean daily temperature; MAP (mm) is the mean annual precipitation; DD (days) is the mean degree days greater than 5 °C; and DMAT (°C) is mean daily temperature distance, DDdif (days) is the degree days distance, and DMAP (mm) is the mean annual precipitation distance between the site and the provenance.   Table A1).