Semi-domesticated reindeer avoid winter habitats with exotic tree species Pinus

The introduction of exotic tree species can have profound effects on the native environment, including habitat use and movement patterns of animals, as well as becoming a management challenge for other land users. Here, we used GPS data from reindeer ( Rangifer tarandus ) and remote sensing measurements of lichen cover and soil moisture to assess the effects of the exotic lodgepole pine ( Pinus contorta ) on reindeer husbandry by the Indigenous S ´ ami in northern Sweden. We used locational data from 67 reindeer for three winters to analyze their habitat selection at the second-order selection (placement of home range in the landscape) and third-order selection (selection of sites within the home range) in relation to land cover class, terricolous lichen cover as measure of winter forage abundance, topographic features


Introduction
Land use and landscape fragmentation affect how animals are willing or able to move across a landscape (Zeller et al., 2012;Kauffman et al., 2021).Animal response to disturbance and anthropogenic change in the landscape may lead to increased energy expenditure on movement, decreased foraging time, and risks involved in finding alternative movement paths (Shepard et al., 2013;Monteith et al., 2018;Shaw, 2020).Consequently, movements, including animal migrations, are becoming increasingly disrupted around the globe as anthropogenic encroachment increases (Tucker et al., 2018;Doherty et al., 2021;Barker et al., 2022).Equipping animals with GPS devices allows us to monitor and quantify their habitat use in anthropogenic landscapes and develop conservation strategies to preserve movement and migration patterns, and thus population viability (Kays et al., 2015;Kauffman et al., 2021).
Purposeful or accidental introduction of exotic species can modify both native habitat structure and ecological functionality, affecting ecological interactions such as habitat use and movement patterns of native animal species (Strauss et al., 2006;van Wilgen & Richardson, 2014;Livingstone et al., 2020).Management of exotic species thus becomes a challenge for conservation, as formulating and implementing environmental policies require detailed knowledge about their ecological impacts (Larson et al., 2011;Simberloff et al., 2013;Pyšek et al., 2020).Furthermore, the relationship between people and introduced exotic species may be complex, as costs and benefits derived from these species depend on values placed on them by different interest groups (Crowley et al., 2017;Beever et al., 2019).One example is when exotic species have been introduced for positive economic reasons, such as food or timber production, but have unintended negative consequences on other ecosystem services valued differently by other members of society (Estévez et al., 2015).
Trees have been widely introduced beyond their native ranges to increase production of timber and fiber, hereafter 'exotic trees' (Castro-Díez et al., 2019).Although these introductions can have positive economic outcomes in the short-term, exotic trees may have negative longterm impacts on native biodiversity and ecosystem function, which can result in destabilizing feedbacks on the economy (van Wilgen & Richardson, 2014;Pyšek et al., 2020).Management of exotic trees can therefore create conflicting goals between different stakeholder groups, if the services and disservices produced by these forest ecosystems are experienced differently (Dickie et al., 2014;Shackleton et al., 2019).
Lodgepole pine (Pinus contorta) is a coniferous tree species native to western North America.Between the 1960s and 1980s, Swedish forestry companies started introducing P. contorta in Northern Sweden to avert an anticipated scarcity in timber supply (Elfving et al., 2001).Under optimal conditions, the growth rate of P. contorta surpasses that of the native Scots pine (P.sylvestris) by 36 % and can reduce rotation cycles between harvests by 10-15 years (Elfving et al., 2001;Engelmark et al., 2001), thereby increasing timber production.At present, P. contorta is the dominating tree species on approx.520,000 ha, mostly in the northern half of the country, corresponding to ca. 4.6 % of the productive forest of that area (Skogsdata, 2022).At the time being, most of the stands dominated by P. contorta are between 40 and 60 years old (Skogsdata, 2022).
So far, the ecological effects of non-native P. contorta have mainly been investigated from a biodiversity perspective.In Sweden, stands with P. contorta affect e.g., species richness and abundance of several taxa, including cryptogams, vascular plants, invertebrates, and birds in Sweden (Engelmark et al., 2001;Nilsson et al., 2008;Bäcklund et al., 2015Bäcklund et al., ,2016;;Löfroth et al., 2022).Higher stem densities, longer, lower, and thicker branches, denser canopies, and higher levels of needle litter limit light availability at the forest floor in stands with P. contorta compared to P. sylvestris (Bäcklund et al., 2018).This affects species diversity and community composition particularly for light-demanding species, including terricolous lichens (Cladonia spp.) (Nilsson et al., 2008;Bäcklund et al., 2015Bäcklund et al., , 2018)).Terricolous lichens constitute the main forage resource for semi-domesticated reindeer (Rangifer t. tarandus) during winter (Åhman & White, 2018).In addition to having negative impacts on winter forage availability, reindeer herders report stands with P. contorta to impede reindeer movements or herding activities, such as gathering the animals for roundups (SSR, 2019;Roos et al., 2022).Thus, the question how P. contorta affects reindeer movement, habitat selection and terricolous lichens as forage resource is of particular concern for reindeer husbandry not only from an economic perspective (e.g., reduced access to winter forage may result in lower body weights in reindeer), but importantly as being the keystone in culture and society of the indigenous Sámi (SSR, 2019).As Indigenous people, the rights of the Sámi to practice their traditional livelihood are affirmed among others by the UN Declaration on the Rights of Indigenous Peoples (2007).Multipleand partly opposingland-use interests overlap spatially in the forest landscape of Northern Sweden (Sandström et al., 2016).Within this context there is an ongoing debate how Indigenous property rights can be understood beyond cultural arguments (Brännström, 2017), underscoring the need for (scientific) evidence of impacts by one land use on the other within land use planning.
The reindeer husbandry area in Sweden covers the northern half of the country, sharing a vast area with commercial forestry with a long history of disputes between the two forms of land use (Kivinen et al., 2010, Sandström et al., 2016).Compared to southern Sweden, the occurrence of P. contorta is higher in the three northern counties, largely overlapping with the reindeer husbandry area (Fig. 1a).Norrbotten, as the northernmost and largest county has a share of 2.5 % of P. contorta stands on productive forest area, and a share of 3.2 % in Västerbotten and 6.6 % in Jämtland (Skogsdata, 2022).
The aim of this study was to assess the impact of stands with P. contorta on reindeer behavior and lichen occurrence to provide information for sustainable forest management.Despite the extensive documentation of the conflict between reindeer husbandry and forestry about P. contorta (Roos et al., 2022), the habitat selection of reindeer in relation to P. contorta has so far not been assessed and quantified.We used reindeer locational data collected during three winters between 2017/18 and 2019/2020 to analyse reindeer habitat selection in relation to exotic P. contorta versus native P. sylvestris.We predicted that (i) habitat selection by reindeer differs between stands with P. contorta and stands of the native P. sylvestris and (ii) that reindeer select areas with a higher lichen abundance.Further, we used remote sensing data to quantify (iii) terricolous lichen cover between different forest types and (iv) explored on what soil moisture classes these different forest types occur to assess if stands with P. contorta occupy habitats potentially suitable for terricolous lichens.

