Reindeer control over shrubification in subarctic wetlands: spatial analysis based on unoccupied aerial vehicle imagery

Herbivores can exert a controlling effect on the reproduction and growth of shrubs, thereby counter‐acting the climate‐driven encroachment of shrubs in the Arctic and the potential consequences. This control is particularly evident in the case of abundant herbivores, such as reindeer (Rangifer tarandus tarandus), whose grazing patterns are affected by management. Here, we tested how different reindeer grazing practices on the border between Finland and Norway impact the occurrence of willow (Salix spp.) dominated patches, their above‐ground biomass (AGB) and the ability of willows to form dense thickets. We used a combination of multispectral and RGB imagery obtained from unoccupied aerial vehicles field data and an ensemble of machine‐learning models, which allowed us to model the occurrence of plant community types (Overall accuracy = 0.80), AGB fractions (maximum R2 = 0.90) and topsoil moisture (maximum R2 = 0.89). With this combination of approaches, we show that willows are kept in a browsing‐trap under spring and early summer grazing by reindeer, growing mostly small and scattered in the landscape. In contrast, willows under the winter grazing regime formed dense stands, particularly within riparian areas. We confirm this pattern using a random forest willow habitat distribution model based on topographical parameters. The model shows that willow biomass correlated with parameters of optimal habitat quality only in the winter grazing regime and did not respond to the same parameters under spring and summer grazing of reindeer.


Introduction
With warmer temperatures, shrubs are becoming increasingly dominant across northern peatlands (Forbes et al., 2010;Mekonnen et al., 2021).In arctic and subarctic habitats, this phenomenon, referred to as "shrubification", is mostly driven by willows (Salix spp.).The impacts of shrubification on tundra wetlands are manifold and complex.Shrubs play a key role in peatland greenhouse gas balance and carbon accumulation, increasing peatland net CO 2 sink, particularly in warmer conditions (Ward et al., 2013).Furthermore, shrubs produce high amounts of fine roots, thereby increasing belowground carbon accumulation (Malhotra et al., 2020).However, shrubification of peatlands could also feedback in the opposite direction: the decomposition rate of aboveground litter is enhanced in the presence of shrubs (Ward et al., 2013), and shrubs, as well as graminoids, may enable a warming-induced loss of ancient carbon from deeper peat via priming (Walker et al., 2016).Furthermore, by competing for space, light and nutrients with the peat-producing Sphagnum mosses, the long-term carbon accumulation as peat may be reduced with shrubification (Larmola et al., 2013).Shrubification may also have an impact on soil moisture regimes through microclimate alterations such as reduced albedo, rainfall interception and increased evapotranspiration (Myers-Smith et al., 2011;Te Beest et al., 2016;Zwieback et al., 2019).In turn, altered moisture dynamics could trigger further changes in ecosystem function and trait shifts (Bjorkman et al., 2018).However, with the narrow understanding of the spatial nature of willow encroachment dynamics, there has been limited extrapolation of the potential impacts of willow encroachment at the landscape scale.
Previously, the distribution and habitat preferences of willows in the Fennoscandian tundra have been linked to local topography, soil moisture regime (Pajunen et al., 2010;Tape et al., 2012) and recently also to herbivory (Bråthen et al., 2017;Kitti et al., 2009).On wet patches, such as flood-prone river plains and wetlands, both large and small herbivores have been shown to selectively forage on the nutrient-rich leaves and the emerging shoots of willows (den Herder et al., 2008;Ravolainen et al., 2011), controlling the growth and reproduction of willows (Christie et al., 2014;Kolari et al., 2019).Moreover, it has been previously shown how reindeer grazing could potentially lead to increases in summer albedo and higher soil moisture by reducing shrub abundance (Te Beest et al., 2016).Across the Arctic, the herbivory exerted by the abundant reindeer (Rangifer tarandus tarandus) could have a substantial impact, and therefore the grazing patterns of wild and semi-domesticated animals should be included in predictions of willow biomass change (Yu et al., 2017).In Fennoscandia, these patterns are tightly linked to the Centuries-old reindeer husbandry regimes, with selective foraging of willows taking place particularly on the grazing grounds in spring and early summer (Kitti et al., 2009).
To improve understanding of the influences of reindeer grazing patterns on the observed trends in shrubification and 'arctic greening' (Berner et al., 2020;Myers-Smith et al., 2020;Verma et al., 2020), novel tools have been provided by the recent emergence of unoccupied aerial vehicles (UAVs) and lightweight sensors.These allow for the study and characterization of the structure, composition and functions of tundra ecosystems (Assmann et al., 2020;Beamish et al., 2020) at very high ground sampling distances (GSDs).In contrast with the medium and coarse GSDs supplied by satellites, the centimeter-scale imagery provided by UAV platforms reflects the spatial nature of ecosystem structures, processes and functions.For instance, Siewert and Olofsson (2021) used changes in the normalized difference vegetation index (NDVI) as a surrogate of herbivory effects on arctic vegetation, whereas Eischeid et al. (2021) used four optical bands, the normalized vegetation index (NDVI), a digital surface model (DSM) and textural metrics derived from a UAV to map herbivory-induced disturbances on arctic vegetation.These advances indicate the potential of the combination of UAVs, optical sensors and machine learning algorithms to provide accurate spatially explicit representations of the effects of reindeer herbivory on Salix encroachment in tundra wetlands.
In this article, we aim to model the contrasting distribution and effects of spring-summer and winter reindeer herbivory on Salix biomass in oroarctic fens using multispectral and RGB imagery from sensors mounted on a UAV in combination with field observations.To assess reindeer impact, we use a natural experiment that captures the impacts of two different long-term reindeer management regimes at the Finnish-Norwegian border.
Here, the Finnish side is grazed in spring and early summer, whereas the Norwegian side has only been grazed during the winter months.We hypothesize that compared to winter-only grazing, the spring and early summer grazing by reindeer leads to a decrease in the occurrence of willow-dominated patches, their above-ground biomass (AGB), and the ability of willows to form dense thickets.Further, we test whether Salix abundance affects topsoil moisture levels.Finally, we used species distribution models (SDMs) to verify and visualize the controlling effect of reindeer grazing on shrub expansion.

