Multivariate characterisation of environmental conditions for reindeer husbandry in Sweden

Pastoralism using semi-domesticated reindeer (Rangifer t. tarandus) is a traditional livelihood in northern Fennoscandia. The area used for reindeer herding in Sweden covers as much as half of the country’s area. Variation in the productivity of reindeer husbandry is clearly affected by many biotic and abiotic factors. The aim of this investigation was to identify factor combinations which describe the spatial variation in conditions that plausibly determine productivity in reindeer herding. Initially, 37 variables representing geographical location, climate, weather episodes related to ice crust formation and insect harassment, topography, vegetation, forage abundance and qualities, and fragmentation of the ranges were derived, using prior ecological knowledge and spatially explicit data. The variables were mapped in a raster of 1958 squares of 100 km each, covering the entire reindeer herding area in Sweden. Reductions of variables were performed with multivariate analyses in steps, ultimately retaining 15 variables. The first five principal components of these variables explained 84% of the total variation. The first component, related to major western mountain/eastern lowland gradients, already accounts for 49% of the variation. The following components explained variation ranging from 10% to 5.4%, and revealed spatial patterns in summer versus winter forage, climatic conditions and ice crust formation, abundance of forests and winter forage, and northward slopes together with valuable forest areas, respectively. A tentative zone division of the Swedish reindeer herding area into seven zones was made, based upon cluster analysis and spatial distribution of component scores. Extending this approach and method seems useful also in the understanding and management of other natural resources and national parks, especially with an ongoing global climate change perspective.


Introduction
Subspecies of Rangifer (reindeer and caribou) occur naturally in arctic, subarctic and some temperate regions of the northern hemisphere.These latitudes are characterised by high spatial and temporal diversity in environmental conditions, with varying season lengths depending on location and year.The physiological and behavioural adaptations of Rangifer typically include marked seasonal rhythms in several physiological functions (e.g., van Oort et al., 2000) and various energy-conserving capabilities.The adaptations also include considerable ability to utilise body depots as a buffer between seasons (Tyler & Blix, 1990).Furthermore, Rangifer show spectacular migratory behaviours to exploit beneficial seasonal habitats.
Abstract: Pastoralism using semi-domesticated reindeer (Rangifer t. tarandus) is a traditional livelihood in northern Fennoscandia.The area used for reindeer herding in Sweden covers as much as half of the country's area.Variation in the productivity of reindeer husbandry is clearly affected by many biotic and abiotic factors.The aim of this investigation was to identify factor combinations which describe the spatial variation in conditions that plausibly determine productivity in reindeer herding.Initially, 37 variables representing geographical location, climate, weather episodes related to ice crust formation and insect harassment, topography, vegetation, forage abundance and qualities, and fragmentation of the ranges were derived, using prior ecological knowledge and spatially explicit data.The variables were mapped in a raster of 1958 squares of 100 km 2 each, covering the entire reindeer herding area in Sweden.Reductions of variables were performed with multivariate analyses in steps, ultimately retaining 15 variables.The first five principal components of these variables explained 84% of the total variation.The first component, related to major western mountain/eastern lowland gradients, already accounts for 49% of the variation.The following components explained variation ranging from 10% to 5.4%, and revealed spatial patterns in summer versus winter forage, climatic conditions and ice crust formation, abundance of forests and winter forage, and northward slopes together with valuable forest areas, respectively.A tentative zone division of the Swedish reindeer herding area into seven zones was made, based upon cluster analysis and spatial distribution of component scores.Extending this approach and method seems useful also in the understanding and management of other natural resources and national parks, especially with an ongoing global climate change perspective.
Key words: Cluster analysis, GIS, landscape ecology, PCA, Rangifer tarandus tarandus.Rangifer,27 (1): 5-23 rely on the animals' adaptations to the natural conditions of the respective areas and are generally parts of indigenous land uses evolved from different kinds of previous subsistence uses of the land (e.g., Lundmark, 1998;Forbes et al., 2006).
The Fennoscandian reindeer herding area covers most of the land north of lat.61˚30'N in Norway and Sweden, and north of lat.64˚30 'N in Finland (Fig. 1A;County Administration Boards, 2000).The Swedish part comprises approximately 200 000 km 2 , which is about 50% of the land area of Sweden.The area is divided into 51 reindeer herding districts, used by associated herding communities, and consists of more than 900 reindeer husbandry enterprises altogether.The ownership structure, slaughter strategies, and economy vary significantly among the herding communities (Statistics Sweden, 1999;Swedish Board of Agriculture, 2005).As part of the management, animal densities are regulated by harvest to a level well below the ecological carrying capacity, thereby keeping the herd in a state of rapid growth with the harvestable surplus near its maximum level of sustained yield (generally corresponding to 1/4 to 1/2 of the winter herd, Statistics Sweden, 1999).In spite of the management measures taken by the herders, which vary between districts, both per area and per capita reindeer productivity vary considerably between years as well as between different parts of the area.The number of yearly slaughtered animals for the entire Swedish reindeer herding area has varied between 29 000 and 112 000 during the last 30 years, according to Statistics Sweden (1999;2005).
In Sweden, the reindeer herding area (Fig. 1A) includes strong differences in climate and weather patterns due to its latitudinal position and a marked west-east gradient from oceanic to continental climate types.A topographical gradient is prevalent from the alpine mountain summer range in the west to forested lowland winter areas in the centre and eastern parts.The latitudinal extent, in combination with marked longitudinal climatic and topographic differences, induces a large variation in the area in such aspects as  2000).Solid arrows in the map indicate the present spring migrations (autumn migrations are in the reverse direction) and broken arrows indicate historical migrations according to Aarseth (1996).Rangifer, 27 (1), 2007 season length, vegetation, and snow conditions, with profound effects on herbivores.Extensive road nets, development of infrastructure and other industrial land uses heavily influence the coastal and adjacent inland areas in particular and potentially induce constraints on reindeer behaviour and movement (e.g., Wolfe et al., 2000;Nellemann et al., 2003;Helle & Hyppönen, 2004), but habituation to disturbances has also been reported (e.g., Klein, 1991;Colman et al., 2001).Within these meta-scale gradients a considerable variation in regional and local conditions is clearly conceivable, but patterns are not easily distinguished from the overriding meta-patterns.The variation in reindeer productivity is likely attributable to this complex web of many correlated and interacting factors.These interact across scales ranging from multiyear and regional levels down to the animals' daily ranges.An understanding of the area-specific spatial and temporal pattern of these would greatly improve the management possibilities.However, so far it has been difficult to uncover them in sufficient detail due to the multitude of conceivable interrelated factors.Tangible and agreed-upon understandings of these patterns are not available.Investigations on compound factors on a larger scale are increasingly of interest when facing a global climate change.According to the Intergovernmental Panel on Climate Change (IPCC, 2001b) and Arctic Climate Impact Assessment (ACIA, 2004), the impact of climate change will be especially pronounced in the arctic with consequently severe effects on ecosystems and societal structures of indigenous people.
In this study we used multivariate ordination techniques to characterise the broad-scale variation in environmental conditions that potentially contributes to regional productivity in reindeer husbandry in Sweden.Based on causal relationships between external factors and reindeer productivity compiled from ecological knowledge of Rangifer in literature, 37 variables were derived from meteorological and geographical databases and reduced to 15 variables by principal component analyses.These results were used for spatial description of major characteristics and a tentative zonation of the herding area.

