Top‐Down Regulation by a Reindeer Herding System Limits Climate‐Driven Arctic Vegetation Change at a Regional Scale

Warming‐driven growth of tall woody vegetation in the Arctic has the potential to accelerate climate change through multiple positive feedbacks. Local‐scale evidence suggests that large herbivores limit this vegetation shift, but there is uncertainty at larger, regional scales whether current herbivory pressure is a major top‐down control on ecosystem structure and functioning. Across a 67,000 km2 region of the Yamal Peninsula in West Siberia, we integrated satellite remote sensing with a novel data set mapping the migrations of herds comprising 151,000 domesticated reindeer. Where reindeer numbers varied over space, higher reindeer herbivory pressure was consistently linked with lower coverage of tall woody vegetation. Within areas dominated by this vegetation type, productivity and climate were increasingly decoupled where reindeer density was higher. Our spaceborne fingerprint detection suggests that large herbivores, at current population densities, counteract Arctic vegetation responses to climate change over large spatial scales.

The expansion of woody vegetation, especially tall shrubs, with climate change in the Arctic is well-documented (Myers-Smith et al., 2011). Studies have shown that, under observed and experimental warming, reindeer herbivory counteracts this vegetation shift, weakening or even negating the positive relationship between annual mean temperature and shrub growth and preventing the shrubline from advancing (Bråthen et al., 2017;Olofsson & Post, 2018;Olofsson et al., 2009;Post & Pedersen, 2008;Väisänen et al., 2014). Since the transition to more shrub-dominated tundra leads to lower albedo and increased permafrost thaw (Parker et al., 2021;Sturm et al., 2005), reindeer herbivory in the rapidly warming Arctic may have important Earth system implications by mitigating these positive climate change feedbacks.
Despite the strong evidence that reindeer herbivory acts as a significant top-down control when reindeer are concentrated over small areas, there is uncertainty over whether reindeer, at current population densities, have a substantial impact at larger spatial scales across the Arctic. Remote sensing is emerging as a tool for analyzing Arctic herbivory impacts, and has been used, for example, to detect a reduction in the normalized difference vegetation index (NDVI), a spectral reflectance metric for vegetation "greenness" (Raynolds et al., 2012), over a 770 km 2 area in the year after vole and lemming peaks (Olofsson et al., 2012). Several studies have found that reindeer reduce NDVI on summer pastures, but these effects were typically weak or inconsistent when comparing across differences in herd sizes or animal densities over space (Campeau et al., 2019;Newton et al., 2014;Rickbeil et al., 2015). Conversely, a study of migratory tundra caribou (same species as reindeer) herds in North America found that an increase in summer NDVI led to declines in caribou populations, likely related to the expansion of shrub species with strong antibrowsing defenses (Fauchald et al., 2017). Their results showed no evidence of a top-down effect of caribou on the vegetation, but instead implied a bottom-up effect initiated by a warming-driven increase in plant biomass. Although testing the relationship between reindeer abundances and NDVI provides an indication of how reindeer grazing is linked with vegetation productivity, analyses relying on NDVI alone cannot detect differences in vegetation cover, the main impact of herbivory found in many field studies. A lack of spatially explicit reindeer data has thus far hindered large-scale comparisons of vegetation cover across variations in grazing pressure over space or time (e.g., Rees et al., 2003;Verma et al., 2020;Walker et al., 2009).
Here, we analyze the influence of reindeer on vegetation on the Yamal Peninsula in West Siberia, where indigenous Nenets nomads manage the world's largest domesticated reindeer herds (∼225,000 animals in total) and migrate up to 1,200 km annually Stammler, 2005). Local observations on Yamal have found that reindeer promote a transformation from shrub-to graminoid-dominated tundra and linked higher grazing intensity to decreased height and growth of willow shrubs Pajunen, 2009). The interactions between reindeer and vegetation in this system are dictated by the style of reindeer herding practised by the tundra Nenets, which affects how the reindeer move across and use the landscape. As described in the Russian literature, this herding style belongs to the so-called Samoyedic type-a tradition common in the tundra zone of Europe and West Siberia, where it was also adopted by several other ethnic groups, including the Khanty and Komi-Izhma (Habeck, 2007;Istomin & Dwyer, 2021;Jernsletten & Klokov, 2002). The defining features of Samoyedic reindeer herding, common to all the local variants, are the large herds of up to several thousand reindeer, extended migrations that enable the rotation of seasonal pastures, and herders' constant control over their herds (i.e., 24 hr a day, 7 days a week). The greatest attention is required during the summer, especially during periods of intense bloodsucking insect harassment in July. The emphasis on controlling animals dictates the nomadic lifestyle of the reindeer pastoralists of the Yamal Peninsula, who live year-round in the tundra. Thus, analysis of the influence of reindeer on the ecosystem of Yamal must incorporate an understanding of the role of humans as actors directing and modifying reindeer movements in response to multiple factors (ecological, weather, social, etc.) (Dwyer & Istomin, 2008). These features of herding in Yamal distinguish the reindeer herbivory patterns from those of both wild reindeer populations and free grazing semi-domestic reindeer in Fennoscandia.
Within the Yamal Peninsula, reindeer herding collectives are either private households (most local reindeer herders) or brigades of reindeer herding enterprises-the "heirs" of Soviet-era state-owned reindeer enterprises (sovkhozy). Both types of collectives have followed relatively stable patterns of seasonal movement, inherited largely from the routes of past generations of Nenets nomads, described since the 19th century (Evladov, 1992). Dozens of households migrate along the same "corridors" with consistent movement timings from year-to-year. Although the winter pasture locations of a certain Nenets household may vary within their corridor, the nomadic route and pasture areas visited in the summer months are more stable. In fact, the local identity of most Nenets is defined mainly by their summer pastures (Volkovitskiy & Terekhina, 2020).
Given the intensity of reindeer herbivory on the Yamal Peninsula and the long-term stability of the seasonal migration patterns of the herds, we investigated whether variations in vegetation productivity and composition over space were linked to differences in the timing and density of reindeer pressure. Central to our analysis is a new data set we present that maps the annual June-August migrations of all the herds on the southern and central parts of the peninsula encompassing the Yar-Sale and Panayevsk tundra (the territory of the Yarsalinskii and Panayevskii sovkhozy, Figure 1), which has the highest summer density of reindeer on Yamal. Prior depictions of herd migrations on the Yamal Peninsula have either focused on the specific routes of a small subset of herds or the general movement patterns within the different sovkhozy (Degteva & Nellemann, 2013;Golovnev et al., 2014;Stammler, 2005), making our data set the first comprehensive map of reindeer numbers and summer ranges on an extensive region of 67,000 km 2 covering more than half of the Yamal Peninsula.
Our detailed, spatially explicit data set provides the opportunity to test hypotheses linking reindeer herbivory to vegetation responses at the regional scale. Alternative ecological approaches involving experimental manipulations are not feasible at such broad scales, and process-based models similar to those used in climate change attribution studies are rarely appropriate for biotic systems due to the need to incorporate myriad extrinsic factors with complex interactions (Parmesan et al., 2013). Thus, our approach adopts a framework along the lines of "associative pattern attribution," a method that evaluates whether multiple lines of evidence point to a coherent and parsimonious conclusion in accordance with a predicted biotic response (Hegerl et al., 2009). Our lines of evidence are derived from comparisons over space between areas with varying reindeer pressure. The analysis consists of three separate spatial comparisons: (a) between the Yar-Sale/Panayevsk tundra and the neighboring Taz Peninsula, which has considerably fewer reindeer but otherwise resembles Yar-Sale/Panayevsk geographically, geologically, and climatically (see Text S1 in Supporting Information S1 for a detailed comparison); (b) between areas in Yar-Sale/Panayevsk that are within summer herd ranges and areas outside of the migrations; and (c) within the summer ranges, considering differences in the density of animals on different pasture land. In all comparisons, we controlled for landscape position by dividing the areas into upland, slope, and lowland landform units (see Section 2, Figure S2 in Supporting Information S1) and separated the grazing areas within Yar-Sale/Panayevsk depending on the month in which they were used ( Figure 1).
We formulated two main hypotheses informed by observations of reindeer impacts at the local scale that we tested at the regional scale using each of the three comparisons. First, we hypothesized that reindeer herbivory would alter the distribution of landcover types, specifically linking lower coverage of tall woody vegetation with higher reindeer numbers. Our second hypothesis was that vegetation productivity in areas with higher reindeer numbers tundra, whereas the Taz Peninsula to the east has many fewer reindeer. The locations of towns and villages serving as reindeer herding centers, Panayevsk, Yar-Sale, Seyakha, and Nyda, are also shown. We present a new data set of polygons depicting the herd locations in June, July, and August, colored to represent the number of reindeer in each herd in 2020, for both the private and enterprise herds with summer pastures on the Yar-Sale and Panayevsk tundra. The imagery base map is from Esri and Earthstar Geographics, and the water bodies base map is derived from the MOD44W Land Water Mask product (Carroll et al., 2017). would be less responsive to climate variability, measured by the vegetation sensitivity index (VSI, after Seddon et al. (2016)). Herbivory in the Arctic tends to have the strongest effect on the species expected to be favored by warming (Olofsson & Post, 2018), so we anticipated that areas with higher herbivory pressure would have lower VSI, indicating that reindeer were constraining vegetation responses to climate change. Whereas differences in the proportions of landcover types present can be interpreted as the cumulative impacts of reindeer herbivory over time, differences in the VSI represent herbivory's effects on the responsiveness of vegetation to recent climate variability. Within each of the comparisons testing our main hypotheses, we expected landscape position and seasonal timing to mediate the magnitude and the significance of the vegetation responses to variation in reindeer pressure. Despite the notably flat terrain of Yamal and Taz, which is mostly less than 70 m above sea level, landform is an important determinant of both vegetation cover in the region (Walker et al., 2009) and the amount of time reindeer spend on a given area. We hypothesized that differences due to herbivory would be greatest in the uplands, which are preferred by herders for migrating and campsites because they provide greater insect relief for the reindeer, dry ground for the campsite, and an elevated view of the landscape (Degteva & Nellemann, 2013;Skarin et al., 2020). We also hypothesized that the pastures used earlier in the growing season would exhibit stronger responses to herbivory. Reindeer tend to prefer deciduous shrubs for their diet in the early growing season, at which point they are most vulnerable to herbivory, before the amount of nitrogen and phosphorus in their leaves declines rapidly at the end of June . By comparing these hypothesized patterns to remote sensing observations, we can address whether the spatial variation in vegetation composition and productivity over the vast 67,000 km 2 Yar-Sale/Panayevsk tundra in Yamal has the fingerprints of reindeer herbivory.