Study area
The Jávrrešduottar area lies at the border between Norway and Finland (Fig. 1) within the Fennoscandian oroarctic tundra zone (Virtanen et al., 2016).The landscape in the area is characterized by a mosaic of heaths with Betula nana and Cladonia lichens, and poor to moderately rich fens in the depressions dominated by sedges (e.g., Carex rostrata, Carex rotundata and Eriophorum angustifolium) and bryophytes (e.g., Sphagnum riparium, S. lindbergii, S. teres, Paludella squarrosa and Sarmentypnum sarmentosum) (Kolari et al., 2019).The main willow species in the study area, Salix lapponum, forms dense stands in more favorable areas accompanied by Salix glauca and S. phylicifolia.
Reindeer (Rangifer tarandus ssp.tarandus) husbandry is an important form of land use in the Jávrrešduottar area, both in Finland and Norway.In 1852, reindeer movement across country borders was prohibited and since then, the two sides of the border fence have been subjected to different grazing management practices.In the mid-1950s, a fence was constructed along the border to prevent reindeer movement across countries.Since then, the border area in Norway is only used as reindeer winter grazing grounds, when the ground is mostly snow-covered and reindeer feed almost exclusively on Cladonia spp.lichens.Under this regime, willows in the Norwegian side of the border remain largely undisturbed by reindeer.The border area in Finland has been used as spring and summer pasture since the 1930s (Näkkäläjärvi, 2013), and reindeer forage on willows, dwarf birches, forbs and graminoids in the fens (Kitti et al., 2009), having a stronger impact on fen vegetation.
Within the Jávrrešduottar area, we selected two study sites due to the representativity of the fens (in terms of size and species composition) found in them.The sites are located in the Eastern and Western flanks of the Má deroaivi hill along the border fence (Fig. 1B) and are referred to as Má deroaivi East and Má deroaivi West hereinafter.

Overview
In this study, we used a combination of field data and multispectral and RGB images derived from an UAV to assess the landscape scale effects of long-term reindeer herbivory on Salix in two oroarctic fens located at the border between Finland and Norway.Firstly, we masked out the heath tundra areas of the landscapes under study, which allowed us to focus exclusively on oroarctic mires (Additional explanations on the masking method and results are provided in Supplementary Materials).We modeled differences in occurrence, extent, and AGB of Salix and topsoil 689 moisture across the border fence utilizing an Extreme Gradient Boosting Algorithm (XGBoost).Based on the Salix extent maps, we applied a set of landscape indices to characterize the spatial patterns of Salix growth through the landscape.Finally, we modeled and mapped Salix potential habitat suitability in the study areas and analyzed against the actual primary productivity of Salix using spatially binned generalized additive models (GAMs).Figure 2 outlines a summary of the methodological steps undertaken to achieve the objectives of this study.

In situ vegetation and soil moisture data collection
We identified six plant community types in the field, based on plant functional types (Chapin III et al., 1996;Mekonnen et al., 2021) and microtopography: (1) Sedge wetland, (2) Sphagnum spp.lawn dominated by Sphagnum lindbergii, (3) Salix patches, (4) Betula nana patches, (5) hummock tops with ericaceous dwarf shrubs and Sphagnum fuscum, and (6) springs with Philonotis spp.At each study site, 10 sampling plots were randomly assigned to each class in the field, five plots within the spring-summer grazing area and five within the winter grazing area.Five sample plots were assigned to the "Spring side with Philonotis spp." class, as it was found only at Má deroaivi East, within the spring-summer grazing area.This resulted in a total of 105 sampling plots (Fig. 3).Within each plot, a 50 × 50 cm pointintercept frame with 15 points was used to record the number of species hits, and a Delta-T WET Sensor (Delta-T Devices Ltd., London, UK) was used to measure topsoil moisture at 10 cm depth (% volumetric water content).Three topsoil moisture measurements were taken and averaged within each quadrat.Finally, above-ground vascular plant biomass (AGB) was harvested.The AGB was dried at 40°C for 72 h, sorted based on species or species groups and by fractions of woody and green and then weighted.The location of each plot was recorded using a RTK GNSS (Topcon Hiper V).We collected field data between 17th and 21st of July 2021.To rule out spatial autocorrelation among the field plots, we tested each sample set (Má deroaivi East and Má deroaivi West) using Moran's I with a k-nearest neighbors connectivity.The analysis was undertaken using the spdep package (Bivand et al., 2015) in R Studio v 1.4.1106.

