Introduction

In the southern parts of arctic environments, shrub communities are one of the most prominent types of vegetation and therefore play an important role in ecosystem functioning (Chernov, 1985; Bliss and Matveyeva, 1992; Blanken and Rouse, 1994). According to recent aerial photo surveys, the abundance of tundra shrubs in northern Alaska has considerably increased in recent decades apparently due to climate change (Sturm et al., 2001; Tape et al., 2006). There is also evidence showing that the experimental warming of air temperature favors growth of shrubs at the expense of other tundra plants (Chapin et al., 1995; Walker et al., 2006). The ongoing and predicted climatic warming in the Arctic (ACIA, 2005; IPCC, 2007) is likely to affect growth patterns of shrubs directly and indirectly through interactions with other plants. However, in the current predictions, other environmental and biotic growth determinants of shrub growth have not sufficiently been taken into account, and they still remain poorly understood (see also Post and Pedersen, 2008).

A reduction of woody plant communities along increasing altitude and latitude is determined by summer temperatures and the length of the growing season. The radial growth of trees and shrubs correlate with mean summer temperatures (Walker, 1987; Knorre et al., 2006), even though other factors such as delayed snowmelt in the spring may interfere with this relationship (Kirdyanov et al., 2003; Schmidt et al., 2006). Warm air and soil temperature directly enhance the metabolism rate and the photosynthetic activity of plants, and have indirect effects through changes in soil resources such as water and nutrient availability and the depth of annual thaw (Nadelhoffer et al., 1991; Jarvis and Linder, 2000). These direct and indirect effects of temperature on vegetation are highly species and site specific (Hobbie and Chapin, 1998; Bret-Harte et al., 2002) and in most cases still remain inadequately known (Jarvis and Linder, 2000). In the Eurasian Arctic, the most common species in the shrub communities are grey-leaved willows, Salix glauca and S. lanata, the growth patterns of which remain unstudied (see, however, Walker, 1987).

Permafrost in tundra areas affects plants not only through changes in available rooting space but also through changes in hydrological, chemical, and thermal conditions of soils. Willows and other deciduous shrubs may effectively utilize areas with discontinuous or deep permafrost, where they have plenty of rooting space, available nutrients, and good drainage (Schickhoff et al., 2002; Lloyd et al., 2003; Walker et al., 2003; Schuur et al., 2007). Therefore, if the observed warming of permafrost (Osterkamp and Romanovsky, 1999; Pavlov and Moskalenko, 2002) will lead to degradation of permafrost and deepening of the active layer as predicted (ACIA 2005; Anisimov and Reneva, 2006), then shrub abundance is likely to increase. However, in predicting vegetation changes in relation to permafrost, local conditions including vegetation type, small-scale topography, snow conditions, and soil characteristics should be taken into account as they are likely to complicate the relationship (Smith and Riseborough, 1996; Walker et al., 2003).

The vertical structure of arctic willow vegetation along toposequences remains largely undocumented (see, however, Edlund and Egginton, 1984; Schickhoff et al., 2002). Despite the latitude of the overall area, wind and snow patterns create microclimate variability along slopes which is likely to affect the growth of shrubs (Seppälä, 2004). In addition, soil conditions and nutrient availability vary along toposequences thereby affecting the growth of plants. In their turn, shrubs may also change the conditions of slopes, for instance by increasing the trapping of snow (Liston et al., 2002; Sturm et al., 2005) and enhancing the turnover of nutrients (Epstein et al., 2004). On a broader geographical scale, the distance from the sea is likely to affect the growth of shrubs due to strong winds often associated with areas near arctic seas (Seppälä, 2004).

In the Eurasian Arctic, the reindeer (Rangifer tarandus) is the most important herbivore affecting the growth of willow (Chernov, 1985). Willows are of special importance for reindeer at the beginning of the growing season when herbs and graminoids are scarce and reindeer suffer from protein and mineral deficiency after the winter (Vakhtina, 1964). Reindeer have a tendency to suppress the growth of their forage shrubs and fragment their distribution in a landscape (den Herder et al., 2004, 2008; Manseau et al., 1996; Pajunen et al., 2008). This may have effects on other herbivores such as willow grouse (Lagopus lagopus) that feeds on willow during the winter (Chernov, 1985; Ims et al., 2007; den Herder et al., 2008). The effects of reindeer grazing on the growth of willow are likely to depend on abiotic factors such as climate and soil nutrient availability (Eskelinen, 2008; Post and Pedersen, 2008), and yet, these interactions remain to a large extent unknown.

