Relationships between biotic and abiotic range characteristics and productivity of reindeer husbandry in Sweden

Reindeer husbandry is a form of pastoralism where vast areas are used as forage ranges throughout the year. The productivity of the reindeer industry in Sweden is affected by a multitude of factors on different geographical and temporal scales. Our aim was to find combinations of factors characterizing the environmental conditions for reindeer husbandry in the 51 herding districts in Sweden, which correlate strongly with variations in productivity both between herding districts in general and between years within districts. Productivities were described by estimated herd growth rates and carcass condition of slaughtered females and calves. These dependent variables were related to the environmental independent variables using linear regression models and structural equation modelling. The independent variables were either considered as stable (e.g. topography, vegetation and infrastructure) or temporally changing (e.g. season lengths, weather events, disturbances and animal slaughter strategies). The most relevant independent variables were included in a cluster analysis to suggest a grouping of herding districts based on similarities in environmental conditions. Considerably larger variation in productivity was found between herding districts than between years. Different variables were found to be important for between-district and within-district variations, respectively. Season lengths and animal densities were found significant at both levels of variation. Other variables found to be relevant were ruggedness, snow condition, harassing insect activity, supplementary feeding, calf slaughter ratio and previousyear animal condition. Snow precipitation, ice-crust formation and forage quality were presumed to be relevant for reindeer productivity, but were not found to have a large impact on productivity. These factors, however, may have been counteracted by husbandry measures, statistically incorporated in animal density variables, or by being strongly correlated with other, more significant variables. Several of the variables that were found to be important for productivity are correlated with climate and weather and therefore foreseen to be altered in a climatic change perspective.


Introduction
The reindeer and the industry Reindeer (Rangifer tarandus) is a migratory cervid well adapted to arctic climatic conditions.Domesticated reindeer (Rangifer t. tarandus) are managed in reindeer herding systems in the northern circumpolar areas.In Sweden, the reindeer herding area covers approximately 200 000 km 2 (County Administration Boards, 2000).The reindeer herding industry in Sweden includes about 940 reindeer enterprises with approximately 4600 reindeer owners, organised in 51 herding communities (SJV, 2005b).The grazing lands are divided into herding districts with legally defined year-round and winter ranges.In Sweden, the average winter stock of semi-domesticated reindeer has been approximately 225 000 animals and has fluctuated between 150 000 and 300 000 in cycles of 25 to 30 years over the past 120 years.The an-Rangifer, 29 (1): 1-24 nual harvest during the last 12 years has been on average 55 000 animals, where calf slaughter accounts for approximately 60% (Statistics Sweden, 1999;SJV, 2005a) officially registered for subsidies purposes.Thus, the economy of reindeer husbandry is strained and the majority of reindeer owners are usually obliged to have additional sources of income (Statistics Sweden, 1999;Riseth, 2006).
Large carnivores preying on reindeer are wolverines, lynx, wolves, golden eagles and brown bears (e.g., Boman, 1995;Nybakk et al., 1999;Pedersen et al., 1999;Nybakk et al., 2002;Persson, 2003).The red fox has also increasingly affected reindeer as predator due to its recent northbound dispersion (Elmhagen, 2003).The animal losses due to depredation, excluding the red fox, have officially been estimated to be 20 000 to 30 000 animals per year (Swedish Government Official Reports (SOU), 1999).These figures are probably underestimated as fresh estimations based on kill rates and recent official investigations suggest a current yearly mean predation loss of approximately 40 000-45 000 reindeer (Ö.Danell, unpubl.).Of these animals, around 25% to 30% were assumed to be adult females, which in turn may have brought up calves.The fragmenting effects on ranges of large carnivores and the cost in animal condition have not yet been estimated.
The 51 reindeer herding districts in Sweden differ in environmental conditions for reindeer husbandry, structure and organisation, and vary in productivity (Statistics Sweden, 1999;SJV, 2005a;b).The differences include a large number of interacting factors, many of which are associated with the properties of the ranges.The majority (n = 33) of the herding districts (called mountain herding districts) have oblong northwest-southeast shapes, following the direction of large rivers and fitting natural migration routes between mountainous summer areas in the west and winter ranges in the coastal and inland areas in the east/southeast.
The remaining districts (n = 18) have more equilateral extents and practice a more stationary type of reindeer husbandry (forest and concession herding districts).The summer ranges of the mountain herding districts are in some cases partly on the Norwegian side of the western national border.In earlier times, the reindeer migrations extended to the Atlantic coast in the west, but are now restricted by the Swedish-Norwegian country border and political conventions (Svensk-Norska Renbeteskommissionen, 2001).The herding districts, and thereby the herds, in Norway and Sweden are hence separated.The mountain ridge dividing the countries also acts as a climate border with more continental climate in Sweden compared to Norway, and the policies regarding large carnivore are very different in the two countries.To keep the herds in rapid growth, the animal densities are regulated by harvest to levels well below the ecological carrying capacities of the available land.The slaughter strategies differ between herding districts as e.g.many of the herders practice calf slaughter.The majority of the slaughtered calves are males, since higher recruitments of female calves are needed.The different management strategies and measures affect reindeer herding productivity in addition to the underlying environmental factors.
The aim of the present study was to quantify the relationship between environmental variables and reindeer productivity for the total Swedish reindeer herding area.A second aim was to suggest a division of herding districts based on the factors found relevant for productivity.The analyses make use of annual slaughter and herd statistics from the entire Swedish reindeer herding area for the derivation of productivity measurements.These include herd growth, carcass weights, fatness, and conformation classifications.The produc-tivity measurements were related to range data previously derived by Lundqvist et al. (2007) on the basis of literature on biotic and abiotic factors suggested to affect reindeer productivity (Lundqvist, 2003).The environmental factors include season lengths, vegetation, topographical ruggedness, snow precipitation, ice-crust formation, forestry clear-cuts, and infrastructure such as roads.In addition, the present study includes management strategies such as animal densities and slaughter strategies, and animal losses due to depredation.