UAV data collection and processing
We carried two UAV surveys at each study site using a Sensefly Ebee UAV in mid-July 2021 equipped with a Parrot Sequoia 1.2 megapixel monochromatic multispectral sensor and a SenseFly S.O.D.A camera.The final GSD was 11 cm per pixel for multispectral images and 3.5 cm per pixel for RGB images.The multispectral images were collected in four spectral bands: Green (530- 570 nm), red (640-680 nm), red edge (730-740 nm), and near-infrared (770-810 nm).Before each multispectral flight, we radiometrically calibrated the Parrot Sequoia sensor using an Airinov calibration panel.We defined the flight areas at each study site using a buffer of c. 500 m along the grazing fence.The total surveyed area was 155 ha at Má deroaivi West and 94 ha at Má deroaivi East.We processed both the multispectral and RGB imagery sets in eMotion 3 ® using a post-process kinematic correction, aimed at accurately geolocating the images.RINEX observation and navigation files from the FinnRef CORS network were used to increase the positional accuracy (Zhang, Aldana-Jague, et al., 2019).We then generated RGB and multispectral orthomosaics in Pix4D v.4.3.31® , resulting in one RGB and four spectral orthomosaics per study site.We used four ground control points (GCPs) and four additional natural features to estimate the positional accuracy of the mosaics.

Vegetation indices and spectral separability
Using the Python console in ArcMap 10.5, we calculated an initial set of 15 vegetation indices (VIs; Table 1) to capture a wide range of ecosystem biophysical properties, structures, and functions (Huete, 2012) while simultaneously avoiding spectral saturation and sensitivity to background factors associated with individual VIs (Huang et al., 2021).The selection of VIs in this study was driven by their relationships to different vegetation properties, such as leaf area, AGB, water content, and the content of chlorophyll and other pigments.VIs were complemented with textural indices based on the gray-level cooccurrence matrix (GLCM).We used the glcm package (Zvoleff, 2020) in R Studio v 1.4.1106 to calculate GLCM dissimilarity (Table 1).GLCM dissimilarity based on the red-edge band mosaics estimating the gray level distance between pairs of pixels and quantifying similarities and dissimilarities between pairs of pixels (Park & Guldmann, 2020).
To avoid potential multicollinearity among variables, we used the M-statistic (Torres-Sánchez et al., 2013) as a tool to evaluate the spectral discrimination capacity of the VIs.The M-statistic computes the difference between the mean values of a VI corresponding to two landcover classes, normalized by the difference of the corresponding standard deviations, as follows: We applied the M-statistic to all the possible pairs of classes formed by the six plant community types in our study.We extracted the values of each of the 15 VIs within all sampling quadrats and selected only those that yielded the highest M value for each pair of plant community types.

Canopy height models
While spectral information and VIs are useful to estimate certain vegetation properties, they may show shortcomings in estimating AGB (Passalacqua et al., 2021).To complement the spectral dataset, we implemented a Structure-from-Motion (SfM) algorithm to model vegetation height in the study sites.We used Pix4Dmapper v.4.3.31® to process all the RGB images collected via the SenseFly S.O.D.A camera, utilizing a combination of SfM  G: green band (530-570 nm), NIR: near-infrared band (770-810 nm), R: red band (640-680 nm), rededge: Red edge band (730-740 nm), P(i, j): frequency that two pixels occur in the image, with gray levels i and j.
algorithm with the Multi-View stereo photogrammetry to detect common features in the images and reconstruct 3D point clouds (Smith et al., 2016;Villoslada et al., 2021).
We then imported the dense 3D point clouds into Cloud-Compare 2.11.3,where we applied a Cloth Simulation Filter to categorize 3D points as ground and non-ground.
We computed a digital terrain model (DTM) and a DSM based on the classified point clouds in CloudCompare.We subtracted the DTM from the DSM to generate the vegetation height model and used the vegetation heights measured at the plots located within Salix and Betula nana patches to assess the accuracy of the model.

