Spatial and temporal dynamics of Mexican spotted owl habitat in the southwestern US

Understanding habitat dynamics is essential for effective conservation as landscapes rapidly change. In a companion paper in this issue, Shirk et al. (2023) introduced an automated habitat monitoring system using Google Earth Engine and applied this framework to develop a dynamic model of Mexican spotted owl (Strix occidentalis lucida) habitat across the southwestern US from 1986 to 2020. We explored the application of this dynamic model of Mexican spotted owl habitat in the context of the species’ ecology. We evaluated environmental correlates of Mexican spotted owl habitat, assessed potential spatial non-stationarity in habitat selection, estimated long-term trends in habitat by quantifying changes in habitat amount and quality between 1986 and 2020, and evaluated the extent to which habitat changes over the past 35 years have been driven by wildfire. Topography and climate appeared to outweigh reflectance-based (vegetation) metrics in describing Mexican spotted owl habitat and habitat selection was non-stationary across modeling sub-regions. Total habitat area for Mexican spotted owls declined by ~ 21% since 1986 (0.6% annually), but trends varied spatially and some even reversed over the past decade. Wildfire was responsible for between 8 and 35% of total habitat loss, depending on the sub-region considered. The automated habitat monitoring system allowed trend estimation and accurate assessment of current habitat status for Mexican spotted owls; maps were accurate, spatially detailed, and current. The ability to continually produce accurate maps for large land areas for threatened species such as the Mexican spotted owl facilitates science-based land management on public lands in the southwestern US.