Data
All variables included in the study are shown in Table 1.The variables were divided into dependent variables describing productivity and independent variables describing the abiotic and biotic factors that affect reindeer productivity.The independent variables were further divided into stable variables (not changing between years in the time range of this study) and yearly shifting variables.
To derive data on reindeer productivity, we used annual herding district and slaughter statistics from 1994 to 2004 (SJV, 2005a).Statistics on inventories of lynx, wolverines and wolves from the years 1994 to 2004 (SP, 2005) were used to estimate the reindeer losses due to predation, which were incorporated in the herd growth calculations.Digital maps on topography, vegetation, roads and snow precipitation were obtained from the National Land Survey of Sweden.Weather data were obtained from the Swedish Meteorological and Hydrological Institute.Digital maps on herding districts and seasonal ranges were obtained from the Reindeer Husbandry Database REN2000.Spatial calculations and analyses were made with ArcGIS geographic information system software (ESRI, 2005) and all statistical analyses were done using SAS statistical software (SAS Institute Inc., 2004).

Productivity variables
Reindeer productivity can be described by e.g.herd production of meat per area unit, herd growth (finite growth rate), and animal body condition (carcass weights, fat depots and conformation).Meat production is a combined quantification of herd growth rate and carcass weights, and is strongly influenced by slaughter strategies and management measures.This measure was therefore omitted in this study.Herd growth rate is partly influenced by slaughter strategies and management, but judged to be less influenced by management.Body condition is relatively unaffected by management, but related to environmental conditions and resources (e.g., Helle & Kojola, 1994;Kumpula et al., 1998;Eilertsen et al., 2001;Pettorelli et al., 2005;Reimers et al., 2005).Hence, to obtain a diverse set of productivity variables that are affected as little as possible by management measures, we used the finite growth rate together with body condition measurements at slaughter (carcass weight, fatness and conformation).
The data used to derive productivity variables included distribution of reindeer herding enterprises and owners, herd sizes, herd and age gender structure, and slaughter statistics including animal category (female, male, calf), carcass weights, fatness and conformation (EUROP) classifications (Swedish Board of Agriculture's Code of Statutes (SJVFS), 2002).In addition, calf slaughter percentage, animal densities and the herders' own outtakes, i.e. the number of slaughtered animals for personal consumption which are not registered in the slaughter database, were estimated from the statistics on slaughter and herding districts.
For each reindeer herding district, the finite annual herd growth rate (rG) was computed from data on animal numbers and estimated depredation losses adjusted for compensatory mortality using assumed natural survivals where adult animals are assumed to survive in a higher degree than calves (Ö.Danell, unpubl.): where N t = winter herd size in year t H t = number of slaughtered reindeer in year t E = estimated yearly own outtake (numbers of animals) P t = estimated number of reindeer killed by predators in year t S = assumed finite survival (proportion of N) each year (adults 0.98 and calves 0.70, weighted mean based on herd age structure).
The estimated own outtake is 190 kg meat per reindeer husbandry enterprise and year (Statistics Sweden, 1999).With 947 enterprises in Sweden, the total own outtake per year becomes approximately 180 000 kg per year.
With 4600 reindeer owners in Sweden, this corresponds to 39.1 kg per owner and year.Since the number of enterprises per herding district and the number of reindeer owners per enterprise vary substantially between regions, the mean of estimated yearly outtake in number of animals per herding district was calculated as: where n 1 = number of enterprises in each herding district n 2 = number of reindeer owners in each herding district W = average slaughter weight in kg Depredation losses for each herding district and year were estimated from inventories of lynx, wolverines and wolves, used in the economic compensation calculations for preda-tor-killed reindeer (SP, 2005).Depredation estimates were based on recent compilations of predation kill rates reported in literature and official investigations relevant for Scandinavian conditions, which suggest a current yearly mean loss due to predation of approximately 40 000 to 45 000 reindeer, as mentioned above, excluding secondary effects on herd productivity via altered herd structure.The inventories of lynx, wolverines and wolves are registered as rejuvenations, regular presence and temporary presence, which represent different predation rates.When applied on individual herding districts and years, average predation rates were derived from the overall depredation and similar predation rates were used per individual of all predator species.The estimations also included golden eagle and brown bear, but the estimated losses due to these species are proportional to the area of each herding districts, due to insufficient inventories, and were therefore unaltered between years, as in SP ( 2005).
Animal condition was measured as carcass weights, and fatness and conformation classifications.These variables were derived for adult females (FemW, FemFat, FemForm) and calves (CalfW, CalfFat, CalfForm) for each herding district and year (Table 1).The fatness and conformation data are ordinal classifications done by the slaughterhouses.Prior to the statistical analyses, the fatness and conformation classifications were transformed to quasinormally distributed variables by deriving the expected value of the underlying variables of each class using a threshold model.The score of class i was obtained as within each animal category, where ϕ is the normal density function evaluated at thresholds i-1 and i, and Ф is the cumulative normal distribution function evaluated at the same thresholds.

Environmental variables
The environmental variables were, as mentioned, defined as either stable or temporally altering.However, some variables, which actually vary between years, were defined as stable, as data were not available in sufficiently fine scales.In order to enable the use of stable and temporally varying variables in the same analyses, small normally distributed random errors were added to all stable variables to obtain a better approximation to the normal distribution.Hence singularities, such as non-invertible matrices, in the calculations were also avoided.These errors were set to a percentage of each variable's standard deviation based on assumed measurement errors and rigidity of the variables.The topographical variables were given errors of 1%, infrastructure variables 2% and the vegetation variables 3% of their standard deviation.The errors are hence minute and do not change the means of the stable variables, but enable simultaneous inclusion of all variables in the analyses.
The independent variables related to range conditions were derived in previous studies (see Table 1 and Lundqvist et al. (2007) for further details).Variables relevant during the vegetation-growing season were derived for the spring, summer and autumn ranges (SSA), and the variables relevant during the snow season were derived for the winter ranges (W).To measure the quality instead of the amount of each environmental factor or resource, the means of range variables for each herding district were used instead of sums of range variables.In some cases seasonal ranges overlap within or between herding districts.In such cases, the shared areas were included in the observations of several seasons or districts.
The temporally varying variables were derived for each of the years 1994 to 2004 (SJV, 2005a;SMHI, 2005).The season length variables were modified in comparison with the corresponding variables in Lundqvist et al. (2007).The growing season length (GrowthS) was defined as the number of days with average temperature above 5 °C and only applied on the SSA ranges.The winter season length (WinterS) was set to the number of days with an average temperature below 0 °C and applied on the winter ranges.
Seven new varying independent variables were introduced in this study: Calf slaughter, animal densities in summer and winter, animal condition lag, temperatures in August, number of days with deep snow, and pre-slaughter supplementary feeding.Calf slaughter ratio (CalfR) was derived by dividing the number of slaughtered calves by the total number of slaughtered animals.Animal densities on winter ranges (DensW) were calculated using the yearly winter herd size divided by the district's winter range area.Shared areas were evenly divided between the herding districts involved, to compensate for the overlap.The reindeer densities on the SSA ranges (DensSSAr) were derived using the sum of the winter densities (N t ) and calculated herd growth (rG × N t ) divided by the reachability sum (indicated by r in variable name), i.e.ReachSSA × area, for each SSA range.To account for the effect of previous-year animal condition on current-year production, an arbitrary construct (WFF_1: Weight, Fatness and conFormation with oneyear lag) was derived according to the formula: where n = the numbers of years (= 10) included in the variable.
To avoid mean animal weights being affected by different adult vs. calf ratios, the means of ratio adjusted calf and female weights were used.Carcass fatness and conformation scores were included as total means of both females and calves for each herding district and year.These means were weighted by 0.2 after having tested alternatives between 0 and 1, where 0.2 seemed like a well balanced alternative as it showed highest variation in a principal component analysis of all alternative weightings.Lower or higher coefficients resulted in a too high correlation with any of its parents.High temperatures during summer and autumn are reckoned to be negative for reindeer condition, probably due to increased insect harassment, low content of nutrients in forage, and intense thermoregulation (e.g., Hagemoen & Reimers, 2002;Kumpula & Colpaert, 2003).A variable, the temperature average in August, was hence included in the study.The number of days per year with snow depths over 0.75 m was included to measure the effect of deep snow on reindeer productivity, as digging and walking through deep snow is energetically costly.The variable was derived by counting the number of days with a snow depth deeper than 0.75 m for the 21 weather stations in the reindeer herding area that measure snow depth.The variables August temperatures (TempAug) and days of deep snow (SnowDepth) were spatially interpolated over the reindeer herding area using ordinary kriging technique with prediction output, according to earlier variables extracted (Lundqvist, 2007;Lundqvist et al., 2007).Kriging was also chosen as it aims at minimizing the variance of the errors of the weighted linear combinations while it is unbiased as the mean of the errors are zero.
Supplementary feeding before slaughter is performed in some areas to decrease radioactive isotopes in reindeer due to the fallout of the Chernobyl accident in 1986 (County Administration Boards, 2008).Since the bio-logical half-life of caesium 137 in ruminants is not more than a few weeks (Åhman, 1996), the supplementary feeding prior to slaughter are enough to achieve animal meat according to food standards.However, this supplementary feeding is assumed to effect carcass body weight and fat content, hence the estimated percentage of fed animals prior to slaughter of each herding district and year (SupFeed) was included in the study.This applies mostly to the southern central herding districts in the counties of Jämtland (Z) and Västerbotten (AC).

Statistical analyses
Initial stepwise linear model analyses (LM) of each productivity variable were performed to test the effects of environmental variables on each productivity variable.First, all environmental variables were included simultaneously with Year as class variable in order to account for the average temporal variation in the variables and when estimating effects of stable and altering variables on district productivities.In a second set of analyses, only the temporally varying environmental variables were included with Herding district as class variable.This was done in order to estimate effects of yearly (interannual) variations in these variables on the productivity, while adjusting for constant characteristics of the herding districts.
In order to test hypotheses about the causality structure of the system, we used Structural Equation Modelling (SEM) (see e.g.Hair et al., 1998;Everitt & Dunn, 2001;Hung et al., 2007), where path diagrams describing hypothetical causal trees were constructed, tested and adjusted.SEM comprises latent variables, i.e. variables not measured directly but estimated from other measured variables.Ordinary SEM is, however, unable to handle class variables such as Herding district and Year.This could be done with multi-level SEM, but in this case the observations are too few to support an appropriate multilevel SEM approach.Therefore, SEM was applied on both between-district and interannual within-district covariance matrices.A multivariate one-way analysis of variance was used to estimate the covariance matrices for the two levels, whose correlation counterparts were used as input to the SEM analysis.The number of observations was originally 561 as we had data from 51 herding district over 11 years.Due to the inclusion of previous year animal numbers in rG, there were 10 complete years included all productivity variables, resulting in 510 observations.Hence, the number of degrees of freedom in the SEM analyses was 50 in the between-district model and 457 in the between-year model, the latter being adjusted for two missing values.
Through confirmatory factor analyses, two measurement models for analyses of differences between districts and between years, respectively, were set up.Measurement models are used to investigate the fit of the constructed latent variables in relation to the indicator (measured) variables.The latent variables were subsequently revised in order to better fit the indicator variables, guided by Lagrange's and Wald's tests in the 'Calis' procedure in the SAS statistical software (SAS Institute Inc., 2004).Thereafter, the measurement models were modified to represent the theoretical models of interest, i.e. the proposed models, which were defined after examinations of the results in previous statistical analyses.The proposed models were also tested and adjusted (correlation paths and variables dropped or moved) based on indications from Lagrange's and Wald's tests, hence retaining only significant path correlations in the final models.
The herding district means of the identified main variables affecting productivity variation between herding districts were finally used to cluster the herding districts into a relevant number of groups.Cluster analyses with average linkage and Ward's method were used according to the work sequence of Lundqvist & Danell (2007), and the resulting groups were shown in a map.

Results
The linear model (LM) analyses showed highly significant effects of years (Table 2) and districts (Table 3), respectively, on the dependent productivity variables.
In the analyses with Year as class variable, the environmental variables affected each productivity variable differently (Table 2).Variation between herding districts in herd growth (rG) was mainly explained by topography, growing season length, snow depth, harassing insect activity and calf slaughter ratio.Variation between districts in calf fat and conformation variables were mainly explained by forestry intensity (ClearCut), road density, insect activity, animal densities, animal condition the previous year (WFF_1), and supplementary feeding.The female fat and conformation variables were affected mainly by August temperatures, SSA animal density, animal condition previous year, and supplementary feeding.Snow pre- cipitation, reachability, winter season length, and ice-crust formation were not found to significantly affect the dependent productivity variables.
No general distinctions between the herd growth, the calf variables and the female variables were seen in the analyses with Herding district as class variable (Table 3).Herd growth (rG) was mainly affected by animal density, insect activity and calf slaughter ratio.Growing season length, snow depth, insect activity, SSA animal density and supplementary feeding mainly affected the weight variables.Season lengths, August temperatures, winter animal density and supplementary feeding mainly affected the fatness variables.The calf and female conformation variables were affected by different variables.Insect activity and temperatures in August affected calf conformation while snow depth and supplementary feeding affected the females.Neither annual variation in ice-crust occurrences nor animal condition lag seemed to affect the productivity variables.
In the SEM analyses, regarding the variation between districts, the productivity variables were, through tests of measurement models, used to describe the three endogenous latent variables Herd production, Calf condi tion and Female condition (Fig. 1).The primary forcing explanatory variables describing Season lengths and Topography were related to the two exogenous latent variables SSA Primary forces and Winter primary forces, which were allowed to affect the four endogenous latent variables: Summer forage, Climate stress, Winter forage and Snow conditions.These four endogenous latent variables were expressed by 1 -3 manifest variables, as seen in Fig. 1.Furthermore, three exogenous latent variables related to management (Slaughter strategy, Animal density and Animal feeding) were constructed.Both winter and growing season animal densities expressed the Animal density variable.In the proposed model, the variable Animal condition lag (WFF_1) was left out since the SEM between herding districts are in a sense based on variable averages over the years and hence become highly autocorrelated.
In the initial model, all exogenous latent variables were allowed to correlate.Female condition was allowed to affect Calf condition, which was allowed to affect Herd production.Slaughter strategy was allowed to affect only Herd produc tion, since calf slaughter ratio is not supposed to affect animal condition.Animal feed ing was only allowed to affect animal condition variables since the feeding are performed on animals chosen for slaughter only and therefore not affecting reproducing animals.The path coefficients for latent variables describing only one manifest variable were set to unity.Moreover, the path coefficients between the manifest variable showing the strongest relation with its latent variable in the preliminary runs of the model were also set to unity in order to stabilize the model and obtain a standardized scale of all estimated correlations and covariances.
In the results of the proposed between-district model, several correlations between exogenous latent variables, path coefficients and manifest variables were suggested to be set to null by the Wald's test or to be moved to another position by the Lagrange's test.This was done in a stepwise procedure, and subsequently the manifest variables ReachSSA and Win terF were nullified as suggested by Wald's test.These independent variables had rather small Fig. 2. Estimated parameters of the final structural equation model describing variation between herding districts.Significant (P ≤ 0.10) path coefficients and covariances between manifest variables (white), exogenous latent variables (green), and endogenous latent variables (orange) are shown as solid arrows while dashed lines were used for the less significant path between Slaughter strat egy and Herd Production.Fixed coefficients are marked with (x) and correlations between exogenous latent variables are shown as curved double-head arrows.
impact in the LM.Furthermore, the latent variables Summer forage, SSA Primary forces and Climate stress were joined into one (SSA condi tion).The latent variables Winter forage, Winter Primary forces and Snow condition were similarly joined into Winter condition.Several covariances between the exogenous latent variables were set to zero.The causal paths between the latent variables of seasonal conditions (Winter condition and SSA condition) and Female condition, as well as, between Female condition and Calf condition were insignificant and therefore set to null.The causal paths between Animal density and Calf condition, as well as Herd production, were also dropped.Slaughter strategy was suggested by Wald's test to be dropped, probably due to its high correlation with SSA condition.Due to its significant impact in the LM analysis the variable was retained to indicate its effect.The final model, describing the variation structure between herding districts with remaining significant variables (P ≤ 0.10 in this model), non-zero covariances and path coefficients, is shown in Fig. 2. The latent variables Slaughter strategy, SSA condition and Animal feeding were positively correlated with each other, and their path coefficients to Herd production and Animal conditions were positive, except from SSA condi tion to Calf condition.The path coefficient from Winter condition to Calf condition and from Ani mal density to Female condition were negative, sug-gesting negative impacts of these variables on animal condition.
T h r o u g h o u t the model reduction, the difference of chisquare (χ 2 ) for the revised model and the theoretical unconstrained model was too high compared to the difference in degrees of freedom between the models, which suggests that the model was unsuccessful in accounting for the relationships between the latent constructs according to advises in e.g.Hair et al. (1998).Other indices to validate the model are shown in Table 4, both for the model between districts and years.Lower, or close to zero, values of the root mean square error of approximation (RMSEA) index, and higher, or close to unity, values of the comparative fit index (CFI) and the non-normed fit index (NNFI) indices suggest improved model fit.All indices suggest poor fit of the models in this study.
In the measurement models describing variation within herding districts (i.e. with focus on interannual variation), only the eleven varying causal variables were used, described by the seven exogenous latent variables Slaughter stra tegy, Animal condition lag, Season length, Summer condition, Snow condition, Herd density and Feeding.All exogenous latent variables were allowed to correlate initially.The productivity variables were used to describe three exogenous latent variables as in the previous model, capturing the variation between herding districts.The measurement model did, however, suggest one latent variable to describe each condition type, i.e. one for each of the variables weight, fatness and conformation.Therefore, four latent Table 4. Goodness-of-fit indices for the final structural equation models describing between-district and between-year variations.For definitions of fit indices, see SAS Statistical Software (SAS Institute Inc., 2004).

Index
Between herding districts (Fig. 2) Between years (Fig. 4 For the proposed interannual model (Fig. 3), Wald's test suggested to subsequently drop the manifest variables TempAug and IceCrust.Furthermore, the latent variables Condition lag and Slaugh ter strategy did not significantly affect the productivity variables and therefore dropped.Lagrange's test suggested to merge the latent variables Season length and Summer condition, which was done accordingly into a new latent variable called Climate conditions, as seen in Fig. 4. All correlations between exogenous latent variables except between Snow condition and Cli mate conditions and Animal density, respectively, were not significant and hence dropped.This was Fig. 3. Proposed structural equation model used to explore the relationship structure between environmental, constructed latent and productivity variables, within reindeer herding districts (focus on interannual variation).Measured manifest variables are showed in white, exogenous explanatory latent variables in green and endogenous productivity latent variables in orange.Fixed correlations are marked with (x) and correlation between exogenous latent variables are drawn in grey.done as previously in a stepwise procedure to detect the consequences of the changes in the model by every adjustment done.
The path coefficients between the latent animal condition variables and their manifest variables in the final model (Fig. 4) were all positive, but especially high from Animal fat ness, suggesting a high correlation between females and calves in this variable.The path coefficients from the latent variable Climate conditions to Animal weight and Animal fatness were positive.The path coefficients from the latent variable Animal density to Animal weight were negative, but positive to Herd production.The Animal density did not seem to affect Ani mal conformation or Animal fatness.Snow condition did, however, affect Animal weight negatively.Animal feeding significantly affected all animal conditions positively.Animal weight positively affected Animal conformation and Animal fatness.Animal fatness was the only animal condition latent variable affecting Herd production and did so positively.Animal density and Snow condition correlated negatively.The fit indices of the interannual model suggested poor fit according to standards given in literature, but somewhat better than the between-district model.
Cluster analyses (CA) were used to group the 51 reindeer herding district based upon the most relevant productivity determinants.Among all variables identified to affect reindeer productivity, only those distinguishing dif- ferences between herding districts (Fig. 2) were extracted for the analysis, i.e.GrowthS, Temp Aug, InsectAct, TRISSA, WinterS, SnowDepth, TRIW, SnowPrec and IceCrust.Variables that were significant only in the interannual analyses or were related to management or animal density were hence excluded.The results of the two clustering methods (average linkage and Ward's method) differed slightly.Average linkage suggested larger span between the sizes of small and big groups.Ward's method suggested groups of more equal size, although the main pattern was similar.Ward's method was therefore chosen as the main method.The cubic cluster criterion (Sarle, 1983) suggested seven groups and the suggested division is shown in Fig. 5a.

Discussion
The results show that the correlations between the explanatory environmental variables and the productivity variables were stronger between herding districts than between years.Moreover, different factors were relevant for variation in productivity between herding districts and between years.The main productivity gradient (not shown), and herd growth in particular, followed a geographical gradient determined by latitude and topography, which also determines many of the explanatory variables included in this study.The most important factors found to affect variation between herding districts in reindeer productivity were season length, topography, animal density, insect harassment, snow conditions, supplementary feeding and calf slaughter.When analysing the interannual variation, season length, animal density, insect harassment, supplementary feeding and snow depth fell out as the most relevant factors determining productivity.Thus, all these variables are expected to be important determining factors for productivity in reindeer husbandry.Several other factors were also supposed to be relevant, but were not found to be highly determining factors in the analyses.However, there are explanations for this that hinder us from discarding them as irrelevant, such as being strongly correlated with other relevant factors on this scale, being counteracted by reindeer husbandry measures, being poor in data quality, or strongly affecting husbandry economically, i.e. causing higher costs, rather than directly affecting production.The clustering of herding districts, based on the factors identified as being most relevant for variation in productivity between herding districts, was to a large extent similar to the grouping of herding districts made by Lundqvist & Danell (2007).It is natural that they differ from the present administrative division of the reindeer herding districts, which are based on county borders rather than biogeographical variables.Other methods besides CA that could be suitable for linking observations and defining clusters could be e.g.Random Forest (Breiman, 2001), which are available as a CRAN package for the R-project (R Development Core Team, 2008).
The LM and SEM methods suggested different variables for being the least important in determining the productivity variables.This is inherent in the methods.In our study, LM was chosen as guidance prior to the structural equation modelling, as the concept and weaknesses of LM is well known (e.g., Whittingham et al., 2006).The variables were also analysed using canonical correlation analysis with congruent results, but not included in this article.The shortcomings of LM were, however, circumvented using SEM.If we had relied on solely using SEM, some correlations might have been undetected due to the complexity of analysing and extracting the results.The primary advantage of SEM in this context was the possibility to investigate more complex causalities and including latent variables to capture not directly measurable variables.The combined results from these methods were helpful to understand the complexity of the system and the difference between the two statistical methods LM and SEM.It was desirable to perform a multi-level SEM to capture both the characteristics of the herding districts and the interannual variations, as different factors seemed important for the two levels.This approach was, nevertheless, unsuccessful as the number of observations was too small for a realistic multi-level SEM, as well as, different latent variables were suggested for the different levels.Other methods to accompany SEM could be mixed models with random effects, but due to this lack of degrees of freedom it was not chosen in this study.Hence, two different measurement models in SEM preceding the two proposed models were chosen to obtain a satisfactory model for each level, i.e. analyses between herding districts and between years.
The SEM approach was unsuccessful in accounting for the complete relationship between the latent constructs.However, a SEM model could be applied to identify significant causal paths and search for adjustments that improve the fit indices, in order to get an insight into the studied system.No index that evaluates the entire significance of the SEM model has been agreed upon, so it is suggested in SEM literature to use several indices to evaluate a model.In biology, and especially in ecology, relationships are often weak with rather low R 2 values due to the noise from a multitude of uncontrolled factors (e.g., Barry & Elith, 2006).The relationships in the studied system were also expected to be non-linear in some cases, and therefore a non-linear approach could be more successful in achieving a better fit.The degrees of freedom in the data for such an analysis are however too few.The fairly weak fit indices found in this study were therefore quite expected.The results from the SEM analyses did, however, back up the results from the preceding statistical analyses regard-ing the positive and negative path coefficients between the independent and the dependent variables.We therefore considered the model estimates in this study to be informative.
The differences in animal condition and productivity were larger between herding districts than between years.The correlation values in the between herding districts and between year correlation matrices, together with stronger path coefficients in the SEM describing variation between herding districts compared to analyses within herding districts, also supports this conclusion.The results are reassuring for the possibility of using animal condition measures for adaptive management of e.g.stocking rate within herding districts (Olofsson et al., 2008).
In the SEM-analyses, calf slaughter did not fall out as an important management measure to increase herd production, as suggested by e.g.Rönnegård (2003), in contrast to the results of the LM.This result may however have been obscured by a geographically uneven distribution of management systems, where more southern herding districts with longer growing season also have adopted calf slaughter to a larger extent, i.e. creating correlations with e.g.season lengths.The path coefficients did, however, show a weak correlation between calf slaughter (Slaughter strategy) and productivity (Herd production) in the between-district model.
The results from the LM analyses suggested that animal condition was persistent over time, as this was indicated by significance for the condition lag variable in all animal condition variables.They also suggested that there is a consistency in the productivity measures used.Furthermore, the results from the present study showed that pre-slaughter animal feeding positively affects carcass condition, which hence is a possible measure in reindeer husbandry to increase production and economic result, provided that feeding costs are low.Snow condition correlated negatively with herd growth and calf conditions, which was expected due to the energetic costs of foraging, possibly followed by competition for resources.
Years with higher animal densities achieved higher production.This seemed, however, to incur a cost in animal condition.This density dependence seems to be most apparent in the SSA ranges, which confirms a view where the snow-free season is most important for the productivity (e.g., Klein, 1968;Reimers, 1997).High animal densities may have stimulated higher predator densities (bottom-up effect) and hence increase the predatory losses.Depredation rates were included in the herd growth variable (rG) and could thereby explain the increased herd growth.
Herding districts with long growing seasons generally showed better herd growths and female conditions than those with shorter growing season, as suggested by the LM analyses.Years with longer growing season also showed increased animal conditions in both calves and female adults.Years with higher insect activity and August temperatures showed higher herd growth and calf conditions but poorer female condition.This might be an effect of higher calf survival, indicated by the higher herd growth that negatively affects female condition since they invest in their calves at a cost of their own body reserves (e.g., Reimers,

Group
Herding districts  5a for the location of each herding district and Fig. 5b for the geographical distribution of the suggested grouping of herding districts.
According to the results, herding districts with more heterogeneous topography had higher animal productivity.This has been suggested in previous studies where ruggedness has been found positive for reindeer (Nellemann & Thomsen, 1994;Nellemann & Fry, 1995;Nellemann, 1996;Mårell & Edenius, 2006).High ruggedness in SSA areas seemed however to negatively affect animal fatness, which might be explained by higher energetic costs for diurnal movement between forage areas in valleys and insect-free areas on hills, or lower access to grazing areas from the insectfree areas.
In the linear model, intensity of forestry (ClearCut) was suggested to be negatively correlated with reindeer productivity.This variable was non-significant in the SEM analyses and hence discarded there.The findings that forage, road density and snow precipitation would be insignificant for reindeer production is surprising, but this may suggest that neither forage nor roads are determining factors for reindeer productivity on the herding district scale, or that the variation might be explained by the animal density variables.Effects of poor forage conditions might also have been mitigated through e.g.adaptive husbandry adjustments and measures, such as a better utilization of grazing resources through controlled migrations, supplementary feeding (Kumpula et al., 2002), and more optimised herd structures and slaughter strategies (Lenvik, 1990).Furthermore, forestry is correlated to primary production, which follows the latitudinal and topographical main gradients.The forestry variable could instead have been integrated in the animal density variables since forestry indisputably makes less area available for grazing, and thereby indirectly increasing animal density, especially during snow season when lichen is an important forage for reindeer (Kumpula, 2001).Forestry is clearly negative for lichen abundance (Dettki & Esseen, 1998).Such forestry and forage range variable integration was however not done in this study.
Ice-crust formations and other severe weather events can be locally devastating for the reindeer industry, but effects of ice crust were not found very significant in this investigation.The reindeer herders may under severe ice-crust conditions choose to use alternative ranges or initiate supplementary feeding.Some herders might also increase their culling to decrease the grazing pressure on winter ranges.Supplementary feeding can be very costly and have a large impact on the herding enterprises as such, but this is not seen in the production figures.Therefore, these weather events can be difficult to identify as important using correlation analyses with productivity.
The clustering of the herding districts agreed to a large extent with the grouping of herding districts suggested by Lundqvist & Danell (2007).The similarities were larger than the dissimilarities and the main groups in Lundqvist & Danell were found in this grouping as well, with a few exceptions.Some of the outliers found by Lundqvist & Danell were also distinguished here.The outliers found in this study should preferably be grouped together with adjacent herding districts for administrative purposes, but should still be considered as unique in their surroundings in matters concerning reindeer husbandry conditions and productivity.Based on the results of Lundqvist & Danell, together with this analysis, we ultimately suggest an administrative division into seven groups as shown in Table 5.
Multivariate analyses of production systems such as reindeer husbandry are essential to get an insight into such a complex system of causalities and multi-causal effects.This analysis has sorted out some of these multidimensional effects and sheds new light on factors that were considered important for reindeer productivity on a herding district scale.All in all, rein-deer husbandry is an industry that handles a variety of determining factors on a daily basis, which are handled to overcome incidents that negatively affect the reindeer and the grazing resources.Large-scale weather cycle patterns, such as the North Atlantic Oscillation (NAO), El Niño Southern Oscillation (ENSO) and the Arctic Oscillation (AO) affect the climate in Fennoscandia and therefore are expected to affect the ecosystems (e.g., Stenseth et al., 2002) and subsequently we should expect them to affect reindeer ecology.They were not included in this analysis due to the short time-scale of the study, but to achieve a sustainable longterm productivity in a variable environment, it is important that the reindeer husbandry is prepared for changes.This should include close watch of forage and animal conditions, and to have a buffer to mitigate negative impacts of disturbances and severe weather conditions.This is probably more important than ever when facing a climatic change (ACIA, 2004;IPCC, 2007) and the impact it will have on reindeer husbandry as such, as well as the effect of increasing land competition from other industries and society as a whole.Hence, in such scenarios, an adaptive management approach in the reindeer husbandry will be increasingly important.
In conclusion, the results of this study suggest that the extraction and reduction of variables in Lundqvist et al. (2007) and the reindeer herding district division of Lundqvist & Danell (2007) were purposeful plausible procedures.The results regarding the identification of determining environmental variables are in good agreement with previous research, although there are variables not included in this study, which could infer some bias in the analyses.Hence, this analysis should help to find new hypotheses and adjust existing ones.The results fully support the view that growing season lengths and animal densities primarily determine the reindeer productivity.These variables are however also affected by several other factors, which some could be identified to correlate with productivity in our study.The remaining variables could however not be concluded to be insignificant for reindeer productivity, as they may be hidden in the complex pattern of correlations, or being mitigated by husbandry measures.

Fig. 1 .
Fig. 1.Proposed between-district structural equation model used to explore the relationship structure between the environmental, constructed latent and productivity variables.Measured manifest variables are shown in white, exogenous explanatory latent variables in green, endogenous explanatory latent variables in yellow and endogenous productivity latent variables in orange.

Fig. 4 .
Fig. 4. Estimated parameters of the final structural equation model describing interannual variation within herding districts.Significant (P ≤ 0.10) path coefficients and latent variable covariances between manifest variables (white), exogenous latent variables (green), and endogenous latent variables (orange) are shown as solid arrows.Fixed correlations are marked with(x) and correlation between exogenous latent variables are drawn in grey.

Fig. 5
Fig. 5. a) Clustering of herding districts based the nine most relevant productivity factors identified in the LM and SEM analyses describing variation structure between herding districts.Some herding districts are overlapping and are therefore partly hidden in the map.See Table 5 for district abbreviations.b) The suggested administrative grouping of the herding districts based on Fig. 5a and administrational conveniences.

Table 1 .
Description of variables.

Table 2 .
Estimated regression coefficients of productivity variables on stable and temporally varying explanatory variables, obtained in stepwise linear model analyses with Year included as class variable.Insignificant coefficients are not shown, and CalfR was included in the rG analyses only.

Table 3 .
Estimated regression coefficients of productivity variables on temporally varying explanatory variables, obtained in stepwise linear model analyses with Herding district included as class variable.Insignificant coefficients are not shown, and CalfR was included in the rG analyses only.