Models for the vegetation classification and prediction of AGB and topsoil moisture
For mapping the five plant community types under assessment, we tested two different machine learning algorithms: Random Forest (RF) and Extreme Gradient Boosting (XGBoost).Due to its ability to predict beyond the limits of the training dataset, we used XGBoost to model AGB fractions and topsoil moisture.
We implemented the RF and the XGBoost in R (version 4.1.1;R Core Team, 2021), using the plot-level vegetation, AGB, and topsoil moisture data collected at the study sites as the training/validation datasets.All pixels intersecting each of the sampling polygons were assigned the corresponding plant community type, AGB and topsoil moisture values and used them as training/validation dataset by extracting each intersecting pixel center point.This approach is essential to capture the spectral diversity of both the individual sampling plots and the plant communities under assessment.To avoid imbalanced sampling, we restricted the number of pixels per sampling plot to 20.We used a split-sample approach with a random 70%/ 30% split for training and validation, resulting in 1470 and 630 pixel points, respectively, both for the RF and the XGBoost methods.We used the VIs selected by the Mstatistic and the vegetation height models as co-predictors.The RF model was built using the R packages raster (Hijmans et al., 2015) and ranger (Wright et al., 2019).The number of trees (ntree) within the RF ensemble was set to 600 and the number of randomly selected predictors at each split (mtry) was set to 4. We built the XGBoost model using the raster and xgboost packages (Chen et al., 2019).We optimized the XGBoost hyperparameters [learning rate, number of trees, minimum number of samples required at a leaf node, maximum depth, the number of features for the best split and gamma (γ)] (Pham et al., 2020) using the mlr package (Bischl et al., 2016), using a fivefold cross-validation and 100 optimization rounds.We assessed the contribution of each input variable in the classification performance using the mean decrease in accuracy in the RF (Villoslada et al., 2020) and the Gain metric within the XGBoost (Chen et al., 2019).
At each study site, we generated five models (Plant community distribution, total vascular AGB, total woody AGB, total vascular green AGB and topsoil moisture), and validated them using the validateMap function within the RStoolbox package (Leutner et al., 2017).For each of the plant community classification maps generated, we selected the best-performing model based two validation metrics: overall accuracy and kappa.We further evaluated per-class accuracies using sensitivity (i.e., the rate of true positives) and specificity (i.e., the true negative rate in a confusion matrix; Nhu et al., 2020).We validated the AGB fractions and topsoil moisture regression models and compared the prediction accuracies using the coefficient of determination (R 2 ) and the root mean square error (RMSE).
Finally, we used spatial overlaps between the plant community distribution maps and the XGBoost regression models to assess differences in total vascular AGB and topsoil moisture across plant community types.To detect significant differences between plant communities, we used a Welch ANOVA with Games-Howell pairwise post-hoc tests.

Effects of differential grazing on Salix biomass and topsoil moisture
To characterize the effects of reindeer grazing on Salix growth across the fence, we overlaid the plant community types and AGB and topsoil moisture raster layers obtained in the previous steps and extracted AGB and topsoil moisture values corresponding to each Salix pixel.To avoid sampling bias, we accounted for local variation in topography (Chen et al., 2021) clustering four topographic variables [DTM, Topographic Wetness Index (TWI), Topographic Position Index (TPI), and distance to water bodies] into three clusters at each study site using a k-means algorithm implemented in QGIS v3.24.3.TWI integrates catchment area, flow width and slope gradient in the calculation (Kopecký et al., 2021) and TPI measures the relative position of a raster cell within a buffer of neighboring cells (De Reu et al., 2013).We calculated the DTM, TPI and distance to water bodies in QGIS v 3.24.3 and the TWI in SAGA GIS v 7.8.2.For subsequent analyses, we randomly selected Salix pixels proportional to the size of the topography clusters on each side of the grazing fence, starting with all Salix pixels falling within the smallest sampling cluster at each study site.We then performed a Welch one-way ANOVA to detect significant differences in Salix AGB fractions and topsoil moisture between winter grazing and springsummer grazing regimes.Before the analysis, we checked the variables for normality and transformed them when necessary, via the bestNormalize package in R 4.1.2(Peterson & Peterson, 2017).