Factors with relevance for reindeer productivity
Reindeer feed on a broad spectrum of vascular plants, mushrooms and lichen (e.g., Llano, 1956;Skuncke, 1958;Klein, 1990;Warenberg et al., 1997).The snow-free season provides the nutritional surplus, which can be converted to growth and production (e.g., White, 1983), and the snow season represents a period of scarcity of forage where reindeer often lose weight.The forage needs to be lush with low levels of cellulose and high levels of nitrogen (Skogland, 1984;Hofmann, 1989;Danell et al., 1994).The availability of suitable seasonal ranges with high-quality forage is hence important.
Topographical features, such as altitude variation, ruggedness, slope and aspect, interact with climate and weather features and affect primary production (e.g., Raven et al., 1998;Krebs, 2001) and several other conditions essential for large herbivores.Snow patches, often found in troughs and on north-facing slopes and remaining over the summer season, are used by reindeer for insect relief (Anderson & Nilssen, 1998;Hagemoen & Reimers, 2002).Ruggedness is considered to extend and diversify the melting time of snow, with a similar effect on plant phenology (Nellemann & Thomsen, 1994).A variation in slopes and aspects seems positively related to body weight in red deer (Mysterud et al., 2001) and is therefore suggested to be favourable also for reindeer.
Deep snow has been reported to negatively influence the percentage of pregnant females in reindeer (Ropstad et al., 2001;Kumpula & Colpaert, 2003).Increased energy costs to obtain forage in deep snow, and delayed melting of snow cover on calving grounds, may explain these results (Thing, 1977;Crête & Payette, 1990).Ice crust and deep firn are critical since they make the necessary ground lichen forage inaccessible (Solberg et al., 2001;Putkonen & Roe, 2003), and force reindeer to switch to arboreal lichens (Johnson et al., 2001).The ice crust/firn formation is determined by complex interactions between altitude, air and ground temperature, snow density, wind velocity, humidity, air pressure, rain-and snowfall, and radiation fluxes (e.g., Colbeck, 1989;Craven & Allison, 1998).
Weather conditions during summer directly affect the nutrient composition and digestibility of green forage, where heat and drought in particular may decrease forage quality (e.g., van Soest, 1994;Lenart et al., 2002).Summer climate and weather also affect the abundance of harassing insects, which reduce the effective summer grazing time for the reindeer, with negative effects on their autumn body condition and the productivity of the reindeer (Downes et al., 1986;Mörschel & Klein, 1997;Colman et al., 2003).

Conceptual basis
Based on the scientific knowledge on Rangifer ecology and local knowledge often shared by herders in the context of reindeer husbandry, a conceptual framework for the environmental influence on reindeer productivity was suggested (Fig. 2).In this framework, broad-scale environmental conditions connect to productivity in reindeer husbandry primarily via the animals' energy and nutrient budgets.These in turn affect growth and body mass and thereby indirectly affect vital life history traits, which determine the demography of the herd and thus the productivity in terms of harvestable surplus.
The environmentally caused variation in energy and nutrient budgets was conceptually attributed to variation in four major components: balance between the winter and snow-free seasons, reachability (accessibility) of forage in cost-benefit terms, available daily foraging time affecting the forage intake rate, and the energetic costs as a consequence of external circumstances.Behind these, a parsimonious set of exogenous driving factors includes latitudinal location, topographical characteristics, climate features, and weather variation.A minimum set of relevant secondary factors may include season lengths, vegetation cover and composition with respect to forage vegetation, snow conditions, and insect harassment, with particular reference to their impact on reindeer foraging.Anthropogenic factors with broad-scale effects add to this and may be attributed to forestry and infrastructural developments in general.

