Stand Dynamics, Humus Type and Water Balance Explain Aspen Long Term Productivity across Canada OPEN ACCESS

: This study examined the relative importance of soil, stand development and climate hypotheses in driving productivity for a species that is widely distributed in North America. Inventory plots, 3548 of such, either dominated by aspen or made up of species mixture of which aspen occurs in dominant canopy position were sampled along a longitudinal gradient from Quebec to British Columbia. Site index (SI), was used as a measure of productivity, and soil, climate and stand attributes were correlated with site index in order to determine their effects on productivity. Results show a decline in productivity with high moisture deficit. Soil humus correlates significantly with SI but does not sufficiently capture differential rates of litter deposition and decomposition effects over the long-term. Consequently, aspen composition, stand ageing, and stand structural changes dominate variability in productivity. Within the context where deciduous cover has being increasing, there are implications for forest productivity.


Introduction
The inherent nutritive potential of a soil plays a major role in the nutrition and growth of trees [1] hereafter described as the soil hypothesis.Tree growth is affected by biotic factors such as the activity of soil microbes [2], which is influenced by litter quality and therefore changes in stand structure and composition [3].Stand structure and compositional dynamics represent developmental changes that are driven by disturbances, tolerance, and competition [4][5][6].Consequently, areas with long fire cycles tend to develop complex stand structures and are characterized by higher species richness [7], thus stand development hypothesis.Also, climate directly impacts the rate at which a tree acquires and assimilates nutrients [8], and therefore climate hypothesis.The contribution of each factor to the overall health and growth of a tree vary both spatially and temporally.
The period 1983-2012 was the warmest 30-year period in the last 1400 years for the northern hemisphere.This is because atmospheric CO2 has increased by 40% since pre-industrial times [9].Warming is expected to continue beyond 2100 and northern latitudes are predicted to warm more rapidly and also the contrast between drier and wetter regions is expected to be more pronounced [9].These changes have consequences for forest productivity as lengthening of growing seasons leads to increased CO2 uptake.Conversely, higher soil temperature and increased biomass return to the forest floor would lead to an increase in decomposition activities and therefore higher respiratory losses [8].Effective resource monitoring is essential in better understanding how forest ecosystems are responding to global changes.Developing tools capable of tracking these changes is therefore paramount to achieving this goal.
Recent studies report that climate warming will result in an expansion of deciduous forests and reduction in conifer cover [10].A climate-driven increase in boreal fire activity is also expected to promote the expansion of aspen stands [11] in the future.Aspen (Populus tremuloides Michaux) is the most widely distributed broadleaf species in North America [12] and therefore ecologically important.The wood is light and soft and used for lumber, oriented strand-board, and pulp [13,14].Across Canada, the species may occur as single species stand, or in mixed species stands, commonly in combination with Balsam fir (Abies balsamea (L.) Mill.), Black spruce (Picea mariana Mill.BSP), White spruce (P.glauca Moench Voss), Yellow or White birch (Betula alleghaniensis Britt.and Betula papyrifera Marsh., respectively).In most of interior and western Canada moisture balance was reported to be the major factor influencing aspen growth and mortality [15][16][17], while warming was identified as the major driver of aspen growth in aspen dominated stands in eastern Canada [18].However, within mixed species stands, direct effects of warming on productivity were minimal [19].Aspen growth could also be significantly impacted by forest tent caterpillar (Malacosoma disstria Hbn.) defoliation to the extent that the direct effects of climate could become difficult to detect at shorter time scales [20].In better tracking the response of aspen forests to environmental changes over the longer term, models need to reflect a wide range of environmental conditions that are relevant for the species growth and productivity.
Site index (SI), defined as mean height of trees at a reference age (normally 50 years), is a measure of long-term productivity of forest stands.It is widely applied within even aged single species stands [21], but also applicable to mixed species forest stands [19].The direct relationship between stand height growth and volume growth, and the non-dependence of height growth on stocking [21], makes SI a useful and widely applied measure of forest productivity.The goals of this work were to apply the concept of SI in measuring productivity within aspen forests and identify the major factors driving productivity along climatic gradients and also along gradients of varying stand structure, composition, age, and soil attributes.This study sought to (1) determine the relative importance of the soil, developmental and climate hypotheses in driving aspen productivity across Canada; (2) parameterize a productivity model that could be used in monitoring the productivity of aspen forests across Canada; based on identified factors and (3) identify potential areas of emphasis owing to peculiar sensitivity to the identified drivers.
Aspen is a fast growing pioneer, intolerant hardwood species with preference for moist-textured soils [12].It is therefore hypothesized that both warming and moisture balance would be a significant factor controlling aspen growth.Being widely distributed, aspen is found on sites where substrate nutritional attributes vary geographically [22].It is therefore expected that aspen growth would vary geographically with soil nutritional status explaining the majority of this spatial variability.Given that aspen mixed species forest stands are widespread across Canada [6] and that productivity is influenced by competition and disturbances [19], it is hypothesized that developmental changes would dominate the dynamics of site productivity within mixed species stands.However, the relative contributions of each of these three factors to the overall health and productivity of aspen forests, is unknown.