Salix growth patterns, distribution range and AGB
We transformed the plant community type map into binary raster masks displaying Salix occurrence only.To describe Salix growth patterns throughout the landscape in response to different grazing strategies, we used the R package "landscapemetrics" (Hesselbarth et al., 2019).The package quantifies the configuration and spatial arrangement of patches of contiguous cells belonging to the same class in a raster layer.We selected four metrics that characterize different elements of landscape heterogeneity, patch size and density, namely LPI, APA, PD and LSI.
Assuming that Salix within the winter grazing occupy the full range of their habitat, we established a Salix distribution model (SDM) based on presence-absence data of willows on the Norwegian side of the fence.We derived Salix presence and pseudo-absence using the mask obtained in the previous step.We generated a set of 500 random presence points at Má deroaivi East and 580 at Má deroaivi West, proportional to the area covered by Salix.Additionally, 5000 random pseudo-absence points were generated at each site.Both sets of presence and pseudo-absence points were generated within the winter grazing areas.To avoid double-sampling individual willows, we set a minimum distance of 2 m between points.We created a set of Salix distribution co-predictors, based on the DTM and the landcover maps: DTM, distance to water bodies, TWI and TPI.Local topography and moisture regime have been highlighted as key drivers of willow distribution in the Fennoscandian tundra (Bråthen et al., 2017;Pajunen et al., 2010).
Following the workflow proposed by Valavi et al. (2021), we tested two approaches to model Salix distribution using presence and pseudo-absence data: Equalsample RF and Down-sampled RF.In the equal-sample RF, we fitted the model on 50 datasets, each with the same number of presences and pseudo-absences, averaging the scores of the 50 models for the final prediction.In the down-sampled RF, the bootstrapping of each tree used the same number of presence and pseudo-absence samples.The SDMs were trained using the presence and pseudo-absence dataset generated on the Norwegian side and subsequently predicted over the whole extent of the study areas as raster maps where each cell represents the probability score (0-1) of Salix occurrence.This score is considered a proxy for species habitat suitability (Valavi et al., 2021).We evaluated the performance of each model using a separate presence and pseudo-absence validation dataset (200 presence and 200 pseudo-absence points at each study site) and the area under the receiver operating characteristic curve (AUC).The SDMs were built using the R package randomForest (Liaw & Wiener, 2002) and validated using the R package pROC (Robin et al., 2021).We evaluated the contribution of each co-predictor to the models using the mean decrease in Gini coefficient (Valavi et al., 2021).
To explore the control exerted by reindeer on Salix abundance and growth, we aggregated both total AGB and the average SDM scores in hexagonal bin (hexbin) maps.We generated hexbin maps of 25 m diameter hexagons in QGIS v 3.24.3.Within each hexagon containing Salix, the sum of total Salix AGB of all raster cells within the hexagon was computed, as well as the average SDM score.We then used GAMs with a restricted maximum likelihood approach (Simpson, 2018) to characterize the relationships between Salix productivity and habitat suitability.GAMs were calculated separately for spring-summer grazing and winter grazing.

