Satellite-based habitat monitoring reveals long-term dynamics of deer habitat in response to forest disturbances

. Disturbances play a key role in driving forest ecosystem dynamics, but how dis- turbances shape wildlife habitat across space and time often remains unclear. A major reason for this is a lack of information about changes in habitat suitability across large areas and longer time periods. Here, we use a novel approach based on Landsat satellite image time series to map seasonal habitat suitability annually from 1986 to 2017. Our approach involves charac- terizing forest disturbance dynamics using Landsat-based metrics, harmonizing these metrics through a temporal segmentation algorithm, and then using them together with GPS telemetry data in habitat models. We apply this framework to assess how natural forest disturbances and post-disturbance salvage logging affect habitat suitability for two ungulates, roe deer ( Capreo- lus capreolus ) and red deer ( Cervus elaphus ), over 32 yr in a Central European forest landscape. We found that red and roe deer differed in their response to forest disturbances. Habitat suit- ability for red deer consistently improved after disturbances, whereas the suitability of disturbed sites was more variable for roe deer depending on season (lower during winter than summer) and disturbance agent (lower in windthrow vs. bark-beetle-affected stands). Salvage logging altered the suitability of bark beetle-affected stands for deer, having negative effects on red deer and mixed effects on roe deer, but generally did not have clear effects on habitat suitability in windthrows. Our results highlight long-lasting legacy effects of forest disturbances on deer habitat. For example, bark beetle disturbances improved red deer habitat suitability for at least 25 yr. The duration of disturbance impacts generally increased with elevation. Method- ologically, our approach proved effective for improving the robustness of habitat reconstruc-tions from Landsat time series: integrating multiyear telemetry data into single, multi-temporal habitat models improved model transferability in time. Likewise, temporally segmenting the Landsat-based metrics increased the temporal consistency of our habitat suitability maps. As the frequency of natural forest disturbances is increasing across the globe, their impacts on wildlife habitat should be considered in wildlife and forest management. Our approach offers a widely applicable method for monitoring habitat suitability changes caused by landscape dynamics such as forest disturbance.

Abstract. Disturbances play a key role in driving forest ecosystem dynamics, but how disturbances shape wildlife habitat across space and time often remains unclear. A major reason for this is a lack of information about changes in habitat suitability across large areas and longer time periods. Here, we use a novel approach based on Landsat satellite image time series to map seasonal habitat suitability annually from 1986 to 2017. Our approach involves characterizing forest disturbance dynamics using Landsat-based metrics, harmonizing these metrics through a temporal segmentation algorithm, and then using them together with GPS telemetry data in habitat models. We apply this framework to assess how natural forest disturbances and post-disturbance salvage logging affect habitat suitability for two ungulates, roe deer (Capreolus capreolus) and red deer (Cervus elaphus), over 32 yr in a Central European forest landscape. We found that red and roe deer differed in their response to forest disturbances. Habitat suitability for red deer consistently improved after disturbances, whereas the suitability of disturbed sites was more variable for roe deer depending on season (lower during winter than summer) and disturbance agent (lower in windthrow vs. bark-beetle-affected stands). Salvage logging altered the suitability of bark beetle-affected stands for deer, having negative effects on red deer and mixed effects on roe deer, but generally did not have clear effects on habitat suitability in windthrows. Our results highlight long-lasting legacy effects of forest disturbances on deer habitat. For example, bark beetle disturbances improved red deer habitat suitability for at least 25 yr. The duration of disturbance impacts generally increased with elevation. Methodologically, our approach proved effective for improving the robustness of habitat reconstructions from Landsat time series: integrating multiyear telemetry data into single, multi-temporal habitat models improved model transferability in time. Likewise, temporally segmenting the Landsat-based metrics increased the temporal consistency of our habitat suitability maps. As the frequency of natural forest disturbances is increasing across the globe, their impacts on wildlife habitat should be considered in wildlife and forest management. Our approach offers a widely applicable method for monitoring habitat suitability changes caused by landscape dynamics such as forest disturbance.