Study area
All reindeer grazing land used by Swedish reindeer husbandry within the borders of Sweden was included in the study.To achieve a sufficient number of observations the area was divided into 1958 squares (Fig. 1A), each square covering 100 km 2 and to which variable values were assigned.Incomplete grid squares along borders and areas prescribed for exclusive use by Norwegian reindeer husbandry according to the reindeer grazing convention between Sweden and Norway (Swedish Code of Statutes (SFS), 1972) were excluded.Areas available for Swedish reindeer husbandry in Norway according to the same convention and an area in the northernmost part of Sweden were left out due to incomplete data.The investigated squares cover 86% of the total area encompassed by the outer borders of the Swedish reindeer herding area, which includes the convention areas disposed in Norway.

Data
The analyses were based on spatially explicit data on topography, vegetation cover, climate, weather, forest age and the road network.Temporally variable data were included either as 10-or 30-year averages (climate and weather-related variables) or as the recent situation.
The topographic information in a grid format with 50-m spatial resolution (DEM50) was acquired from National Land Survey of Sweden (2002a).The climate information was collected from the National Atlas of Sweden (National Land Survey of Sweden, 2002b).Daily weather information was obtained from the Swedish Meteorological and Hydrological Institute (SMHI, 2004).Reindeer forage classes (Table 1) were obtained from the reindeer husbandry database REN2000 (County Administration Boards, 2000), which uses data from satellite imagery and the vegetation maps from National Land Survey of Sweden and the Swedish Environmental Protection Agency.Vegetation data with 25-m spatial resolution were obtained from the National Land Survey of Sweden, GSD -Land Cover Data (2004).Data on old forests was obtained from a survey of northern Sweden made available by the Swedish Environmental Protection Agency (Jacobson et al., 2002).Road data to account for road fragmentation was obtained from GSD -Blue Map (National Land Survey of Sweden, 1998).All cartographic data was handled with ArcGIS™ software (ESRI, 2005).
Altogether 37 variables (Table 2), representing the conceptual background factors for the energy and nutrient retrieval shown in Fig. 2, were derived and assigned to each grid square from the available data.In order to obtain a better fit to the normal distribution, the variables were, when necessary, Box-Cox transformed (Box & Cox, 1964), with the exponent rounded to the nearest 0.5.The variables were standardized to zero means and unity standard deviations.For several factors a number of variants were developed to achieve a varied selection of variables with the intention of retaining those that described as much variation as possible, while being multivariate normally distributed.Some variables are considered non-linear in comparison to productivity but when searching for correlations between variables through PCA, non-linearity against productivity is not an issue.In this study we only investigated the product-moment correlations between the independent variables without considering the order of causality amongst them.Productivity is a dependent variable in this context and was not included (not available in this scale).Latitude (Lat) and longitude (Long) for the centre point of each grid square were derived in order to account for the general north-south and east-west gradients of the data.
Eight variables characterizing topographical features (elevation, slope, proportion of north-facing slopes, and five ruggedness indices) were derived from the DEM50 data.Water bodies larger than 2500 m 2 were excluded in these calculations in order to avoid bias due to water surfaces.Slope was derived as the mean slope (without considering the direction) within each square, where the slope was the angle of the steepest slope of a plane fitted to a 3 x 3 grid cell area around each 50 m x 50 m elevation grid cell.The north-facing slope variable (SlopeN) was arbitrarily classified as the proportion of slopes with bearings between 330 and 30 degrees within each square, based on the slope directions for each elevation grid cell.The five topographic ruggedness indices were variable variants measuring the amount of small-scale variability in local topography.RugTRI was derived from the terrain ruggedness index, as originally done by Riley et al. (1999), but with a slight modification.The values assigned to each square were the means of the elevation grid point TRI raw values, prior to the classification of the values used by Riley et al.These raw values were calculated by taking the square root of the sum of squared differences in elevation between a centre cell (x 22 ) and its 8 neighbour cells (x ij ; where i, j = 1, 2, 3; (i, j) ≠ (2, 2)).To reduce skewness of the distribution of RugTRI, the values were transformed to natural logarithms after adding 1 in order to avoid logarithms of zeros: Other ruggedness indices included the mean grid point value of each square after the classification according to Riley et al. (RugTRIclass), the average of the elevation range (RugFocRng), the elevation standard deviation (RugFocStd), and the mean elevation difference of the centre cells and their eight surrounding grid cells (RugFocMean).
Climate variables included were growth season length (GrowthSea), snow season length (SnowSea) and snow precipitation (SnowPrec), all derived by isoline interpolation from the climate data.The values used were mean values for the 30-year period from 1961 to 1990, with a spatial grain resolution of 1000 m (1 km 2 per grid cell).GrowthSea was computed as the number of days with average daily temperatures above +5 ˚C and SnowSea as the average number of days between the first and last day of snow cover.SnowPrec was derived as the product of mean annual precipitation and the mean proportion of annual snow precipitation.
Spatial variation in probabilities of icing or ice crust formation was captured with three variables derived from available daily weather observations of temperature and precipitation over a 10-year period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) at 58 weather stations (Fig. 1A; SMHI, 2004).An ice crust probability index was first derived based on four different constellations of daily maximum, minimum and average temperatures and precipitations, as described in Table 3.If any of these events occurred at a weather station on a certain day, the day was regarded as having an ice crust formation incidence equal to 1, otherwise 0. Three ice crust probability grids were constructed by interpolating the average number of such days per year at each weather station by using ordinary kriging (IceOrdKr), and with elevation as covariate (Todini & Ferraresi, 1996;Goovaerts, 1998;Goovaerts, 2000) cokriging (IceOrdCokr) and disjunctive cokriging (IceDisCokr), since the ice crust showed negative correlation with elevation.Each 100-km 2 square was given the average value of 10-year average probability values of its 2500 200-m grid cells.The square values ranged from 7 to 35 days of high ice crust probability per year.
Daily data on wind speed and temperature, measured every third hour from 58 meteorological stations over a 10-year period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002), were used to calculate cumulative warble fly activity indices during the snowfree season.Two insect harassment probability grids were based on a model for warble fly activity suggested by Mörschel (1999).This model projects the probable insect activity level using wind speed (m/s) and temperature (°C) data in a logistic regression setting: The two insect activity grids were extrapolated for the entire herding area from the weather station values using ordinary cokriging with elevation as a covariate.MI12h was based on the accumulated sum of all measurements available for each 12-hour daytime period between 6 a.m. and 6 p.m., and MI24h was based on 24 hours a day.
Vegetation and forage-related variables included seven variables capturing forest area, forest age structure and lichen abundance as indication of winter forage availability, four variables related to forage during the snow-free season spring, summer and autumn (SSA), as well as two variables capturing variation in fragmentation of linear structures on habitats, and six indices combining the latter two aspects in a cost-benefit concept called reachability (Lundqvist, 2007).
Table 3. Description of weather events assumed to create a high probability of ice crust formation.