Spatial accuracy, co-predictors selection and spatial autocorrelation
The comparison of eight known GCPs at the study sites with their location on the multispectral orthomosaics yielded a mean AE SD horizontal error of 9 AE 2.4 cm.The vertical error of the DSMs was 5.5 AE 4.2 cm (comparison of the RTK GNSS receiver-recorded elevation at the sampling quadrats and the elevations generated by the SfM model).The M values for vegetation and textural indices ranged between 0.00 and 2.62 (Table S2).Based on the results of the M-statistic, we selected the following six vegetation and textural indices: GLCM-dis, CVI, Datt4, GRDI, GRVI and NDVI.None of the filed plots sets exhibited spatial autocorrelation (Moran's I P-value [ 0.05).

Distribution of plant community types, AGB and topsoil moisture
Within the study areas, Salix appeared along water courses (Fig. 3), where also the highest AGB was found (Figs. 4 and 5).This pattern was particularly distinct in the winter grazing regime.Salix also occurred within sedge fens and Sphagnum patches.
While XGBoost yielded higher overall accuracy and kappa values than RF at Má deroaivi East, RF outperformed XGBoost at Má deroaivi West (Table 2).The models were able to forecast Salix abundance with a sensitivity and specificity of 1.00 and 0.89 at Má deroaivi East, and 0.69 and 0.97 at Má deroaivi West.Specificity was higher than sensitivity for all plant community types at both study sites (Table S1).The lowest sensitivities recorded corresponded to Sphagnum lawns (0.46), sedge wetland (0.44) and spring (0.33) at Má deroaivi East.In the classification of plant communities, the vegetation index Datt4 was the most important co-predictor at both study sites (Figure S1).
The XGBoost model provided robust descriptions of AGB fractions and topsoil moisture patterns (R 2 0.66-0.90;Table 3), with the prediction accuracies being generally higher at Má deroaivi East.At both study sites, accuracies gradually increased from total vascular AGB to total vascular green AGB.In general canopy height and GNDVI stood out as key co-predictors in the topsoil moisture and AGB fractions models, although copredictor contributions varied depending on the AGB fraction (Figures S2 and S3).
There was a clear decreasing trend in AGB from communities dominated by deciduous shrubs (Salix and Betula nana, with mean total AGB values of 595.12 and 385.35 g m -2 , respectively) to sedge-dominated wetland communities and ultimately vascular biomass in Sphagnum moss lawns (mean total AGB of 113.79 g m -2 ) (Fig. 6).On the other hand, sedge-dominated wetlands and Sphagnum lawns occupy areas characterized by very high soil moisture, while willows were located at intermediate moisture levels (62.8% on average).Hummock tops represented the driest habitats, and dwarf birch occupied the intermediate sides of the hummocks.
The Welch ANOVA confirmed the difference across plant communities for total AGB (F 4, 28384 = 10 112; P \ 0.0001) and topsoil moisture (F 4, 27092 = 8977; P \ 0.0001).According to Games-Howell pairwise post-hoc test, total AGB of all plant communities differed from each other.For the topsoil moisture, all communities except for the sedge wetland and Sphagnum lawns differed significantly from each other.

Effects of reindeer summer grazing on Salix biomass and topsoil moisture
Above-ground biomass (all fractions and total) was significantly higher under the winter grazing regime as compared to the spring-summer grazing regime (Table 4; Fig. 7).In contrast, soil water content within the Salixdominated pixels was significantly higher within the spring-summer grazing regime (mean value of 92.17%) than within the winter grazing regime (mean value of 71.93%), as shown in Table and Figure 7.
Notably, the distribution of AGB was multimodal within the winter grazing regime (Fig. 7), as opposed to the unimodal or nearly unimodal distribution in the spring-summer grazing regime.Conversely, the shape of the distribution in the spring-summer grazing regime  shows the largest concentration of willows around the highest topsoil moisture values.

Landscape patterns of Salix growth
The results of the landscape metrics assessment showed clear trends regarding willow spatial distribution within the spring-summer and winter grazing regimes (Table 5).
The largest patch index (LPI) was higher in the winter grazing than in the spring-summer grazing regime (with values of 500 vs. 2289 m 2 ; Table 5).The winter grazing regime also showed a larger average patch area for willows (APA) and lower patch density (PD) and landscape shape index (LSI).These results indicate a higher degree of fragmentation and landscape heterogeneity in the spring-summer grazing regime, characterized by a higher density of small, scattered willows.On the contrary, in the winter grazing regime, willows appeared clustered and showed a tendency to form thickets with closed, dense canopies (Fig. 8).

SDMs and total Salix AGB
The best-performing SDM was the down-sampled RF, which yielded AUC values of 0.95 for Má deroaivi West and 0.94 for Má deroaivi East.In both Má deroaivi West and East, the best-performing predictors for willow biomass were elevation (DEM) and the distance to water bodies (Figure S4).The hexbin aggregation of total Salix AGB and average SDM score (Figure S5), followed by single-variable GAMs (Fig. 9), indicate an evident limitation to the growth of Salix within the spring-summer grazing regime.Within the winter grazing regime, willows occupy the areas of higher habitat suitability following the expected trend, with a positive relationship between total willow AGB and the SDM score.Here, the single-variable GAM was able to explain 35% of the deviance (P \ 0.001).The same pattern was, however, not found in the spring-summer grazing regime.Here, willows barely responded to increases in habitat quality, with the GAM showing a nearly flat trend curve (Fig. 9) and a smaller explanation of deviance (6.44%, P \ 0.001).

Discussion
Plant community distribution models Our workflow produced highly accurate models of plant community distribution, with minimal differences between the accuracies yielded by the XGBoost and the RF algorithms (Table 2).Per-class statistics (Table S1) revealed contrasting trends at both study sites.Sensitivity was the highest (1.00) at Má deroaivi East and relatively high (0.69) at Má deroaivi West for Salix patches.We attribute this difference to the fact that the training data as well as the UAV images were collected over a considerably larger area at Má deroaivi West, encompassing a larger spectral variability and consequently lower classification accuracies (Villoslada et al., 2020).Sphagnum lawns at Má deroaivi West and sedge wetland at Má deroaivi East present relatively low sensitivities (0.58 and 0.47 respectively).This can be explained by the complex spectral characteristics of different Sphagnum species (abundant also in the sedge wetland class), influenced by factors such as moisture content (May et al., 2018), moss layer thickness (Mikola et al., 2018) and differences in pigmentation (Vogelmann & Moss, 1993).
Effects of reindeer summer grazing on Salix AGB and topsoil moisture Above-ground biomass and topsoil moisture models yielded high accuracies, with better predictions at Má deroaivi East, following the same trend as for Salix classification models.The RMSE achieved by our Salix AGB models (Table 3) were generally similar or higher than those in previous studies (Aalto et al., 2021;Greaves et al., 2016Greaves et al., , 2017)).
We found that spring-summer reindeer grazing could control shrub expansion in the Arctic tundra, as highlighted by the significant differences in biomass fractions between the spring-summer grazing regime and the winter-only regime (Table 5).The biomass suppression suggested by the models aligns well with the "browse trap" concept proposed by Bråthen et al. (2017), according to which, under certain herbivore densities, a portion of shrubs are kept inside a "browse trap" unable to grow into taller form.The higher biomass fractions in the Norwegian side (winter grazing) indicate that in the absence of reindeer browsing pressure, Salix shrubs have shifted to taller stages.This idea is further reinforced by the biomass distributions revealed by the violin plots in Figure 7.In the Finnish side of the grazing fence, the AGB violin plots display a nearly unimodal distribution, indicating willow growth under spring-summer grazing is restricted.Conversely, the Salix AGB distribution in the Norwegian side of the fence shows a bimodal or multimodal distribution for all AGB fractions.We suggest that without browsing pressure, this bimodality reflects the ability of Salix willows to achieve taller stages, driven by the habitat suitability.In the Norwegian side of the study area, a share of the landscape is characterized by riparian habitats.Within these productive and nutrient-rich streambanks, and in the absence of grazing disturbance, willows form dense and tall thickets (Liljedahl et al., 2020).
Our results suggest that soils under denser willow thickets and higher willow AGB show lower topsoil moisture levels than soils under short and sparse willows (Table 5).Local hydrological and microtopographical factors may be responsible for these differences.Willows in the Finnish side of the fence are located in moderately rich fens, occurring around ponds or in local depressions, characterized by a shallow peat layer and a high water table levels.Other local microtopography elements such as hummocks play a key role in topsoil moisture regimes in wetlands through sub-meter scale differences in height and variations in soil properties (Zwieback et al., 2019).
Larger willow AGB values in the Norwegian side of the fence may, in turn, lead to drier soils.Mechanisms such as rainfall interception (Zwieback et al., 2019), reduced surface albedo (Te Beest et al., 2016), increased evapotranspiration (Myers- Smith et al., 2011), and direct water uptake by shrub roots could cause a decrease in topsoil moisture under the dense willow thickets.
The topsoil moisture violin plots in Figure 7 revealed further insights regarding the distribution of willow shrubs in the landscape.Specifically, most willows in the springsummer grazing regime are located in areas characterized by high topsoil moisture.These results suggest that wet, less accessible areas are likely to be avoided by reindeer and constitute the only sites within the spring-summer grazing where willows prevail.This potentially indicates a reindeer-induced displacement of willows out of their niche and partially explains the topsoil moisture differences between the winter and spring-summer willow populations.It has been previously suggested that the impacts of reindeer grazing are not uniform throughout the landscape and may largely depend on underlying local environmental conditions and landscape structure.For instance, Skarin et al. (2010Skarin et al. ( , 2020) ) suggest that reindeer may avoid wet, low-lying areas where insect harassment hinders foraging.Moreover, energy expenditure for locomotion is known to influence movement patterns in wild caribou populations (Mallory & Boyce, 2018).The avoidance of energy-demanding areas such as wet patches could further explain the patterns revealed by the topsoil moisture violin plots in the spring-summer grazing regime.

