Contrasting shrub species respond to early summer temperatures leading to correspondence of shrub growth patterns

The Arctic-alpine biome is warming rapidly, resulting in a gradual replacement of low statured species by taller woody species in many tundra ecosystems. In northwest North America, the remotely sensed normalized difference vegetation index (NDVI), suggests an increase in productivity of the Arctic and alpine tundra and a decrease in productivity of boreal forests. However, the responses of contrasting shrub species growing at the same sites to climate drivers remain largely unexplored. Here, we test growth, climate, and NDVI relationships of two contrasting species: the expanding tall deciduous shrub Salix pulchra and the circumarctic evergreen dwarf shrub Cassiope tetragona from an alpine tundra site in the Pika valley in the Kluane Region, southwest Yukon Territories, Canada. We found that annual growth variability of both species at this site is strongly driven by early summer temperatures, despite their contrasting traits and habitats. Shrub growth chronologies for both species were correlated with the regional climate signal and showed spatial correspondence with interannual variation in NDVI in surrounding alpine and Arctic regions. Our results suggest that early summer warming represents a common driver of vegetation change for contrasting shrub species growing in different habitats in the same alpine environments.


Introduction
Arctic-alpine ecosystems are sensitive to climate change (Settele et al 2014) and the rate of warming increases with latitude and elevation (Pepin et al 2015). Consequently, the Arctic-alpine biome in the northern hemisphere is expected to shift north-and upward (Settele et al 2014). Evidence from experimental warming (Walker et al 2006, Elmendorf et al 2012a, repeated vegetation surveys (Elmendorf et al 2012b), repeated photography (Sturm et al 2001, Tape et al 2006, and dendrochronology (Myers-Smith et al 2011, 2015a, suggests that shrub growth is sensitive to summer warming and that shrub cover has increased in response to climate warming throughout the tundra biome. Pollen records indicate greater shrub dominance in tundra during past warm episodes in Late Quaternary North America (Hu et al 2002, Higuera et al 2008 and Russia (Velichko et al 1997). In addition, a growing number of studies indicate shrubline ecotone advance up-and north-ward into tall-shrub free tundra (Myers-Smith and Hik 2018). However, growth responses of contrasting shrub species with different traits and habitats, occupying separate niches, are rarely compared in the same plant communities. Thus, we do not yet know whether all shrub functional groups are responding similarly to the rapidly warming tundra climate.
In general, taller deciduous shrubs are expected to replace lower statured species, and positive effects of ambient and experimental warming were predominantly observed on the abundance of taller deciduous shrubs (Elmendorf et al 2012a(Elmendorf et al , 2012b. Dwarf shrub cover was found to decline in a tundra biome-wide synthesis study on experimental warming (Elmendorf et al 2012a), especially in warmer areas; perhaps as a result of increased competition for light. Evergreen dwarf shrub species have, however, been shown to be climate sensitive (Bär et al 2008, Buizer et al 2012 and have been observed to increase their leaf size and height in response to experimental warming (Hudson et al 2011) and cover in response to ambient warming (Hudson and Henry 2009) at some High Arctic sites.
Satellite observations of vegetation productivity, as measured by the normalized difference vegetation index (NDVI), have shown a greening of large parts of the tundra biome in the northern hemisphere (Jia et al 2003, Goetz et al 2005, Guay et al 2014. Increases in NDVI over northern (Forbes et al 2010, Macias-Fauria et al 2012, northeastern Siberia (Blok et al 2011), and northwest North America (Tape et al 2012) have been related to increases in growth of deciduous shrubs and summer temperatures. Such links have yet to be explored for evergreen dwarf shrub species. Using aerial photography and satellite imagery, an increase in NDVI over the Low Arctic Tuktoyaktuk Coastal Plain, Northern Territories, western Canada, was shown to be related with an increase in shrub canopy cover (Fraser et al 2014). If remote sensing data are documenting tundra shrub expansion, then these data indicate that multiple species are responding synchronously to the changing climate conditions.
Here, we explore the links between NDVI, climate, and annual shrub growth of the tall deciduous shrub Salix pulchra and the evergreen dwarf shrub Cassiope tetragona co-occurring at an alpine tundra site in the Pika valley in the Kluane Region, southwest Yukon Territories, Canada. In this region Salix species have recently been shown to be expanding upward through new recruitment along the shrubline in response to winter warming and their growth has been shown to be sensitive to summer temperatures (Myers-Smith and Hik 2018). Annual growth of C. tetragona has been shown to be driven by summer temperatures at many High Arctic sites (Callaghan et al 1989, Havström et al 1995, Johnstone and Henry 1997, Rayback and Henry 2006, Rozema et al 2009, Weijers et al 2010, but may be more sensitive to fluctuations in nutrient availability at warmer sites in the Subarctic (Havström et al 1993). However, climate sensitivity of contrasting shrub species including a tall deciduous willow and a dwarf evergreen Cassiope species have yet to be directly compared.
Our research questions were: 1. What are the main climate drivers of annual growth variability of the co-occurring evergreen dwarf shrub Cassiope tetragona and the deciduous tall shrub Salix pulchra in the Pika valley of the Kluane Region, Yukon Territory, Canada?
2. Does the climate signal in the growth chronologies of the two species represent the regional climate?
3. How does the interannual variation in growth correspond with interannual variation in productivity as observed through the NDVI?

Study site and sampling
Our research site is located on the Pika valley slopes in the Kluane Region, Yukon Territory, northwest Canada (61.22 • N, 138.28 • W; figure 1(a)). On August 16, 2010, 16 samples of Cassiope tetragona were collected from the east-northeast-facing slope of the valley (figure 1(b)) at elevations of 1640-1673 m above sea level (a.s.l.). Several meters distances were kept between sampling spots to prevent repeated sampling of the same genet. Cassiope tetragona (L.) D. Don. (Ericaceae), or Arctic bell heather, is a multi-branched, clonal, hemi-prostrate, evergreen dwarf shrub with a circumarctic distribution (Eidesen et al 2007, Weijers et al 2017. At this site, C. tetragona is dominant in depressions with long-lasting snow cover and grows up to 30 cm in height. Samples of the tall deciduous shrub Salix pulchra were collected between 19 and 21 August, 2007 from both valley slopes ( figure 1(b)). Samples were taken at the shrubline (∼1812 and 1970 m a.s.l. on the east-northeast and west-southwest-facing slope, respectively) and at more downslope positions with a Salix species cover of approx. 50% (∼1715 m a.s.l.). In total, 17 S. pulchra specimens were sampled: four and three at the shrubline and three and seven downslope on the east-northeast and west-southwest-facing slopes, respectively. A 3-5 cm long disc was taken just above the stem-root boundary of the largest stem of each individual for growth ring analysis. Only distinct Salix patches were sampled, likely representing distinct genets. Salix pulchra Cham. (Salicaceae) is a canopyforming deciduous shrub found in the Siberian and northwest North American tundra (CYSIP 2017). It is the most dominant willow species in the Kluane Region east of Kluane Lake (Myers-Smith and Hik 2018). The species reaches heights between 20 and 80 cm at our site.

Climate
We calculated monthly precipitation sums (mm), and mean, mean minimum, and mean maximum monthly temperatures ( • C) as valid, i.e. lapse rate adjusted, for our research site for each month over the period 1901-2012 in ClimateWNA v5.40 (Wang et al 2017). In this program monthly normal data from climate mapping systems are downscaled to scale-free point data, which are then used as a baseline, together with monthly anomaly data from the Climate Research Unit Time Series (CRU TS3.23; Harris et al 2014), for the calculation of historical climate variables for specific locations and elevations (Wang et al 2016).
The mean annual temperature at our study site is −4.35 • C and mean annual precipitation sum is 502 mm (figure 2(a)). July is both the warmest and the wettest month with a mean temperature of 9.28 • C and mean precipitation sum of 69 mm. January is the coldest month (mean temperature −17.25 • C). We defined November-March as the winter months, April-May as spring, June-August as summer, and September-October as autumn, and calculated seasonal means and trends therein over two periods: 1950-2012 (figures 2(b) and (c)) and 1902-2012 (figure S1), as trends over these periods may have been relevant for shrub growth (cf. figure S2). Seasonal mean temperatures have risen significantly in winter (r = 0.29, p < 0.05), spring (r = 0.38, p < 0.01), and summer (r = 0.33, p < 0.01) over the period 1950-2012. Over the period 1902-2012 mean temperatures have risen significantly in spring (r = 0.26, p < 0.01) and summer (r = 0.38, p < 0.001; figure S1). Trends in seasonal precipitation sums were not significant.

Climate-growth analyses: linear mixed models
We measured annual growth of C. tetragona as shoot length increments and that of S. pulchra as ring width increments. See supplementary information available at stacks.iop.org/ERL/13/034005/mmedia for details on annual growth measurement. We used linear mixed model analyses to test the influence of monthly and seasonal climatic parameters on shrub growth. Climate-growth models were compared over the periods 1902-2009 and 1949-2009 (C. tetragona), or 1949-2006 (S. pulchra). The R-package nlme (Pinheiro et al 2017) was used for the mixed model analyses, with maximum likelihood estimation for model comparison and restricted maximum likelihood estimation for slope estimates (Crawley 2007). Before the analyses, the climate and shrub shoot length and ring-width chronologies were normalized at the individual level through subtraction of the mean, followed by a division by the standard deviation. Annual shoot lengths or ring widths of individual shrubs were included in the models as the response variable, and climate variables were included as fixed effects. A random intercept for year and an autocorrelation structure (AR1, autoregressive process of order one) were additionally included in the models. Conditional pseudo-R 2 values were calculated for each model with the r.squaredGLMM function of the MuMIn package (Nakagawa and Schielzeth 2013). We tested 94 climate-growth models including temperature means, mean maximums, mean minimums, and precipitation sums from 17 individual months (previous June to current October), the four seasons, and early (June-July) and late (July-August) summer as fixed effects, besides a null model for both species. In a first step, climate models that performed better than the null model were selected based on the Akaike Information Criterion (ΔAIC>2). Subsequently, Akaike weights for the selected models were calculated for model comparison. Akaike weights are a relative weight of evidence for each model and can be interpreted as the probability a certain model is the best model, given the selected set of models, for the observed data (Johnson and Omland 2004). High numbers of models in a comparison analysis increases the chances that a model has a lower AIC-value than the accompanying null model. By including 94 models for each species, models might thus come up by chance having a ΔAIC-value greater than 2. However, it is unlikely that such models will have an AIC value much lower than the related null model or have a high Akaike weight.

Climate-growth analyses: site chronologies
We tested climate-growth relationships between monthly climate parameters and standardized site chronologies. See supplementary information for details on standardization, and chronology construction and statistics (table S1 and figure S2).
We calculated Pearson's correlation coefficients and response function coefficients between each standardized site chronology and monthly mean temperatures as well as monthly precipitation sums, with monthly parameters from previous June to current October. This was done over the period 1907-2009 for C. tetragona and 1971-2006 for S. pulchra. Parts of the chronologies with n < 5, were excluded from the analyses, as were the most recent growth years. The latter were excluded due to possible incomplete growth at the time of harvest. In addition, Pearson's correlation coefficients and response coefficients were calculated over 25 year moving windows over the period 1907-2009, with a two-year window offset with the same monthly parameters, separately for mean monthly temperatures and monthly precipitation sums. This was only done for C. tetragona ; n ≥ 5), due to the relatively shortness of the S. pulchra chronology (1971-2006; n ≥ 5). Correlation and response coefficient significance was determined through a 1000 bootstrapped iterations. Response function analysis takes multicollinearity between climate parameters into account, through the regression of growth chronologies against principal components of climate parameters (Zang and Biondi 2015). The correlation and response analyses were performed with the R-package treeclim (Zang and Biondi 2015) in R version 3.4.1 (R Core Team 2017).
Spatial relationships in the area between 50 • -75 • N and 170 • -120 • W over the period 1981-2006 were calculated between the standardized site chronologies, June-July NDVI (0.22 • spatial resolution), and June-July mean maximum temperatures from the 0.5 • × 0.5 • monthly gridded meteorological dataset CRU TS4.00 (Harris et al 2014). Furthermore, the spatial relationships in the same area and over the same period between local mean June−July maximum temperatures (ClimateWNA v5.40), mean June-July NDVI, and mean June−July maximum temperatures (CRU TS4.00) were calculated. These calculations were carried out in KNMI Climate Explorer (Trouet and van Oldenborgh 2013), which uses a Monte Carlo approach for the determination of confidence intervals. The months June and July were chosen for the spatial analyses, because of their influence on growth of both species (see below). NDVI data for the period 1981-2006 were obtained from the GIMMS dataset (Tucker et al 2005) via KNMI Climate Explorer.

Climate-growth analyses
Variation in annual growth of both species was best explained by early summer temperature models (tables 1 and S2), with evidence for the mean June-July temperature model and the mean maximum June-July temperature as best model for C. tetragona and the mean maximum June-July temperature model as best model for S. pulchra. In addition, we found some support for mean (maximum) summer temperaturegrowth models for C. tetragona and mean (maximum) July temperature-growth models for S. pulchra (table  1). Besides summer climate models some current spring as well previous year late summer and early autumn temperature models were among the selected models for C. tetragona.
A similar picture arises from the correlation and response function analyses, with positive correlations between mean monthly temperatures of April−August of the growing season, as well as previous August and September and C. tetragona growth ( figure 3(a)). Mean June, July, August temperatures of the current year, and previous August remain significant predictors of C. tetragona shoot length growth, taken collinearity between the monthly climate parameters into account ( figure 3(b)). Mean June and July temperatures were again identified as the main predictors of S. pulchra radial growth (figures 3(c) and (d)).
The influence of mean June and July temperatures has been stable throughout the growth record of C. tetragona ( figure 4). There has been a shift, however, in the influence of August temperatures from the current year, which have influenced growth until the late 1950s. Thereafter, August temperatures started to influence C. tetragona growth in the next year ( figure 4(a)). The negative moving correlations found between C. tetragona growth and monthly precipitation sums, mainly of summer months ( figure 4(b)), are likely a result of collinearity between the monthly climate variables (cf. figure 4(d)). The influence of June and July temperatures may have been shifting back and forth between these two months (figure 4(c)).

Spatial analyses
Both standardized chronologies as well as the Pika valley mean June-July maximum temperatures are strongly correlated to mean June-July maximum temperatures over a large area in northwest North America (figures 5(a)-(c)). Furthermore, the C. tetragona chronology reflects mean June-July NDVI-values from a large area, mainly from the tundra and boreal treeline vegetation at higher elevations and/or latitudes such as at the Stikine plateau, Mackenzie mountains, Coast mountains, Babine range, Muskwa ranges, and in the area north of the Brooks ranges ( figure 5(d)). The standardized S. pulchra chronology shares the spatial correlations with June-July NDVI over the alpine tundra southeast to our research site, as found for C. tetragona, but lacks those to the north and northwest (figure 5(e)). Neither chronology is related to the June-July NDVI-values and temperatures of the boreal forests of the Alaskan interior. The spatial relationship between the lapse rate corrected mean June-July maximum temperatures as valid for the Pika valley and NDVI (figure 5(f)) is similar as that for the C. tetragona chronology, but it lacks the relationships with June-July NDVI north of the Brooks Range.

Discussion
We have shown that Cassiope tetragona and Salix pulchra growth at the Pika valley in the Kluane Region, Yukon Territory, Canada is driven by early (June-July) growing season temperatures. Despite contrasting plant traits of dwarf evergreen versus tall deciduous shrubs and different habitats related to snow distribution, consistent climate sensitivity to early growing season temperatures was observed. Early summer temperatures, which coincide with the time of year with maximum insolation, are warming across northwest North America and shrub growth chronologies for Table 1. Results of the mixed model analyses with annual shoot lengths of Cassiope tetragona or annual ring widths of Salix pulchra included in the models as the response variable and climate variables as fixed effects, calculated over the period 1949-2009 (C. tetragona) or 1949-2006 (S. pulchra). Selected models are models with AIC values of at least two lower than the corresponding null model. R 2 -values are conditional pseudo-R 2 values. Akaike weights are the relative weight of evidence for each model. T: mean temperature; T max : mean maximum temperature; T min : mean minimum temperature; P: precipitation sum.

Species
Selected both species were correlated with the regional climate and showed spatial correspondence with interannual variation in NDVI in surrounding alpine and Arctic regions. Our results suggest a common driver of vegetation change for contrasting shrub species growing in different habitats in the same alpine environments. In contrast to S. pulchra, we found some influence of late summer and early autumn temperatures of the previous year on growth of C. tetragona, besides the influence of early summer temperature, as was earlier observed for this species on Svalbard (Weijers et al 2010). This is likely a result of the evergreen nature of C. tetragona leaves and the formation of its primordial leaves at the end of the growing season. Photosynthesis in C. tetragona leaves starts before snowmelt is completed (Starr and Oberbauer 2003) and shoot elongation may thus benefit from warm early summer temperatures. Annual growth of C. tetragona has previously been shown to be driven by summer temperatures at several High Arctic sites, through experimental warming in Ellesmere Island, Canada (Hudson et al 2011), Greenland (Campioli et al 2013) and Svalbard (Havström et al 1993, Rozema et al 2009, and dendrochronological analyses in Ellesmere Island, Canada (Havström et al 1995, Johnstone and Henry 1997, Rayback and Henry 2006, North Greenland (Weijers et al 2017), and Svalbard (Callaghan et al 1989, Aanes et al 2002, Weijers et al 2010, 2013b. However, at warmer Subarctic alpine tundra sites, as our site, environmental growth controls of C. tetragona may be less uniform. Growth in Subarctic Abisko, North Sweden, for example, showed a greater response to nitrogen addition than to experimental warming, and may be more nitrogen-than temperature-limited (Havström et al 1993). Still, after longer timespans (22 years) neither fertilization nor warming was found to affect C. tetragona growth at the same experiment (Campioli et al 2012).   The C. tetragona shrubs in our study were from downslope depressions with late snowmelt. Winter snow depth at Arctic-alpine localities is, unlike at flat lowland sites, relatively independent of the winter precipitation sum, as excess snow is removed by wind and redistributed according to topography (Erickson et al 2005). Snowmelt date and growing season length at these C. tetragona localities are thus likely largely determined by temperature instead of precipitation, which may explain the strong relationship found between early summer temperatures and growth. In addition, mean growth rate at our site (5.21 mm year −1 ) is close to the 5.05 mm year −1 reported for High Arctic Svalbard , despite the warmer conditions at our site (mean July temperatures of 9.28 • C versus 6.43 • C). This suggests relatively harsh conditions at the downslope snowbeds on the eastnortheast-facing slope at our site, which may explain the dominance there of C. tetragona, which is generally a more High Arctic species. Relatively harsh conditions may persist at these places due to a shortened growing season in snowbeds, cooling of soils by (upslope) meltwater, and low amounts of direct sunlight due to the east-northeast slope aspect. During mornings, sunlight is blocked by the opposing mountain.
Early growing season temperatures have earlier been identified as the most important factor influencing S. pulchra growth in the northeastern Siberian tundra (Blok et al 2011) and may stimulate its growth as leaf flush in S. pulchra takes place within days after snowmelt (Borner et al 2008). Moreover, Salix spp. growth was found to be driven predominantly by summer temperatures, and their recruitment by winter temperatures, throughout the Kluane Region (Myers-Smith and Hik 2018). In contrast, secondary growth of S. pulchra showed little response to fertilization and warming treatments when in competition with Betula nana (Bret-Harte et al 2002) and its abundance declined under both treatments (Bret-Harte et al 2001) in the moist tussock tundra at Toolik Lake in the northern foothills of the Brooks Range, Alaska. However, at a nearby site, June temperatures were found to be important for radial growth of S. pulchra (Ackerman et al 2017). Despite some variation among studies, June and July temperatures appear to be consistent drivers of variation in S. pulchra shrub growth.
Annual C. tetragona and S. pulchra growth at our site correspond with mean maximum June-July temperatures for a large part of northwest North America.
C. tetragona shoot length growth furthermore tracked June-July NDVI-values of the vegetation at higher elevations and latitudes in this region and radial growth of S. pulchra corresponded with mean maximum NDVI values of a large area southeast of the study site. This suggests that annual growth variability of C. tetragona and S. pulchra respond to climate drivers with wide spatial extents and that productivity in these alpine regions may correspond across the whole region. Our findings correspond with NDVI observations of a greening Alaskan and Yukon tundra since 1982 (Verbyla 2008, Beck and Goetz 2011, Guay et al 2014, Ju and Masek 2016. Similar inter-correlations between annual shrub growth, summer temperatures, and NDVI have earlier been found for the deciduous shrubs S. lanata and Alnus fruticosa in northern Siberia (Forbes et al 2010, Macias-Fauria et al 2012 and for S. pulchra in northeastern Siberia (Blok et al 2011).
The climate signal in S. pulchra is somewhat less strong than that in C. tetragona shrubs and its growth corresponded with NDVI variability over a smaller area than C. tetragona growth. This might be because the evergreen dwarf shrub C. tetragona is not grazed upon, in contrast to many deciduous Arctic-alpine shrub species (Mallik et al 2011). At our site, the deciduous S. pulchra shrubs may have been prone to ptarmigan browsing new buds in spring, non-cyclic insect herbivores, stem herbivory by small mammals including marmots and rare browsing by moose or other large herbivores (Myers-Smith and Hik 2018). Expansion of deciduous shrubs has earlier been shown to be suppressed by reindeer grazing in the Scandinavian alpine tundra, in contrast to that of evergreen shrubs (Vowles et al 2017).
Although Arctic tundra greening has mainly been attributed to an expansion of erect and tall deciduous shrub cover (Elmendorf et al 2012a, 2012b, Fraser et al 2014, we found that annual growth of the hemi-prostrate evergreen dwarf shrub C. tetragona largely corresponds with NDVI as a proxy for productivity over large parts of the northwest North American tundra. C. tetragona might be able to expand its cover through a densification of existing shrub patches in the warmer parts of the tundra biome in certain environmental niches with long lasting snow cover, where relatively harsh conditions may persist during the growing season. C. tetragona is capable of forming dense mats, which might inhibit the recruitment of other taller shrubs in such locations.

Conclusions
We found a strong positive growth response of the cooccurring evergreen dwarf shrub Cassiope tetragona and the deciduous tall shrub Salix pulchra to early summer warming in the Pika valley of the Kluane Region, Yukon Territory, Canada. Despite differences in plant species traits and habitats, the two contrasting species demonstrated surprisingly consistent growth responses to climate drivers. Moreover, our findings show that the annual growth variability of these species are likely representative of the annual variability in tundra vegetation productivity of large parts of northwest North America. Early summer warming has likely enhanced growth rates of entire shrub communities in this region, resulting in the densification of shrub stands and a greening of the Arctic-alpine tundra of Alaska and the Yukon, as reported in the literature.