Study Area
The study area is located within the Yamal Peninsula and neighboring Taz Peninsula in the Yamal-Nenets Autonomous Okrug, West Siberia, Russia ( Figure 1). The area within the Yamal Peninsula was confined to the territory of the Yar-Sale/Panayevsk tundra, where herd migrations are typically long and linear compared with the more circular migrations in the northern Seyakha tundra (Degteva & Nellemann, 2013). The area on the Taz Peninsula was delineated based on the extent included in the Circumpolar Arctic Vegetation Map (Walker et al., 2005).
A detailed comparison of the Yar-Sale/Panayevsk tundra and the Taz Peninsula based on topography (Porter et al., 2018), geology (Drozdov et al., 2005), and temperature (Karger et al., 2017(Karger et al., , 2018 is presented in Text S1 in Supporting Information S1. Taz was chosen as a region for comparison with Yar-Sale/Panayevsk that differs mainly in the abundance of reindeer while having a similar physical environment. Overall, the whole study area is low and flat and traverses primarily bioclimate subzones D and E, although the northwest extent of Yar-Sale/Panayevsk is in subzone C (Walker et al., 2005). The land was unglaciated during the last glaciation and consists of sandy and clayey marine, alluvial, and lacustrine sediments (Walker et al., 2009). The region is underlain by continuous, ice-rich permafrost, and the vegetation consists of a mosaic of tall woody shrub tundra (dominated by Salix spp., Betula), low stature shrub tundra (Betula, ericoids), and mires dominated by sedges (Carex spp., Eriophorum spp.) (Botch et al., 1971;Forbes et al., 2009). The geomorphology is shaped by big meandering rivers, thermokarst lakes, and frequent cryogenic landslides (Walker et al., 2009).

Mapping of Reindeer Herds and Migrations
We recorded the number of reindeer in each herd on Yar-Sale/Panayevsk and mapped the extent of the area on which they graze annually in June, July, and August using official statistics, semi-structured interviews, and participatory mapping with indigenous experts. The official statistics for the Yar-Sale herding enterprise were based on corral (round-up) calculations of animals in the enterprise herds in 2019-2020. Interview participants included "zootechnicians," experts who are consulted on reindeer herding matters such as breeding and seasonal movements, from the Yar-Sale and Panayevsk herding enterprises as well as employees of the Kharp obshchina, a community structure facilitating the cooperation of private reindeer herders. The Kharp obshchina is the largest in Yamal, accounting for the ownership of more than 60,000 reindeer. All applicable international, national and/ or institutional guidelines for the social research were followed. Interviewees consented to be interviewed and have approved the information and the sharing of it publicly.
We used printed GGTS (State GIS Center) maps (1:200,000) to mark the herd locations during the interviews. After separately recording the reindeer numbers for each herd, we checked their totals against the 2019-2020 reconstructed official estimate of all reindeer on Yar-Sale/Panayevsk (155,000-160,000 animals). We counted 37,747 reindeer in Yar-Sale brigades, 19,700 in Panayevsk brigades, and 94,000 in private herds, giving 151,447 total reindeer. Since we mapped all major herds, we presume that the reindeer constituting the small discrepancy between our count and the reconstructed estimate are spread over the whole region.
After digitizing the mapped routes into polygons for each herd's location in each month, we created monthly maps giving the number of reindeer across the Yar-Sale/Panayevsk tundra. Animal densities were calculated by dividing the number of reindeer by the amount of land area (excluding water bodies) within their respective range for each month based on the 20 m landcover classification of Bartsch et al. (2019). We dealt with overlapping polygons by adding the reindeer numbers and densities, but since this inflated the amount of reindeer pressure occurring within small intersecting areas at the edge of the ranges of multiple herds, we discarded any new polygons smaller than 10 km 2 (<2% of the total pasture area in total) in subsequent analyses in which reindeer density was considered.
Many fewer reindeer graze on the Taz Peninsula compared with Yar-Sale/Panayevsk. The most recent published information on reindeer herds grazing in the summer on Taz within our study area estimated their numbers at ∼10,000, although this is reflective of the 1980-1990's (Kvashnin, 2009). Our conversations with local people in Nyda, the reindeer herding center in the region, confirmed that the number of reindeer in Taz is significantly lower than in any part of the Yamal Peninsula.

Landcover Classification
The mapping of landcover types used the landcover classification of Bartsch et al. (2019) (Figure S1c in Supporting Information S1). This landcover product was developed using a supervised classification of Sentinel-1 and Sentinel-2 imagery from 2016 to 2018 and groundtruthed by local experts. It contains 18 terrestrial classes and three aquatic classes. We merged the terrestrial classes into six broad classes, including "sparse vegetation," "low stature shrub tundra," "tall woody vegetation," "meadow," "inundated/landslides," and "barren/rare vegetation." The associations between the original class definitions and the landcover classes we used in our analysis are presented in Table S1 in Supporting Information S1. To facilitate our adaptation of the landcover classes, we compared the Bartsch et al. (2019) landcover map with the 1 km resolution Circumpolar Arctic Vegetation Map (Raynolds et al., 2019). This comparison helped us determine, for example, which classes aligned with our definition of tall shrubs as exceeding 0.4 m height.
"Low stature shrub tundra" includes areas with a mixture of graminoids, forbs, erect dwarf shrubs (0.1-0.4 m), and prostrate dwarf shrubs (<0.1 m) (Myers-Smith et al., 2011). This class is distinguished from "tall woody vegetation," which consists of tall shrub tundra as well as deciduous, coniferous, and mixed forest, although tree-dominated pixels cover less than 4% of the study area and are exclusively concentrated in the far south of the area. The "sparse vegetation" class describes land on mostly sandy soil with few shrubs and potentially dry biological crusts. The "inundated/landslide" class consists of areas within disturbed habitats that are seasonally inundated or within landslide scars. Landslide scars are commonly found on Yamal, and landslides are increasing in frequency with ongoing rapid permafrost thaw (Leibman, 1995;Walker et al., 2009). In our analyses, we discarded pixels in the "meadow" and "barren/rare vegetation" classes because they represented no more than 8% of the total pixels within any comparison. The meadow pixels were primarily confined to riparian areas rather than occurring throughout the landscape where they might be indicative of herbivory-driven transitions from moss-and shrub-dominated tundra to grasslands ).

Landform Delineation
Our landform classification was performed by applying a modified version of the Geomorphons algorithm (Jasiewicz & Stepinski, 2013) to the ArcticDEM v3.0 release (Porter et al., 2018) at 32 m resolution ( Figure S2 in Supporting Information S1). The algorithm assigns each pixel to one of 10 common landform elements (e.g., ridge, pit) in a procedure that automatically chooses the most appropriate spatial scale for landform delineation at each location. We further divided the pixels assigned to the flat landform class into flat uplands and flat lowlands based on the gradient in their neighboring pixels, and then reclassified the remaining landform pixels into "uplands," "slopes," and "lowlands."

Vegetation Sensitivity Index (VSI)
We adapted the VSI from Seddon et al. (2016) to quantify the sensitivity of vegetation productivity to variability in temperature, cloudiness, and water availability. We computed the VSI over the study area at 1 km resolution using monthly time series of remote sensing data from 2000 to 2019 ( Figure S1a in Supporting Information S1). The model used the moderate-resolution imaging spectroradiometer (MODIS) enhanced vegetation index (EVI) MOD13A2 version 6 product (Didan, 2015) as the response variable representing vegetation productivity. The three climate predictor variables were monthly maximum temperature (obtained from CHELSA version 2.1 (Karger et al., 2017(Karger et al., , 2018), ratio of actual to potential evapotranspiration (AET/PET) for water availability (obtained from MODIS MOD16A2GF version 6 (Running et al., 2019)), and cloudiness fraction (obtained from the MODIS MOD35L2 version 6.1 cloud mask product (Ackerman, 2015)). All time series were transformed to z-score anomalies using monthly means and standard deviations, and the months were filtered based on a minimum EVI (0.1) and temperature threshold (0°C) to identify the growing season for each pixel. Principal component regressions were used to identify the relative importance of the predictor variables in driving variations in EVI for each pixel, and these coefficients were multiplied by the log 10 -transformed ratios of EVI variability and variability in each climate time series to arrive at sensitivity. The sensitivity (VSI) is a relative metric, scaled from 0 to 100 over the study area, and was mapped in the MODIS sinusoidal projection.

Statistical Analyses
To address our first hypothesis, that reindeer had influenced landcover, we analyzed the variation in landcover distribution with herbivory pressure. We compared the proportions of pixels in the four selected landcover classes between regions (Yar-Sale/Panayevsk and Taz) and between different pasture areas within Yar-Sale/Panayevsk by calculating odds ratios. These comparisons were performed separately for pixels within "upland," "slope," or "lowland" landforms. We used log-linear modeling to test whether reindeer pressure (categorically differentiated by region or by pasture area), landform, and the interactions between them significantly affected the proportions of sparse vegetation and tall woody vegetation pixels. We tested a suite of possible models each representing different hypotheses regarding which factors were significant predictors using the glm function from the stats package (R Core Team, 2020) with a Poisson variance and log link function. We evaluated the fit of each model using the residual deviance statistic G 2 (von Eye et al., 2012).
For our analysis of landcover across variations in reindeer density within the area of Yar-Sale/Panayevsk used as summer pastures, we quantified landcover anomalies as the difference from the expected ratio of sparse vegetation to tall woody vegetation pixels. The ratio was calculated as the fraction of sparse vegetation pixels, (sparse/ (sparse + tall woody)), and the expected ratio was based on the landforms present, calculated as a weighted average combining the number of pixels and the mean sparse fraction for each landform over the whole pertinent area.
To test whether reindeer herbivory had affected the sensitivity of vegetation productivity to climate as suggested by our second hypothesis, we analyzed variations in VSI with reindeer pressure using our three comparisons. The VSI map, at 1 km resolution, is much coarser than the landcover (20 m) and landform classifications (32 m). To enable analysis at a common spatial scale, each VSI pixel was assigned to a landcover class if that class covered at least half of the 1 km pixel and to a landform class if it covered at least 2/3 of the pixel. The higher threshold for landforms was chosen because of the fewer landform classes and to avoid potentially complex dynamics in more rugged terrain. Pixels that did not meet these thresholds were not used in the subsequent comparisons, but the classified pixels were well-distributed over the area with ample coverage ( Figure S3 in Supporting Information S1).
Comparisons of VSI in which only two samples were compared (i.e., between Yar-Sale/Panayevsk and Taz within a common landform or landcover) used the Wilcoxon rank sum test (wilcox.test from stats (R Core Team, 2020)) to determine significance. Comparisons involving more than two samples or multiple categorical variables used one-way ANOVAs or multi-way ANOVAs with type II sum of squares (anova test from rstatix (Kassambara, 2020)). When the ANOVAs found significant differences, we evaluated the differences between pairs of samples using Tukey's test (TukeyHSD from stats (R Core Team, 2020)). To quantify the relationship between VSI and reindeer density within pasture areas, we performed linear regressions between these two variables (lm from stats (R Core Team, 2020)).

Variation in Landcover With Reindeer Herbivory Pressure
In comparisons of the relative proportions of landcover types between the Yar-Sale/Panayevsk tundra, where reindeer numbers are highest on Yamal, and the Taz Peninsula, where there are considerably fewer reindeer (Figures 2a, 2c, and 2e), and between areas with and without summer pastures within Yar-Sale/Panayevsk (Figures 2b, 2d, and 2f), the areas with higher reindeer pressure consistently had higher coverage of sparse vegetation and lower coverage of tall woody vegetation than areas with fewer reindeer. This pattern occurred within all landform units, but in accordance with our hypothesis, the effect was greatest in the uplands, where reindeer spend the most time. Based on odds ratios, the ratio of tall woody vegetation area to sparse vegetation area was 4.0 times lower on the uplands, 3.6 times lower on the slopes, and 3.2 times lower on the lowlands of Yar-Sale/Panayevsk compared to Taz. Also in accordance with our expectations, the pastures grazed in June, when deciduous shrubs are most vulnerable to herbivory, generally had the most sparse vegetation area and least tall woody vegetation within Yar-Sale/Panayevsk. This result would be surprising in the absence of herbivory, as especially for the enterprise herds, these pastures tend to be located farther south in the low Arctic (Figure 1), where greater coverage of tall shrubs is expected (Myers-Smith et al., 2015). When comparing the June pastures to the no pastures area, the ratio  Bartsch et al. (2019). The distributions are compared by region (a, c, and e) and by the different pasture areas used in the summer months by herds in the Yar-Sale/Panayevsk tundra (b, d, and f). "No pastures" represents areas without any known reindeer herds. The pixels are also divided topographically into landforms: (a, b) uplands, (c, d) slopes, and (e, f) lowlands. The percentage of pixels within each landcover type is written inside the bar segments, and the total number of pixels in the specified area is given above the bar. Note that for each landform class, the sum of the number of pixels in the different pasture areas (b, d, and f) exceeds the number of pixels in the Yar-Sale/Panayevsk region (a, c, and e) because pixels may belong to multiple pasture areas (i.e., grazed in multiple months). of tall woody vegetation area to sparse vegetation area was 2.3 times lower on the uplands, 1.5 times lower on the slopes, and 1.2 times lower on the lowlands. Using a suite of log-linear models formulated to test the dependence of the proportions of landcover classes on categorical predictors, we found that both landforms and comparison areas (region in the first case, Yar-Sale/Panayevsk pasture areas in the second) significantly improved the model fit when their main effects and the interaction effects between them were included (see Table 1 for deviance statistics).
While there were also differences in the relative proportions of low stature shrub tundra and inundated/landslide area, they showed mixed responses to reindeer presence. These differences likely were not caused by variation in herbivory pressure, but potentially due to a higher density of rivers and lakes in Yar-Sale/Panayevsk (11.7% area covered by water, compared to 4.8% in Taz, calculated from the landcover classification), leading to greater coverage of inundated/landslide area than in Taz, and a tendency for reindeer herds to avoid landslide and mire areas (Skarin et al., 2020), resulting in lower coverage of this landcover in pasture areas within Yar-Sale/Panayevsk. We also found that higher reindeer pressure remained linked to less tall woody vegetation and more sparse vegetation when considering variations in animal density across the pasture areas. Within each herd area of a given reindeer density, we calculated how the observed proportions of sparse vegetation and tall woody vegetation differed from the expectation based on the average coverage proportions associated with the landforms present ( Figure 3). A positive value for this difference implies that the area contains more sparse vegetation and less tall woody vegetation than expected. The relationship between reindeer density and this difference was highly nonlinear, as the most pronounced effect occurred once the reindeer density exceeded a high threshold: 21 reindeer/km 2 on June pastures and 47.5 reindeer·months/km 2 when integrating over June-August. In all comparisons, points were mostly absent from the lower right of the plot area, representing high animal density and more tall woody vegetation than expected based on the landforms present. This tendency is highlighted by the contours, which display the results of a 2D kernel density estimation. Especially when only considering the uplands (Figures 3b  and 3d), the outermost contour lines slope upwards (from lower left to upper right in the plot), indicating a similar probability of pasture areas having low reindeer density with more tall woody vegetation than expected and high reindeer density with more sparse vegetation than expected.

Effect of Reindeer Herbivory on Vegetation Sensitivity to Climate
Contrary to our second hypothesis, VSI values were consistently greater in Yar-Sale/Panayevsk than in Taz, albeit by a small margin, despite the much greater herbivory pressure in Yar-Sale/Panayevsk (Figure 4). Paired comparisons of the VSI between regions using the Wilcoxon rank sum test indicated a significant difference for all comparisons within landforms and all comparisons within landcover, except within tall woody vegetation. Note. The suite of models evaluates possible explanations for the differences in the proportions of tall woody vegetation and sparse vegetation landcover pixels (see Figure 2). One set of models evaluates differences in landcover between Yar-Sale/ Panayevsk and Taz ("Regions"), while the other set evaluates differences between June pastures, July pastures, August pastures, and land not used for summer pastures in Yar-Sale/Panayevsk ("Pasture Areas"). For each model, df represents the degrees of freedom and G 2 gives the residual deviance comparing the fitted model against the saturated model with all possible effects including interactions.

Table 1 Deviance Statistics for Log-Linear Models
When dividing Yar-Sale/Panayevsk into pasture areas, the VSI in each area remained slightly higher than in Taz within every landcover type except, notably, within tall woody vegetation. Whereas the mean VSI in tall woody vegetation still was higher in the no pasture area of Yar-Sale/Panayevsk than in Taz (+2.59), all reindeer pasture areas had a lower mean VSI (June: −1.38, July: −1.64, and August: −0.533) (Figure 5a). A one-way ANOVA indicated significant differences in the mean VSI in the comparisons within tall woody vegetation (F = 4.728, p < 0.001), but Tukey's test showed a significant difference only in the comparison of no pasture area with July pastures (p = 0.002).
The impact of reindeer herbivory on the sensitivity of tall woody vegetation to climate is further supported by the negative correlation between VSI and reindeer density within this landcover (Figures 5b and 5c). The estimated slope of the linear regression was −0.708 VSI/(reindeer/km 2 ) (p = 0.0415) for June pastures and −0.268 VSI/ (reindeer·months/km 2 ) (p = 0.0155) for all pasture areas.

Discussion
Our results suggest that reindeer herbivory on Yamal constrains the response of tall woody vegetation to warming and limits shrub expansion, shaping the landscape to consist of less tall woody vegetation and more sparse vegetation than it would have without the influence of animals. Employing an associative pattern attribution approach, Figure 3. Link between reindeer density and the amount by which the ratio of sparse vegetation to tall woody vegetation differs from the ratio that would be expected based on the distribution of landforms in each herd area. Each point represents a continuous herd area with a given reindeer density in (a, b) June or density summed over (c, d) June, July, and August. The difference from the expected ratio was calculated as (sparse fraction − expected sparse fraction), where the sparse fraction is (sparse/(sparse + tall woody)) and expected sparse fraction is a weighted average combining the number of pixels and the mean sparse fraction in each relevant landform class over all June pasture areas (a, b) or all summer pasture areas (c, d) in Yar-Sale/Panayevsk. The number of sparse vegetation and tall woody vegetation pixels on all landforms is included in panels (a, c), while the sparse fractions are only calculated on upland pixels in panels (b, d). The contour lines and shading display the results of 2D kernel density estimation.
we evaluated six lines of evidence, which entailed using three spatial comparisons of herbivory pressure to test each of two hypotheses regarding the responses of landcover distributions and VSI to reindeer herbivory. The results of the landcover comparisons were consistent with our hypotheses-higher reindeer pressure was linked to lower coverage of tall woody vegetation, as well as more sparse vegetation, especially on uplands and early in the growing season. However, the higher vegetation sensitivity to climate observed in Yar-Sale/Panayevsk with respect to the Taz Peninsula was unexpected, was not explained by differences in landforms, general relief, landcover, or substrate, and occurred in both pastures and non-pastures (Figures 4 and 5, Figure S4 in Supporting Information S1). We speculate that this pattern might be driven by differences in growing season climate variability between both regions during the period 2000-2019 over which the VSI is computed. Yet, despite this general pattern in VSI, vegetation sensitivity within Yar-Sale/Panayevsk was consistently lower in the more heavily grazed areas exclusively in the tall woody vegetation. This finding suggests that reindeer herbivory does not have an extensive impact on vegetation sensitivity, but that it still has an observable effect in tall woody vegetation, where vegetation is otherwise most sensitive and where we would expect herbivory to counteract climate effects as drivers of vegetation change (Olofsson et al., 2009;Post & Pedersen, 2008). Therefore, after considering all the evidence simultaneously, top-down control from reindeer herbivory emerged as the most parsimonious, consistent explanation for the observed patterns. Our findings contrast with earlier regional-scale work concluding that bottom-up, climate-driven increases in unpalatable shrubs in North America contributed to declines in caribou abundance (Fauchald et al., 2017). However, our analysis benefitted from finer spatial grain-whereas Fauchald et al. (2017) used the GIMMS NDVI3g product at 1/12° resolution (∼8 km), our coarsest remote sensing product was at 1 km resolution, >60 times finer resolution-which enabled the pixel-based observations to better represent landscape-scale ecosystem structure and functioning. Additionally, the spatial comprehensiveness of the herd size and movement data provided a more accurate picture of the heterogeneity in herbivory pressure across the region, avoiding the implicit assumption often found within prior studies, including those in Yamal (e.g., Verma et al., 2020;Veselkin et al., 2021), that reindeer and changes in their population sizes are spread evenly over a region. The integration of high-resolution remote sensing products with spatially explicit reindeer population data provided the opportunity for multiple tests of our hypotheses, allowing us to carefully examine the competing effects of climate and herbivory on the vegetation on Yamal at a regional scale. This pattern-based methodology could be adapted for future studies testing whether, or the degree to which, large-scale ecosystem responses can be attributed to the influence of animals.
The magnitude of the vegetation response to herbivory detected in this study was significant but relatively modest, observed as small percentage differences in VSI and landcover composition between the reindeer areas and their controls, and primarily in the areas with the strongest grazing pressure. This result evokes prior work which found that reindeer limited tall shrub expansion in Northern Norway but only when their densities exceeded approximately 5 animals/km 2 (Bråthen et al., 2017). The threshold we found (i.e., 21 reindeer/km 2 on June pastures, Figure 3) was higher than that observed in Northern Norway, but this disparity may be related to differences between the Nenets' reindeer herding in Yamal and the practices in Fennoscandia. Although Nenets herders are generally guided by the natural ecology of reindeer, at the regional scale, the seasonal choice of pastures and the timing of movements are determined entirely by herders. In the summer, Nenets on the Yar-Sale/Panayevsk tundra migrate every 3-5 days on average and even while a household remains in one place between movements, the herders try to direct reindeer to new areas as time passes. Thus, we hypothesize that a lower density of Fennoscandian free-grazing reindeer without constant human control could exert greater pressure on vegetation, especially the preferred shrub pastures in June, by remaining longer in certain areas. Within the Yamal Peninsula, the greater vegetation response we detected on June pastures may be not only due to the vulnerability of the vegetation to herbivory in the early growing season, but also because many households pass through the same areas as they migrate before dispersing to their July and August pastures. As a result, it would primarily be the June pastures that exceed the threshold of reindeer density needed to cause significant impacts on the vegetation, despite the characteristic high rates of movement.  Figure 4 for the explanation of the boxes. Only significant differences based on Tukey's tests are annotated: a single asterisk (*) indicates p < 0.05; double asterisks (**) indicate p < 0.01; and triple asterisks (***) indicate p < 0.001. Each point in panels (b and c) represents a 1 km VSI pixel within tall woody vegetation with a given June reindeer density (b) or density summed over June, July, and August (c). The linear least squares fits have a significant negative slope ((b)−0.708 VSI/(reindeer/km 2 ), (c) −0.268 VSI/(reindeer·months/km 2 ), p < 0.05) and are displayed with shaded 95% confidence intervals.
Furthermore, detecting the impacts of herbivory in an uncontrolled setting is complicated by challenges such as rapid vegetation recovery blurring the signal of the herbivory impact (Ravolainen et al., 2011) and uncertainty in the precise positioning and timing of animal movements. Thus, the true ecological magnitude of the herbivory impacts on Yamal may exceed our findings, and given the 67,000 km 2 regional expanse over which we made our observations, there likely are implications for carbon storage and albedo relevant to Earth system functioning (Meredith et al., 2019).
Our observations represent a snapshot and do not inform whether the impact of reindeer on Yamal is changing over time nor the extent to which it may counteract climate-driven shrub growth in the future. Repeated monitoring of vegetation is needed to track and anticipate how the region will respond as the Arctic continues to warm. Our data on summer locations of reindeer herds and their migratory corridors reflect at least the last 30 years, but reindeer herding practices in the area are also changing, as herders adapt to infrastructure development for hydrocarbon extraction ) and climate change, which affects the timing of the freezing and thawing of rivers and has led to an increase in rain-on-snow events (Forbes et al., 2016(Forbes et al., , 2022. The latter cause the severe icing of pastures, leading to high reindeer mortality. There also have been calls from both biologists and administrators to reduce the number of reindeer on Yamal and restrict their migrations in response to a recurring argument that the land is overgrazed, based primarily on a decline in lichen availability (Golovnev, 2017;Volkovitskiy & Terekhina, 2020). While the definition and perception of overgrazing is a complex, interdisciplinary topic, our findings can inform these discussions by providing empirical insights into the impact reindeer are having on the vegetation over a regional scale.
The potential for rewilding the Arctic with large herbivores to reduce the abundance of shrubs and slow climate change feedbacks has already drawn considerable attention, but questions remain as to how effectively such projects will provide these ecosystem services Olofsson & Post, 2018). Some of this discussion arises from uncertainty over whether Pleistocene megafauna had a large top-down role in maintaining the graminoid-and forb-rich vegetation, or whether they were a more passive element within a trophic chain in which vegetation cover was determined by climate as a bottom-up driver (e.g., Monteath et al., 2021). Our results support that herbivores, at current population densities, affect vegetation cover in the modern Arctic at a regional and potentially biome-wide scale. The domesticated reindeer herds of the Nenets are the largest in the world today, whereas estimates of past megafauna population sizes are highly uncertain (Monteath et al., 2021). The density of large herbivores across the Arctic historically may have been similar or likely much higher than the current reindeer density on Yamal (Zimov et al., 2012), and the modern perception of this density as high may be a case of "shifting base-line syndrome" (Pauly, 1995) occurring over long periods of time. Thus, especially given the larger and more diverse Pleistocene herbivore assemblage (Olofsson & Post, 2018), it is likely that herbivory in the Arctic during the Pleistocene had even stronger impacts than what we have found here. Further work is necessary to translate the fingerprints of reindeer herbivory we detected into ecosystem variables and their quantitative estimates, but our findings demonstrate that top-down regulation by herbivory must be considered to understand how tundra ecosystems have responded, are responding, and will respond to climate change.

Data Availability Statement
The data used in this study are available in the Oxford University Research Archive via https://doi.org/10.5287/ bodleian:E2x6bnggk under a Creative Commons license (CC BY 4.0) (Spiegel et al., 2022b). V1.0.0 of the code (Spiegel et al., 2022a) used in the pre-processing of the data, analysis, and figure generation is preserved at https://doi.org/10.5281/zenodo.7387883, available under a GNU General Public License v3.0 and developed openly using R version 4.0.0 (R Core Team, 2020).