Landscape patterns of Salix growth and SDMs
Metrics related to the areal characteristics of patches (LPI and APA) were higher within the winter grazing regime, reflecting the ability of willows to form large thickets under limited or no grazing pressure.In comparison, willow thickets appeared fragmented in the spring-summer grazing regime as confirmed by the higher PD index and LSI, representing a lower degree of patch compaction, in line with previous work showing indications that reindeer grazing could lead to willow fragmentation (Bråthen et al., 2017;Henden et al., 2011;Ravolainen et al., 2011).
The SDMs showed the probability of Salix distribution at the study areas, highlighting low-lying depressions and riparian corridors as optimum habitats.The variable importance metrics revealed the pivotal role of soil moisture and topography in shaping the distribution of willows (Figure S4).Following our expectations, the GAMs indicated a clear effect of willow growth suppression by reindeer grazing.In Figure 9, GAMs produced an almost flat trendline in the spring-summer grazing areas, highlighting the inability of willows to reach tall life stages.The opposite trend is observed in the winter grazing regime, where willows were concentrated along water courses (Tape et al., 2012) forming dense and tall canopies.This provides strong evidence of the effects of reindeer herbivory on the distribution and growth of willows, confirming that the reindeer grazing pressure at the study sites is strong enough to limit the dominance of willows (Bråthen et al., 2017).

Conclusion
In this study, we provided a detailed map of plant community types, AGB, and topsoil moisture in two oroarctic fens located between Finland and Norway.The UAVbased model accurately showed the current distribution of Salix shrubs under two contrasting reindeer grazing strategies.In line with our hypothesis, we found that the spring-summer grazing by reindeer led to decreased occurrence and AGB of Salix when compared to winteronly grazing by reindeer.This confirms that the previous observations at the plot level (Kitti et al., 2009;Kolari et al., 2019) hold true at the landscape level.Further, we found that higher willow abundance was associated with lower topsoil moisture levels.
The use of UAV data allowed us to characterize in detail the spatial configuration of willow distribution and abundance.We found strong evidence of the effects of reindeer herbivory on the distribution of willows on the landscape.Salix shurbs formed large thickets in riparian areas, associated with the negligible browsing pressure on willows under the winter grazing regime.Conversely, a higher density of smaller patches was observed under spring-summer grazing, confirming that the reindeer grazing pressure at the study sites is strong enough to limit the dominance of willows (Bråthen et al., 2017).
Our results demonstrate that the use of highly detailed UAV data improves our understanding on the spatial nature of shrubification dynamics.This is especially important in oroarctic wetlands characterized by a fine-scale distribution of heterogeneous plant communities (Virtanen et al., 2016).Our spatially explicit datasets provide essential knowledge on herbivory controls and its effects on tundra ecosystems, confirming that, due to the omnipresent distribution of herbivores, the realized extent of shrubification might be much lower than climatically allowed.