I conducted a field study in the central and southern Yamal Peninsula in order to find out which environmental factors determine the growth and height of willow shrubs in a tundra region of West Siberia that is extensively grazed by large herds of reindeer. Particularly, I addressed the following questions: (1) What is the relative importance and mutual dependency between climate, soil thaw, slope position, and distance from the sea for the height and growth of willow shrubs in Arctic conditions? (2) How does reindeer grazing affect the height and growth of willow and how do the effects of reindeer interact with the above-mentioned abiotic factors and herbivory by willow grouse?

Study Area

I sampled willow thickets along a latitudinal transect spanning from the central Yamal Peninsula to the northern Polar Urals (Fig. 1). The total length of the transect was over 300 km. The sampling transect includes a gap in the southern Yamal Peninsula, because field work had to be restricted to sites that were accessible during the summers of 2006 and 2007. Therefore, the transect consists of two parts: (1) an approximately 125-km-long section in the more northern area (sites 1–9), and (2) an approximately 30-km-long section in the more southern area (sites 10–13) (Fig. 1). The transects followed the only two available access routes into the region, an all-terrain vehicle trail in the more northern area and a road in the more southern area.

Figure 1

(a and b) The study was conducted along a latitudinal transect across the Yamal Peninsula, northwestern Russia. The isotherms of July temperatures are redrawn from Shiyatov and Mazepa (1995). The arrows indicate the direction of the migration routes of Nenets reindeer herders. (c) The thermal sum indexes (sum of positive monthly mean temperatures May–September, °C) of the weather stations in the vicinity of the study transects are presented during the years 2002–2007.