Event Method
Rational 1 3 days (average temp > 1 ˚C ) followed by 2 days (average temp < -1 ˚C ) Snow melts and gets packed with increasing density when temperature is above 0 °C.Since snow has high temperature capacitivity due to its high water content, high temperatures are needed for a longer time.The same is the case for low temperatures.Probably a thicker ice crust is formed during this event than during the following three events.
2 Day when max temp > 2 ˚C min temp < -2 ˚C max temp -min temp > 10 ˚C When irradiation and albedo are high, usually during days with clear skies, the snow surface melts and freezes, which causes ice crust formation, likely leading to a shallow ice crust.

3
Day when max temp > 2 ˚C min temp < -2 ˚C precipitation > 3 mm Wet precipitation that freezes directly or shortly after reaching the ground and creates an ice crust, which also can lead to firn layers in the snow pack.

4
Day when average temp > 1 ˚C precipitation > 3 mm, followed by day with average temp < -1 ˚C As above but instead of using minimum and maximum temperatures of the same day, average temperatures of two following days are used to include days of decreasing temperatures with precipitation on the first warmer day.Rangifer, 27 (1), 2007 Forest structure variables included, with relevance for the winter grazing resources, were the amount of total forest (Forest), lichen-rich forest (ForestL), clear-cut areas (ClCut) and proportion of uncut forest (UncutFor) against total forest including clear-cut areas.Two variants of mature forest indices (MatFor) and (MatFor-Marsh) capturing variation in availability of connected mature forest areas were derived from classified satellite imagery (Jacobson et al., 2002), where stands older than approx.50-70 years were defined as old forest.
MatFor showed the relative occurrence of connected old-growth forests, quantified as the proportion of old forest within a 400-m radius from each cell.MatFor-Marsh also included marshes as a connecting vegetation type.A more direct winter forage variable (ForageW) was derived by summing all areas with coniferous and birch forests of heath type with lichen, as well as dry and extremely dry heaths.These vegetation types are all classified as "very good" winter habitats in the REN2000 database.
The four variables related to forage quality during the snow-free seasons (ForageSSA, ForageSSAvg, Forage-SSA3w and ForageSSA4w) were developed from spatial data from vegetation databases and classifications of vegetation types (Table 1).The assigned vegetation values were arbitrary values based on the ranked forage value classification of reindeer forage in REN2000.All forage variables were expressed as area-weighted sums of the classified vegetation type values within the 100-km 2 squares.
Two indices describing density of fragmenting linear structures such as trails, roads and railroads were constructed by measuring the length of these linear structures in each 100-km 2 square.LinStrLgh included all structures with equal weights, while LinStrWgh arbitrarily weighted the linear structures according to size (Table 4), similarly as in the following reachability indices.The weights are multipliers to the linear structure lengths of each structure type.
Snow-free season forage abundance, quality and costs to reach the forage were combined in six reachability indices (Lundqvist, 2007).The reachability indices used the SSA forage values shown in Table 1 divided by the cost of movement, reflecting energy demands and animal behaviour when moving over terrain or linear structures to reach the forage from evenly distributed sample points 2 km apart.Hence, the forage values of all grid cells are divided by the cost of reaching them from predetermined starting points, traversing both ground and fragmenting structures such as roads, and thereafter summed in each observation square.The cost of movement (friction value) in terrain was set to unity per km, while friction values on water bodies were arbitrarily set to friction value 10 per km and different linear structures in terrain were given friction values as listed in Table 4. Edge effects, capturing the gradual forage quality decline close to forage patch edges, were arbitrarily set to 75 m according to an estimated average of line of sight.All friction values were chosen according to calibration studies on summer ranges of the Handölsdalen reindeer herding district (Lundqvist, 2007), which covered an area of 3400 km 2 .Differential costs of movement due to topography and vegetation were hence not included.The values assigned to each 100-km 2 square were the sums of the grid reachability values.Two indices used a 3-class forage classification including fragmenting linear structures; the first index (Re3wEd) included edge effects in opposition to the second index (Re3w-Noed).The following two indices used the 4-class forage classification and included fragmenting linear structures, with (Re4wEd) and without (Re4wNoed) edge effects respectively, and the last two indices used the 4-class forage classification excluding linear structures with (Re4wNoroEd) and without (Re4wNoro-Noed) edge effects, respectively.