Study Area
The study area is a longitudinal transect extending from Qué bec (from longitude 64° W) to British Columbia (to longitude 125° W) and from latitude 42° N in Ontario to 60° N in British Columbia (BC), Figure 1.The transect is located within three boreal ecological zones; plots in Quebec, Ontario and eastern Manitoba were within the boreal shield ecozone, those in mid-western Manitoba, Saskatchewan and Alberta were within the boreal plains while plots in BC were within the boreal cordillera ecological zone.The Boreal cordillera is dominated by mountains and valleys with wide lowlands and relatively sparse populations, mainly the Yukon population.The Boreal plains is the low lying interior of Canada, agriculture is an important activity here however, forestry remains the main industry with about 80% of the area forested.The boreal shield derives its name from the fact that geographically, at this point the boreal zone meets the Canadian shield rock.
The study area is comprised of 32-ecological regions (areas characterized by a particular distribution of vegetation types) and 126-ecological districts (land units characterized by particular relief, soil deposits, and drainage regimes).Drainage is largely imperfect and the landscape is generally gently sloped (<10%).Total annual precipitation varies along the study transect from about 600 mm in the west to about 1300 mm in the eastern part of the transect (Figure 1).Annual snow fall is about 160 mm in the boreal cordillera, 90 mm in the plains with an eastward increasing trend within the boreal shield to about 400 mm.Annual mean temperature varies from about −1 °C in the western part of the transect to 2.5 °C in the east.

Tree Measurement Data
In this study, inventory data was obtained from the provincial jurisdictions in all six provinces along the study transect.Quebec's Ministère des Ressources Naturelles (MRN) started a decadal province-wide measurement program in 1970, where circular plots (400 m 2 ) were positioned in forested lands through a stratified sampling scheme.Data used in this study comes from its third inventory program (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002).In Ontario, the Ontario Forest Research Institute has a forest Growth and Yield Program where 6400 m 2 main plots are sampled from managed and natural forests within which are nested three 400 m 2 growth plots and nine 25 m 2 shrub/regeneration plots.Data from the growth plots were used in this study.In Manitoba, data was collected from 500 plots by Louisiana Pacific Canada Ltd., Swan Valley Forest Resources Division.In Saskatchewan, trees were sampled within plots of varying sizes from 100-1000 m 2 depending on the stand density and height.Historically, data was collected by the Saskatchewan Ministry of Environment but responsibility was transferred to licensed holders.Data from different sources were assembled by Weyehaeuser into one database for Saskatchewan.In Alberta, 200-3400 m 2 plots were sampled and re-measured every five to ten years from 1960 by the Ministry of Environment and Sustainable Resource Development.Finally, data from British Columbia was obtained by the Ministry of Forests, Lands and Natural Resource Operations from plots monitored every five to fifteen years from 1920.The area of plots varied from 400-1000 m 2 , depending on stand density.
From a total of 8560 inventory plots within the transect area, 3548 plots had aspen in dominant or co-dominant canopy positions.Plots were identified as either pure aspen stands (19% of such) or mixture species stands (81% of such).A pure aspen plot was defined as having at least 75% of total plot basal area contributed by aspen.Plot density varied along the transect (Figure S1).Consequently, the number of trees selected for top height estimation varied by plot.By definition, top height refers to the 100 largest trees per hectare, the four largest trees per 400 m 2 , or the 30 largest trees per 3000 m 2 .For each plot, data was collected on trees with a diameter at breast height (DBH) > 9 cm.For age determination, trees were cored at breast height (1.3 m) or stump height (0.30 m) corresponding to breast height age or total age, respectively.Since the study focused on aspen, in stands where aspen proportion was relatively low, it was expected that plot mean of aspen age would be different from mean age estimated from all trees.