The mean July temperature at the northernmost study site is 6 °C and increases approximately 1 °C/100 km in a southward direction. Along the western coast of the peninsula, the July isotherms curve southwards due to the cooling effect of the sea (Shiyatov and Mazepa, 1995) (Fig. 1). The mean annual temperature in the more northern area (Marre-Sale weather station) is −8.4 °C and in the more southern area (Salekhard weather station) −6.0 °C. The mean annual precipitation is 289.0 mm and 494.8 mm, and the mean snow depth is 21.4 cm and 35.6 cm, respectively, at more northern and more southern areas ( http://meteo.infospace.ru). The thermal sum indexes (sum of positive monthly mean temperatures May–September, °C) of the study areas during the years 2002–2007 are presented in Figure 1c.

The study areas have undergone Pleistocene glaciation, although details of glaciation history are unknown (Svendsen et al., 1999). In the western Yamal Peninsula, the oldest exposed sediments are laminated clayed silts and silty clays (Astakhov et al., 1996), which are overlain by the Kara till which is up to 20 m thick (Gataullin, 1988). These sediments are mantled by Baydarata sandy silt. This can be partly regarded as loess, analogous to the Yedoma formation of East Siberia and Alaska (Astakhov et al., 1996). This sedimentary structure is modified by various Holocene formations, such as limnic fills of thermokarst sinkholes and fluvial sediments. Podsols are normally nonexistent, but some signs of leaching occur on well-drained hills. These sites can have fairly acidic soils and frost boils also occur.

The entire sampling area is underlain by continuous permafrost; the annual depth of thaw varies between 0.3 and 1.7 m depending on temperature regimes at variable levels above and below ground. Therefore the depth of annual thaw is also affected by soil and snow conditions and vegetation (Pavlov and Moskalenko, 2002). In the Russian classification scheme (Matveyeva, 1998), Sites 1–9 belong to the typical tundra zone, whereas sites 10–13 belong to the southern tundra (Fig. 1). In another classification scheme, both of the tundra zones belong to the Low Arctic vegetation zone (Walker et al., 2005, and sources cited there). Typically, the abundance and height of shrubs is greater in the southern than in the typical tundra zone (Walker et al., 2005).

The number of semi-domesticated reindeer has increased dramatically during the past decades on the Yamal Peninsula (Baskin, 2005). The number of reindeer in the whole Yamal region (500,000 km²) was approximately 630,000 animals in 2007 ( http://www.uralpolit.ru/regions/jamal/28-06-2007/page_63081.html). Nenets nomads migrate with their animals through the study areas in several herding units (brigades) twice a summer, moving north in mid-June and back south in mid-August). When a reindeer herd moves through an area the number of animals (approximately 3000 animals occupying 100 km²) may for several days approach 30 animals km⁻² in a location (Stammler, 2005). The direction of the migration routes of Nenets reindeer herders is indicated in Figure 1.

Data Collection

The data was collected in the more northern area (sites 1–9) during the period 18 July–1 August 2006 and in the more southern area (sites 10–13) during the period 19–25 July 2007. The criteria for selection of sampling locations were easily met along the only two available access routes into the region. In both of the areas, sampling was restricted to locations that had a slope of a watershed growing willow shrubs within walking distance from a road/all-terrain vehicle trail. In the more northern area, the distance between sampling locations ranged from 5 to 27 km, whereas in the more southern area, the distance between sampling locations ranged from 7 to 15 km. The study site was the nearest slope of a watershed from a road/all terrain vehicle trail, where willow shrubs were growing.

At each study site three distinct topographical parts based on vertical position and the steepness of the slope were selected. (1) The uppermost position (hereafter called top) was at the lowest flat part of the top, under which the angle of the slope began to grow. This part of the slope was characterized by tundra heath vegetation types (Andreyashkina and Peshkova, 1995). The species composition of these tundra heath vegetation types varied depending on the moisture regime, but typically bryophytes, lichens, graminoids, and dwarf shrubs dominated. Low shrubs (mainly Salix glauca and Betula nana) occurred sporadically. The cover of willow was on average 16%. (2) The mid-slope position, (hereafter called middle), represented the steepest part of the slope. This part of the slope was characterized by shrub vegetation types (Andreyashkina and Peshkova, 1995) with ubiquitous willow cover (average cover of willow 52%). At this slope position, the understorey of vegetation was relatively sparse constituting mainly of horsetails (Equisetum arvense) and forbs (e.g., Polemonium acutiflorum, Rubus arcticus, Stellaria peduncularis) and sporadically occurring bryophytes (e.g. Sanionia uncinata). Towards the more southern study sites, willows formed a thicket with Betula nana. (3) The down slope position (hereafter called bottom) was near the end of the steep part of the slope, but on level or nearly level ground. At this part of the slope the vegetation type was either a willow thicket, as in the middle slope, or a mire (Andreyashkina and Peshkova, 1995) that was dominated by tall willows in the upper layer of vegetation and by Sphagnum mosses or sedges (Carex spp.) in the understorey of vegetation. The cover of willow in the down slope was on average 58%. This part of the slope was restricted by a river, lake, or mire in the down slope.

In each part of the study site (top, middle and bottom), 1–3 study plots 10 × 10 m in size were chosen by using a randomization procedure. To ensure that willow vegetation was present at the top sites, the sampling was restricted to sites with an overall cover of willow of at least 10%. In the middle and bottom sites, the cover of willow was always high and sampling was performed at all sites.

In each study plot, the height of 10 random willow ramets per species (Salix glauca or S. lanata) representing the prevalent uppermost canopy was measured from the ground. The top shoots with at least four years of growth increments of the measured individuals were collected for later retrospective growth measurements.

The top shoots collected in 2006 were dried and measured retrospectively in the laboratory during winter 2007. The top shoots collected in 2007 were measured in the field. In the retrospective measurements, the length of shoots was measured starting from the leading shoot and proceeding annually back to the shoot that had been growing in 2003. The living shoot and total length of the shoots were measured separately and the frost- or herbivore-damaged shoots were marked. For the analyses of length of shoots, only the measurements from non-damaged shoots were used. If the damage had occurred in later years and the actual growth for a specific year could be seen from a dead, but undamaged shoot, the result was included in the analyses. The shoot diameter was measured using cross measurements in 2003. The number of feeding marks by willow grouse (Lagopus lagopus) was counted annually between years 2003 and 2005.

The number of reindeer faeces was counted along a side of one randomly chosen plot (10 × 2 m) at each part of a study site (top, middle, bottom), and the mean value of these three transects was used in the analyses as a representative figure for the grazing intensity for that site. Faeces from previous years could not be differentiated. In addition, the depth of thaw (five measurements) was measured within each study plot.

Statistical Analyses

Prior to statistical analyses, the mean value of each response variable—the height of willow shrubs and the length and diameter of willow shoots—was calculated for each plot in order to avoid pseudoreplication. The explanatory variables were latitude (distance from the Equator [km]), distance from the sea (km), position on the slope (“top,” “middle,” “bottom”), depth of thaw (cm), the mean number of reindeer faeces on the study slope, the number of willow grouse feeding marks per ramet, and the identity of Salix species (S. lanata or S. glauca). In addition, a year (2003, 2004, 2005, or 2006) was included as an explanatory variable in the model for the length of willow shoots.

Linear mixed effects (LME) models (Pinheiro and Bates, 2000; Crawley, 2002) were used to investigate the relationships between response and explanatory variables. Since the plots within a given site were expected to be more similar to each other than to random plots drawn from the entire data set, a study site (1–13) was set as a random factor, whereas all the explanatory variables mentioned above were set as fixed factors. The quadratic terms of all the continuous explanatory variables were included as explanatory variables in the initial models. The modeling procedure followed the suggestion for multiple regression models by Crawley (2002). All the possible one-way interactions were included in the initial models, whereas two-way interactions were excluded. To obtain minimum adequate models (Crawley, 2002), non-significant terms were removed, and final models contained only significant explanatory variables or their interaction terms. In order to obtain satisfactory residual plots of the models, the data of length and diameter of willow shoots were log-transformed prior to the modeling. To graphically illustrate the results, linear or curvilinear regression models with only one variable were used to draw the prediction lines for each variable. In order to graphically illustrate the interaction between two continuous explanatory variables, some of the variables were split into three groups by subtracting the smallest value from the largest and dividing the figure by three. In this way, the response of a continuous variable could be plotted at different levels of a grouped variable. The grouping classes are described in the legends for Figures 2–Figure 34. The statistical analyses were carried out using the R statistical environment (R Development Core Team, 2005).

Figure 2

The height of willow ramets (cm) in relation to (a) the distance from the sea (km), (b) the number of reindeer faeces, and (c) the position on the slope. The lines (a)–(b) represent a fitted value for each relationship. The relationships between all these illustrated explanatory variables (and their interactions) and willow height were significant (p < 0.05) in the LME models (see text for details). In Figure 2b, the black triangle marks and the solid line indicate the plots that are characterized by deep thaw (thaw depth above or equal to 76.8 cm), and the open circles and the dashed line indicate the plots with middle thaw depth (53.3 cm ≤ thaw depth ≤ 76.7 cm). The open triangles and the dotted line indicate the plots with low thaw (thaw depth less than or equal to 53.2 cm). In the boxplot Figure 2c, the middle line represents the median, and the other two lines the lower and upper quartiles. The extreme values are marked as separate dots if they go beyond 1.5 times the interquartile, which is the maximum length of the line. The abbreviation B indicates bottom; M middle; and T the top part of the slope.

Figure 3

The length of willow shoots (cm) in relation to (a) the distance from the Equator (km), (b) the distance from the sea (km), (c) the depth of thaw (cm), and (d) the number of reindeer faeces. The lines (a)–(d) represent fitted value for each relationship. The relationships between all these illustrated explanatory variables (and their interactions) and the length of willow shoots were significant (p < 0.05) in the LME models (see text for details). In Figures 3a and 3c, the triangle plots and the solid line indicate Salix lanata and the open circles and the dashed line indicate S. glauca. The dot-dash line indicates the quadratic term of distance from the Equator. In Figure 3b, the black triangle marks and the solid line indicate the plots that have deep thaw (thaw depth above or equal to 76.8 cm), and the open circles and the dashed line indicate the plots with middle thaw depth (53.3 cm ≤ thaw depth ≤ 76.7 cm). The open triangles and the dotted line indicate the plots with low thaw (thaw depth less than or equal to 53.2 cm). The dot-dash line indicates the quadratic effect of distance from the sea. In the Figure 3d, the black triangle marks and the solid line indicate the plots at the bottom of the slope, the open circles and the dashed line plots at the middle of the slope, and the open triangles and the dotted line plots at the top of the slope.

Figure 4

The diameter of willow shoots (mm) in relation to (a) the distance from the Equator (km), (b) the depth of thaw (cm), (c) the number of reindeer faeces, and (d) the number of willow grouse feeding marks. The lines (a)–(d) represent the fitted value for each relationship. The relationships between all these illustrated explanatory variables (and their interactions) and the diameter of willow shoots were significant (p < 0.05) in the LME models (see text for details). In Figure 4c, the black triangles and the solid line indicate the plots close to the sea (23.7 km or less from the sea) and the open circles and the dashed line indicate the plots with middle distance from the sea (23.8 km ≤ distance from the sea ≤ 41.2 km). The open triangles and the dotted line indicate the plots far away from the sea (41.3 km or more from the sea). In Figure 4d, the black triangle marks and the solid line indicate high reindeer pressure (the mean number of faeces greater than or equal to 4.9), the circles and the dashed line indicate medium reindeer pressure (2.7 ≤ the mean number of faeces on the slope ≤ 4.8), and the open triangles and the dotted line indicate low reindeer pressure (the mean number of faeces on the slope less than or equal to 2.6).

Results

The height of willow ramets linearly increased as a response to increasing distance from the sea (sea: F_1,46  =  4.9, p  =  0.03; Fig. 2a). The effect of reindeer grazing on the height of willow ramets depended on the depth of thaw (grazing × thaw: F_1,46  =  4.6, p  =  0.04). In the areas with deep or middle thaw, the increase in the grazing intensity decreased the height of willow ramets, whereas at low thaw depths, the increase in the grazing intensity increased the height of willow ramets (Fig. 2b). Willow ramets were highest in the middle of the slope and lowest at the top of the slope (slope: F_2,46  =  12.8, p < 0.0001, Fig. 2c). The effects of explanatory variables did not interact with the identity of Salix species, but the difference of average ramet length of Salix lanata (74.1 ± 32.6 cm) and S. glauca (45.2 ± 22.3 cm) was consistent (species: F_1,46  =  19.6, p < 0.0001).

The length of willow shoots curvilinearly decreased northwards (latitude ˆ 2: F_1,8  =  6.4, p  =  0.04; Fig. 3a). The effect of latitude on the length of willow shoots also varied between Salix species (latitude × species: F_1,226  =  10.0, p  =  0.002): the length shoot of Salix glauca declined more strongly northwards than the length shoot of S. lanata (Fig. 3a). There was a significant concave curvilinear relationship between the distance from the sea and the length of willow shoots (sea ˆ 2: F_1,226  =  7.0, p  =  0.009). From the coast until approximately 30 km inland, the curve showed a decline in the length shoot as a response to increasing distance from the sea, whereas from 30 km onwards, the curve showed an increase in length shoot as a response to increasing distance from the sea (Fig. 3b). The effect of distance from the sea on the length of willow shoots depended on the depth of thaw (sea × thaw: F_1,226  =  5.7, p  =  0.02). At shallower thaw depths, where the shoot length was smaller than at the greater thaw depths, there was a strong increase in shoot length as a response to increasing distance from the sea (Fig. 3b). In the areas of middle thaw depth, the length of willow shoots slightly increased as a response to increasing distance from the sea, and in the areas with deep thaw, the shoot length remained nearly constant irrespective of the distance from the sea (Fig. 3b). The length of willow shoots linearly increased as a response to increasing depth of thaw. For Salix glauca this effect was slightly stronger than for S. lanata (thaw × species: F_1,226  =  3.8, p  =  0.05) (Fig. 3c). The length of willow shoots was negatively related to reindeer density, and this effect depended on the position on the slope (grazing × slope: F_2,226  =  6.0, p  =  0.003). At the middle and the bottom of the slope, where the shoot length was greater than at the top of the slope, the negative effect of grazing was stronger than at the top of the slope (Fig. 3d). In 2005, the length of willow shoots was smaller (2.8 ± 1.3 cm) than in the other investigated years (2003: 3.8 ± 1.3 cm; 2004: 3.8 ± 1.4 cm; 2006: 3.9 ± 1.5 cm) (year: F_3,226  =  22.5, p  =  0.0001).

The diameter of willow shoots in 2003 decreased northwards (latitude: F_1,9  =  451.0, p < 0.0001) (Fig. 4a) and increased as a response to increasing depth of thaw (thaw: F_1,46  =  24.3, p < 0.0001) (Fig. 4b). The increasing intensity of reindeer grazing decreased the diameter of willow shoots and this effect depended on the distance from the sea (grazing × sea: F_1,46  =  22.5, p < 0.0001). In the areas far away from the sea, where the diameter of willow shoots was greatest, the response of shoot diameter to increasing grazing pressure was more clear than in the areas closer to the sea. The effect of willow grouse on the shoot diameter of willow depended on the intensity of reindeer grazing (grazing × willow grouse: F_1,46  =  6.8.3, p  =  0.01). In the areas of high or low intensity of reindeer grazing, shoot diameter increased with increasing number of willow grouse feeding marks, whereas in the areas of medium intensity of grazing, there was a negative association between shoot diameter and willow grouse feeding (Fig. 4d). The diameter of willow shoots was smaller at the top of the slope (2.3 ± 0.7 mm) than at the other two slope positions (middle: 3.0 ± 1.1 mm; bottom: 3.1 ± 1.0 mm) (slope: F_2,46  =  27.6, p < 0.0001). The effects of explanatory variables did not interact with the identity of Salix species, but the difference between shoot diameter of Salix lanata (3.4 ± 1.0 cm) and S. glauca (2.3 ± 0.7 cm) was constant (species: F_1,46  =  167.3, p < 0.0001).

Discussion