Study area and period
The reindeer husbandry area in Sweden is divided into 51 reindeer herding communities (RHC).A RHC is a defined geographical area, as well as an economic and administrative unit that organizes reindeer husbandry within its borders.Within each RHC borders, semidomesticated reindeer migrate between seasonal grazing areas (Skarin et al., 2022).Winter grazing areas are located in the boreal forest, which simultaneously is used intensively by commercial forestry.Scots pine (P.sylvestris) dominates on xeric and oligotrophic sites, as well as on mires, while Norway spruce (Picea abies) occupies mesic-moist soils (Esseen et al., 1997).P. contorta is planted on varied sites, but mostly in areas suited for P. sylvestris.Deciduous trees include birches (Betula pendula, B. pubescens), aspen (Populus tremula), rowan (Sorbus aucuparia) and willow (Salix spp.).The forested landscape is relatively flat, interspersed with mires and lakes, ranging from the coastline of the Baltic Sea to hills with ca.700 m in elevation further inland.During the winter months (December-March), average temperature ranges from − 9 • C further inland to − 5 • C at the coast.A permanent snow cover usually forms in late October to early November and disappears in late April to early May, depending on weather conditions and distance from the sea (Swedish Meteorological Institute, https://www.smhi.se).Snow depth reaches its maximum in March, with an average of 80 cm for the period 1961-1990 (Swedish Meteorological Institute).
Our study area covered the winter grazing area of Vilhelmina norra RHC in Västerbotten county (Fig. 1a, b).The herding community is divided into two or more winter groups (siidas).The Marsfjäll-siida uses the southern and western part of the winter grazing area, while the Vardofjäll-siida uses the northern parts.Except for strategic moments of animal handling, such as round-ups in corrals to divide herds into winter groups, reindeer move and graze freely in the landscape.Herders guard the herds to minimize dispersal of the reindeer into grazing areas belonging to other winter groups or RHCs, and to minimize disturbances and protect them from predators, mainly lynx (Lynx lynx), wolverine (Gulo gulo), and occasionally wolves (Canis lupus).Our study period covered the winters 2017/18, 2018/19 and 2019/2020.During the winters 2017/18 and 2018/19, reindeer remained west of the mostly unfenced railroad cutting through the eastern part of the winter grazing area (Fig. 1b).Herders seek to keep reindeer away from the railroad to prevent train accidents.Due to difficult snow conditions, reindeer however also used the area east of the railroad in 2019/2020 (Fig. S1).We thus only included the western part of the winter grazing area (8870 km 2 ) in the analysis of the two first winters, while we included the whole winter grazing area in the analysis of the third winter (10,830 km 2 ).We added a buffer of 2.5 km around the winter grazing area to encompass GPS-positions outside borders of the RHC.

GPS data
Vilhelmina norra RHC provided us with positional data from female reindeer (n = 67) fitted with GPS-collars (Followit Lindesberg AB, Telespor AS).Positions were recorded at 05:00 and 17:00 every day during the three winter periods.Prior to analysis, we screened GPS-positions for erroneous relocations and spike positions, following Bjørneraas et al. (2010).We also removed positions when collars were clearly not on an animal (e.g., in houses), when the animals were transported on the road with vehicles, kept in corrals, or when herders moved the animals between herding infrastructure and the movements did not represent an animal's habitat selection.Further, we excluded time series of collars with >15 % missing data.In total, we included data from 67 reindeer, resulting in 12,158 positions after data cleaning (Table 1).Weather conditions (e.g., snow arrival) during late autumn and early winter determine when herders gather their animals and move to the winter grazing area (Horstkotte et al., 2014, Rosqvist et al., 2022).Likewise, the start of migration back to calving grounds during spring varied between years and herding groups.To analyse a common period throughout all winters we limited our study period from 1st of November to 30th of April.

Habitat variables
We reclassified the Swedish Land Cover Map (NMD, Environmental Protection Agency, 2018a, 10 m × 10 m resolution) that originally consists of 25 land cover classes into land cover classes that reflect the major vegetation types of the study area (Fig. 1c, Table 2, Table S1).For each year, we updated the map with information about forest harvests and aging of forest stands.Areas harvested were reclassified as clear cuts (Swedish Forestry Agency, https://www.skogsstyrelsen.se/skogligagrunddata), and clear cuts were reclassified as young forests when becoming older than five years.Over the study period, differences in landscape composition remained low.For instance, the proportion of clear cuts varied from 2.8 % of the study area in the winter 2017/2018 to 2.5 % in the winter 2018/19.The NMD did not discriminate between P. sylvestris and P. contorta but referred to both as pine forest.Therefore, we retrieved geospatial information on stands with P. contorta from databases provided by the forest companies Holmen, SCA and Sveaskog.We extracted areas with information about age, height, and proportion of P. contorta and reclassified the overlapping pixels in the original NMD.In our study area, P. contorta occurred only in a few cases on private owned land (Skogsdata, 2010).We classified

Table 1
Number of reindeer and number of GPS locations (min, median, max per individual) used to assess habitat selection in the winter grazing area of Vilhelmina norra reindeer herding community (Sweden) between winters 2017/18 and 2019/20.The population range size for each winter was estimated using a Brownian Bridge movement model.stands with P. contorta into four different classes.Stands where the proportion of P. contorta in the tree species composition was <20 % were combined with the land cover class "mixed coniferous forests" in the original NMD.Stands where P. contorta comprised between 20 % and 60 % of the tree species composition were separated from stands with >60% of P. contorta, i.e., where P. contorta was the dominating tree species.Furthermore, we divided P. contorta-stands according to the tree height into two height classes.Tree height was provided in the intervals 0 m, 0.5-3 m, 3-5 m, and further intervals of 5 m (Environmental Protection Agency, 2018b).We chose 3 m as a cut-off value to divide height classes of P. contorta.Below that threshold, we considered P. contorta too low to present a significant obstacle for reindeer or to develop a structure becoming disadvantageous for terricolous lichens (hereafter: stands ≤ 3 m).Above that threshold, stands with P. contorta are likely to negatively affect reindeer movements and terricolous lichens due to the development of a dense structure (hereafter: stands >3 m).In total, we used twelve land cover classes in the analysis of reindeer habitat selection (Table 2).Stands where the proportion of P. contorta is >20 % cover 3 % of the productive forest lands of Vilhelmina Norra RHCs winter grazing area (Table 2).About 77 % of the area covered by P. contorta, and 81 % of all individual P. contorta stands that make up this area, belong to the class with P. contorta as the dominant species (i.e., >60 %).Stands with >20 % of P. contorta cover 3.5 % of the productive forest lands west of the railway (as defined by reindeer use during the winters 2017/18 and 2018/19) (Fig. S1).
Of the four classes defined for P. contorta, the class with P. contorta as the dominant tree species and >3 m in height is the most abundant, exceeding in extent all three other P. contorta-classes combined (Table 2, Fig. S2).P. contorta stands have an average size of 15.8 ha (minimum: 0.07 ha, maximum: 144.2 ha), while stands of the other forest classes have an average size of 8.1 ha with substantial variation (minimum: 0.05 ha, maximum: 528.4 ha).
We quantified topographic variation in the study area using elevation, terrain ruggedness index (TRI) and topographic position index (TPI).Following De Reu et al. ( 2013), we classified TPI into valleys, flat areas, lower and upper slopes, and ridges.We calculated all topographic variables from a digital elevation model (DEM) at a resolution of 10 m × 10 m (Swedish Forestry Agency, https://www.skogsstyrelsen.se/skogligagrunddata).Furthermore, we calculated the minimum Euclidean distance to major public roads (>5 m wide, www.lantmateriet.se)for each GPS position, as reindeer avoid roads with regular traffic (Skarin & Åhman, 2014).We tested for collinearity between any of the variables with the R-package performance (Lüdecke et al., 2021), and variance inflation factors did not exceed a critical threshold ≥3 (Zuur et al., 2009).

Lichen data
To test the significance of terricolous lichen cover for habitat selection by reindeer, as well as its relationship to forest classes, we used a raster map of predicted lichen cover.We predicted lichen cover with data from the Swedish National Forest Inventory, collected during 2014-2019 (Fridman et al., 2014;Adler et al., in preparation).We combined data on terricolous lichen cover with the Normalized Difference Vegetation Index (NDVI), a measure of vegetation productivity, derived from Sentinel 2 satellite data.We obtained additional explanatory variables from remote sensing data (LiDAR), including tree height, soil type, altitude, slope, and aspect from the Swedish Land Survey Agency (http://www.lantmateriet.se/en/).We modelled lichen cover with Generalized Additive Models.We checked for multicollinearity between the explanatory variables using the function vifstep in the Rpackage usmd (Naimi et al., 2014).The explained deviance of the resulting model was 57.4 %, with tree height, NDVI and Sentinel 2 bands 2 and 11 as the most important predictors.In general, the model underestimated lichen cover, as the training set included only few sample plots with high lichen cover (for detail see https://sametinget.se/170594).Based on the model, we produced a raster map of lichen cover for the whole study area with a resolution of 10 m × 10 m.

Statistical analysis of habitat selection
We analysed habitat selection by reindeer using habitat selection functions (HSF).The HSF estimates the probability of a particular habitat being used disproportionally relative to its availability compared to a reference category (Johnson, 1980;Manly et al., 2002;Fieberg et al., 2021).As habitat selection is a hierarchical process, different levels of space and time need to be considered to understand the observed behaviour (Senft et al., 1987;Northrup et al., 2022).To capture reindeer selection at both regional and landscape scales, we analysed selection of home range placement within the winter grazing area (second order of selection sensu Johnson, 1980), and selection of sites within that home range (third order of selection sensu Johnson, 1980).Decisions by reindeer herders may influence reindeer habitat selection and behaviour, such as when to start migration between seasonal ranges or selection of grazing areas within a seasonal range (Skarin et al., 2022).Habitat selection within seasonal grazing areas and movements between them, however, usually reflect reindeer decisions.
To quantify these relationships, we applied habitat selection functions with a use-availability design, i.e., generalized linear mixed-effects

Table 2
Land cover classes used as explanatory variables used to assess reindeer habitat selection and their contribution to the landscape composition in the winter grazing area of Vilhelmina norra reindeer herding community (Sweden) between winters 2017/18 and 2019/20.models (GLMM) with a binomial distribution (Fieberg et al., 2021).
Positive model coefficients for a given class variable indicate a higher likelihood of use (i.e., selection) compared to the reference class, while negative coefficients indicate a lower likelihood of use (i.e., avoidance) of that habitat class compared to the reference class (Fieberg et al., 2021).For continuous variables, positive coefficients indicate a higher use as the variable increases, while negative coefficients indicate that use decreases as the variable decreases (Fieberg et al., 2021).We calculated the relative selection strength, i.e., the likelihood to be chosen relative to the reference level for each variable by exponentiating the model coefficient of the GLMM (Fieberg et al., 2021).We compared habitat components at reindeer GPS locations, i.e., used locations, with available locations.We analysed habitat selection at the second order by comparing used locations to available locations randomly spread within Vilhelmina norra's winter grazing area each winter.We analysed habitat selection at the third order by comparing used locations to available locations randomly spread within the home range of each individual.At both scales of selection, we generated ten times as many available locations relative to used locations (Signer et al., 2019).We included year and animal ID as two different random effects to account for variance between animals using the R-package lme4 (Bates et al., 2015).We defined individual home ranges as the estimated 99% utilisation distribution using Brownian Bridge Movement Model with the kernelbbfunction in the adehabitatHR package in R (Calenge, 2006, Horne et al., 2007, Kranstauber et al., 2012), using a position error of 20 m.For the land cover classes, we used the forest class of mature P. sylvestris (dominated by trees > 5 m in the NMD) as the reference category, as oligotrophic P. sylvestris forests best reflect the ecological niche of terricolous lichens (Esseen et al., 1997).For the second categorical variable, TPI, we used flat areas as reference category to compare reindeer habitat selection against the remaining TPI classes (Table 2).
We used a log10-transformation for the distance to roads as animals' responses to roads probably decreases at larger distances (Nielsen et al., 2009).We centred all other continuous variables to be in similar range, using the scale-function in R. We used a significance level of α < 0.05 for the variables in the HSF.
To identify the most parsimonious model, we started with a model including all variables (Table 2) and progressively removed variables until we retained only land cover class and lichen cover.We selected the best-fit model based on the lowest value of the Akaike Information Criterion (AIC, Burnham and Anderson, 2004) (Table S2a, b).To evaluate the prediction success of our most parsimonious model, we used kfold cross validation (Boyce et al., 2002).We randomly divided the data set into a training set consisting of 80% of the data and a test set consisting of the remaining 20%.To assess the consistency between the predicted probability of selection derived from the training model and the test model, we binned the corresponding model outcomes into ten equal-sized classes.Spearman's Rank correlations at a significance level of α < 0.05 between the frequencies for the two sets indicate that the model is applicable to both the training and the test set.To evaluate model fit, we repeated the procedure of random data-partitioning and correlation of bin frequencies ten times and calculated the average correlation coefficient and significance (α < 0.05) of the ten resulting correlations.To assess the variance explained by the model, we used the r.squaredGLMM-function in the package MuMIn (Barton, 2020).

Differences between forest classes relative to lichen cover and soil moisture
To test for differences in lichen cover between forest classes, we used the habitat classes dominated by different types of forests (Table 2) and the data on lichen cover described in section 2.2.3.To investigate the relationship between forest classes, including clear cuts, and soil moisture, we extracted data on soil moisture, classified as dry, mesic, and moist soils from the Environmental Protection Agency (2018c).We resampled the original raster of 2 m × 2 m to a resolution of 10 m to fit the raster of habitat classes, using the r.resample function in GRASS (GRASS Development Team, 2017).

Statistical analysis of forest classes relative to lichen cover and soil moisture
We tested for differences in lichen cover between forest classes by calculating the area-weighted mean lichen cover over all stands representing the forest classes (Table 2), using the 8-neighbour rule to define a stand.We excluded clear cuts from this analysis, as they were masked in the original lichen map.Due to the non-normal distribution of the data, we applied a Kruskal-Wallis one way ANOVA on ranks.To test for significant differences of pairwise comparisons between forest classes, we used a Dunn-test for multiple comparison with a Bonferronicorrection. To test for an association between the three levels of soil moisture and forest classes, we extracted the most represented soil moisture class for each stand and applied a χ 2 -test and a correspondence analysis to represent that relationship (Quinn & Keough, 2002).Standardized residuals of the χ 2 -test exceeding ± 1.96 indicate a significant positive or negative association between forest class and soil type (Agresti, 2003).For this analysis, we combined stands ≤3 m and stands >3 m of both dominance classes of P. contorta, as they represent similar forest types that only differ in their successional age.As we tested for the association between of forest classes and soil moisture, we removed all lakes and mires from the soil moisture dataset.We used R version 4.1 (R core team, 2021) for all data handling and statistical analyses and QGIS version 3.16.3for spatial analyses.

Habitat selection by reindeer
At both scales of selection, the most parsimonious model included all variables (Table 2), as the AIC value increased when variables were removed (Table S2a,b).Reindeer avoided stands with P. contorta >3 m in height relative to stands with mature P. sylvestris (Fig. 2a).Within these stands, the relative selection strength indicated that reindeer were 57 % less likely to choose stands dominated by P. contorta (>60 % coverage) than the reference class at the second order of selection, and 61 % less likely at the third order (Tables S3 and S4).Stands with 20− 60 % P. contorta were 36 % and 42 % less likely to be chosen at the second and third order, respectively (Tables S3 and S4).Stands with P. contorta where trees were ≤3 m did not affect reindeer habitat selection relative to mature P. sylvestris (Fig. 2a).
Notably, the number of used locations was low for the four classes of P. contorta, due to their comparatively low abundance in the study area (Table 3).
At both selection scales, reindeer selected for young forests of native conifers (Fig. 2a, Table S3 and Table S4).Reindeer preferred upper slopes and ridges relative to flat areas at both orders of selection (Fig. 2b) and more rugged terrain but preferred lower elevation (Fig. 2c).Reindeer preferred areas with higher cover of terricolous lichens and areas at larger distances from roads at both scales of selection (Fig. 2c).
The average Spearman correlation coefficients of the k-fold cross validation for the HSF was 0.95 (p = 0.001) for both the second order and third order of selection, indicating a strong fit between the training and the test data set.The variation in habitat selection explained by the models remained low at both orders of selection (conditional R 2 = 0.14 at the second order, 0.13 at the third order).

Differences between forest classes relative to lichen cover and soil moisture
Area-weighted mean lichen cover differed between forest classes (Kruskal-Wallis χ 2 = 79675, df = 7, p < 0.001).Area-weighted mean lichen cover was higher in stands of mature P. sylvestris compared to all T. Horstkotte et al. other forest classes (Dunn-test, Table S5).All other mature forests classes had a lower mean of area-weighted lichen cover than any class of young forests, including the two P. contorta classes lower than 3 m (Fig. 3).However, area-weighted mean lichen cover did not differ between stands where P. contorta was present with >60 % and where trees were ≤3 m compared to stands of lower abundance of P. contorta with trees >3 m (Fig. 3).Stands with >60 % P. contorta and trees >3 m in height had lower area-weighted mean lichen cover than any other forest class, except for mixed forests where neither deciduous nor coniferous trees dominate (Fig. 3).Averaging the estimates on total of lichen cover for each forest class across the whole study area showed that about 49 % of all lichens occurred in mature stands with P. sylvestris, which constitute about 30 % of the total winter grazing area (Table 2), excluding non-forested areas.
Forest classes differed in soil moisture (χ 2 = 128,804, df = 12, p < 0.001; Table S6).Both dominance classes of P. contorta, as well as young forests were closer associated to dry soils (Fig. 4).In our study area, stands of mature P. sylvestris were associated stronger with moist soils compared to dry soils, and so did mixed coniferous forest (Fig. 4).Clear cuts were intermediate between dry and mesic soils (Fig. 4).

Habitat selection by reindeer
Our study on reindeer behavioural response to stands with P. contorta within a managed boreal forest landscape highlighted that reindeer differentiate between forest stands with P. contorta and P. sylvestris in relation to forests' successional stage.Specifically, we found that reindeer used stands with of P. contorta with trees >3 m less than stands dominated by native P. sylvestris at both scales of selection, i.e., placement of the home range in the landscape and in selection of habitat patch within the home range.Our results suggest that stands with P. contorta reduce the grazing grounds available to reindeer in winter, at least at progressed successional stages.This implies that the quality and usability of the winter grazing area for reindeer are reduced at the landscape level.
Our results confirm reindeer herders' reports about reindeer avoidance of stands with P. contorta for stands with tree heights >3 m.The low occurrence of terricolous lichens in stands with P. contorta in this study may further explain why reindeer avoided stands with P. contorta trees >3 m.Contrastingly, we did not find that reindeer avoid stands with P. contorta trees ≤3 m, relative to stands dominated by mature

Table 3
Number and % of GPS positions and available locations in each land cover class at the second and third order of selection.The ratio of the percentage of GPS positions to the percentage of available locations is an indicator of greater use than expected (>1) or less than expected (<1).P. sylvestris.Most likely, these younger stands did not hinder reindeer movements or have not yet had a negative effect on lichen cover.However, we expect that stands with P. contorta that today are <3 m in height will become avoided once they have grown higher than 3 m, as we have found for stands that are above that height threshold today.These stands will affect the suitability of the future forested landscape as grazing area for reindeer.Hence, the area of avoided habitat will continue to increase, if not stands of mature P. contorta will be harvested before this and replanted with native species.We found that reindeer selected areas with higher lichen cover, i.e., a higher availability of forage resources, a well-known driver of habitat selection by Rangifer in winter (e.g., Helle, 1984, Johnson et al., 2002, Hins et al., 2009).Furthermore, we found that reindeer selected young forest stands dominated by a mixture of native conifers relative to stands dominated by mature P. sylvestris.While increased visibility for predator detection may contribute to selection of young forests by reindeer (Anderson & Johnson, 2014), our results pinpoint the possibility of young forests being important habitats additionally to old growth forests during winter, as the relatively open canopy and low basal area may contribute to high lichen growth rate after final harvest (Horstkotte & Moen, 2019).However, old-growth forests can fulfil different functions as reindeer grazing grounds, as they can provide suitable foraging conditions under a range of snow and winter weather conditions (Roturier & Roué, 2009).

Differences between forest classes relative to lichen cover and soil moisture
We found that stands with of P. contorta primarily occur on dry soils, the prime habitat where terricolous lichens have competitive advantages over vascular plants (Cornelissen et al., 2001).This indicates that potential habitat for terricolous lichensin other words, valuable winter grazing areas for reindeeris lost when stands with P. sylvestris are turned into stands with P. contorta.According to our results these stands do not support high lichen abundance.We therefore anticipate that this type of habitat that most likely would be preferred by reindeer may ultimately turn into avoided habitat, which would result in a net reduction of grazing grounds.However, estimation of lichen cover before planting of P. contorta was beyond the scope of our study, and we recommend further research to assess the effect of stands with P. contorta on previously lichen-dominated sites in more detail.In our study area, mean lichen cover was lowest in stands with a mixture of coniferous and deciduous trees.We therefore assume that overall low occurrence of lichens is the most likely reason why reindeer avoided stands with a mixture of deciduous and coniferous trees.
The association we found between P. sylvestris and moist soils is most likely an effect of the species' wide ecological amplitude that allows it to grow also on forested wetlands (Esseen et al., 1997).Old-growth forests of P. sylvestris and Picea abies on productive sites can be important habitat for epiphytic lichens, and thus are vital when deep and/or hard snow prevents reindeer from digging for terricolous lichens (Berg et al., 2011;Esseen et al., 2016).Under such events, rugged terrain with smallscale topographic variability can offer sites of suitable grazing conditions (Inga, 2007).This explains why we found reindeer to prefer rugged terrain and selected for ridges and upper slopes.
We found that nearly 50% of all lichens occurred on stands  dominated by mature P. sylvestris, comprising 30% of the study area.
Since the introduction of clear cutting in Northern Sweden in the 1950s, stands of P. sylvestris where lichens cover exceeds 50% of the ground vegetation have decreased by 71%, leading to a high grazing pressure on the remaining areas (Sandström et al., 2016).Even though the proportion of stands with P. contorta in the study area and the entire reindeer husbandry area may seem low, the fact that they can occur on dry soils where terricolous lichens thrive may reduce the available winter grazing grounds for reindeer.Stands with P. contorta therefore affect reindeer husbandry disproportionally to their proportion in the landscape.

Limitations
Our study focused on the selection of different habitats by reindeer at the second and third order, rather than the selection of forage resources per se (e.g., terricolous lichens) compared to other forage resources, i.e., the fourth order sensu Johnson (1980).We acknowledge that the concept of habitat comprises the sum of all environmental variables at a given location that affect animal behaviour across scales.As any habitat selection model, our model is therefore a simplification of the reality that drives behavioural patterns of reindeer.Within our study area, the proportion of P. contorta was low at the landscape level, resulting few reindeer GPS positions and corresponding available locations within stands with P. contorta.This may introduce uncertainties into the estimation of habitat selection of comparatively uncommon P. contortaclasses relative to the other land cover classes.Uncommon cover classes are challenging to estimate in the process of the animals' habitat selection, as important drivers of habitat selection might be missing (Northrup et al., 2022).Our analysis on reindeer habitat selection was based on GPS-positions with a fix rate of 12 h.Higher temporal resolution on reindeer movement data would have provided more detailed insights into reindeer movement behavior and habitat selection (Northrup et al., 2016), such as residence time or activity levels within stands of different trees species composition (Rautiainen et al., 2022).The temporal resolution reflects the prevailing herding practices in Vilhelmina Norra RHC, where disturbances levels and battery life length were balanced with information needed in the daily work with the reindeer.

Management implications
Within the managed forest landscape of Northern Sweden, few options exist today for herders to rotate between grazing areas to allow lichen recovery in temporarily ungrazed areas (Axelsson-Linkowski et al., 2020).Planting P. contorta as part of forest management practice thus add to other external pressures on reindeer husbandry, including competition for space with other forms of land use such as infrastructure, wind power development and mining (Skarin et al., 2018;Rosqvist et al., 2022), but also with negative effects of predation (Åhman et al., 2022) and climate change (Löf, 2013;Rasmus et al., 2022).In our study, reindeer avoided areas close to major roads, demonstrating and confirming the impact that such disturbances have on reindeer habitat use (Anttonen et al., 2011).Thus, the percentage of the study area covered by P. contortaan impediment to reindeerneeds to be understood in the context of a landscape already heavily affected and as contributing to the direct and behavioral losses of grazing grounds.In our study area, stands dominated by P. contorta > 3 m in height were the most common of the four P. contorta classes.As these stands will be ready for harvest within the coming decades, there is a window of opportunity to restore these sites to functional grazing areas that support higher lichen abundances.Contrasting the assumption that P. contorta would depend on forest fires to regenerate, it is spreading naturally in Sweden, where natural forest fires are prevented by forest management (Jacobson & Hannerz, 2020).Notably, even though planting P. contorta is prohibited within 1 km from national parks or nature reserves, concerns have been raised that the species will spread beyond its current distribution limits, and individuals have already been reported in the mountain birch forest (Engelmark et al., 2001;Kullman, 2016).These high levels of uncertainty involved in understanding the effects of exotic species requires a context-specific planning approach at the regional and local level.While P. contorta is still actively planted by several of the major forest companies in the reindeer husbandry area, they follow regulations by the Forest Stewardship Council (FSC), a certification scheme for responsible forest management (FSC, 2020).The FSC regulations involve consultations with the herding communities including time perspectives of 5-7 years planning of the forest activities at the landscape scale.This includes not planting P. contorta on or removing it from areas important for reindeer herding.Planning at the landscape level also means to consider the connectivity of reindeer grazing areas to respond to different grazing conditions within and between different winters (Horstkotte et al., 2014).In particular, the interaction between winter weather and forest structure influences snow conditions and accessibility of forage resources below the snowpack, and thus influences which forest types are most suitable for grazing in any given winter (Horstkotte & Roturier, 2013;Kater & Baxter, 2022).According to reindeer herders in Vilhelmina Norra RHC, winter conditions differed between the three years of this study.Deep snow and frequent thawfreeze cycles during the winter 2019/2020 resulted in ice formation, while grazing conditions were better in winters 2017/18 and 2018/ 2019.To assess the impact of such conditions on habitat selection will require an analysis over more years and studies of fine-scale behaviour of reindeer (Rautiainen et al., 2022).

Conclusions
Our findings suggest that stands with P. contorta reduce the area available for reindeer grazing during winter, due to the animals' avoidance of these sites and a lower occurrence of foraging resources compared to forest with the native P. sylvestris.This implies that removing stands with P. contorta and restoring them with native P. sylvestris would be a way forward to improve the functionality of the landscape as winter grazing grounds for reindeer.Yet, such mitigation measure could come at the expense of some production loss for forestry.Such trade-offs are indispensable in forests managed for multiple interests.The cumulative effects of climate change and land use change by forestry, wind power development and mining on the winter grazing areas have become the most serious threats to reindeer husbandry in accordance with Sámi self-determination (Larsson Blind, 2022).This underscores the urgent need for improved decision-making procedures within landscape planning to overcome present misalignments of spatial and temporal reference points between traditional reindeer husbandry and extractive forms of land use, as well as a power imbalance to influence decision-making in land use planning that currently compromise reindeer herders' options to negotiate or reaching consensus (Widmark and Sandstrom, 2012;Widmark, 2019).Here, awareness of social--ecological mismatches in risk perception and valuation of exotic species, including their unintended side-effects, are inevitable in negotiations to enable sustainable multi-functional landscape forest landscapes.This implies recognition of movements between different grazing areas beyond a socio-cultural marker of Sámi reindeer husbandry, but as key feature of the ecology of reindeer and the pronounced seasonality of northern forest ecosystems.

Fig. 1 .
Fig. 1. (A) The reindeer husbandry area (dark grey) and the three northern most communities in Sweden (1 -Norrbotten, 2 -Västerbotten, 3 -Jämtland).The herding community in this study highlighted in white.(B) Seasonal grazing areas in Vilhelmina norra RHC and major infrastructure.(C) A detail from the study area, showing habitat classes, major roads and reindeer GPS positions.For clarity, stands with P. contorta are only shown according to the species composition.

Fig. 2 .
Fig. 2. Habitat selection of reindeer at the second order (blue, filled circle) and third order (red, open circle) of selection.Negative coefficients indicate lesser use compared to the reference class for categorical variables (panel a and b), or lesser use with decreasing values for continuous variables (panel c).95 % confidence intervals not crossing the zero-line indicate significant relationships.(a) Land cover classes, mature P. sylvestris as reference class.(b) Topographic position index, flat areas as reference class.(c) Continuous variables.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3 .
Fig.3.Area weighted mean lichen cover for the different forest classes, indicated by the white diamond symbol.Classes with the same letter do not differ significantly from each other.Note the square root scale of the y-axis to include outliers (grey).

Fig. 4 .
Fig. 4. Correspondence analysis plot of the association between forest classes (blue) and soil type (red).The closer a forest class is to a soil moisture class, the closer they are associated with each other.