Greater effect of increasing shrub height on winter versus summer soil temperature

Shrub expansion is increasingly observed in arctic and subarctic environments. The development of shrub structure may significantly impact the abiotic environment at the local scale. Our objective was to reconstruct the development of the vertical structure of Betula glandulosa Michx. and to evaluate its effects on winter and summer soil temperature and on snow depth. Stratified sampling of the shrub revealed that shrub biomass distribution followed a similar pattern in stands of contrasting heights. Woody biomass was maximal in the lower stratum and relatively stable in the intermediate strata, while the foliar biomass tracked the vertical development of the shrub structure. Dendrochronological analysis revealed that shrub stands are relatively young; most of the dominant stems started their development after 1990. Shrub height was positively associated with both the dominant stem age and its vertical growth rate. Temperature differences among sites were greater during winter (ca 10 °C) than during summer (ca 2 °C), while the sum of freezing degree-days varied from 680 °C to 2125 °C. Shrub height was the most plausible variable explaining snow depth, winter ground level temperature and the sum of freezing degree-days. However, woody biomass in the 30–40 cm strata best explained summer ground level temperature. Our results suggest that the development of a shrub structure will have far-reaching consequences on the abiotic environment of subarctic ecosystems.


Introduction
Arctic and subarctic ecosystems have been subjected to rapid and intense climatic change over the last few decades, triggering permafrost degradation, reduction in sea-ice cover and trophic mismatch, to name a few (IPCC 2013). In terrestrial ecosystems, shrub expansion appears to be generalised within the circumpolar area according to numerous studies using either aerial photos or satellite image analyses (Russia: . Long-term vegetation models predict that erect shrubs could replace graminoids, herbs and prostrate shrubs and affect 50% to 70% of all arctic communities by 2050 (Pearson et al 2013).
It is generally accepted that both low temperatures (Parsons et al 1994, Chapin et al 1995, Van Wijk et al 2004, Wu et al 2011, Elmendorf et al 2012 and nutrient availability (Chapin et al 1995, Chapin and Shaver 1996, Zamin and Grogan 2012, Heskel et al 2013 restrict plant growth in arctic and subarctic environments. Shrub species were shown to respond positively to warmer temperatures and increased nutrient availability (Forbes et al 2010, Tape et al 2012, Frost andEpstein 2014) and precipitation (Blok et al 2011b). Under experimental conditions, increases in foliar, woody and/or root biomass were observed for different shrub species (Parsons et al 1994, Chapin et al 1995, Chapin and Shaver 1996, Wu et al 2011, DeMarco et al 2014, Paradis et al 2014. Sturm et al (2005) were the first to suggest that shrub expansion was likely associated with a positive feedback loop involving soil temperature, nutrient availability and snowy precipitations. Shrub cover decreases both the redistribution (Liston et al 2002, Pomeroy et al 2006) and sublimation of snow by wind (Liston et al 2002), resulting in greater snow accumulation under shrubs than in the adjacent open tundra (Sturm et al 2001b, Pomeroy et al 2006, Myers-Smith and Hik 2013. Since snow has an insulating effect on soil temperature during winter (Hardy et al 2001, Zhang et al 2008, Myers-Smith andHik 2013), winter soil temperatures under shrub cover are warmer than in the open tundra (Sturm et al 2001b, Myers-Smith andHik 2013). This phenomenon could increase microbial activity and, therefore, nutrient availability for plant growth during the subsequent growth season (Nadelhoffer et al 1991, Schimel et al 2004. Since shrub species are known to efficiently use newly available nutrients for growth (Chapin and Shaver 1996), an increase in nutrient availability resulting from a greater shrub cover could exacerbate shrub development (Sturm et al 2005). Although one might argue that the shade provided by a denser shrub cover during summer might reduce soil temperature (Blok et al 2010, Epstein et al 2013, Myers-Smith and Hik 2013 and microbial activity, it is generally accepted that the warming effect of shrub cover on soil temperature during the winter is greater than its cooling effect during the summer (Lawrence and Swenson 2011, Loranty et al 2011, Bonfils et al 2012, Epstein et al 2013. The aim of this study was to evaluate the effect of the vertical structure of Betula glandulosa Michx on soil temperature and snow cover in order to understand how shrub expansion could influence subarctic terrestrial ecosystems in the near future. To do so, we characterised and reconstructed the development of B. glandulosa vertical structure (height, age, biomass, stem length, wood volume) in 22 sites in the Umiujaq region of Nunavik (subarctic Québec) and described the abiotic environment (summer and winter soil temperature, annual freezing and thawing degree-days and snow depth) in each of the stands. We hypothesised that B. glandulosa vertical structure influences the abiotic environment mainly during winter through its impact on snow cover. We predicted (i) that snow depth and winter soil temperature would be higher in sites with taller shrubs, (ii) that annual freezing degree-days would be lower in sites with taller shrubs, (iii) that summer soil temperature and thawing degree-days would not be influenced by B. glandulosa vertical structure, and (iv) that shrub height would be the most significant variable in explaining winter soil temperature, annual freezing degree-days and snow depth.

Study species
Betula glandulosa is an erect deciduous shrub species that can grow up to 2.5 m (DeGroot et al 1997). Its geographical distribution ranges from Greenland to Alaska and all the way south to California (DeGroot et al 1997). It reproduces either sexually via abundant viable seed production (Hermanutz et al 1989) or asexually via clonal growth Hermanutz 1988, 1993). It can tolerate a wide array of abiotic conditions and can be found on well-drained sites as well as in wetlands. In Eastern Canada, it has been shown to expand rapidly over the last decades in response to climate change

Sampling design
A total of 22 stands dominated by B. glandulosa were selected for the study. Stands (>150 m 2 ) were selected to cover the height range of B. glandulosa in the region (ca 20 to 140 cm) and were located either in the coastal area (7), in the Tasiapik Valley (8) or on the plateau between the two (7). In each stand, B. glandulosa height was calculated as the average of ten randomly located height measurements. Total shrub cover (all shrub species; nomenclature used: VASCAN; Brouillet et al 2010+) was also estimated using the following cover classes: 0%-1%, 1%-5%, 5%-15%, 15%-25%, 25%-50%, 50%-75%, 75%-90%, 90-%100%. Soil temperature was recorded at ground level every 6 h with iButton Thermochron® data loggers (±0.5°C, model DS1921G-F5#, Dallas Semiconductor Corporation, Dallas, TX, USA) from September 2013 to August 2014. A long pole (at least 10 cm taller than the shrub stand) was placed in the vicinity of each data logger, making it possible to locate it during winter. In late February 2014, 16 randomly distributed measures of snow depth were taken in each stand, including one measure near the location of the data logger (ca 25 cm), since most of the poles were still visible above the snow cover (18 sites). For the remaining sites (4), we used GPS coordinates to locate the data loggers. The 15 other measurements were taken within 3 m around the data logger location to make sure snow measurements were taken within the shrub stands. Snow cover was not disturbed during the measurements since it was hard enough to support the weight of a person.

Betula glandulosa vertical structure characterisation
In August 2013, we sampled the vertical biomass of B. glandulosa within four 0.5 m 2 quadrats, randomly located in each stand but at least 4 m apart from the data logger. We sampled B. glandulosa leaf and woody biomass within 10 cm vertical strata (starting from the top of the shrub). In the laboratory, we sorted, dried (at 70°C for 72 h), and weighed leafy and woody biomass with a high precision scale (A&D, Model ER-182A, ±0.0001 g) for each stratified sample. Afterwards, we randomly selected one quadrat per stand and measured the length and diameter of each woody ramification (>1.5 cm long) with an electronic caliper (±0.001 cm). Wood volume (cm 3 ) was then calculated with the following equation: where V=volume, r=radius and l=length of the ramification.

Betula glandulosa vertical structure development
In August 2013, we also collected stems to reconstruct B. glandulosa vertical structure development. Five dominant B. glandulosa stems per stand were randomly sampled. We sampled the stems at heights (vertical distance from the ground) of 0, 10, 20 cm and so on. Back in the laboratory, samples were dried at room temperature for six weeks, then boiled in water for at least two hours, and sliced with a microtome (ca 20 μm, WLS core-microtome). Slices were coloured with a safranine solution (1% w/v in distilled water) and fixed on permanent microscope slides. Growth rings were counted along two radii for each sample, since dwarf shrubs may form incomplete growth rings. We used the age of the stems (samples collected at height 0 cm) to determine the minimal age of the stands and used the age of the samples collected at the different heights to reconstruct the development (vertical growth rate in cm yr -1 ) of the vertical structure of the stand. In order to obtain an average growth rate per stand, the growth rate of the 5 stems sampled in each stand was averaged.

Statistical analyses
All statistical analyses were performed using the software R (version 3.0.3, R Development Core Team, Vienna). Linear regressions were conducted to evaluate if B. glandulosa biomass (total, woody and foliar) was a function of stand height and to evaluate the relationship between dependent variables (soil temperature, degree-days and snow cover depth) and shrub structure properties (height, age, biomass). We also used linear regressions to evaluate the relationships between stand height and moss and lichen covers.
We built ecologically relevant statistical models that could explain the following climatic variables: ground-level winter (December-March) and summer (July-August) temperatures, snow depth and annual soil freezing and thawing degree-days. The explanatory variables available to build the linear regression models were all descriptors of the shrub vertical structure and were classified into five categories: height of the stand, minimal age of the stand, woody volume, stem length and biomass. For the latter three categories, data were available for each of the 10 cm-strata. In order to decrease the number of models, we selected only one explanatory variable from each of these categories, based on the Pearson correlation coefficient between these variables and the dependent variables. For example, we built the models to explain snow depth using the woody biomass in the 30-40 cm strata because it was this variable that showed the highest correlation coefficient with snow cover. We did not include interactions between the explanatory variables in the linear models as the degrees of freedom were restricted due to the number of sites.
We used a model selection approach (AICcmodavg package, Mazerolle 2014) based on Akaike's information criterion to determine which model(s) was the most plausible to explain the variable of interest. We tested each model for colinearity using Variance Inflation Factors (VIF; car package, Fox and Weisberg 2011) in order to remove variables from the model when VIF>10. Because no model had a VIF>10, we decided to leave all models as initially built even though some variables were correlated due to the nature of the data set. We followed the recommendations of Burnham and Anderson (2002) who suggest using the multimodel inference if the best model has a wAICc<0.90 in order to obtain the model-averaged estimate of the parameters, the unconditional standard error and the 95% confidence interval. We concluded that one particular variable had an effect on the dependent variables only if the confidence interval excluded 0.

Betula glandulosa stand characteristics
The height of the B. glandulosa stands varied from 21 to 121 cm. Apart from the dominant B. glandulosa (>50% cover in all sites), other erect shrub species found at the sites included Salix planifolia Pursh. and S. glauca L., although their abundances were lower (<25% of total shrub cover; table 1). Ground layer vegetation was dominated by lichens, mosses and herbaceous species as well as by low shrub species such as Empetrum nigrum subsp. hermaphroditum (hagerup) Böcher, Rhododendron groenlandicum (Oeder) Kron and Judd, R. lapponicum (L.) Wahlenb., Vaccinium uliginosum L., V. vitis-idaea L. and Salix uva-ursi Pursh. Stand height of B. glandulosa had a significant positive effect on moss cover ( F 1,20 =42.00, p<0.001) but a negative one on lichen cover ( F 1,20 =44.13, p<0.001).
The abiotic environment varied greatly among the 22 stands. Winter (December-March) and summer (June-August) temperatures at ground level ranged from −14.9°C to −4.2°C and from 11.2°C to 13.3°C, respectively. Annual freezing degree-days varied from 680°C to 2125°C while annual thawing degree-days varied from 1241°C to 1563°C. Snow depth, measured in February 2014, ranged from 22 to 101 cm.
For all sites, and regardless of the height or age of the B. glandulosa stems, woody biomass was always greatest in the 0-10 cm stratum (figures 2(a)-(c), but see figure S1 for all 22 sites). For B. glandulosa <50 cm, woody biomass decreased sharply above the 0-10 cm stratum ( figure 2(a)). For taller shrubs, woody biomass tended to stabilise in the strata above the 0-10 stratum before declining sharply in the 2-3 higher strata (figures 2(b) and (c)). In contrast, foliar biomass appeared to track the development of B. glandulosa vertical structure, moving upwards as the shrubs grew in height (figures 2(a)-(c)). For taller shrubs, leaves were nearly absent from the lower strata.
3.3. Betula glandulosa age and vertical growth Average stem age for the different sites ranged from 13 to 31 yr. Shrub height (21-121 cm) was positively associated with age ( F 1,20 =30.9, p<0.001; figure 3(a)). The vertical growth rates of individual stems ranged between 1.04 and 8.94 cm yr −1 (figures 3(b)-(d), but see figure S2 for all 22 sites). Once averaged, growth rate per stand varied between 1.80 and 4.67 cm yr −1 . Furthermore, average vertical growth rates per stand were positively correlated with shrub height ( F 1.20 =18.37, p<0.001; figure 3(e)), lower shrubs having lower vertical growth rates than taller shrubs.
3.4. Betula glandulosa structure impacts on soil temperature and snow depth None of the most plausible models that explained variation among sites in winter and summer temperatures at ground level, in annual freezing and thawing degree-days and in snow cover depth had a wAICc >0.90 (table 2). We therefore computed the weighted average of the estimates of the variables included in all models instead of relying solely on the estimates of the best model.
Betula glandulosa height was the sole variable in the most plausible model to explain variations among sites in winter temperature at ground level, in freezing degree-days and in snow depth (table 2). The weighted average of estimates procedure gave slightly different results. For winter temperature, the weighted average of the estimates of all variables revealed that both B. glandulosa height (estimate: 0.06, unconditional SE: 0.02, 95% confidence interval: 0.02, 0.10; figure 4(a)) and age (estimate: 0.023, unconditional SE: 0.12, 95% confidence interval: 0.01, 0.46; figure 4(b)) had a significant positive effect on winter temperature. Similar results were obtained for snow cover depth since B. glandulosa height (est.: 1.94, uncond. SE: 0.94, 95% conf. int.: 0.10, 3.78; figure 4(c)) and age (est.: 0.40, uncond. SE: 0.18, 95% conf. int.: 0.05, 0.75; figure 4(d)) were the only two variables that had a significant effect on this dependent variable. However, the weighted average of the estimates of all variables revealed that B. glandulosa height remained the sole variable with a significant impact on freezing degreedays, (est.: −8.63, uncond. SE: 2.61, 95% conf. int.: −13.75, −3.52, figure 4(e)).
For summer temperature, the model including only the woody biomass H30-40 was the most plausible to explain variations among sites (table 2). This result was supported by the weighted average of the estimates of all variables that revealed that only the woody biomass H30-40 had a significant impact on summer temperature at ground level (est.: −0.01, uncond. SE: 0.00, 95% conf. int.: −0.01, −0.00; figure 4(f)). For thawing degree-days, age was the sole variable in the most plausible model to explain the variation between sites. However, the weighted average of Table 1. Vegetation cover and characteristics of Betula glandulosa in 22 sites sampled during summer 2013, near Umiujaq, Nunavik. Vegetation cover was visually evaluated by cover classes in all shrub stands (>150 m 2 ).

Discussion
According to Sturm et al (2005), recent observed shrub expansion is likely associated with a positive feedback loop involving warmer soil temperature, increased nutrient availability and increased snow depth under shrub cover. In this study, our objective was to evaluate how shrub vertical structure influences soil temperature and snow cover. We were able to demonstrate that winter temperature at ground level, snow depth and annual freezing degree-days were all significantly influenced by B. glandulosa height, while summer temperature at ground level was influenced solely by the woody biomass in the H30-40 stratum (but the effect size was much lower than for winter ground level temperature). Moreover, none of the variables measured significantly explained the variation in annual thawing degree-days. Thus, our results support the hypothesis that an erect shrub cover has a The positive association between moss abundance and shrub height suggests that shrub canopy development and closure can lead to the replacement of lichens by mosses, a phenomenon that is likely to persist in subsequent years. The vertical structure of the studied B. glandulosa stands was relatively young, ranging from 13 to 31 years. Although we cannot ascertain it, we believe that The stem analysis conducted in the different sites allowed us to reconstruct the development of B. glandulosa vertical structure. The vertical development appears to be linear for most stems (see slopes on figures 3(b)-(d) and figure S1), a somewhat unexpected result. Such regular increase in height in all studied sites revealed that B. glandulosa individuals are able to grow vertically throughout the landscape. Such results, corroborated by personal observations of rapid B. glandulosa vertical development at several sites in the region, could have significant implications for the abiotic environment in the Umiujaq region, including an increase of winter soil temperature and an exacerbated rate of permafrost degradation.
Moreover, B. glandulosa vertical structure was similar between stands of similar height. Overall, woody biomass was greater in the lower stratum. Woody biomass in this stratum supports most of the  aboveground biomass and provides resistance to gravity and wind (King and Loucks 1978). It is likely the result of a longer period of wood accumulation in the lower strata. In fact, while primary growth leads to axis elongation during the first year of wood formation, secondary growth, which occurs over several years afterwards, leads to stem diameter increase (Raven et al 2014). Foliar biomass did not follow woody biomass distribution patterns, as it tracks the development of the vertical structure. We found few or no leaves in the lower strata of the taller shrub stands, where photosynthetically active radiation is believed to be limited as the increased foliage density in the upper half of the shrub's vertical structure absorbs a significant percentage of the incoming solar radiation  to show that snow depth is positively associated to the height of B. glandulosa stands, although shrub height explained only 27.1% of the variation observed in snow depth between sites. The percentage of variance explained is however higher (39.3%) if we removed a single outlier. Moreover, minimal stand age also explained 25% of the variability in snow depth. The relatively low percentage of variation explained by shrub height and minimal stand age suggests that other variables, such as local topography, have a strong influence on snow accumulation (Sturm et al 2001b, Pomeroy et al 2006. Although our experimentation did not allow us to draw any conclusion about snow cover duration, preliminary temperature data suggest that snow accumulates rapidly under shrub cover at the beginning of the winter season, a likely consequence of snow redistribution by wind. The same data also suggest that snow tends to melt rapidly in spring, corroborating the results of studies that have shown that snow albedo was influenced by shrub cover during spring, resulting in warmer air temperature and faster snowmelt ( Shrub height also had a positive effect on the winter temperature at ground level and a negative effect on the sum of freezing degree-days at the soil-snow interface. While winter temperatures were best explained by shrub height and by minimal stand age, the sum of freezing degree-days was best explained by shrub height only. The range difference in winter ground level temperature (10.7°C) found among our study sites is comparable to the one observed by Sturm et al (2001b) but is higher than the ca 5°C difference observed by Lawrence and Swenson (2011)  Overall, our results show that shrub height is the most significant predictor of the measured winter climatic parameters (temperature, freezing degree-days, snow depth). Although age is also significant for 2 out of 3 variables, we believe that these relationships with age might not stand as well in other arctic or subarctic regions where the development of the vertical structure is hindered by harsher climatic conditions, leading to shorter shrub stands, regardless of the age of the individuals. The inclusion of other variables (biomass, total stem length and wood volume) did not increase the performance of our models. However, these variables might have a greater impact on snow physical properties (Domine et al 2015). For example, reduced snow compaction under shrub cover increases the thermal resistance of snow (Sturm et al 2001b), increasing its insulation capacity.
The effect of shrub vertical structure on average summer temperature was barely significant and calls for caution in extrapolating these results. Variability of ground level temperature in summer was better explained by the woody biomass in the 30-40 cm vertical stratum, a somewhat unexpected result. It is possible that woody biomass in this stratum allows for a greater segregation of the different sites since shorter shrubs had no woody biomass in this stratum (height<30 cm) while taller ones had an important one. Regardless of this particular result, the variation of ca. 2°C in average summer temperature between shorter and taller stands was comparable to the one observed between shrub-dominated and shrubless sites in other studies (Blok et al 2010, Myers-Smith andHik 2013). When compared to the 10.7°C observed among sites during winter, these results suggest that the impact of the shrub vertical structure is more important in winter than in summer.
One limitation of our study is the difficulty to discriminate the relative effect of shrub height and moss cover on soil temperature due to the fact that both variables are strongly associated. Changes in the abundance of moss species are believed to be an important driver of water and energy balance between soil and atmosphere in tundra ecosystems (McFadden et al 2003, Beringer et al 2005, Blok et al 2011c. In fact, moss species' low thermal conductivity can have a cooling effect on soil in summer (Gornall et al 2007), a phenomenon likely to favour the persistence of permafrost (Zimov et al 2006). When combined with the shade provided by dense shrub canopy, both phenomena could result in lower summer soil temperature and reduced permafrost degradation (Blok et al 2010). However, our results strongly suggest that the observed effect of shrub vertical structure development, regardless of moss cover, is greater in winter than in summer. While annual freezing degree-days decreased sharply from 2125°C to 680°C with increasing stand height, thawing degree-days varied only from 1241°C to 1563°C among sites. When considering that an increase of 1°C or 2°C in mean annual temperature is sufficient to increase permafrost degradation , the expansion of shrub cover, and its overall impact on soil thermal regime, is expected to significantly accelerate permafrost degradation.

Conclusion
Our study revealed that shrub stand height had a significant impact on snow depth, winter soil temperature and freezing degree-days in the Umiujaq region. Moreover, it highlighted that the vertical structure of B. glandulosa stands had developed rapidly over the last three decades. Although our study provided useful information on the relationship between B. glandulosa stand structure and snow depth, we did not evaluate how B. glandulosa biomass distribution influenced the quality of the snow. Future studies should then focus on the relationship between woody biomass structure and the characteristics of the snow pack (density, isolation coefficient, etc).