Introduction
Habitat loss and fragmentation are critical threats to biodiversity (Pardini et al. 2017;Tilman et al. 2017;Fletcher et al. 2018). For forest-dependent species, habitat loss and fragmentation are often driven by multiple factors, including timber harvesting (LaManna and Martin 2017) and land-use conversion (Gaston et al. 2003;Cushman et al. 2017). But, wildfire (van Wees et al. 2021) and drought-induced tree mortality (Allen et al. 2015), which are now both strongly linked to climate change (Diffenbaugh et al. 2015;Abatzoglou and Williams 2016), are rapidly transforming forest structure and composition (McDowell et al. 2020). These changes are likely to lead to large-scale habitat alteration and loss for forest-dependent species (Roberts et al. 2019;Kelly et al. 2020;Ward et al. 2020).
Patterns of fire effects on species habitat (Fontaine and Kennedy 2012;Jones and Tingley 2022) and tree mortality following disturbances are far from uniform (Guarín and Taylor 2005;Stevens et al. 2017). Quantifying and mapping such patterns concomitantly with their effect on habitat loss is complex, partly because of how rapidly landscape changes are occurring. Given the complexity and magnitude of recent disturbance impacts and the pace at which these are occurring in the American southwest (e.g., fire alone has affected 4.1 million ha in Arizona and New Mexico in the last three decades; Singleton et al. 2019) it is clear that new methods are needed to assess such landscape scale impacts on threatened species.
Species that occur at low densities, occupy naturally fragmented habitat, and specialize on rare resources are often considered to be at greater risk of extinction because of habitat loss (Davies et al. 2004). One such species is the Mexican spotted owl (Strix occidentalis lucida, Fig. 1; hereafter "MSO") whose range consists of naturally fragmented forest and canyonland habitat, characterized by disjunct mountain ranges sometimes referred to as "sky islands" across the southwestern United States and northern Mexico (Gutiérrez et al. 1995). However, the Mexican spotted owl is also a rare old-forest specialist that was listed as "Threatened" under the United States Endangered Species Act (ESA) in 1993 because of actual and potential loss to its old-forest habitat from (i) "even-aged forest management" that resulted in loss of forest structure and (ii) large, stand-replacing wildfires (USFWS 1993(USFWS , 2012. Consequently, these imminent threats necessitate the application of monitoring tools capable of assessing habitat conditions across a broad and highly variable range. In a companion paper in this issue, Shirk et al. (2023) describes an automated habitat monitoring system for modeling habitat over time using geographically extensive and temporally dynamic modeling in a cloud-computing framework and applied that framework to MSO in the southwestern US. Briefly, Shirk et al. (2023) fitted random forest species distribution models to MSO nest/roost data using a suite of 40 environmental covariates on Google Earth Engine (GEE). Here, we analyzed outputs from this automated habitat monitoring system for MSO and provided an ecological interpretation of their models. We analyzed the 35-year time series (1986-2020) of habitat models produced by Shirk et al. (2023) to (1) evaluate MSO habitat associations with vegetation, topography, and climate, their relative importance, and the extent to which these associations fluctuated across the modeling region, a phenomenon often called non-stationarity in habitat selection (Osborne et al. 2007;Miller and Hanham 2011); (2) evaluate long-term trends in MSO habitat to aid recovery planning; and (3) attribute the extent to which habitat changes have been driven by recent wildfire. We interpret the conservation implications of this dynamic modeling effort in light of MSO ecology in the context of the future role of the automated habitat monitoring system for MSO to assist land management planning.

Development of MSO occurrence database
We compiled a database of known Mexican spotted owl nest/roost locations from 1989 to 2020 recorded on national forests in Arizona and New Mexico (Fig. 2). This database included nest/roost sites from eight independent data sets, including two previous demographic studies (Seamans et al. 1999), species recovery planning efforts (USFWS 1995;White et al. 1995;May et al. 1996) Service project-level surveys, and opportunistic nest/ roost observations reported by cooperating observers (Table S1). Relevant to our interest in conservation of owls, the U.S. Fish and Wildlife Service has partitioned the range of the Mexican spotted owl in the southwestern USA into Ecological Management Units (EMUs; USFWS 2012), five of which were within our study area (see below). Thus, we focused our analysis on these five subregions.
We began with a database of nest/roost locations (n = 7455), but applied decision rules for quality control to yield a refined database of high-quality locations. In this case, we applied five decision rules for quality control to maximize data quality and standardization, reduce pseudoreplication, and reduce the likelihood that observations might represent behaviors other than nesting or roosting. First, we eliminated roost location(s) when surveyors detected a nest location during the same survey. This situation arose frequently in data originating from demographic studies (Seamans et al. 1999) because one member of the owl pair or young were often detected on the nest while another owl was detected roosting elsewhere; we considered nest locations to contain more biological information than roost locations. Second, we used only one roost location when two were recorded for a pair at the exact same location (e.g., pair roost). Third, we selected only observations that occurred during the breeding season (1 March through 31 July). Observations after this period are less likely to represent breeding habitats as owls begin to range more widely in late summer (Ganey 1988;Ganey and Balda 1989). Fourth, we selected only one nest location per territory per year. Spotted owls typically do not renest if their first nest fails (Gutiérrez et al. 1995) so multiple within-year nest locations likely reflect recording inaccuracies of the same nest location. If only two locations existed (n = 2), we selected the earlier one because early season locations are more likely to be determined from direct observations of owls on the nest before adults and offspring begin to wander more widely in their home range. If more than two locations existed (n > 2), we selected the most frequently recorded nest location. Fifth, we eliminated any detections recorded as roosts that occurred during twilight or nocturnal hours (adjusted by geographic location within the region) because owls may have already moved from their day roost for nocturnal foraging activities by the time surveyors recorded these detections (Ganey 1988;Delaney et al. 1999).
Applying our decision rules reduced the original raw data set of n = 7455 nest/roost locations to n = 2913 nest/roost locations (Table S1). Sample sizes and temporal sampling frames in this final dataset varied considerably among the five EMUs that we modeled as sub-regions as follows: Basin andRange East, n = 312 (1990-1993); Basin andRange West, n = 293 (1990-2018);Colorado Plateau, n = 13 (1990-2002; Southern Rockies, n = 33 (1989Rockies, n = 33 ( -2020; Upper Gila Mountains, n = 2253 (1990-2019). We then re-projected all locations to the WGS84 coordinate system for compatibility with available high-resolution satellite imagery and performed visual inspection of all locations overlaid on that imagery to ensure locations were not mapped in obvious non-forested areas (e.g., agricultural lands) because of data recording errors. This final dataset was then used to train the dynamic models described in Shirk et al. (2023).

Potential non-stationarity in MSO habitat
We evaluated the potential for non-stationarity in MSO habitat selection across EMUs in two ways: (i) variable importance and (ii) predictability. Potential spatial non-stationarity in habitat selection was allowed by Shirk et al. (2023) by fitting a different species distribution model for each of the 5 EMUs.
First, we identified the factors that shaped MSO habitat by computing the variable importance for the 40 predictor variables (22 reflectance-based, 9 climatic, 9 topographic) used to fit the EMU-specific species distribution models in Shirk et al. (2023). These three groups of covariates were selected because previous work has suggested the interplay between them shapes MSO habitat in the arid Southwest (Peery et al. 2012;Wan et al. 2017). We then assessed the degree to which these factors varied across EMUs. Variable importance describes the degree to which each predictor contributes to model performance (Thuiller et al. 2009). We computed the mean variable importance across 10 bootstrapped replicate model runs within the GEE environment (Gorelick et al. 2017) to summarize variable importance for each EMU. Apparent variable importance of individual predictors in random forest models can be incorrect in the presence of high collinearity (Genuer et al. 2010;Stijven et al. 2011). We found that correlations were generally high within the three different variable types (reflectance, climatic, topographic) but low among variable types (Fig. S1). We therefore avoided interpreting and comparing the relative importance of individual predictors and instead focused on comparing relative importance among pooled variable types; if correlations between types are low, comparing relative importance among types of variables in the models should remain valid (see below). We plotted response functions for all variables to assess their biological meaning and potential differences across EMUs. These marginal relationships may share explanatory power among correlated variables, but the relative magnitudes of the marginal effects can be interpreted meaningfully.
We evaluated the relative contributions of climate, topography, and vegetation (reflectance) on MSO habitat to variation among EMUs by summing variable importance scores within each predictor class. We then compared the cumulative variable importance for each class to a null expectation of variable importance. Our null expectation was simply the fraction of variables contributed by each class; this represented the expectation that each variable would contribute to overall variable importance in proportion to its abundance.
Second, we evaluated non-stationarity by projecting models that were trained within each EMU to the other EMUs and computed predictive transferability using AUC. We computed AUCs across 10 bootstrapped replicate model runs with available sites conserved across the bootstrapped replicates. Nonstationarity in habitat selection was evident if models built for one region (EMU) poorly predict habitat for other regions.

Evaluating trends in MSO habitat
We evaluated trends in MSO habitat from the 35-year time series of annual model projections produced by Shirk et al. (2023). Each annual model was a gradient surface, where the continuous value in each pixel represented relative probability of nest/roost occurrence. Accordingly, to compute equivalent habitat area from a continuous habitat probability surface, the probabilities can be summed within a specified area; in a correctly calibrated random forest model, the sum of predicted probabilities precisely equals the number of cells that have experienced habitat loss . Therefore, we computed the area of habitat for each year 1986-2020 by summing the predicted probability surfaces within each EMU. This can be thought of as a weighted sum; the area of all pixels in each EMU are summed but are weighted by the probability that they contain habitat. We evaluated evidence for habitat trends, and the slope/direction of trends by fitting linear trend lines to the annual habitat area within each EMU, and for the whole study area combined. We evaluated trends over the entire 35-year study period  in addition to the most recent 10 years (2011-2020) because the latter corresponds with the time period required to assume a stable or increasing habitat trend in the MSO recovery plan (USFWS 2012).

Drivers of habitat change: the role of wildfire
We evaluated the extent to which wildfire had contributed to MSO habitat declines by linking the habitat probability model to a spatially extensive fire perimeter dataset that spanned the modeling period (Monitoring Trends in Burn Severity program, www. mtbs. gov). Ideally, we would evaluate multiple drivers of habitat change, including insects, disease, drought, harvesting, and other drivers in addition to wildfire; however, accurate and/or complete spatial databases of other drivers do not currently exist.
We associated declines in MSO habitat with wildfire by computing the total annual change in habitat area for a given year (for each EMU), and then the proportion of that total change that occurred within fire perimeters (sensu Wan et al. 2020). We used a two-year moving window to compute annual change in habitat for year t for a given EMU; we subtracted habitat area from year t + 1 from year t − 1 to account for fires that occurred after the mapped habitat date for a given year. So, for example, habitat change for the year 2000 (habitat 2000 ) was computed as habitat 2001 − habitat 1999 . Similarly, the change attributable to fire in year t was then computed as the habitat within fire perimeters in year t + 1 minus the habitat within fire perimeters in year t − 1 such that habitat change caused by fire for the year 2000 (fire 2000 ) was computed as fire 2001 − fire 1999 . We could then compute the proportion of habitat change owing to fire by division (e.g., fire 2000 /habitat 2000 ), which we computed for all years 1986-2020.

Results
We observed a relatively high degree of spatial nonstationarity among modeling regions with respect to the environmental factors that shaped MSO habitat (Fig. S2); different factors appeared to govern MSO habitat across modeling subregions. For example, while monsoon season precipitation appeared to most strongly describe MSO habitat in the Basin & Range West and Upper Gila EMUs, its importance appeared very low in the Colorado Plateau; similar patterns can be found for numerous other variables (Fig. S2); however, we caution interpretation of individual variable importance in the presence of collinearity (see Methods and Fig. S2). However, when variable types were pooled, topographic and climate variables appeared to dominate variable importance across all EMUs (Fig. S2). Topographic and climate variables explained more of the total variable importance than a null expectation (Fig. 3). In contrast, reflectancebased variables explained an average of 49% of total variable importance, which was lower than expected given they constituted 55% of predictor variables (Fig. 3).
Spatial non-stationarity was also apparent when examining variable response curves across EMUs (see Fig. 4 for illustrative variables; see also Figs. S2-S4 for full reporting). In many cases the direction of the effect of a predictor variable on habitat probability was similar across EMUs, but the shape of response curves showed considerable variation ( Fig. 4a-h). In general, habitat probability was higher in areas characterized by greater monsoon season precipitation, (Fig. 4b), more topographic complexity (Fig. 4g), and less cumulative heat represented by degree days > 5 °C (Fig. 4h). Other relationships tended to be unimodal with a peak in habitat probability at intermediate values including mean winter temperature (Fig. 4c), winter precipitation (Fig. 4e), and precipitation falling as snow (Fig. 4f). Areas characterized by drainages (lower topographic position index) were associated with higher habitat probability (Fig. 4a); however, the relationship between topographic position index and habitat probability was bimodal, where flatter areas (topographic position index ≈ 0) were less associated with habitat. Areas with warmer minimum winter temperature generally had lower habitat probability, but this relationship was opposite for the Southern Rockies EMU (Fig. 4d).
Non-stationarity was apparent because of quantitative variation in predictive performance among modeling subregions (Fig. 5). Models that were trained in each region unsurprisingly performed best Fig. 3 Relative importance of predictors related to climate, reflectance, and topography for each EMU. The height of the bars is the sum of variable importance scores for all predictors within a given variable class (e.g., see Fig. S2). Error bars are the pooled standard deviation across 10 bootstrapped replicates. The horizontal black bars show a null expectation of variable importance for each class based on the abundance of variables in each category applied to all other EMUs. Colors represent the region in which the model was trained, and the x-axis represents the region into which the model was projected. The height of the bars represents the average AUC across 10 random forest model runs as part of k-fold cross-validation (Shirk et al. 2023). Error bars show the mean ± 1 standard error averaged across the ten model runs when they were projected onto that same region, with AUCs ranging from 0.993 to 0.998. Model performance decreased when models were projected outside the regions in which they were trained, although performance still remained reasonably high (all AUCs > 0.6). In contrast, the model trained in the Basin and Range West EMU always performed reasonably well when projected to other regions, with AUCs ranging from 0.85 to 0.95 (Fig. 5).
The model estimated that in 1986, there was approximately 1 M ha (999,688 ha) of MSO nest/ roost habitat in Arizona and New Mexico. That amount declined by 21% over the study period (0.6% annual decline) to 788,578 ha by the year 2020 (Fig. 6a). Declines in habitat from 1986 to 2020 were observed in all EMUs (Fig. 6b). The largest percent decline was observed in the Basin and Range West EMU (34% decline, ~ 0.98% per year) whereas the smallest percent decline was observed in the Southern Rockies EMU (7% decline, 0.2% per year).
The pattern and rate of decline was not constant over the study period 1986-2020 or among EMUs (Fig. 7). Overall rate of decline across the study area and within all individual EMUs (except for the Upper Gila EMU) was greater over the 35-year study period than it was in the most recent 10-year period (Fig. 7), when the rate of decline slowed. Conversely, in the Basin and Range West and Southern Rockies Fig. 6 Trends in Mexican spotted owl habitat from 1986 to 2020 a across the entire study area (all 5 Ecological Management Units, or EMUs) and b within each EMU. The y-axis in each panel shows 1000s of hectares in habitat (i.e., a value of 500 represents 500,000 ha)

Fig. 7
Linear trend coefficients estimated from long-term Mexican spotted owl habitat data. The y-axis represents the model coefficient from a linear fit (habitat area ~ year), which is interpreted as the annual change in habitat in 1000s of hectares. Circles represent the coefficient using all 35-years of habitat data, while triangles represent the coefficient when only using the most recent 10-year period (2011-2020). Error bars show 95% confidence intervals around the point estimate for the linear coefficient. Error bars are larger for the 10-year period because of smaller sample sizes EMUs, the most recent 10-year period showed a positive trend in MSO habitat (Fig. 7) even though overall trends from the 35-year period were negative.
Wildfire accounted for 17% of the loss in habitat area across the study area over the period 1986-2020 (Fig. 8a). However, there was relatively high spatial variation in this effect across EMUs (Fig. 8b). In both the Southern Rockies and Basin and Range West EMUs, wildfire accounted for 35% of habitat loss over the study period. In the remaining EMUs, wildfire accounted for less of the overall loss, ranging from 8 to 15%.

Discussion
In this paper, we assessed the ecological implications of the dynamic species distribution model developed by Shirk et al. (2023). We drew four main conclusions from our assessment: (i) MSO habitat selection appears to be shaped and limited by different environmental features across large portions of its range in the southwestern US (spatial non-stationarity), (ii) topography and climate disproportionately explain MSO habitat compared to vegetation (reflectance), (iii) MSO habitat area has declined over the past 35 years, but the rate of decline has slowed and in some cases reversed, and (iv) the degree to which wildfire has contributed to habitat declines varies widely across subregions.

Implications of spatial non-stationarity and habitat associations
Non-stationarity in habitat selection suggests that the same conservation actions may influence a species in different ways depending on the location within its range if there are different limiting ecological factors among locations (e.g., see Rissman and Wardropper 2020). Moreover, non-stationarity implies that different thresholds may be needed to trigger management actions across space (Nie and Schultz 2012;Schultz et al. 2013;Schwartz et al. 2015). Our work showed apparent strong non-stationarity in habitat selection by MSO in the southwestern US, indicating that environmental conditions that predict nest/ roost habitat in one EMU may not adequately predict nest/roost habitat elsewhere (as previously shown by Wan et al. 2019a). Interestingly, models trained in the Basin and Range West EMU always performed relatively well when applied to other sub-regions (Fig. 5). We think this is because the Basin and Range West EMU encompasses a broader range of environmental conditions suitable for nesting and roosting activities that is well-represented across the other EMUs, making predictions more transferrable. In contrast, Fig. 8 The role of wildfire in Mexican spotted owl habitat loss for a the entire study region and b each individual Ecological Management Unit (EMU). Darker shaded lines shows the total habitat area lost in each year 1986-2020, while the lighter shaded lines show the total area lost to wildfire environmental conditions in other EMUs may be more uniform and therefore models perform poorly when projected to 'novel' conditions. Ecological drivers of MSO habitat varied across EMUs, but topographic and climatic predictors consistently outperformed reflectance-based predictors ( Fig. S2 and Fig. 3). The importance of topography and climate in shaping MSO habitat in the arid and naturally fragmented southwestern US has been shown previously (Timm et al. 2016;Wan et al. 2017Wan et al. , 2019a. However, Wan et al. (2017) showed that canopy cover was the most important predictor of multi-scale MSO habitat in both the Sacramento mountains (located in Basin & Range East EMU) and the Mogollon Plateau (located in Upper Gila EMU); similar results were found by Timm et al. (2016). Individually, climate and topographic variables were most important to MSO habitat in the present study, but a broad range of reflectance variables were also useful at explaining variation in selection (Fig. S2). So, while no single reflectance covariate appeared to carry much weight, collectively they explained a considerable fraction of model variability and therefore also form a key model component.
Several factors can potentially explain differences in results between our study and previous studies investigating ecological drivers of MSO habitat. First, Shirk et al. (2023) used unclassified Landsat reflectance images as predictors representing the vegetation component of MSO habitat. Unclassified imagery is less susceptible to classification error and therefore is likely to improve performance in habitat models (McKelvey and Noon 2001). Yet, the biological interpretation of reflectance is less clear, whereas canopy cover is a metric with clear biological meaning for MSO; previous studies may have identified this and other vegetation predictors as relatively more important for this reason. Additional explanations for this discrepancy could include collinearity between covariates (Genuer et al. 2010;Stijven et al. 2011) and that previous studies used multi-scale optimization during model fitting Wan et al. 2017). Second, different datasets were used in Shirk et al. (2023) compared to previous studies; sampling variation could have given rise to variation in those variables deemed most important for describing each unique dataset. Third, different analytical approaches were employed; Timm et al. (2016) and Wan et al. (2017) both used standard logistic regression approaches and multi-scale optimization to model habitat whereas Shirk et al. (2023) used random forest without scale optimization. Random forest has been shown to outperform logistic regression in predictive habitat modeling (Cushman and Wasserman 2018) and can identify non-linear response functions that could lead to different species-environment relationships. Regardless of differences between our efforts to quantify environmental features associated with MSO habitat and previous studies, it is clear from all existing research that MSO habitat cannot be adequately described without topography and climate.

Spatial and temporal habitat dynamics
Compared with factors shaping MSO habitat, relatively little is understood about how MSO habitat has changed in recent decades. Yet, understanding recent habitat trends is a requirement for ESA delisting outlined in the MSO recovery plan (USFWS 2012). Recovery of threatened and endangered species like the MSO requires effectively monitoring trends in habitat and populations (Robinson et al. 2018). The need to understand trends in-and current conditions of-MSO habitat is underscored by rapid changes to the landscapes they inhabit Hessburg et al. 2019). In recent years, MSO habitats in Arizona and New Mexico have seen large increases in the area burned by wildfires (Singleton et al. 2019); future projections of fire area and severity with climate change suggest considerable increases in fire impacts on MSO habitat in the coming decades (Littell et al. 2018;Wan et al. 2019b).
We found that MSO habitat declined since the mid-1980s, but declines have slowed and, in some cases, reversed in the past decade (Fig. 7). Currently, we do not know why habitat declines have slowed and/or reversed; this is an area we are exploring in ongoing research. Yet, the spatial variation in habitat trajectories has important conservation implications because different trends could be partly the result of different management actions, in addition to other cofactors (e.g., spatial variation in wildfire occurrence). Further research will be needed to determine what management drivers (if any) are associated with trends. While habitat gains might indeed be attributed to recruitment of conditions suitable for MSO nesting and roosting through time, map error and uncertainty could also contribute to unrealistic gains in suitable habitat (Witt et al. 2022). Accurate mapping of growth and recruitment of old-forest conditions with Landsat data can be affected by several artifacts, including canopy shadows and 'false recruitment' produced by spectral signals that are similar to oldforest conditions but are actually from shrub and forb regreening that can occur rapidly after a high-severity fire (Davis et al. 2015). Future updates to the models will seek to improve techniques that account for this potential source of bias. Nevertheless, a key takeaway is that the trajectories of MSO habitat as well as the factors shaping habitat across this large bioregion are highly variable.
Evaluating drivers of landscape change is a rapidly advancing research area . While many factors that contribute to landscape change are difficult to measure, fire records in the continental US are relatively complete (Eidenshink et al. 2007;Finco et al. 2012), which facilitated our ability to link wildfire with changes in habitat. We demonstrated that 17% of the total decline in MSO habitat across the southwestern US from 1986 to 2020 could be associated with wildfire, but the contribution was as high as 35% in the Basin and Range West and Southern Rockies EMUs. Therefore, this finding supports the conclusion that wildfire may pose a major threat to the persistence of MSO habitat in this region (USFWS 2012). Moreover, other projections have indicated that area burned in the MSO range will increase 13-fold by the 2080s, with increased risk of severe fire (Wan et al. 2019b). High-severity fire was associated with relatively large declines in predicted habitat quality (Wan et al. 2020), while extensive high-severity fire have resulted in persistent population declines (Jones et al. 2016 and foraging avoidance (Jones et al. 2020b;Kramer et al. 2021) in the California subspecies (S. o. occidentalis). Clearly, increasing extent and severity of wildfire in the southwestern US poses a conservation threat to MSO Wan et al. 2018;Jones et al. 2020a).
Other important natural and anthropogenic drivers contributed to observed declines that we were unable to evaluate fully at this time. Although wildfire explained a meaningful amount of overall habitat decline (8-35%), most variation in habitat loss remains to be explained. For this reason, we think that relatively complete and accurate spatial layers of other landscape drivers of change (e.g., insects, disease, drought, timber harvest, fuels management) are needed to conduct a more comprehensive analysis of the causes of habitat decline. Moreover, interactions among change factors (e.g., drought and insect outbreaks) also complicate linking individual factors to declines (e.g., salvage logging following wildfire). As automated habitat monitoring systems are applied in other systems, a robust approach for attributing causes of habitat change will be needed to allow conservation efforts to react to new knowledge about what is causing declines.
Attempting to link potential drivers to habitat change will likely provide opportunities for rapid learning. For example, we observed a major, synchronous loss in habitat over the 2001-2002 period across three EMUs (Fig. 8), suggesting the role of a regional climatic event on habitat loss. Indeed, the southwestern US was affected by a major drought over at least the period 1999-2004 (Piechota et al. 2004). In particular, previous research has shown that extremely warm and dry conditions in our study area during 2002 triggered sudden and widespread tree mortality (Ganey and Vojta 2011;Kane et al. 2014;Ganey et al. 2021), which coincided with the period we observed abrupt and large-scale habitat loss. Little is currently known about how spotted owls and their habitat are affected by drought 'events', but warmer and drier conditions associated with climate change are predicted to reduce owl populations in Arizona and New Mexico (Peery et al. 2012). Owl habitat that contains high tree densities may be more susceptible to drought-related tree mortality in the southwestern US (Bradford and Bell 2016). If this is indeed the case, active management that reduces drought-related tree mortality might benefit owls. For example, thinning of smaller trees and increased use of prescribed fire, both of which support forest resilience to drought and wildfire, have the potential to buffer owl habitat against these events in the future (Stephens et al. 2018;North et al. 2021). Certainly, more research is needed to identify linkages between less 'discrete' events such as drought (and even climate change) and species' habitat change.

Additional limitations and caveats
The automated habitat monitoring system, as applied to MSO (Shirk et al. 2023), was imperfect, but no habitat monitoring system has been. So, we provide some important caveats.
Sampling of MSO nest/roost locations was not systematic or balanced in space and time, and therefore limits inferences in some ways. Sample size varied strongly across modeling sub-regions (see Methods in Shirk et al. 2023), which means our confidence in models also varies spatially. Even though model performance metrics (e.g., AUCs; Fig. 5) were still good even in low-sample conditions, additional samples would likely have a strong influence on the spatial distribution of habitat in low-sample regions (in particular, the Colorado Plateau EMU, n = 13 and the Southern Rockies EMU, n = 33); indeed, a critical part of the adaptive management loop is to identify areas that require more sampling (Schwartz et al. 2015). However, random forest is thought to be one of the most robust algorithms for small sample size (Qi 2012). In addition, while the Ecological Management Units, or EMUs, were a useful modeling sub-region because they align with recovery planning efforts (USFWS 2012), they were partly an administrative construct, which often does not reflect biology. Therefore, we think there are likely other, more robust approaches for developing more ecologically appropriate modeling sub-regions based on multivariate environmental similarity that future work could explore to better describe spatial non-stationarity. Future research could explore the potential role of non-representative sampling and other sampling biases in shaping inferences using simulation methods to control these factors (Atzeni et al. 2020;Chiaverini et al. 2021).
Temporally unbalanced MSO nest/roost locations lead to other important caveats related to the interpretation of our results. Temporal clustering of sampling locations (i.e., most of our data comes from the period 1990-1993) might mask changing climate-habitat relationships and climatic thresholds. For example, if data are mostly collected within a cool/wet period and MSO have an unmeasured hot/dry tolerance boundary that is crossed as warming occurs during our study period, climate impacts on habitat trend could be underestimated. Likewise, we have made the implicit assumption that MSO nest/roost habitat selection patterns were stationary through time. It was possible that MSO changed how they select nest/roost habitat through time (temporal non-stationarity), but because our current database was under-represented in more recent years it could not support such an analysis. Ideally, longitudinal owl data could inform the way nesting/roosting habitat is affected by longer-term climate variability. Our ability to estimate long-term trends in MSO habitat, and our specific estimates of habitat loss (e.g., 21% over the study period) is therefore predicated on the validity space-for-time variance substitution (Pickett 1989).
Our approach for estimating habitat area and computing habitat loss is robust but requires careful attention by the reader to avoid misinterpretation. Specifically, we viewed habitat as a continuous gradient, which, while providing certain advantages, also carries caveats about how habitat 'loss' should be interpreted. Studies of habitat loss and fragmentation increasingly represent habitat as a gradient rather than a habitat vs. non-habitat binary, because of information loss that occurs using a binary approach (Fardila et al. 2017). Our procedure for computing overall habitat change, or 'loss', implicitly includes types of habitat 'degradation', which may reduce the relative suitability of habitat but does not represent a binary 'loss'. If the predicted habitat probability declines from 0.5 to 0.4 (or even from 0.2 to 0.1) across a group of pixels, this more subtle change can result in habitat 'loss' according to our method. Moreover, multiple ecological scenarios can produce the same total habitat area (or loss). For example, ten pixels with a probability of 0.1 and one pixel with a probability of 1 would produce the same habitat area. This gradient-based approach for quantifying habitat is philosophically distinct from a binary threshold-type approach (e.g., McGarigal and Cushman 2005;) because habitat is not a categorical entity and habitat loss is not a binary event in most cases. Thus, our gradient-based measure of habitat change is more realistic, sensitive, robust, and less prone to classification errors related to commission and omission than typical binary measures of loss based on thresholds of predicted probability. We recognize that although our approach may be more holistic, it does not isolate potential changes to different 'levels' of quality habitat (e.g., 'best' quality habitat, where probabilities are close to 1). As a result, inferences must be adjusted accordingly.

Concluding remarks
One key benefit of the automated habitat monitoring system (Shirk et al. 2023) is that products can be regularly updated to address existing limitations through time. For example, for modeling sub-regions in which we currently have limited MSO nest/roost data, future targeted surveys could improve sample sizes to improve confidence in inferences and mapped habitat. Additional data that spans a greater temporal window could also inform models that examine temporal nonstationarity. Future iterations could also hone modeling approaches, apply different algorithms, and add new covariates based on new ecological information. Therefore, automated habitat monitoring systems not only facilitate tracking of near real-time changes in habitat (Shirk et al. 2023), but they also can facilitate the use of the scientific process to update, refine, and improve models as new information becomes available (Cushman and McKelvey 2010), which is the key for predictive models to be powerful tools in the adaptive management cycle (Cushman and McGarigal 2019). From a species conservation perspective, automated habitat monitoring systems provide a tool for rapid identification of recovery goals related to habitat and disturbances. Updating, calibrating, and improving broad-scale predictive models, such as those evaluated here, will provide continually improving information to managers, and hopefully, better conservation outcomes for threatened species like the MSO. David Anderson, Jack Triepke, and Ryan Heaslip contributed to important discussions that shaped the direction of this study. Over the years, numerous field technicians contributed to the collection of data used to develop these models and help us learn about how to better conserve the owl; we deeply appreciate and value their contributions.
Funding This work was supported in part by the USDA Forest Service Southwestern Region.

Conflict of interest
The authors have not disclosed any competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.