Climate Data
Climate data was generated from the BioSIM software [23] using the "climatic annual" model.This model computes annual climate variables from a specified number of weather stations nearest to the point of interest and corrects for differences in elevation and aspect of the weather stations to the point of interest.It also accounts for the distance of the weather stations to the point of interest by assigning greater interpolation weights to nearer stations.Two climate normals, 1961-1990 and 1981-2010 were used in generating climate data.
Climate variables included annual means of lowest daily minimum temperature, daily minimum temperature, daily mean temperature, daily maximum temperature, highest daily maximum temperature, annual total precipitation (mm), annual snowfall (mm), annual mean of daily mean dew point temperature (°C), annual mean of daily mean relative humidity (mm), mean wind speed, annual growing degree-day above 5 °C, frost days (number of days in the year when the daily minimum temperature is <0), frost free days (longest uninterrupted period without frost in the year), growing season days (the period between the last three consecutive days with frost in the spring and the first three consecutive days with frost in the fall), annual total potential evapotranspiration (mm), and aridity index, which is accumulated monthly water deficit (Thornwaite potential evapotranspiration-monthly precipitation, in mm; zero if negative).Other variables included annual mean of daily vapor pressure deficit (kpa), annual total radiation (MJ/m 2 ), number of days with precipitation >1 mm, consecutive days without precipitation (precipitation <1 mm) and climate moisture index (total annual precipitation-potential annual evapotranspiration) [15].

Soil Data
As part of the inventory program, data on soil texture, drainage class, slope class, surface deposit type, humus type, and depth of organic layer were collected.In dealing with differences in classification systems across the provinces, a particular system was selected and others re-classified along the selected system.For the purposes of this study, the Saskatchewan drainage classification system was adopted; classes 1-7 denote respectively; very rapid, rapid, well, moderately well, imperfect, poor, and very poor drainage classes.Slope percentages were grouped into 10 classes where 1 is level (<0.6%) and 10 is very steep slope (>100%).Humus types were also regrouped into five classes; Mull (MU), Moder (MD), Mor (MO), Anmoor (AN), and peatland humus (OT).The peatland humus type is a combination of the organic humus type and the peaty humus types.Both soil humus and texture are indicative of soil nutritional status [3,24].

Stand Attributes
Sampled stands vary in stand age, diameter structure and species composition.Diameter structure was measured with the Shannon diversity index computed in Equation (1) [25]: where pi is the basal area proportion of the diameter class i relative to the total stand basal area, and s is the number of diameter classes.

Model Selection
Across the provinces covered by the study transect, there are several site index models potentially applicable in this study.A literature search allowed consideration of five models.A model is considered because it is used in growth and yield estimation within a particular province or because of its potentially wide applicability.Considered were i) the model by Pothier and Savard [26] currently applied in growth and yield estimation in Quebec, ii) the model by Plonski [27,28] widely applied in Ontario, iii) the model by Huang et al. [29] used in growth and yield modelling (GYPSY model) in Alberta and validated for parts of central and western Canada [30], iv) the model by Nigh et al. [31], otherwise called the site tools software used widely in British Columbia and recently validated for central and western Canada [30], and finally, 5th model) the model by Garcia [30], which was calibrated using data spanning from Manitoba to British Columbia and thus valid for central and western Canada.
Models were first applied to data from the particular region for which they were calibrated.The site index values produced in these cases were called reference SIs.For a given region with a reference SI, the other models were applied to obtain four alternative SIs per region.The reference and alternative values were used in computing root mean square error for each alternative model and by region.The model selected was the one which overall, produced the minimum root mean square error.As the model of Nigh et al. [31] produced the minimum root mean square error, it was selected for site index estimation.