Analyses
To reduce the number of variables, the correlation matrix of the 37 variables was inspected and analysed with principal component analyses (PCA) in steps and cluster analyses (CA) with average linkage, using SIMCA-P (Umetrics, 2004) and SAS software (SAS Institute Inc., 2005).Principal components (PCs) were extracted according to the following commonly suggested stopping criteria (e.g., Hair et al., 1998;McGarigal et al., 2000): size of successive eigenvalues of the correlation matrix (the Kaiser-Guttman criterion), inspection of a Cattell scree plot, and comparison with an average scree plot based on 10 000 simulations of PCA for independent random variables.Correlated variables and variables describing similar factors (indices) were analysed as subsets in further PCAs.One variable per biologically relevant factor was retained by choosing those that described appreciable variation and had a normal-like and continuous distribution over the observation squares.To distinguish similar variables showing resembling weights in the PCA, its underlying data resolution and ecological knowledge of the variable's effect on reindeer biology guided the decisions of which variable to retain.Fifteen variables were retained and as a verification of sufficiency, scatter plots of pair-wise component scores from the PCs based on 37 variables, versus the reduced set of 15 variables, were finally examined to investigate the variation alteration caused by the reduction of variables.As an aid in the interpretation of the extracted PCs, the spatial distributions of the PC scores were visualized with score maps.
To suggest interpretations of the PCs, individual variable loadings were evaluated against the critical value 0.13, which is the doubled (as suggested by Hair et al., 1998) critical value for correlation coefficients at experiment-wise testing of 15 estimates with 1958 observations at the significance level 0.05 (Rohlf & Sokal, 1995;Sokal & Rohlf, 1995).
Finally, a zone division into 7 zones was suggested based on CA (Ward's method, standardized distance) on the 1958 observations and the 15 retained variables, and the characteristics of each zone were compiled in a table.

Results
The 15 retained variables are indicated in bold in Table 2.The latitude and longitude variables were excluded in the final synthesis in order to avoid confusion with other related gradients and artefacts due to the correlation caused by the SW-NE extension of the study area.This exclusion then also allows for a better biological interpretation, which was what we strived for.Each variable were, as expected, spatially autocorrelated according to 'Moran I' spatial autocorrelation test (e.g., Legendre, 1993;ESRI, 2005).Elevation and north-facing slopes (SlopeN) were moderately correlated (r = 0.34) and were both retained in the final study.All five ruggedness indices and slope were highly correlated (0.86 ≤ r ≤ 1.00), and among them we retained the raw TRI variable (RugTRI) due to its high resolution and ability to capture saddle areas and evenly sloping planes.Growth and snow season were negatively correlated and were retained together with snow precipitation.Among the highly correlated ice crust variables (0.91 ≤ r ≤ 0.98) we retained the ordinary cokriging index (IceOrdCokr, mean = 19.51,S.E.= 5.69) since it showed higher variation among observed squares than variables based on other interpolation techniques.The two insect harassment indices MI12h and MI24h behaved identically (r = 1.00), but the index covering 12 hours a day (MI12h, mean = 102.8,S.E.= 24.43)showed higher accumulated communality in the PCA and was therefore retained.ClCut, describing the amount of clear-cuts, was excluded due to its large proportion of zeros caused by many observations outside the boreal forest area.The variable UncutFor better represented the amount of forestry impacts since it measures existing forest against clear-cuts and was therefore retained.The mature forest variables (MatFor and MatForMarsh) were subsequently left out due to strong correlations with Forest (r = 0.80 and 0.84), as well as due to their low spatial resolution (400 m grain).The road density variables LinStrLgh and LinStrWgh were highly correlated (r = 0.88).The weighted road density variable LinStrWgh was chosen over the unweighted variable due to the different impact of the various linear structures, e.g.track vs. railroads, on reindeer ecology.Among the vegetation indices, ForageSSA4w was retained due to its higher vegetation class resolution in comparison with the alternative indices and because it had a better fit to the normal distribution.All reachability indices were highly correlated (0.57 ≤ r ≤ 1.00).Re4wEd was retained as it showed large variation over the area and was assumed to capture a larger part of utilisable forage variation in the landscape.
A dendrogram from a CA of the 15 variables is shown in Fig. 3, where the variables clearly appeared in two subclusters.The first subcluster included variables with high values towards the north and east, and the second subcluster included variables with increasing values towards the south and east.
The eigenvalues of the first five PCs of the retained variable set were 7.3, 2.0, 1.5, 0.96 and 0.81, explaining 49%, 13%, 10%, 6.4% and 5.4% of the total variation, respectively (Fig. 4).A pronounced 'elbow' in the scree plot occurred already at the second PC and a second 'elbow' appeared at the fourth PC.The sample scree plot intersected the average of 10 000 randomly generated scree plots just before the fourth PC.These criteria supported extraction of four PCs.The fifth eigenvalue was 0.81, which is above the limit of 0.7 recommended by Jolliffe (1972).This criterion supported the extraction of five PCs.
Scatter plots of pair-wise PC scores (Fig. 5), based on 15 variables versus 37 variables, indicated good agreement for the first three PCs (r = 0.96, 0.87 and 0.88, respectively).Lower correlations were found for the fourth and fifth PCs (r = 0.54 and 0.37), suggesting altered compositions of these components.Reducing further to 10 variables (not shown) resulted in lower score correlations for 2 nd to 5 th PC between the reduced 10-variable set and the original 37 variables set (r = 0.98, 0.01, 0.04, 0.44 and 0.37, respectively).The grey dotted lines represent the factor loading value (±0.32), which accounts for 10% of the variance in the component (e.g., McGarigal et al., 2000).Loadings with absolute values exceeding 0.13 are considered statistically significant.
Fig. 6 shows that all 15 variables had strong loadings in PC 1; they were divided into two sets of contrasting variables conceptually connected to summer and winter conditions.The mapped scores reveal a clear east-west (i.e., lowland-mountain) gradient captured by the component.This macro-scale variation explains the first elbow in the scree plot.
The second PC contained four variables with strong loadings, Re4wEd and ForageSSA4w in opposition to ForageW and ForestL, i.e. summer versus winter forage resources.The score map identified rich summer forage, particularly in the central western and the easternmost parts of the area and along several river valleys.Winter forage resources were particularly identified in the central northern and southernmost parts.A more scattered distribution was apparent in the north.
In the third PC, the absolute loadings were all higher than the critical value of 0.13 (Fig. 6).The component contrasted two mixed groups of variables with highest emphasis on a combination of variables related to winter and high elevation terrain versus forest-and forage-related variables.The score map emphasised the highest elevated eastern and western parts, which were not particularly emphasised in the previous components.Four of the eight positively loaded variables fit known climatic patterns, e.g., warmer winters and longer growth season in the south.
Table 5. Characteristics of the tentative zones based on averages of each zone compared to the total average of the retained 15 variables.The zone division is shown in Fig. 7.

Zone Characteristics
A Very short growth season with very high snow precipitation.The ice-crust probability is low as well as the harassing insect activity.Highly elevated and rugged with high abundance of winter forage, but low lichen abundance.Very low amounts of roads, forest and clear-cuts.
B Short growth season with high snow precipitation.Fairly elevated and rugged with low probability of harassing insect activity.Fairly good winter forage but very low amounts of lichen and low amounts of clearcuts and roads.

C
Fairly short growth season with low snow precipitation and ice-crust probability.Fairly low ruggedness with low amounts of winter forage but rather high abundance of lichen.Fairly good summer forage and low amounts of roads.

D
Rather low elevation and ruggedness.Fairly low amounts of summer forage.Very good winter forage and very high amounts of lichen.Fair amount of forest, clear-cuts and roads.

E
Fairly long growth season with low snow precipitation.Fairly high probability of ice-crust formation and harassing insect activity.Low elevation and ruggedness.Very low amounts of winter forage and fairly low amounts of lichen.Good summer forage, many roads and a fair amount of forests and clear-cuts.

F
Long growth season with low snow precipitation.Very high probability of ice-crust formation and harassing insect activity.Very low elevation, very low amounts of winter forage and fairly low amounts of lichen.Fairly good summer forage, much forest and clear-cuts and abundant road net.

G
Fairly elevated and low probability of harassing insect activity.Extremely good winter forage with very high amounts of lichen but very low amounts of summer forage.Fair amount of forest.

Fig. 7.
A tentative zone division of the Swedish reindeer husbandry area, based on a generalization of a cluster analysis of 15 variables (n = 1958).The characteristics of each tentative zone are described in Table 5.
Rangifer, 27 (1), 2007 In the fourth PC, eight variables had loadings higher than 0.13, with the highest emphasis on north-sloping and forested areas.High scores coincided particularly with high-altitude forests, which stretch from southern mountainous areas to the areas near mountains in the north.
PC 5 contrasted one variable (SlopeN) with strong loading versus a mixed group of particularly winter-and forest-related variables.The score map for PC 5 appeared to identify a large-scale gradient between the northeastern part (not emphasised in the other PCs) and the rest of the area.The mean accumulated communality of all variables for the first five PCs was 0.84, and ranged for individual variables from 0.66 (UncutFor) to 0.97 (SlopeN).
A suggested zone division of the Swedish reindeer herding area based on a CA of the 15 retained variables is shown in Fig. 7 and the characteristics of the zones are listed in Table 5.

Discussion
The available scientific knowledge, as well as knowledge shared by herders, suggests a complex causal background for variation in the productivity of reindeer husbandry, e.g.due to landscape properties.To achieve good productivity, each herding district needs to be sufficiently complete as regards availability of forage resources and other vital conditions.The different districts have their own strengths and weaknesses, which are known only in qualitative terms.They also have very different shapes, extensions and strategies for their land use, which reflects the variation in suitability and capabilities for reindeer husbandry, as well as adaptation to the conditions on a whole-year basis (see Fig. 1).This tangled web of conditions is not possible to disintegrate spatially using productivity data from the 51 reindeer herding districts due to the low number of productivity observations, and the large data grain the productivity figures represent.Furthermore, a possible confounding effect between conditions and management measures (including a long-since evolved herding district division and different degrees of traditional subsistence strategies vs. adoption of modern production strategies such as calf slaughter), aiming to benefit from strengths and counteract deficiencies, is present in the system.Therefore, productivity per se is not an unambiguous indication of variation in quality of external conditions.
In this study therefore we attempted to identify the multidimensional spatial pattern of positive and negative external conditions shown to be important in literature.Here, we only focus on landscape effects, and try to quantify these effects, leaving other explanations as 'noise'.When the number of factors and variables is large and spatially correlated, multivariate ordination is commonly used.Therefore, it appeared to be a suitable methodological choice to identify patterns of underlying compound variation in conditions.Finding causal relationships including many variables in this coarse scale for the large study area is not an option.Using PCA when variables are correlated and for variable reduction had the advantage of handling all variables simultaneously and with no pre-ranking order.The opportunity of using a large set of observation squares (n = 1958) gave a considerable statistical strength to our analyses.Spatial autocorrelation for each variable was apparent.However, the estimated correlations between variables are not affected by autocorrelation, because the biases of the variance and covariance estimates are compensated when the quotient for the correlation is calculated.Furthermore, the aim of this study was to investigate the total area for variable reduction, and hence no dependent variable was included.
Several alternative variables were derived to capture the variation of the same factor, and therefore showed high correlation and collinearity.Hence only one of the variables describing the same factor was retained after reduction, as seen in Table 1.This suggests that the effect of collinearity of the original variable set was one of the factors explaining the fairly small difference between the scores of the different variables sets shown in Fig. 5.A specific scale-related bias is inherited in the interpolated variables using data from the 58 weather stations, where interpolation errors depend on the kriging techniques used.These inherited biases are dependant on distances to weather stations.When so they are transferred to the PCs and may have affected the zonation.In this large-scale approach and based on relative comparisons, the biases were however thought to be inferior, especially since the correlations between different interpolation techniques were remarkably high.
Although the initial variable set was large, it is still fairly limited and other variables could have been included.For example, predator densities were excluded due to their temporary nature and the difficulties of fitting these into the observation scale.Furthermore, we lack such data explicit in time and space, since available predator data are based on predation compensation figures, which are available only on the herding district scale.Snow depth was another excluded factor as not all weather stations record snow depth.Further, snow depths depend on several meteorological, topographical and vegetational factors, which lead us to conclude that an interpolation of snow depths would be poor because snow depth is often determined by snow drift and wind, and not only precipitation as such.Rangifer, 27 (1), 2007 Some of the variables, mainly climatologic and meteorological, are temporally variable.These, however, were included only as mean values in order to limit the initial set of variables, although variability measures could have been additional options.The variation of some variables, e.g.insect harassment and ice crust probabilities, were expected to be strongly correlated with the interannual means.For other variables, such as season lengths, the variation is expected to vary independently of the mean of the variables.Including temporal variation for the latter would have increased the initial set by about 1/4 and the reduced set by up to 1/3 more variables and was therefore disregarded in order to gain parsimoniousness.Large variation in temporal varying factors is supposed to have impact on reindeer productivity.Because the complex web of interacting factors is so intricate, the temporal variation itself cannot, however, be declared exclusively positive or negative for reindeer productivity.
The major gradients of conditions, i.e. the east-west gradient connected to major altitudinal and oceaniccontinental climate gradients and the north-south one connected to latitudinal climate gradients (Fig. 4a), were expected a priori.As suggested in Fig. 3, these major trends showed a positive correlation between the topographical variables and snow season length, snow precipitation as well as a negative correlation with forestry intensity (i.e.positively correlated with UncutFor).Elevated and rugged areas inherit low temperatures with more snowfall, and are areas of little interest to forestry due to its low productivity.The second large cluster of Fig. 3 is connected to warmer climate and higher productivity, which are found in the south and east.Winter forage (ForageW) was a somewhat ambiguous variable since it was correlated with the variables of the first large cluster in Fig. 3 containing snow season, and also highly correlated with lichen-rich forest (ForestL) of the other large cluster.These variables were mainly apparent in the central area between the two extremes in the topographical gradient.This exemplifies a deficiency with clustering techniques where multidimensional correlations are summarized in one single graph.However, these gradients, which are commonly exploited by migratory reindeer herding (Fig. 1B), fell out strongly as the outstandingly dominant PC.
Above these major gradients PC 2-PC 5 show considerable irregularities, which have not previously been phrased and documented other than in the form of incoherent local knowledge.Considering PC 1 as a way to remove self-evident meta-scale gradients, PC 2-PC 5 seem quite informative; of the remaining component variation after fitting PC 1, PC 2-PC 5 explained 35%.
The second PC clearly identified a variation between more prominent green forage resources in the southern half of the western mountain chain (the Scandes) and in the northeastern inland at fairly low altitudes.This is in contrast to other areas with higher abundance of winter forage, and describes a more mixed appearance above the main pattern displayed with PC 1.
Similarly, PC 3 revealed more difficult winter conditions elsewhere than in the barren alpine areas in the west, captured in PC 1. Harsh snow conditions were apparent in the southeastern part combined with longer growth season and a rather rugged landscape with high degree of infrastructure.In addition, areas with high amount of forage and forest were apparent in the central and northeastern parts.
PC 4 in particular identified contrasting richer forest resources with high lichen abundance in a band from the central inland in the north towards the southwest part of the Scandes (likely extending over the Norwegian border).The occurrence of northward slopes, which could differentiate a delayed snowmelt, was a pronounced variable, as it also was in PC 5.
Besides identifying contrasting high occurrence in the northernmost part, PC 5, together with PC 2, seems to recognize the northwest-southeast direction of river valleys.The difference of PC 4 and PC 5 between the complete and reduced sets of variables is probably an effect of concealing variation when including several highly correlated variables, as was the case in the complete set, which might have transferred this variation to later PCs.
In a global climate change perspective, where the global mean temperature is expected to increase between 1.4 -5.8 ˚C between the years 1990 -2100 (IPCC, 2001b), the Nordic region is believed to achieve a higher temperature increase than the global average with considerable changes in precipitation and hydrology, especially during winters (IPCC, 2001a;ACIA, 2004;SMHI, 2007).Suggested hydrological changes are decreased spring floods, increased autumn, winter and annual runoff, and increased frequency of highflow events during autumn (Andreasson et al., 2004).The variables correlated with warmer climate, such as growth season length, forage plant physiology, harassing insect activity, forestry and infrastructure, may therefore possibly increase and affect reindeer ecology (e.g., Gunn & Skogland, 1997;Lenart et al., 2002).The shorter winter season and increased precipitation, most likely rain, may affect the conditions for reindeer husbandry in several aspects, e.g.winter areas being more prone to severe thawing and icing events.Altogether, both positive and negative effects on a whole year basis could be foreseen.Therefore the total impact on the reindeer industry of a global warming is difficult to evaluate and predict, but Rangifer, 27 (1), 2007 changes in the compound condition and usage of natural forage and compensating areas should be expected.This, however, needs further analyses based on the data we used in this study.
As previously mentioned, the data grain in these analyses is not directly compatible with reindeer productivity measures, which therefore cannot readily be used for verification of the results.It was necessary to adopt this finer grain in order to enable fitting of the conceptual set of factors suggested in literature, and also to uncover variation in a relevant sub-seasonal scale in reindeer husbandry.Instead the tentative zone division, based on CA, provided an informal test of the results by comparing with the actual use patterns in practical reindeer herding, thereby also indirectly verifying the interpretations of retained PCs.These use patterns have evolved as traditional land use over a period of centuries and are probably based on the wild reindeer's use of these areas.The long-term intensity of use has been stable for more than a century (besides smaller variations due to regular reindeer fluctuations of about ± 25% in cycles of 20-30 years; Statistics Sweden, 1999), and may therefore be considered a settled practice.The husbandry depends upon diversified ranges for different seasons and back-up ranges when severe weather events occur, e.g.ice-crust formation.Lack of such ranges has effect on productivity and may cause existing ranges being overgrazed when such weather episodes occur.As is summarised in Table 5, particularly zone E, and in part C and F are suitable for summer use while zones A, D, G and in part B contain suitable winter forage, but are not necessarily suitable for winter use due to meteorological reasons.Stationary forest reindeer herding communities are in fact found in, or stretches between, D and the northern parts of E and F, and hence achieve complementing ranges throughout the year.To achieve suitable all-season conditions in the rest of the herding area, exceptionally long migrations (200-350 km per direction) are performed by migratory mountain herding communities, in fact moving through the stationary districts largely along the D-E boundary (Fig. 1B).Furthermore, the score distributions indicate deficits of suitable summer ranges in the northernmost part of the reindeer husbandry area.The solution in practice for this area has been to seek summer grazing in Norway.Formerly the migrations extended right to the Norwegian coastline, but in the southern part they are currently hindered by fences along the national border due to political conventions between the countries.Some Swedish reindeer herding districts in the far north have the opportunity to use some areas in Norway according to the convention.
In conclusion, this study allowed for a first screening of areas used for reindeer herding in Sweden in a multi-variate perspective based on documented animals' demands.The area characteristics seem to qualitatively agree with the current use of the Swedish reindeer husbandry area on larger scales.Thus, a new administrative zone division of the land resources and conditions, supporting, for example, comparisons between individual herding communities, might be based both on results from the multivariate analyses and on the current use patterns.These divisions will likely not follow the county borders on which the administrative management decisions are based today.The extensions of the suggested zones may also be a matter of revision in a climate change perspective and when new data and analyses become available.Further zonation needs, however, to be done in a resolution where the compound conditions for each herding district and season are merged.This is possible with clustering techniques based on the reduced set of variables derived in this paper applied on the seasonal ranges of the herding districts.The results in this paper open up possibilities to analyse conditions together with productivity measures.Furthermore, a largescale multidimensional approach such as this could be applicable on other resources and species as well, thereby providing helpful information on landscape characteristics, e.g. in nature conservation planning and administration.

Fig. 1 .
Fig. 1.A) The Fennoscandian reindeer herding area with the imposed raster of 1958 squares, each covering 100 km 2 , and the distribution of the meteorological stations used.Uncovered areas within the herding area are incomplete edge squares, excluded areas for Swedish reindeer husbandry use or areas with missing/incomplete data.B) The 51 Swedish reindeer herding districts and main reindeer migration routes (County Administration Boards,2000).Solid arrows in the map indicate the present spring migrations (autumn migrations are in the reverse direction) and broken arrows indicate historical migrations according toAarseth (1996).