Figure 1 .
Figure 1.(A) Location of the study sites between Norway and Finland, in the Jávrrešduottar area, (B) map showing Má deroaivi West and Má deroaivi East study sites at the border between Finland and Norway, (C) UAV rgb mosaic of Má deroaivi West, and (D) UAV rgb mosaic of Má deroaivi East.(Proj: ETRS-TM35FIN).Source of background maps: OpenStreetMap.UAV, unoccupied aerial vehicle.

Figure 2 .
Figure 2. Interconnected workflows used to assess the controlling effect of reindeer grazing on Salix expansion.ANOVA, analysis of variance; UAV, unoccupied aerial vehicle; XGBoost, Extreme Gradient Boosting.The icons in the figure were created by Freepik and Buandesign and downloaded from www.flaticon.com.

Figure 3 .
Figure 3. Random forest (A) and Extreme Gradient Boosting (B) model predictions for the occurrence and distribution of oroarctic fen plant communities at Má deroaivi West (A) and Má deroaivi East (B).(C) Closeup section showing the drastic change in the occurrence of Salix patches across the fence.The black line represents the reindeer fence that separates the areas with a decadal history of only winter grazing by reindeer (North) and spring-summer grazing (South).The areas depicted with a hillshade correspond with heath tundra.The dots represent the location of the field sampling plots.

ª
2023 The Authors.Remote Sensing in Ecology and Conservation published by John Wiley & Sons Ltd on behalf of Zoological Society of London.

694 ª 2023
The Authors.Remote Sensing in Ecology and Conservation published by John Wiley & Sons Ltd on behalf of Zoological Society of London.

Figure 6 .
Figure6.Violin plots showing average values of total AGB and topsoil moisture for each selected community type across both study sites.The symmetric curves on the sides represent the kernel density estimate of pixels corresponding to each of the communities under study.AGB fractions and topsoil moisture were modeled using an XGBoost regression algorithm.Communities are ordered in a decreasing manner, from the highest mean total AGB and topsoil moisture to the lowest.The black dots represent the mean.AGB, above-ground biomass.

Figure 7 .
Figure 7. Violin plots showing total AGB (top left), woody AGB (top right), green AGB (down left) and topsoil moisture (down right) differences across the border in willow-dominated areas.The black dot represents the mean, whereas the horizontal black line inside the boxplot represents the median.The symmetric curves on the sides represent the kernel density estimate of pixels corresponding to Salix individuals.AGB fractions and topsoil moisture were modeled using an XGBoost regression algorithm.AGB, above-ground biomass.

Figure 8 .
Figure 8. Close-up section at Má deroaivi East showing the spatial distribution of Salix across the grazing fence.The spring-summer grazing regime (Southern side of the fence) shows a more fragmented pattern of small willow patches.The distribution of Salix patches was modeled using unoccupied aerial vehicle-derived images and an XGBoost algorithm.

Table 1 .
List of vegetation and textural indices selected in the present study to classify plant communities and predict above-ground biomass and topsoil moisture.

Table 2 .
Comparison of classification accuracies (Overall accuracy and kappa) between RF and XGBoost.
ª 2023 The Authors.Remote Sensing in Ecology and Conservation published by John Wiley & Sons Ltd on behalf of Zoological Society of London.

Table 3 .
Root mean squared error (RMSE) for the predicted total vascular aboveground biomass, total woody aboveground biomass, total vascular aboveground green biomass and topsoil moisture in the XGBoost model, expressed in the same units as the modeled maps.
ª 2023 The Authors.Remote Sensing in Ecology and Conservation published by John Wiley & Sons Ltd on behalf of Zoological Society of London.

Table 5 .
Four landscape metrics calculated using the Salix coverage on both sides of the border fence.
ª 2023 The Authors.Remote Sensing in Ecology and Conservation published by John Wiley & Sons Ltd on behalf of Zoological Society of London.

Table S2 .
M-statistic values for spectral separability between plant communities.Figure S1.Variable importance metrics for the Random Forest classification in Má deroaivi West (left) and the Extreme Gradient Boosting classification in Má deroaivi East (right).Figure S2.Variable importance metrics for the Extreme Gradient Boosting regression in Má deroaivi East Figure S3.Variable importance metrics for the Extreme Gradient Boosting regression in Má deroaivi West Figure S4.Mean decrease in Gini values for all co-predictors used in the species distribution models Figure S5.Spatial prediction of RF-based species distribution models for Salix spp. at Má deroaivi West (A) and Má deroaivi East (B).
706 ª 2023 The Authors.Remote Sensing in Ecology and Conservation published by John Wiley & Sons Ltd on behalf of Zoological Society of London.