INTRODUCTION
Forest disturbances significantly alter the amount and quality of habitat available to wildlife (Thom and Seidl 2016) and, through their long-lasting effects on forest structure, shape habitat conditions over long time periods (Schieck and Song 2006). Yet, information on how forest disturbances modify habitat suitability for wildlife, particularly across large areas and long time periods, is not commonly available to wildlife management and conservation (Clare et al. 2019). At the same time, generalized predictions about wildlife responses to forest disturbances are difficult, as different disturbance types vary in their impacts on forests. For example, natural forest disturbances (e.g., insect outbreaks and windthrows) create biological legacies promoting structural complexity (Johnstone et al. 2016), while forest management interventions such as clear-cutting or salvage logging homogenize forest structure (Thorn et al. 2017). As natural forest disturbances are becoming more frequent across large parts of globe (Seidl et al. 2017, Senf et al. 2018), a better understanding of forest disturbance impacts on wildlife habitat is important to develop adequate management and conservation strategies.
Habitat suitability models allow the assessment of wildlife responses to forest disturbances (Berland et al. 2008) and provide spatially explicit habitat suitability maps that are critical for informing wildlife management and conservation (Guisan and Thuiller 2005). Typically, however, habitat models only provide a snapshot of habitat for single points in time (Franklin 2010), limiting the usefulness of such models for monitoring the effects of landscape dynamics on species' habitat. Integrating satellite time series into habitat models has great potential for overcoming this limitation (Randin et al. 2020). Satellite images offer detailed characterizations of wildlife habitat (Bellis et al. 2008, Lahoz-Monfort et al. 2010, and long-term satellite records, such as the Landsat archive, provide consistent information on landscape change, including forest disturbances . Linking habitat models to Landsat time series thus could allow a wall-to-wall mapping of habitat suitability several decades back in time. Yet, whereas Landsat time series have been widely and successfully used for monitoring forest disturbances (Banskota et al. 2014), their application for monitoring disturbance-related dynamics in wildlife habitats has been scarce (but see Kearney et al. [2019] for a recent example).
Reconstructing habitat dynamics based on Landsat time series, however, entails challenges particularly with regards to ensuring the consistency of habitat maps over time. As available wildlife records often only cover parts of multi-decadal Landsat time series, model transfers in time are typically necessary, which pose challenges for habitat models (Tuanmu et al. 2011, Yates et al. 2018. At the same time, sampling schemes for collecting wildlife data often vary over time. To tackle these issues, species records collected over time can be integrated into single habitat models, where species records are matched with time-varying predictors (Nogu es-Bravo 2009). Such a multi-temporal calibration approach of habitat modeling (hereafter: multi-temporal habitat models) makes model predictions independent of year-to-year variations in sampling bias and can improve model transferability (Maiorano et al. 2013). While satellite time series generally facilitate building multi-temporal habitat models (Randin et al. 2020), time series of satellite-based metric also often exhibit year-to-year fluctuations unrelated to actual changes in habitat, caused for example by variations in satellite image availability or cloud cover. Thus, to avoid mapping spurious habitat change, satellite time series must be harmonized across time. In this context, satellite time series algorithms have been developed to separate change signals from noise (Zhu 2017), but their applicability for monitoring habitat dynamics remains untested.
Multi-decadal time series of habitat suitability maps also have considerable potential to provide new insights into how forest disturbance affect wildlife habitat, particularly regarding long-term effects (Kearney et al. 2019). Understanding such impacts of forest disturbances on large herbivore species, such as ungulates, is an important concern for forest and wildlife management (Vospernik and Reimoser 2008). Forest disturbances are a key driver of ungulate habitat dynamics (Fisher and Wilkinson 2005) and, in turn, ungulate herbivory has important effects on forest ecosystem dynamics and functions by altering forest structure and successional trajectories (Côt e et al. 2004). Forest openings created by disturbance typically provide good foraging conditions to ungulates due to a high abundance of plant biomass on the ground (Kuijper et al. 2009). On the other hand, forest openings offer low levels of cover, which increases the exposure to thermal extremes and predation risk (Mysterud and Østbye 1999). The studies available indicate variation of ungulate responses to forest disturbance across species and study systems (see, e.g., Moser et al. 2008, Kuijper et al. 2009, Lamont et al. 2019. However, these studies are overwhelmingly limited to single types of disturbance, as well as initiation stages immediately after disturbance events. How ungulate habitat changes as disturbed areas recover, and whether habitat suitability varies between different natural disturbances and post-disturbance management interventions remain, therefore, largely unresolved questions. Here, we use a novel approach linking habitat models with Landsat satellite time series for monitoring wildlife habitat dynamics consistently and continuously over multiple decades. Specifically, we assess how habitat suitability for red deer (Cervus elaphus) and roe deer (Capreolus capreolus) changed in response to forest disturbance in the Bohemian Forest Ecosystem, a Central European montane forest landscape. We characterized forest disturbance and recovery dynamics using time series of Landsat-based spectral-temporal metrics, which we harmonized across time using a temporal segmentation algorithm (Kennedy et al. 2010). Then, we integrated these Landsat-based metrics with GPS telemetry data collected over 13 yr (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) in multi-temporal habitat models. By projecting habitat models across the full length of our satellite time series, we mapped changes in summer and winter habitat suitability annually for all years between 1986 and 2017. Finally, we used this time series of habitat suitability maps to investigate how deer respond to forest disturbance caused by bark beetle outbreaks, windthrow, and salvage logging. Specifically, we sought to answer the following research questions: (1) Does our approach of multi-temporal data integration and harmonization improve the transferability and consistency of habitat models across time? (2) Do red and roe deer respond differently to forest disturbances by bark beetle, windthrow, and salvage logging? (3) How does habitat suitability evolve as disturbed sites recover?

Study area
We conducted our study in the Bohemian Forest Ecosystem situated along the border of Austria, Czechia, and Germany (Fig. 1). Two protected areas form the center of the study site: the Bavarian Forest National Park (240 km 2 ) and the Sumava National Park (680 km 2 ). Elevations range from~400 to 1,450 m. The area's forests are typical for Central European mountainous forests, dominated by Norway spruce (Picea abies) and accompanied by mountain ash (Sorbus aucuparia) at higher elevations (>1,100 m). Lower elevations are characterized by mixed montane forests composed mainly of European beech (Fagus sylvatica), Norway spruce, and silver fir (Abies alba). Since the 1990s, largescale forest disturbances caused by bark beetle outbreaks and windthrow have restructured large parts of the Bohemian Forest (Oeser et al. 2017). Inside the core zones of the national parks, disturbed sites were left to recover naturally, but salvage logging was applied in the management zones to prevent an uncontrolled spread of bark beetle (Thorn et al. 2017).

GPS telemetry data for roe deer and red deer
We used GPS telemetry records containing a total of 111 red deer and 164 roe deer individuals. Red deer data were collected between 2002 and 2014, and roe deer data between 2004 and 2013. As some of the GPS-collared red deer were kept in winter enclosures (Heurich et al. 2015), we removed all red deer locations falling into these enclosures and only used locations from red deer outside the enclosures for further analysis. We standardized the sampling frequency among individuals to four observations per day, retaining one observation each from four temporal windows: nighttime (22:00-04:00), dawn (04:00-10:00), daytime (10:00-16:00), and dusk (16:00-22:00). We further split observations into summer (1 April-31 October) and winter (1 November-31 March) locations and created a random subset of 100 observations per individual and season. This allowed using data from most deer individuals, while ensuring that individuals monitored over longer time periods than others were not overrepresented in our habitat models. We excluded individuals with fewer than 100 observations per season (after standardizing the sampling frequency). The final data sets used for building habitat models consisted of 10,999 summer and 7,997 winter observations from 110 and 80 red deer individuals, respectively, and 14,723 summer and 14,925 winter observations from 148 and 150 roe deer individuals, respectively.

Characterizing habitat change with time series of Landsat-based metrics
We used Landsat satellite image time series to characterize changes in deer habitat caused by forest disturbance and post-disturbance recovery. Our approach consisted of three main steps (Fig. 2).
FIG. 1. Map of the study area. Forest disturbances between 1986 and 2014 were mapped based on Landsat time series (Oeser et al. 2017). Polygons for red and roe deer show minimum convex polygons around GPS telemetry observations. In our first main step, we derived spectral-temporal metrics from Landsat satellite images. For each Landsat image, we derived the tasseled cap indices brightness, greenness, and wetness (Crist and Cicone 1984), which are sensitive to forest types, structure, and disturbance (Hansen et al. 2001, Dymond et al. 2002, Healey et al. 2005. Calculating the tasseled cap indices for all available images from a given year results in distributions of index values for each image pixel, since spectral indices vary along characteristic seasonal profiles, which are indicative of land-cover and vegetation characteristics (Pasquarella et al. 2016). To summarize the intra-annual phenological variation of the three tasseled cap indices, we then calculated the 10th, 50th (i.e., the median), and 90th percentiles at the pixel level, resulting in a total of nine spectral-temporal metrics. In an earlier study, we found that spectral-temporal metrics based on the tasseled cap indices allow for detailed characterizations of wildlife habitat, including deer habitat ). Here, we extend our approach of linking habitat models with Landsat imagery to map long-term habitat dynamics. Therefore, we created annual time series of the Landsat-based metrics using all images recorded by the sensors 4 TM, 5 TM, 7 ETM+, and 8 OLI between 1985 and 2018 covering our study area (1,455 images with 30-m spatial resolution; see Appendix S1 for a detailed description of the image processing). Due to high levels of cloud cover in our study area, we did not calculate metrics for each year separately, but instead used three-year moving windows, including all images from the year before and after for calculating metrics for a given year. This resulted in annual time series for all nine metrics from 1986 to 2017 (hereafter: raw metrics). While the use of three-year windows effectively reduced the temporal resolution of our time series, it ensured the robust calculation of metrics, as year-to-year variations in the number and seasonal timing of observations otherwise would have strongly affected the calculation of metrics and hence the consistency of our habitat monitoring.
To further improve the consistency of the Landsatbased metrics, we removed ephemeral year-to-year variations in the raw metric time series. To this end, we applied the satellite time series segmentation algorithm LandTrendr, which has been developed for detecting and characterizing land surface changes with Landsat time series (Kennedy et al. 2010). LandTrendr partitions time series of spectral indices by first identifying breakpoints and subsequently fitting straight-line segments between these breakpoints (see Kennedy et al. 2010 for details). This results in linear trajectories of stable, decreasing or increasing values, which strongly reduces noise while capturing both abrupt and gradual changes (Vogeler et al. 2018). We derived trajectories from LandTrendr for all nine metric time series (hereafter: smoothed metrics). We performed all satellite image processing in the FIG. 2. Schematic overview of our workflow. (A) We derived Landsat-based metrics that we harmonized across time using time series segmentation. (B) GPS telemetry records collected over 13 yr were temporally matched with the smoothed metrics to train multi-temporal habitat models. These models were then used to predict habitat suitability across all years . (C) We subsequently linked the annual habitat suitability maps to independent forest disturbance maps to assess habitat suitability changes in response to forest disturbances.
Google Earth Engine (Gorelick et al. 2017, Kennedy et al. 2018) and provide the code via a link in Appendix S1.

Multi-temporal habitat modeling
In our second main step, we built multi-temporal habitat models for red and roe deer. For calibrating models, we used the subsampled telemetry data sets containing observations from 2002-2014 for red deer and 2003-2012 for roe deer. We temporally matched the telemetry observations with the Landsat-based metrics (Fig. 2) and built four habitat models for each deer species, generating separate models using the raw metrics (without temporal segmentation) and smoothed metrics (with temporal segmentation), as well as for summer (1 April-31 October) and winter (1 November-31 March). The latter was important since habitat use of red and roe deer differs between seasons in our study area (Heurich et al. 2015). We built all habitat models using the Max-Ent algorithm (Phillips et al. 2017). As background points (pseudoabsences), we sampled five random points in our study area for every presence observation and assigned the observation's year to these background points. To limit model complexity and to improve transferability, we fitted MaxEnt models using only hinge features ). Finally, we projected each habitat model across the entire time series of Landsatbased metrics to create habitat suitability maps for each year between 1986 and 2017. For predicting models, we used the cloglog output of MaxEnt (Phillips et al. 2017), which we use as an index of relative habitat suitability.
To compare the temporal consistency of habitat suitability maps based on the smoothed and raw metrics, we analyzed the year-to-year variability of habitat suitability in undisturbed forests. In undisturbed forests, habitat suitability should vary little between consecutive years, thus serving as a benchmark for the temporal stability of our habitat models. As a measure of year-to-year variability, we calculated the mean absolute difference in habitat suitability between consecutive years at the pixel level. To evaluate whether variability in the habitat suitability maps was influenced by satellite image availability, we correlated the variability in habitat suitability with the average number of clear Landsat observations at a given location.
To test whether the integration of deer observations from multiple years improved the temporal transferability of habitat models, we built seasonal habitat models using deer observations from each individual year, as well as all possible windows of three and five consecutive years. Then, to validate these models, we created habitat suitability predictions for each of the remaining annual data sets (i.e., presence and background points from deer individuals not used for model training). As our data set only allowed transferring models built from five years of deer data up to eight and five years in time for red and roe deer, respectively, we limited model transfers of all models to this time range to ensure comparability. We assessed model transferability by calculating two performance metrics from the model predictions: (1) the area under the Receiver Operating Characteristic curve (AUC) and (2) the Continuous Boyce Index (CBI; Hirzel et al. 2006). These metrics provide complementary measures for two different aspects of habitat model performance. While the AUC measures the discriminatory ability of models (i.e., their ability to discriminate presence from background points), the CBI can be used as an indicator of model calibration (i.e., how well the predicted habitat suitability values match with the frequency of presence locations; see Phillips and Elith [2010] for more details on both validation metrics).
Assessing habitat suitability changes in response to forest disturbance We linked our time series of habitat suitability maps to Landsat-based forest disturbance maps from our study area (Oeser et al. 2017). These disturbance maps capture forest disturbances occurring between 1986 and 2014 at the level of Landsat pixels (30-m resolution), and further differentiate four disturbance types: bark beetle infestations, windthrows, salvage-logged windthrows and other logging, which includes salvage logging of bark beetleinfested stands inside the national parks. Based on our time series of habitat suitability maps, we built chronosequences of habitat suitability for salvaged and unmanaged bark beetle and windthrow disturbances by grouping values based on the number of years since a pixel had been disturbed. Because bark beetle outbreaks overwhelmingly occurred after 1992 and windthrow disturbances after 2006 at our study site, we capped chronosequences at 25 yr for bark beetle disturbances, and at 11 yr for windthrows.
The disturbance types tend to occur at different elevations in our study area (e.g., windthrows predominantly occurring at higher elevations; see Appendix S2: Fig. S1). As elevation has strong effects on the distribution of deer in our study area (Heurich et al. 2015) and further affects post-disturbance forest recovery ), we thus separately compared salvaged and non-intervention sites at low (600-900 m), medium (900-1,100 m), and high elevations (1,200-1,400 m) and further summarized habitat suitability changes for all bark beetle and windthrow disturbances across elevation-bins of 50 m.

Consistency and transferability of habitat models
Using the smoothed metrics obtained from Land-Trendr improved the temporal consistency of habitat suitability maps, as the year-to-year variability of habitat suitability in undisturbed forests (i.e., stable habitat) was lower compared to habitat suitability maps based on the raw metrics (average reduction of 52% across all habitat models; Fig. 3A). The temporal consistency of habitat suitability time series generally decreased with fewer available Landsat observations at a given location, but this negative effect was dampened when using the smoothed metrics ( Fig. 3A; correlation coefficient for smoothed and raw metrics: À0.72 and À0.47, respectively).
The transferability of habitat models in time improved on average through the integration of deer data from multiple years for model calibration ( Fig. 3B; mean AUC and CBI for 5-yr vs. 1-yr models: 0.74 and 0.91 vs. 0.69 and 0.81). In addition, models based on the smoothed metrics showed better average performance during model transfers than those based on the raw metrics ( Fig. 3B; mean AUC and CBI for five-year models based on the smoothed vs. raw metrics: 0.73 and 0.91 vs. 0.70 and 0.84).

Habitat suitability changes after forest disturbances
Our habitat suitability maps revealed pronounced changes in deer habitat following forest disturbances (Fig. 4). Overall, the widespread forest disturbances in our study area led to improved habitat conditions for both deer species over our monitoring period (see Appendix S3: Fig. S1 for a summary of habitat suitability changes across all forest areas). However, deer responses to forest disturbances varied between deer species, seasons, and disturbance types (Fig. 5). While red deer benefitted from forest disturbances throughout the year, roe deer showed a seasonally changing response to disturbances. Habitat suitability for red deer initially increased steeply after disturbances but started to gradually decrease after a few years (peak of median winter and summer habitat suitability across all disturbance types after 4 yr). For roe deer, disturbances also led to increases in summer habitat suitability (peaking after 5 yr) but caused decreases in winter habitat suitability before recovering again after a few years (reaching a minimum after 3 yr).
Red and roe deer responded differently to bark beetle and windthrow disturbances. For the first few years after disturbance, habitat suitability for red deer was similar at bark beetle and windthrow disturbances (Fig. 5). Conversely, for roe deer, habitat suitability was higher after bark beetle than after windthrow disturbances, which caused substantial decreases in suitability for this species during winter (Fig. 5). The habitat suitability of disturbed sites also varied with the post-disturbance management strategy (i.e., salvage logging vs. no intervention). Salvage logging, however, had much clearer effects on deer habitat suitability in bark beetle-affected areas than in windthrows. Salvage-logged bark beetle disturbances were generally less suitable for red deer and less suitable for roe deer during summer when they were recently disturbed (i.e., <10 yr). Yet, seasonal habitat suitability for roe deer also improved due to salvage logging: During winter, recently disturbed sites were more suitable for roe deer after salvage logging than after no intervention. With advancing regeneration over time, habitat suitability tended to return to pre-disturbance conditions. Changes in habitat suitability, however, were generally long-lasting. For instance, bark beetle disturbances altered red deer habitat suitability for a minimum of around 25 yr (Fig. 5). While we could assess habitat suitability changes after windthrows for only 11 yr postdisturbance (due to windthrows occurring only after 2006 in our study area), our habitat suitability time series indicate that the duration of disturbance impacts on deer habitat differ between windthrow and bark beetle disturbances. For example, after its initial increase, red deer habitat suitability declined faster at windthrows than at bark beetle disturbances (Fig. 5). Regarding differences in post-disturbance management, salvagelogged and non-intervention sites showed overall similar temporal patterns, but habitat suitability generally returned more quickly to pre-disturbance conditions at salvage-logged sites. Finally, we also observed clear effects of the elevation gradient in our study area on disturbance-related habitat suitability changes (Fig. 6). Habitat suitability for both deer species remained changed for longer periods with increasing elevation. Disturbed sites were most suitable for red deer in the spruce-dominated forest zone (>1,100 m), but at lower elevations for roe deer.

DISCUSSION
Wildlife responses to forest disturbance often remain poorly understood and spatially explicit information on disturbance-related habitat suitability dynamics is rarely available to wildlife management and conservation. Here we provide new insights into how natural disturbances and salvage logging affect habitat suitability for two of the most widespread European ungulates, and highlight the great, but currently underused potential of Landsat time series for monitoring forest disturbance impacts on wildlife habitat across long time periods. Specifically, our study yielded three main findings: First, red and roe deer showed distinct, and partly diverging, responses to forest disturbances, indicating that disturbance impacts on ungulate habitat vary between species with different foraging strategies. Second, disturbance impacts on deer habitat were generally long-lasting, but legacy effects varied along the elevation gradient in our study area. This highlights that variation in site conditions, and hence in forest recovery, can cause considerable heterogeneity in ungulate habitat quality after disturbances. Finally, the integration of species records in multi-temporal habitat models, as well as the temporal segmentation of the Landsat-based metrics, improved the transferability and stability of habitat models across time. This underlines the potential of our approach for consistent monitoring of habitat dynamics directly from satellite imagery.

Different responses by red and roe deer to forest disturbances
The seasonally changing response to forest disturbances by roe deer, with lower habitat suitability of forest openings during winter, can likely be explained by two factors: First, during spring and summer roe deer in the Bohemian Forest feed heavily on forbs and grasses, which are widespread at disturbed sites during the vegetation season, but increase forage intake through browsing during winter (Baran cekov a et al. 2010). Second, disturbed sites provide low levels of cover for roe deer during winter. Roe deer are sensitive to thermal exposure and snow accumulation, which inhibits forage accessibility and movement (Mysterud and Østbye 1999). To escape harsh winter conditions, roe deer in the Bohemian Forest migrate to lower elevations and use more densely vegetated habitats (Cagnacci et al. 2011, Ewald et al. 2014b). Our results therefore underscore the importance of cover as a habitat element for roe deer during winter, and corroborate findings that red deer use forest openings more intensively than roe deer (Kuijper et al. 2009).
Roe deer also showed more variable responses with regards to the disturbance type. Whereas habitat suitability for red deer increased similarly after windthrow and bark-beetle disturbances, windthrows were less suitable for roe deer, even decreasing habitat suitability during winter. A likely contributing factor are differences in FIG. 6. Variation of habitat suitability changes along the elevation gradient in our study area, combining salvage-logged and non-intervention sites. Habitat suitability changes are summarized by calculating the average habitat suitability across elevation bins of 50 m for each year on the x-axis. deadwood structure, as windthrow gaps are dominated by lying deadwood, while bark beetle-affected stands initially remain as standing deadwood. Fallen logs impede both the movement and habitat visibility for deer, which can increase predation risk (Kuijper et al. 2013). Thus, roe deer might avoid windthrows to reduce predation risk by lynx, an ambush predator and the main predator of roe deer in the Bohemian Forest (Heurich et al. 2015). In addition, differences in foraging conditions are likely important: Windthrows tend to have lower tree regeneration densities than bark beetle-disturbed sites (Jon a sov a et al. 2010), and full-light conditions in windthrow gaps can lead to a reduction of forage quality (Moser et al. 2008). Our results therefore indicate that concentrate selectors, such as roe deer, show temporally and spatially variable responses to forest disturbances depending on the availability of high-quality forage, while less selective feeders, such as red deer, consistently benefit from the increase in forage quantities following disturbance.
Even though salvage logging of naturally disturbed forests is widely applied, our understanding of how it affects ungulate habitat quality remains largely incomplete (Hebblewhite et al. 2009). We found that salvage logging of bark beetle-affected stands reduced habitat suitability for red deer. This could be due to salvage logging removing ground vegetation that has survived under dead canopy, and thus reducing forage quantities (Jon a sov a and Prach 2008). On the other hand, recently disturbed sites that were salvage-logged were more suitable as winter habitat for roe deer, potentially because the removal of deadwood allows for easier movements and lowers predation risk by lynx (Kuijper et al. 2013). In contrast to this finding, Hebblewhite et al. (2009) found that elk (Cervus canadensis) avoid postfire salvage-logged sites because of an increased predation risk by wolves. The use of clearcut areas by ungulates could therefore also depend on the hunting strategy of their main predators (i.e., stalking vs. cursorial predators; Schmidt and Kuijper 2015).

Evolution of deer habitat with forest recovery
Our approach using satellite time series allowed tracking disturbance impacts on deer habitat for up to 25 yr, thus offering new insights into long-term effects of disturbances on ungulate habitat. Our results show that variations in site conditions, in our case along the elevation gradient of our study site, play an important role in modifying disturbance impacts on ungulate habitat. Disturbed sites in high-elevation, spruce-dominated forests in our study area regenerate more slowly, offer higher forage quantities (e.g., grasses and ferns), but lower forage quality (Ewald et al. 2014a). Forest disturbances thereby amplified the variation of habitat quality along the elevation gradient for both species, thus likely contributing to higher concentrations of red and roe deer at high and low elevations, respectively (Heurich et al. 2015).
The length of disturbance impacts on deer habitat was also affected by salvage logging, which caused habitat suitability to return more quickly to pre-disturbance conditions. Forest recovery progresses faster and more homogeneously at salvage-logged than at non-intervention sites in our study area . This demonstrates that the removal of biological legacies through salvage logging can negatively affect ungulates by shortening the access to highly attractive foraging habitat. Overall, the temporal patterns of habitat suitability we observed here match well with expected changes in forage availability with stand development (Smolko et al. 2018). This indicates that forage availability likely remains the most important factor determining habitat quality for deer in forest openings even as disturbed sites are recovering (Dupke et al. 2016).

Potential and limitations of our approach
Methodologically, we demonstrated that linking habitat suitability models to time series of Landsat-based metrics offers an effective approach for monitoring wildlife habitat dynamics continuously across multiple decades. Despite the potential of Landsat imagery for characterizing wildlife habitat (Bellis et al. 2008, Lahoz-Monfort et al. 2010), very few studies have leveraged the unique temporal depth of the Landsat archive for mapping long-term habitat dynamics (Romero-Muñoz et al. 2018, Clare et al. 2019, Kearney et al. 2019. Information on habitat dynamics is key for informing wildlife management and conservation, but creating consistent time series of habitat suitability maps from satellite-based variables poses challenges due to high noise levels in satellite time series (Tuanmu et al. 2011). Here, we showed that combining all available Landsat imagery to produce time series of spectral-temporal metrics and temporally segmenting these time series with the LandTrendr algorithm allows for robust characterizations of habitat dynamics directly from Landsat imagery. Habitat suitability maps based on smoothed Landsat metrics were more consistent in time, showing less year-to-year variation in undisturbed forests, and were less affected by low levels of satellite data availability. Similarly, habitat models based on the smoothed metrics were more transferable in time. Because our approach relies solely on freely and globally available Landsat imagery, it is potentially widely applicable for informing wildlife management and conservation (Lahoz-Monfort et al. 2010).
We also found that integrating deer observations from multiple years in multi-temporal habitat models improved model transferability in time. In a recent review, Randin et al. (2020) emphasized the potential of multi-temporal habitat models calibrated using satellite time series. Yet, thus far, only few studies (e.g., Sieber et al. 2015) have made use of this potential. Multi-temporal calibration allows using more observations for model training, and further captures species' responses to changing habitat conditions over time (e.g., by characterizing deer habitat use at different stages of forest regeneration in our case; Nogu es-Bravo 2009). Our model validation shows that this improves the transferability of Landsat-based habitat models, as previously shown for species distribution models calibrated with species data from different climate periods (Maiorano et al. 2013). As Landsat imagery is available retrospectively up until the 1970s, it offers unique possibilities to integrate species records for reconstructing habitat dynamics.
While our approach provides a cost-effective and transferable method for monitoring wildlife habitat dynamics, some limitations need to be mentioned. First, although integrating species records in multi-temporal habitat models offers important advantages, it also leads to species responses being averaged over time. Animals, however, can adjust habitat selection when the availability of resources is varying across space or time (through so-called functional responses; Mysterud and Ims 1998). In our case, deer might change their selection for forest openings as more forest area was disturbed in our study area over time.
Methods that allow incorporating functional responses in habitat models are an active field of research (Paton and Matthiopoulos 2016), and could help improve the reliability of long-term habitat monitoring in the future. Second, by only using spectral information at the pixel level, our approach cannot capture effects related to the spatial configuration or size of disturbed patches, which can affect habitat suitability for deer (Mass e and Côt e 2012). Finally, when relying on spectral-temporal metrics alone, additional sources of information are often needed to better understand the drivers of mapped habitat changes (such as the forest disturbance maps from a previous study we used here). Patterns in our time series of habitat suitability maps generally matched very well with the independent forest disturbance maps, but sometimes indicated habitat suitability changes earlier than the disturbance date assigned in the disturbance maps. These temporal mismatches between both data sets are likely caused by the three year-moving windows we had to use for producing the Landsat-based metrics in order to deal with low satellite observation densities, as well as errors and inconsistencies in the satellite-based disturbance maps (Oeser et al. 2017). In areas with better Landsat data availability, metric time series with annual temporal resolutions can be produced without using moving window averages to improve to temporal accuracy of the habitat monitoring. In addition, information from independent forest disturbance maps can also be incorporated directly into habitat models to ensure the temporal consistency between disturbance maps and habitat time series (e.g., by using the number of years since disturbance as a model predictor; Kearney et al. 2019).

Implications for wildlife and forest management
Overall, widespread disturbances in our study area have substantially improved habitat conditions for red and roe deer (see Appendix S3: Fig. S1), also implying increases in the carrying capacity of the landscape for deer. The frequency of natural disturbances in European Forests has been increasing over the last decades (Senf et al. 2018), and likely will continue to increase as a consequence of global warming (Seidl et al. 2017). Thus, considering disturbance impacts on deer habitat will be increasingly important in the management of ungulate populations (e.g., when setting hunting quotas). Deer overabundance in many parts of Europe and North America already has considerable negative economic and ecological effects on forests (Côt e et al. 2004). As disturbed sites provide attractive foraging habitat, they lead to a concentration of deer, as well as higher levels of browsing pressure in forest openings (Kuijper et al. 2009). More frequent disturbances thus have the potential to trigger an exacerbating feedback on forest structure: Deer densities could increase because of improved resource availability (Gaillard et al. 2003), which in turn increases browsing damages that impede forest regeneration (Tremblay et al. 2007). The relationship between more frequent forest disturbances and ungulate populations dynamics, however, remains poorly understood. Linking observed trends in ungulate populations to satellite-derived time series of habitat dynamics opens new possibilities to explore possible links across space and time.