Fig. 2 .
Fig. 2. Conceptual overview of factors suggested to affect reindeer productivity.

Fig. 3 .
Fig. 3. Dendrogram from cluster analysis with the average linkage method of the 15 retained variables (normalized values).

Fig. 4 .
Fig. 4. Scree plot showing the variance described by the first 10 PCs in a PCA including the 15 retained variables (solid line).Dotted line shows average scree plot of 10 000 simulations on 15 uncorrelated variables, and 1958 observations.

Fig. 5 .
Fig.5.Score correlations (n = 1958)  between the original variable set of 37 variables and the reduced set of 15 variables for the first five PCs.

Fig. 6 .
Fig. 6.Geographical distribution of PC scores of 1958 squares based on the 15 retained variables in the first five principal components.Maximum and minimum score values of the PCs are in italics inside the legends.The loadings of the 15 variables in each PC are shown on the vertical axis.The grey dotted lines represent the factor loading value (±0.32), which accounts for 10% of the variance in the component (e.g.,McGarigal et al., 2000).Loadings with absolute values exceeding 0.13 are considered statistically significant.

Table 1 .
Vegetative classification and forage quality classification according to REN2000 (County Administration Boards 2000).Forage values are alternative weightings of vegetation types in the computations.

Table 2 .
Characteristics of the initial set of 37 variables.The variables are grouped according to their exogenous background as conceptualised in Fig.2.The 15 retained variables after variable reduction are in bold.

Table 4 .
Weights used for weighting linear structures in the fragmentation indices and as friction values in the reachability indices.