Estimating Site Index
The model of Nigh et al. [31] is as follows: where H is plot dominant height, SI is site index at reference age of 50-years, A BH is breast height age and a 0 , a 1 , a 2 are species specific model parameters.The height of aspen trees was available in the database from which plot dominant height was computed.
Existing site index models, especially those that produce anamorphic SIs, are not particularly suitable for juvenile forest stands [32].Previous work addressed this problem either by excluding juvenile stands from site index modelling or using growth intercept approach to estimate site index [33].Growth intercept approach computes SI from the early (<50 years) average height growth.Here, we applied the growth intercept method in obtaining SI values for juvenile (<30 years, Supplementary information section 2.0) forest stands.Stem analysis data from Anyomi et al. [34] was used for this purpose, and the resultant SI was denoted as the periodic growth index ( iH S ).Growth intercepts were computed for each breast height age and correlated with iH S .Results showed strong correlation between growth intercept and iH S (R 2 = 0.81), Equation (3).( where iH S is periodic growth index, GI is growth intercept (cm/year) and b1, b2 are parameter values.The parameter values of Equation ( 3) can be found in Table S1.Equation (3) was therefore used in estimating SI for juvenile stands.

Modeling Approach
Soil variables were first accounted for in the model (soil hypothesis) and then other variables (climate and stand developmental hypothesis) were introduced stepwisely.An approach in which the mean value of site index ( SI ) across the study transect is modified by retained explanatory variables ( i x ) was applied [34] in modeling productivity, Equations ( 4) and ( 5).
β and Xi q. β represent the linear and quadratic effects of the explanatory variable i x on SI.Categorical variables were converted to dummy variables so that they could be used in the NLMIXED procedure in SAS (see Supplementary information section 3.0).The product of n modifiers has a value close to unity when the variable i x is equal to the average i x observed value, and increases or decreases when moving away from the average.
A variable was retained in the model if it significantly (p < 0.05) related with SI and contributed substantially (>1%) to the total explained variance (or in reducing RMSE).Inclusion or exclusion of a modifier was determined with a likelihood-ratio test [35].
In verifying if some regional-level factors contribute to stand-level productivity, which were not already included in the model, a random effect variable ( s ) was introduced into Equation ( 4) and tested at defined landscape scales, Equation (6).
The parameters of Equation ( 6) were estimated with the NLMIXED procedure to, offer the potential to capture both fixed and random effects that may occur at different spatial scales, in this case stand and regional-levels respectively.Ecological district and ecological region were the regional-level categories considered.These models were rated against the model that included only fixed effects of local variables with a likelihood-ratio test.

Variability in SI
SI significantly correlated with stand variables (Figure 2A-D).SI decreased with stand age (AIC = 19,048.5;r = −0.51)as well as with an increase in proportion of spruce basal area (AIC = 19,275.8;r = −0.28).Aspen dominated stands exhibited higher SI which decreased with a decrease in proportion of aspen basal area (AIC = 19,355.1;r = 0.24).There was also a positive effect of stand horizontal structure measured with the Shannon index (AIC = 19,535.1;r = 0.05) such that SI slightly increased with structural complexity.SI also slightly varied by latitude (ns) and longitude (Figure 2E-F) (AIC = 20107.1;r = 0.05).Aspen stands located from longitude −110° W to −120° W were the most productive with mean SIs of over 25 m.SI slightly increased from longitude −108° W eastwards and also slightly increased with latitude with the most productive sites located above latitude 55° N. A significant positive correlation with mull humus type and a negative correlation with moder were also observed (AIC = 15,996; R 2 = 6%).There was also a negative effect of moisture balance on SI (AIC = 19,934.3;r = −0.22),particularly when moisture deficit was beyond 60 mm (Figure 2G-H).

Major Drivers of SI
Of the climate variables considered in this study, only relative humidity and potential evapotranspiration were not significantly correlated with SI and aspect was the only site variable not significantly related to SI. Soil texture was a major explanatory variable (R 2 = 8%), but there were 15-textural classes and in ensuring parsimony with respect to model parameters, soil humus type was rather retained.Soil effects were first accounted for in the model, humus was therefore the first explanatory variable retained, explaining 6% of variability in SI.Indeed, results suggest that aspen growth is better on mull humus type: mean SI was higher on mull compared to the other humus types.Higher variability in SI was observed on Anmoor humus type (Figure 2G).Stand age was the second most correlated variable after accounting for soil nutritional status and it explained 29% of the residual variance.Stand diameter structure measured with Shannon index and aspen basal area proportion subsequently explained 11% and 4% respectively of residual variance.Finally, aridity index explained 3% of the remaining variability in SI after controlling for soil and stand features.Even though other climatic variables were additionally significant, their contributions to total explained variance were marginal (<1%), hence they were not retained in the model.The selected variables, Equation ( 7), explained 53% of variability in aspen SI. , 1 where SI is predicted site index, S is the mean site index, , 1 are respectively modifiers for stand age, Shannon index, aspen basal area proportion and aridity index.
Ecological district ( D s ) as a random parameter, contributed 12% to the explained variability in stand-level SI, raising the total explained variability in SI to 65%.However, much of the variability captured by this parameter is related to plot density at that scale (Figure S2).
where D s is random effect parameter that varies by ecological district (Table S2), parameter values of Equation ( 8) can be found in Table 1.

Sensitivity to Climate
Sensitivity to aridity index (climatic variable retained in the model, Equations ( 7) and ( 8) as illustrated in Figure 3A,C reveal a pattern similar to the sensitivity to stand dynamics (effects of the changes in stand age, structure and composition; Figure 3B,D).Sites with annual mean temperature of less than −1 °C (Classes 1 & 2) were the most sensitive along the temperature gradient while areas with annual total precipitation ranging from 900 mm to 1300 mm (Classes 3 & 4) were relatively less sensitive to the effects of both climate and stand dynamics (Figure 3C,D

Stand Development Explains Long-Term Productivity
Consistent with expectations, results suggest that humus type is an important driver of aspen productivity.Humus forms vary spatially, mainly by geology and topography, and temporally, by climate and litter type.The nature of sub-soils influences litter quality and also the accessibility to earthworms, such that more acidic soils tend to promote mor humus formation while mull humus develops in more alkaline environments [24].Mull humus type is transformed organic material, mixed with mineral soil through activities of soil fauna and since it leads to more complete decomposition, nutrient availability is high and thereby supports greater productivity [3].Mors on the other hand are primarily products of fungal decomposition and because of incomplete decomposition, nutrients are fairly immobilized.Higher average SI observed on mull humus type is therefore due to its nutritional potential.Understanding the spatial repartitioning of humus forms and their changes through time within the boreal forest, would greatly assist in the long-term management of aspen productivity.That said, analysis also revealed that humus forms do not entirely reflect litter decomposition regime along the studied transect.
Results show that stand age, changes in species composition, and stand structure, correlate significantly with SI.Stand ageing, composition, and structural changes are dynamics that characterize stand development, and ultimately stand succession [36].In this study, SI declined over time; we observed mean SI of 21 m (standard error = 0.46) at mean stand age of 30 years, a mean of 17 m (0.38) at a stand age of 60 years and a mean SI of 14 m (1.40) at a stand age of 90 years.Three processes can explain the age-related decline in productivity, (1) physiological constraint to CO2 sequestration [37]; (2) accumulation of organic material that results in soil cooling and nutrient immobilization [38]; and (3) stand age as a driver of leaf litter quality decline [39].
Because a significant proportion of the studied plots (81%) has aspen trees in a mixture with other tree species, in line with our hypothesis, SI also vary by the composition of aspen.A mean site index of 20 m (0.78) was observed in pure aspen stands (aspen basal area proportion >75%) compared to a mean of 18 m (0.90) in mixed species stands.Compositional effects reflect the differential rates of litter deposition vis a vis rates of decomposition because hardwood leaves decompose faster than conifer leaves [3].Changes in the species composition of the dominant cohort [40] could also directly impact SI since tree species inherently vary in their growth potential.However, given that the effect of stand structure on SI is marginal (Figure 2D), it implies that compositional changes reflect more accurately the decomposition regime driven by changes in litter quality than by the direct impact of dominant cohort displacement on SI.
Finally and also consistent with research hypothesis, results suggest moisture deficit constraints on aspen SI, especially with deficits higher than 60 mm.By reference to an SI of 19 m, aspen height dropped by 5 m with a doubling in moisture deficit from 60-120 mm (Table 1; Figure 2).This effect is equivalent to the loss of 13 years in height growth.This observation adds to a growing body of evidence on the role of moisture in aspen growth and productivity [15][16][17], especially considering that aspen does not tolerate imperfectly drained soils and have only been observed on moderate to well drained soils.Yet, within the ensemble of environmental factors influencing aspen long-term productivity, the contribution of this variable was only 3%.

Aspen Forests in Cooler, Drier Environments are Most Sensitive to Productivity Drivers
Results suggest that, the productivity of aspen forests located in areas with mean annual temperature less than −1 °C was the most sensitive to growth drivers while relatively warmer regions of aspen stands (>2 °C) are less responsive to changes in growth drivers (Figure 3).Consequently, most of northern aspen stands (>latitude 54° N) of BC, Alberta and Saskatchewan were relatively more responsive to growth drivers when compared to those of Manitoba, Ontario and Quebec (Figure 4A).Along the transect from east to west, two major areas were most responsive to growth drivers, north east of Alberta and western Ontario, because of the relatively low precipitation amounts these areas receive annually (<900 mm).These observations were true, when only moisture balance was the explanatory variable in the model as well as when both moisture balance and stand attributes were considered.Overall, stands in western Quebec and eastern Ontario were most resistant (less responsive) to the impacts of moisture and stand dynamics.

Implications within the Context of Global Changes
Results from this study show that aspen composition has a positive influence on productivity and that stand dynamics have greater impact on forest productivity than the direct effects of climate.While stand dynamics may be directly influenced by climate [10], disturbance events are increasingly common within the circumboreal region [41,42] and substantially affect stand attributes [36,43]; this study shows that such dynamics control productivity of aspen forests across Canada.A recent study [11] suggests that the fire regime in the last millennia exceeds the limit observed within the last 10,000 years.Such intense fires promote shifts to an increasingly deciduous-dominated landscape [11,36].Aspen abundance is also favored by logging [44], and considering that the demand for aspen wood is expanding, primarily because it is a relatively low-cost, easy-to-use wood, an increase in logging would promote aspen expansion.With an increase in aspen abundance, this work reports that aspen productivity would substantially increase.Considering that plots cover a wide spatial extent and data was sampled by different provincial officials, differences in plot selection and measurement protocols could potentially affect our results.

Conclusions
Analyses of inventory data along an east-west transect show that variability in aspen long-term growth is mainly driven by species compositional changes and stand ageing effects.While soil humus plays a major role in long-term growth, it does not entirely reflect litter deposition and decomposition regime effects on productivity.Extreme dry events have major negative consequences for aspen growth but, within the ensemble of growth driving factors, its effects are relatively weak.It is concluded that stand development drives aspen productivity along an east-west precipitation gradient.An empirical model has been calibrated which could aid in monitoring forest site productivity of aspen forests in Canada.

Figure 1 .
Figure 1.Map showing the distribution of plots used in this study.Plot densities are higher in the east and the area of individual plots tends to be higher in the west.The deeper the shade of blue, the higher the amount of precipitation (1981-2010 climate normal).

Figure 3 .
Figure 3. Changes in the strength of productivity response to aridity index (A,C) and stand dynamics (the combined effects of stand age, structure and composition) (B,D); measured with co-efficient of determination (A,B) and Akaike information criterion (C,D).The solid red line depicts the variability in the strength of the relationships along the temperature gradient, while the broken dark blue line reveals changes in the strength of the relationship along the precipitation gradient.

Figure 4 .
Figure 4. (A) Northern lying (deeper shades) aspen forests are more sensitive to productivity drivers relative to more southern stands (lighter shades) along a north south temperature gradient; (B) Areas that are more responsive to growth drivers along an east-west precipitation gradient.The darker the shade, the higher the sensitivity of the stand.