Linking monitoring and data analysis to predictions and decisions for the range-wide eastern black rail status assessment

The US Fish and Wildlife Service (USFWS) has initiated a re-envisioned approach for providing decision makers with the best available science and synthesis of that information, called the Species Status Assessment (SSA), for endangered species decision making. The SSA report is a descriptive document that provides decision makers with an assessment of the current and pre dicted future status of a species. These analyses support all manner of decisions under the US Endangered Species Act, such as listing, reclassification, and recovery planning. Novel scientific analysis and predictive modeling in SSAs could be an important part of rooting conservation decisions in current data and cutting edge analytical and modeling techniques. Here, we describe a novel analysis of available data to assess the current condition of eastern black rail Laterallus jamai censis jamai cen sis across its range in a dynamic occupancy analysis. We used the results of the ana lysis to develop a site occupancy projection model where the model parameters (initial occupancy, site persistence, colo nization) were linked to environmental covariates, such as land management and land cover change (sea-level rise, development, etc.). We used the projection model to predict future status under multiple sea-level rise and habitat management scenarios. Occupancy probability and site colonization were low in all analysis units, and site persistence was also low, suggesting low resiliency and redundancy currently. Extinction probability was high for all analysis units in all simulated scenarios except one with significant effort to preserve existing habitat, suggesting low future re siliency and redundancy. With the results of these data analyses and predictive models, the USFWS concluded that protections of the Endangered Species Act were warranted for this subspecies.


INTRODUCTION
The emerging approach of the US Fish and Wildlife Service (hereafter USFWS) to decision making on protected and candidate species under the US Endangered Species Act of 1973 (ESA) uses a Species Status Assessment (SSA) to gather and organize available data and information on a species of interest (Smith et al. 2018). The SSA report is intended to be the scientific document that supports ESA decision making. The document contains a de scription of the biology and ecological needs of the species, an assessment of current availability of needed resources and of populations in the wild (i.e. current status assessment), and predictions on future status of populations and resources (McGowan et al. 2017, Smith et al. 2018, Voorhies et al. 2019. The SSA framework uses the concepts of 'representation', 'resiliency', and 'redundancy' (the '3 Rs') as the lens through which status is evaluated, but each SSA can define the 3 Rs using different metrics that suit the species and available data.
The USFWS, and any management agency dealing with imperiled species, is often challenged with timesensitive decisions that make designing and implementing monitoring programs specifically designed for a given decision difficult or impossible (Smith et al. 2018). Thus, any analyses and predictive modeling of current and future status must use the best available data for inference (McGowan et al. 2017, Smith et al. 2018. When timesensitive decisions are imminent, analysts are often required to use statistical models that can be adapted to the available data (e.g. Skalski et al. 2010). The USFWS received a petition to list a subspecies, the eastern black rail Laterallus jamaicensis jamaicensis, under the ESA, and we analyzed available absence data (MacKenzie et al. 2002), and developed a predictive model for their populations, to support the decision process. Eastern black rails (Eddleman et al. 1994) are small elusive marsh land birds with high habitat specificity. Eastern black rails were historically widespread but apparently have exhibited a substantial range contraction over the last ≥ 50 yr (Watts 2016). The USFWS was petitioned to list the eastern black rail as a threatened or endangered subspecies and in advance of making that determination, the USFWS developed an SSA following the guidelines of Smith et al. (2018). For this, the USFWS required a rapid analysis of current and future status to meet the decision deadlines. Smith et al. (2018) presented an SSA framework that has 3 core components: (1) a review and description of the biology and ecological needs of the species; (2) an analysis and assessment of the current condition of these needs (i.e. availability of re sources) and populations on the landscape (i.e. demo graphic, abundance, and/or trend estimation); (3) projection of future condition involving some type of probabilistic prediction about future status of needed resources and population metrics. The analysis portions of the SSA (i.e. steps 2 and 3) can range in complexity from qualitative descriptions of species range and status, to complex analyses of existing data sets and projection models (Smith et al. 2018). Novel analyses of existing data sets for an SSA can uncover new ecological relationships or quantify existing ideas and hypotheses about ecological relationships. These analyses thus inform the species' needs and current status components, and also form the basis for making future status predictions. Here, we present novel analyses of presence-absence data for eastern black rails from throughout their range and a simulation model to predict future status for the subspecies. Specifically, we present the analysis of data to quantify ecological needs (SSA component 1) and estimate current status (SSA component 2) using a dynamic occupancy analysis (MacKenzie et al. 2003) and then a projection model to predict future status (SSA component 3; Table 1). 210

SSA components
General approaches Our approach Focus Species biological and Literature review, data Literature review, dynamic NA ecological needs analysis, expert elicitation occupancy analysis with assessment covariates and model selection Current condition Literature review, data Dynamic occupancy analysis, Current redundancy and analysis, expert elicitation assessment of current occupancy, resiliency extinction, and colonization probabilities

Future condition
Predictive modeling Occupancy projection modeling Future redundancy and using parameter estimates from resiliency dynamic occupancy analysis under 6 possible future scenarios of sea-level rise, habitat management, and conservation effort Table 1. Organization of the 3 core components of the species status assessment (SSA) framework (Smith et al. 2018), the possible approaches to addressing those components, how we addressed those components, and on which of the so-called '3 Rs' (representation, resiliency, and redundancy) our analysis focused. NA: not applicable SSAs put metrics about the current and future status of populations and species into the '3 Rs' framework (resiliency, redundancy, and representation). Within the 3 Rs, resiliency is defined as the ability of a population to withstand stochastic events, redundancy is a measure of the number of resilient populations to withstand catastrophic events, and representation assesses whether a species has sufficient genetic and phenotypic diversity to maintain adaptability (Smith et al. 2018). We attempted to adapt the dynamic occupancy parameter estimates (MacKenzie et al. 2003, Popescu et al. 2012 into the 3 Rs terminology for decision makers, focusing on resiliency and redundancy, to describe current status ( Table 1). The SSA included some analysis to address representation across the range, but that analysis is outside the scope of this paper. Our projection model to predict future status outputs the proportion of sites re maining occupied each year, which we also attempted to translate into resiliency and redundancy for describing probable future status (Table 1).

STUDY SPECIES
Eastern black rails occupy relatively high ele vations along heavily vegetated wetlands, with moist soils or flooded up to a depth of 3 cm (Eddleman et al. 1994). Because they require dense vegetation cover that allows movement underneath the canopy, vegetation structure matters more than plant species composition (Eddleman et al. 1994, Watts 2016. Eastern black rails are found in a variety of salt, brackish, and freshwater wetland habitats that can be tidally or non-tidally influenced (e.g. Eddleman et al. 1988). During the breeding season, these birds require shallow pools that are 1−3 cm deep, which may be the most optimal for foraging and for chick-rearing, but nest flooding can be a significant impediment to reproductive success (Eddleman et al. 1994). Some elevational variability in the substrate is needed; eastern black rails require elevated refugia with dense cover to survive high water events due to the propensity of juvenile and adult eastern black rails to walk and run rather than fly and the precocial chicks' inability to fly (Eddleman et al. 1994, Flores & Eddleman 1995. The subspecies was historically widely distributed, with museum and other historical records reporting birds observed in interior wetlands throughout the mid-western and north-eastern USA (Eddleman et al. 1994, Watts 2016. However, reliable records have been scarce in those regions for the last 50 yr (Watts 2016). Most states in New England and the Appa la -chians consider the subspecies extirpated (Watts 2016). The remaining range includes coastal plains of the Mid-Atlantic, Southeast, and Southwest with some additional breeding sites in the Great Plains.
Species status assessments often subdivide the range of widely distributed species to analyze the species at ecologically relevant and statistically meaningful scales (Smith et al. 2018). These spatial subdivisions often, though not necessarily, follow dividing lines for representation units for the species and there is no specific guidance on how to subdivide a species' range into spatial subunits (Smith et al. 2018). We de termined that geographical and ecological differences across the range were aspects of the subspecies' biology that needed to be represented (Smith et al. 2018, USFWS 2018. Using a separate encounter-only point similarity analysis, we divided the range into 7 analysis units that served as representation units under the 3 Rs framework (USFWS 2018). These analysis units may also serve as large-scale redundancy components, i.e. if one unit is lost in the future, the subspecies will still be extant in the other units across the range. These spatial units are often referred to as 'analysis units' in species status assessments (Smith et al. 2018), and we used them as the largest scale for analyzing data and assessing current and future status of the eastern black rail. The 7 units in our assessment were 'New England,' 'Appalachians', 'Mid-Atlantic Coastal Plain', 'Southeast Coastal Plain', 'Southwest Coastal Plain', 'Great Plains', and the 'Central Lowlands' (Fig. 1). The subspecies is considered extirpated from the New England, Appa la chians, and Central Lowlands analysis units (USFWS 2018), so our assessment focused on the remaining 4 units.
We evaluated the following factors influencing eastern black rail population status identified in the SSA as potential threats: (1) habitat fragmentation and conversion resulting in the loss of wetland habitats across the range; (2) altered plant communities, primarily due to fire suppression, changing temperatures, sea level rise, and human modification; (3) altered hydrology resulting in impacts to soil moisture, surface water, sediment and nutrient transport, riparian and wetland vegetation communities, and land subsidence; (4) land management such as fire suppression, grazing, haying and mowing, and impoundments; (5) effects of climate change resulting in increased temperatures, decreased precipitation, increased frequency and intensity of wildfire and severe weather such as drought, flooding, or storms, increased sea level; (6) oil and chemical spills and environmental contaminants such as pesticides; (7) disease, specifically West Nile virus; (8) altered food webs resulting from invasive species (fire ants, feral pigs, mongoose, and exotic reptiles) introductions; and (9) human disturbance such as birdwatchers using excessive playback calls of eastern black rail vocalizations (Watts 2016, USWFS 2018. In the SSA we used data analyses and predictive modeling to evaluate some of the effects of these factors on current and future status.

Occupancy analysis to quantify species' needs and assess current redundancy and resiliency
In an SSA context, analysis of existing data can inform and quantify the ecological needs of a species and assess the current condition of the species across the landscape (Smith et al. 2018). Analyzing data to test associations with environmental covariates (e.g. land cover, weather data, etc.) can confirm and quantify hypotheses about the species' ecological needs (Williams et al. 2002). The analyses also have the potential to estimate metrics about the species in recent years (Smith et al. 2018). Demo graphic parameters like abundance, population trends, or related parameters such as occupancy probability can give information about the current condition of the species by estimating redundancy and resiliency if appropriately interpreted (Smith et al. 2018).
The USFWS requested and received data from state and partner agencies and organizations throughout the current and historical eastern black rail range (Fig. 1). These data varied in quality from direct surveys specifically targeting eastern black rails with callback surveys, to historical encounteronly records from museum collections and eBird data (http:// ebird. org/ content/ ebird/). For assessing ecological needs and current conditions, we used only high quality data from re peated presence−absence surveys across the range (MacKenzie et al. 2003, Conway 2011. These surveys were generally conducted according to the protocols of the Standardized North American Marsh Bird Monitoring Protocol (Conway 2011) and modified specifically for eastern black rail (i.e. targeting shallow-water wetlands and using eastern black rail call-backs). The surveys provided presence-absence data for use in occupancy modeling (MacKenzie et al. 2002(MacKenzie et al. , 2003. With these models, we can estimate the probability of presence at a site and test for eastern black rail associations with en vironmental covariates of interest (Table S1 in Supplement 2; all Supplements available at www. int-res. com/ articles/ suppl/ n043p209_ supp/; MacKenzie et al. 2002MacKenzie et al. , 2003. Concurrently, these analyses estimated the probability of detecting an animal, if it is present, because detection is imperfect and resulting occupancy estimates are in part a function of detection (MacKenzie et al. 2002, 2003, Ruiz-Gutiérrez & Zipkin 2011. Eastern black rails are particularly elusive and cryptic marsh land birds (e.g. Conway et al. 2004), and therefore accounting for detection probability can be especially important (Thompson 2013, Roach & Barrett 2015, Gilbert & Ferguson 2019, Tolliver et al. 2019.
We treated the estimated occupancy at survey points as a proxy for direct assessments of redundancy. For example, high occupancy probability estimates from 100s of survey sites would indicate that a population has high redundancy (i.e. many sites likely  into the future occupied), whereas low occupancy could be interpreted to mean low within analysis unit redundancy since fewer sites are likely occupied. Presenceabsence data and occupancy analyses are useful for understanding patterns and process of the animals' presence on the landscape and how much of the available habitat is used. However, there is no way to infer if a population is stable, increasing, or decreasing, from these data (MacKenzie et al. 2002, 2003, Popescu et al. 2012, Chandler et al. 2015, Zipkin et al. 2017. Unfortunately, there were no abundance data, trend data, or demographic data available for analysis in this SSA, and we therefore relied on these widespread but non-demographic and not-trendbased data. While occupancy analyses do not give direct insight into demographic processes of a population, they are increasingly used to inform conservation and management decisions (e.g. Fuller et al.

2016, Sutherland & Linden 2019).
We focused on surveys from each analysis unit that were repeated across years, so that we could use dynamic occupancy models to estimate site colonization and extinction over time (i.e. Southeast Coastal Plain, Southwest Coastal Plain, and the Great Plains; MacKenzie et al. 2003, our Fig. 1). This requires data collection to occur multiple times per season, for multiple seasons (MacKenzie et al. 2003, Kéry & Chandler 2016. A benefit of dynamic occupancy models over standard occupancy analyses is the ability to estimate extinction and colonization dynamics (MacKenzie et al. 2003). Even though we have minimal information on the demographics at a given site, we can still make inferences about site quality from the ex tinction probability estimates (Chandler et al. 2015, Kéry & Chandler 2016. We considered these 3 para meters combined to be useful for making inferences about population resiliency. For example, if the re sults indicated that initial occupancy was high, extinction probability was low, and colonization probability was high, the combination of those 3 parameters shows current high resiliency. That is, since a high proportion of the habitat is occupied and sites that are not occupied or might go extinct have a high probability of being recolonized, resiliency would be high. The occupancy and model selection analysis served 3 purposes: (1) to assess and quantify the subspecies' ecological needs by testing associations with environmental covariates such as presence of wetland land cover; (2) to evaluate the importance of some of the threats identified in the SSA, such as habitat loss with National Land Cover Database (NLCD) data or invasive species (fire ant data); and (3) to assess current redundancy and resiliency by estimating occupancy, extinction, and colonization probability.
We used data from South Carolina (2014−2017, 396 sites) and Florida (2016−2017, 64 sites) to represent the Southeast Coastal Plain analysis unit, from Texas in 2015 and 2016 (309 sites) to represent the Southwest Coastal Plain analysis unit, and data from Kansas to represent the Great Plains analysis unit (2005− 2008, 28 sites). We did not have dynamic occupancy data from the other extant analysis unit, the Mid-Atlantic, so we applied the results from the Southeast Coastal Plain to the Mid-Atlantic Coastal Plain. Those units where the subspecies is presumed to be extirpated (New England, Appalachians, and the Central Lowlands) were not included in the data analysis steps of the SSA. The parameters estimated in these analyses apply to an effective sampling area around each point for the surveys (MacKenzie et al. 2003, Kéry & Chandler 2016. The specific size of the effective sampling may vary among sites and analysis units, but the results effectively apply to 200− 250 m radius circles. We ran unit-specific analyses in order to estimate parameters for each analysis unit to use in the future projection models (see below, Section 4). We used the package 'unmarked' in R (Fiske & Chandler 2011, Kéry & Chandler 2016, R Core Team 2017 to analyze and compare models of dynamic occupancy and link model parameters to environmental covariates. Using expert opinion and prior literature as guidance, we developed 12 candidate models within each analysis unit comprised of varied combinations of precipitation, temperature, invasive species presence, and other covariates to explain variation in occupancy parameters using an analysis of Akaike's information criterion (AIC) ( Table S1). In these models, we included land cover type (NLCD data, Homer et al. 2015) as a potential covariate on initial occupancy, state (e.g. Florida or South Carolina), presence of fire ants, precipitation, mean spring temperature, annual temperature range, or year as potential covariates on extinction and colonization parameters, and year or state as potential covariates on detection probability.
The candidate model set was developed by discussion amongst members of the SSA authorship team (N. F. A., W. A. B., J. O. W., E. R., N. M. R., A. S., C. E. H., J. K. W., C. S., and N. M. R.). The covariates and model structures were chosen to address leading hypotheses about habitat associations for the subspecies and to evaluate the importance of perceived threats such as invasive fire ants (Korzukhin et al. 2001) or cattle grazing (Thornton 2010, Richmond et al. 2012. We endeavored to find covariates for the data analysis that represented possible threats to the subspecies identified in the SSA. However, for some threats, we did not have data that directly addressed a specific threat listed and instead we used surrogate or closely associated metrics.

Occupancy results, current status, and ecological needs assessment
Model selection results and parameter estimates varied by analysis unit ( Table 2). The results indicate low occupancy probabilities (ψ) in all analysis units. Occupancy was 0.25 (SE 0.048) in the Southwest Coastal Plain, 0.13 (0.075) in the Great Plains, and 0.09 (0.007) in the Southeast Coastal Plain. The statistical models appeared to have effectively converged as indicated by the SE values for most parameters; in other words, the coefficient of variation for each parameter (except for colonization [γ] in the Great Plains analysis unit) was <1. The results also indicated high site extinction probabilities (ε) (i.e. low site persistence) with an estimated extinction probability of 0.32 (SE 0.22) in the Great Plains and 0.61 (0.13) in the Southwest Coastal Plain (Table 2). In all analysis units, a null model (one with no covariates) or a simple, year-specific model was the best model or equally as good (Table 2). In the Southeast Coastal Plain there was evidence of year-specific extinction, with 2016 being as low as 0.001 and 2014 being as high as 0.57 (Table 2).
Inclusion of precipitation and mean temperature as covariates did not improve model performance for the Great Plains and Southeast Coastal Plain (Table 2). Model selection revealed weak evidence that wet season precipitation influenced occupancy dynamics in the Great Plains, and there was weak support that fire ants are determinants of seasonal occupancy in the Southeast Coastal Plain (Table S2). There was stronger evidence in the Southwest Coastal Plain for temperature playing a role in occupancy dynamics, as the best performing model had the range of temperatures as a covariate on colonization and extinction but the model weight was 0.54 and the null model fell within 2 AIC units (Table 2; Table S2).
For the Southeast Coastal Plain analysis unit, we analyzed data from Florida separately from that of South Carolina because there were fewer years (only 2) to analyze, much smaller sample sizes, and the years of the surveys did not match with those available for South Carolina. Occupancy probability was higher in Florida (0.17, SE 0.065) but was very low in South Carolina (0.04, SE 0.04) so we calculated a weighted average of the estimates from the 2 states, weighting the average by the sample size in each data set. Detection probability in the Southwest Coastal Plain and the Great Plains was ~0.25 (Table 2), meaning that when the birds are present at a site, there is a 0.25 probability of detecting them. In the Southeast Coastal Plain analysis unit, there was support for a year-specific detection probability (Table 2), and detection ranged from 0.09 to 0.53.

Model description
We used the results of the current dynamic occupancy analysis (see Section 3.2) to create a dynamic site-occupancy projection model for each of the analysis units. Occupancy simulation models have been used in conservation and management, especially with pond-breeding amphibians, although this model structure is uncommon in the avian literature (e.g. Martin et al. 2011, Green & Bailey 2015, Heard et al. 2013, Chandler et al. 2015. Generally, avian population models have more detailed demographic data on productivity and survival of individuals, allowing for the application of age-or stage-structured population viability models (Morris & Doak 2002 our case, however, we have data on site occupancy from multiple years across the subspecies' range, but lack specific demographic rates, so a dynamic site occupancy simulation model was appropriate (Williams et al. 2002).
Our model used a Markovian process to predict the number of sites occupied in the future based on the current number of sites occupied. Our modeling framework is similar to the SPOM model used by Risk et al. (2011) to model California black rail and Virginia rail populations (Moilanen 2004). The future number of sites occupied (N t +1 ) was a set of Bernoulli trials where the number of trials was the number of previously occupied sites (N t ) and the probability of success was the analysis unit-specific persistence probability estimated in the data analysis described above (i.e. 1 − ε i,t ; Moilanen 2004). In this simulation model, site is the same spatial unit (i.e. 200−250 m radius circles) as in the occupancy analysis described above. The process was modeled as: where the number of trials is N t in analysis unit i, and ε (extinction probability) is modeled as a stochastic, beta-distributed variable where the alpha and beta shape para meters were derived from the estimated mean and variance (see Section 3.1) (Morris & Doak 2002). The data analysis results supported yearspecific extinction probabilities in the analysis unit where we had the most survey points and the longest time series (Southeast Coastal Plain: ~400 survey points visited over 4 successive years). Therefore, we modeled a process that used a different base distribution for each year depending on whether it was a good year, an intermediate year, or a bad year for eastern black rails. Even though these year-specific dynamics were not supported in the dynamic occupancy analysis for the Southwest or the Great Plains, we applied this same dynamic to those analysis units in the simulation model because the data from the Southeast strong ly supported year-specific dynamics, and had more years and far more sites sampled than the other analysis units. Also, in preliminary simulations that did not include year-specific effects, populations went extinct incredibly quickly for each of those other regions (Southwest and Great Plains), a result that disagreed with observational data from recent years. The occupancy analysis in the Southeast Coastal Plain indicated that extinction probability was 0.57 in 2014, 0.49 in 2015, and 0.001 in 2016, so we used a function that first determined whether it was a good, intermediate, or bad year with ~0.33 probability of each, then drew the annual persistence probability from the appropriate distribution. For the Southwest Coastal Plain and the Great Plains analysis units, we did not have support for year-dependent models of occupancy. In those ana lysis units, we used the upper bound of the 95% confidence interval (CI) on our estimated extinction probability to represent good years, the mean to represent intermediate years, and the lower bound of the 95% CI to represent bad years. The initial number of sites occupied was specific for each analysis unit ( Table 2) and was the product of multiplying the total number of possible eastern black rail sites (U ) in an analysis unit by the estimated initial occupancy probability (ψ): N t=1 = U i × ψ i . The ψ parameter varied across simulation replicates and was drawn from a beta distribution where the alpha and beta shape parameters were derived from the estimated analysis unit specific mean and variance (Table 2; Morris & Doak 2002). We set initial U very high in each analysis unit so that the number of sites initially occupied would not be the primary driver of short-term analysis unit extinction.
We incorporated a colonization function into the model to allow sites that were not initially occupied or had previously gone extinct, to be colonized. We used a binomial function where the number of Bernoulli trials was the total number of eastern black rail sites available that year, and the probability of success was the estimated colonization probability in that analysis unit. Colonization probability (γ) was modeled as a temporally varying parameter and drawn annually from a beta distribution where the alpha and beta shape parameters were derived from the estimated mean and variance (Morris & Doak 2002). Therefore, the full formulation on the N i,t +1 model was as follows: The primary output metric for our model was the mean proportion of sites still occupied at each time step into the future (± 95% CI), which we calculated by dividing the number of sites occupied at time t by the initial number of sites occupied in each replicate. We conducted all simulation modeling in R, and replicated the scenarios 5000 times each to capture variability in each scenario (see Section 4.2; R Core Team 2017; Supplement 3).

Simulation scenarios
We incorporated functions to account for habitat quality and possible habitat loss over time to mimic some of the perceived threats to the subspecies. Even though our covariate analyses found no predictive effect of any covariates in the data analysis on extinction or colonization parameters, SSA documents, published literature, and experts still believed that there are habitat loss and management-related threats present for eastern black rails (e.g. Richmond et al. 2012). To that end, we incorporated habitat loss and quality functions into our simulation model. The habitat loss function was a simple reduction in U at each time step in the simulation by a randomly drawn percentage (a beta distributed random variable) that was specified under different simulation scenarios to represent habitat loss due to development (urbanization) or sea-level rise.
For alternative scenarios, we took the mean predicted rate of decline from published sources and increased or decreased the rates by ~10−20% over or under the estimated rate of habitat loss, i.e. to account for uncertainty in the estimated habitat loss rates. We used the change in 'developed' land cover from NLCD data to derive an annual rate of change in each analysis unit (SLEUTH 2014, Homer et al. 2015, and we used climate change (Young et al. 2017) and NOAA sea-level rise predictions to estimate probable coastal marsh habitat loss rates (Sweet et al. 2017). The rates were all negative be cause even though some new habitat will be created as sea levels rise, the overall expectation is habitat loss unless habitat conservation actions are implemented. In the Great Plains, we used ground water loss rates instead of sea-level rise data to represent permanent non-urbanization habitat loss in the analysis unit (Supplement 2).
Habitat loss due to development, sea-level rise, fire suppression, and other causes were listed as major concerns in the SSA, so these habitat loss scenarios are intended to represent all of those potential threats. We used the sea-level rise predictions from NOAA and the NLCD wetland habitat loss metrics be cause they fit well with our modeling framework. Other metrics or predictions may have been useful or appropriate but because of courtordered decision deadlines, we did not have the opportunity to ex plore other approaches to modeling or metrics for climate change factors. The strength of these data was their immediate availability but also their uniformity across the subspecies' range with commonality of scale and precision. We did not explore varied conservation scenarios, other than the one habitat conservation scenario, because this work was done in support of a listing decision and those types of scenarios are more appropriate for a recovery planning effort if the subspecies receives ESA protections.
We also incorporated a function to allow for 'poor habitat condition' related to land management, fire, or agricultural practices not compatible with eastern black rail ecological needs. Using available data, we calculated the mean annual proportion of the land exposed to potentially negative cattle, fire, haying, and water management practices in each ana lysis unit (Thornton 2010, Supplement 2). We implemented a function to reduce the persistence probabilities at the proportion of sites exposed to those practices. The realized extinction probabilities were calculated as a weighted average of the sites ex posed to poor land management and sites not ex posed, weighted by the proportions randomly generated each year. The annual extinction probability was therefore modeled as: ( 3) where ε R is the realized extinction probability in analysis unit i at time t, P(ph) is the proportion of the habitat in poor condition, ε b is the baseline extinction probability, and ph is the poor habitat effect (i.e. the increase in site extinction probability caused by poor habitat quality). The P(ph) value was drawn annually from a beta distribution that was based on the mean and variation estimated from available data (Table S1), and the mean was increased or decreased to represent differing land management scenarios (Table 3). We did not have data to inform the magnitude of the ph factor, so we input a mean of 0.05 in crease in extinction probability (which varied annually and was drawn from a beta distribution) and tested the sensitivity of model predictions to changes in the mean ph value.
We had 2 primary concerns when developing scenarios to present in the SSA: (1) statistical and input sensitivity scenarios and (2) threat-focused scenarios. These components were combined in all scenarios. We identified the threats to project in the model (habitat loss and habitat management), found data or resources to guide the base level of each threat (e.g. NOAA data for sea-level rise), and then created high, medium, or low versions of those threats to represent uncertainty in the future trajectory of the threats (e.g. predicted sea-level rise is over-or underestimated). This approach allowed us to evaluate the various threats to the subspecies but also to explore the possibility that threats may be lessened or increased over time, or were under-or overestimated. We also limited the number of scenarios in the SSA to reduce confusion for decision makers in case there were divergent predictions among scenarios. Decision makers can experience decision fatigue or information overload (e.g. Lenz & Lyles 1989). We settled on 5 primary threat scenarios with varying levels of sea-level rise and habitat management and 1 conservation scenario to evaluate whether conservation effort would help the subspecies (Table 3). For the alternative scenarios, we took the mean predicted rate of habitat decline from published sources (Supplement 2) and increased or decreased the rates by ~10−20% in case those predictions over-or underestimated the rate of habitat loss, i.e. to ac count for uncertainty in the estimated habitat loss rates. We also used the model to explore conservation scenarios in each analysis unit to predict how the subspecies would fare under significant habitat conservation and restoration efforts (Table 3).

Projection modeling results
The model predicted high probability of complete extinction by the year 2100 for all analysis units under all of the primary scenarios (Table 3). The Southeast Coastal Plain and the Mid-Atlantic Coastal Plain analysis units had the longest predicted time to complete analysis unit extinction, between 35 and 50 yr, depending on the scenario -with the exception of the habitat conservation scenario (Figs. 2 & 3). The Great Plains had the shortest time to complete analysis unit ex tinction, be tween 15 and 25 yr, depending on the scenario, and the Southwest Coastal Plain analysis unit was in be tween (Figs. 2 & 3). The simulations exhibited high variability across the 5000 replicates (Figs. 2 & 3), but generally, after the first ~25 yr, all primary scenarios exhibited consistent downward trends in the proportion of sites remaining occupied across most replicates.
Most of the predicted occupancy declines were driven by habitat loss rates input into each scenario. The model results exhibited little sensitivity to changes in the habitat quality components in the simulations (i.e. the P(ph) and the ph components) for the range of values that we explored. When in creasing the mean habitat loss rate from 0.01 an nually to 0.02 annually and holding the mean proportion of the range in poor habitat condition constant at 0.02, the proportion of sites remaining occupied at 100 yr de clined from 0.32 to 0.12. Whereas when holding the mean habitat loss rate constant at 0.01 and increasing the mean proportion of habitat in poor condition from 0.02 to 0.04, the proportion of sites remaining occupied at 100 yr declined from 0.32 to 0.30. Under the conservation scenario, our model predicts that habitat loss rates of 0.005, or 0.5% annually, would likely result in fairly stable populations in the coastal analysis units (> 60% of sites still occupied in 50 yr), but still predicts large declines in the Great Plains analysis unit in the next 25 yr (Table 3, Fig. 2

DISCUSSION
Occupancy modeling offers a potentially useful frame work for evaluating the current status of a species in terms of redundancy and resiliency. Our dyna mic occupancy results enabled decision makers to conclude that current redundancy and resiliency were low, given low occupancy and colonization probability for the subspecies across the range and high site extinction probability. The occupancy probability can give insight to decision makers on the proportion of sites that are likely occupied and, with a 218 Fig. 2. Median proportion of habitat remaining (solid line) and the 2.5 and 97.5 percentile (CI; dashed lines) from (a) Scenario 1 and (b) Scenario 6 of the 6 scenarios (see Table 3) used to predict future status in each of 4 analysis units with extant eastern black rail populations sense of how much habitat is available, redundancy of the populations can be inferred. If the site selection process to gather the data was appropriate, this metric can give decision makers a sense of how much of the available habitat is currently utilized. Occupancy probability for eastern black rails throughout their range was ~0.25 or less, suggesting that 75% or more of available habitat is unoccupied. Further suggesting that population redundancy in each analysis unit was limited. We assert that dynamic occupancy modeling offers the capacity to assess resiliency of a population and species by combining the occupancy probability with the colonization and extinction probabilities. Together these 3 parameters allow decision makers to make inferences about the distribution of a species and its current stability (Sutherland et al. 2014, Chandler et al. 2015. Eastern black rail colonization rates were low in all analysis units, suggesting that the ability of the populations, and therefore the subspecies, to withstand catastrophic stochastic events is limited. Populations with low occupancy or high site extinction probability must also have high annual colonization probability to be a resilient population. This was not the case for the eastern black rail; the initial occupancy and colonization rates were low (less than 25%), coupled with high site extinction probability in 2 analysis units, and year-dependent extinction probability in the third. Our results fit well with other studies analyzing the same or similar datasets. For example, Tolliver et al. (2019) estimated mean colonization rates of 0.09−0.17 in Texas, and our mean colonization rate was 0.05−0.22 in the Southwest Coastal Plain. Our occupancy was 0.15− 0.34 and was also similar to mean occupancy estimated by Tolliver et al. (2019;range: 0.22−0.30). Even though our sample sizes were larger and we used more years of data, our results exhibited more variability than those of Tolliver et al. (2019), probably because our data set had more spatial and temporal variability.
Eastern black rail detection probability was low across the range (~0.25 or less on average), which introduces uncertainty into the other model parameters by increasing variance estimates. This also af fects our ability to investigate and estimate relationships with environmental covariates because of large estimated variance on occupancy, extinction, and colonization parameters (Williams et al. 2002). Our estimates of detection were similar to other studies in the same systems. For example, our Southwest Coastal Plain detec tion was 0.15−0.31, whereas the detection rate re ported by Tolliver et al. (2019) was approximately 0.19, depending on the weather conditions during sur veys. Based on other subspecies of rails, the eastern black rail detection rates found in our study were similarly low, and in fact we had greater eastern black rail detection than in other studies before purported declines. For example, methodologically similar studies found that Yuma clapper rails were 19% detectable (Conway et al. 1993), eastern black rails in Florida were 20% detectable (Legare et al. 1999), and Virginia rails were 22% detectable (Glahn 1974). The estimated occupancy, extinction, and colonization probabilities account for detection probability, which is another benefit of using occupancy and dynamic occupancy analyses to assess redundancy and resiliency (MacKenzie et al. 2003, Kéry & Chandler 2016.
The dynamic occupancy modeling was very useful for assessing the current condition of the populations and subspecies on the landscape. We used the informal combination of 3 parameters (initial occupancy, colonization, and extinction) to convey information about redundancy and resiliency to decision makers. Our results showing high site extinction probability combined with low site colonization and initial occupancy, indicate low current resiliency for the eastern black rail. This informal interpretation of the dynamic occupancy analysis could allow for inconsistent as sess ments across analysis units or species. For example, decision makers may have different definitions of what level of site extinction probability is too high, or what combination of extinction and colonization probability is too low. It may be useful to develop a more impartially derived metric, such as the log(γ/ε), where values above 0.0 indicate stable dyna mics and values below 0.0 would likely have unstable dynamics (i.e. extinction is too high and not 219 Fig. 3. Mean year of analysis unit extinction of eastern black rails based on occupancy projection modeling from 2018 through the year 2100 balanced by colonization), thus ensuring more consistency among SSAs and analytical interpretation. Furthermore, communicating the results of these quantitative analyses to decision makers in accessible and understandable ways is imperative for decision making. The ecological needs assessment, which we attempted to support through covariate analysis of occupancy data, did not work as we had hoped. The null model or a year-specific model with no ecological covariates received some or significant support in each analysis unit. We suspect this resulted from a scale mismatch between the spatial scale of the occupancy data (200−250 m radius circles) and the covariate data (e.g. remotely sensed land cover data, or nearby weather station data). Had there been more time to complete the analysis before the listing decision was required, we might have been able to seek out more useful ecological covariates. If important covariates had been supported in the occupancy analysis, we could have drawn inferences on the conditions that increase or decrease occupancy, extinction, or colonization probability. However, there was only, at best, weak support for any covariates in the model selection analyses, as in other studies (i.e. Tolliver et al. 2019). We instead used expert judgment or literature resources to assess the ecological needs of the subspecies and to incorporate those probable but uncertain effects into the simulation model (Smith et al. 2018).
With the results of the dynamic occupancy modeling, we designed a model to predict the future proportion of sites remaining occupied (Moilanen 2004, Risk et al. 2011, Sutherland et al. 2014). This output metric was particularly useful in our case because we do not know how many sites are currently occupied or how many sites are available for occupation. It is important for decision support models, in an SSA context especially, to incorporate all of the important sources of uncertainty so that decision makers can see the range of expected outcomes in the future (McGowan et al. 2011, 2017, Heard et al. 2013, Howell et al. 2020). Thus our model incorporated uncertainty in the initial number of sites occupied and then assessed future conditions with respect to the uncertain starting conditions. This relative metric of future status avoids making specific predictions of future site occupancy, i.e. number of sites occupied each year, but still provides a prediction about the future state of the population and its trajectory.
Our model incorporated initial state uncertainty (as discussed above), stochastic variation from statistical distributions to represent environmental variability, and we also used multiple scenarios to incorporate uncertainty in the future condition of important environmental variables, such as sea-level rise, land management, or habitat conservation. Under all scenarios, except the habitat preservation and conservation scenario, our model predicted significant and rapid decline in each analysis unit in the future, suggesting low resiliency and low redundancy of eastern black rails.
Our model is intended to be useful for informing decisions about eastern black rails; however, all models are limited in their utility and inference capabilities. One major potential limitation of our modeling is the data we used to parameterize these simulations (Williams et al. 2002, Howell et al. 2020. The projection models are entirely dependent on the data used to estimate occupancy and extinction dynamics. If the survey sites that were sampled for eastern black rails were not in optimal eastern black rail habitat we would likely underestimate initial occupancy and colonization probability and overestimate extinction probability (Williams et al. 2002, MacKenzie et al. 2003. Those biases would result in overestimating extinction risk in the future for each analysis unit and the subspecies with the simulation model (Williams et al. 2002). However, the USFWS is tasked with using the best available science and data to inform decisions, and the marsh bird monitoring protocol of Conway (2011) is widely accepted as the best method for surveying marsh bird populations. The data we used modified that survey design to specifically target eastern black rail habitats and populations (with the exception of the Kansas dataset). Our results correspond to similar studies of eastern black rails across the analysis units or have higher detection and occupancy probability, suggesting that our data were robust enough for the analysis and projection model carried out.
Our overall approach took advantage of existing data collected by the States and other management partners within a widely accepted monitoring framework and therefore did not rely on expert elicitation or large-scale extrapolations from highly detailed but small-scale studies and data sets (e.g. McGowan et al. 2017, Voorhies et al. 2019. Other analyses and predictive approaches may have also addressed the information needs and supported the decision process. However, incorporating data from so many management partners and researchers made the process inclusive and collaborative, which can be important for potentially controversial decision processes (e.g. McGowan et al. 2015). The data that we used represent the best available for supporting the SSA and the decision-making process, but those data varied in complexity, representing many researchers' efforts over multiple years. For example, data from the Southeast Coastal Plain analysis unit had large sample sizes, good spatial coverage, and adequate time series for estimating key model parameters, whereas data from the other analysis units were temporally (Southwest Coastal Plain) and spatially (Great Plains) limited. Even so, the parameter estimates appear to be effectively estimated (i.e. the models adequately converged) because the variance estimates are not excessively large.
Our simulation model incorporated significant variability and uncertainty into the projections, which leads to variability and uncertainty in model output and predictions (Williams et al. 2002, McGowan et al. 2011). The predictive model was focused on estimating redundancy and resiliency for the populations and subspecies in the future. If the proportion of sites remaining occupied remained stable or increased over time, the population would exhibit resiliency in the face of stochasticity. Assessing resiliency across the number of occupied sites within a representative analysis unit and across representative analysis units gives insight on future redundancy of the populations and species (Smith et al. 2018). This was not the case for the eastern black rail. The SSA found evidence of historical declines, current low redundancy and resiliency across the range, and predicted future resiliency and redundancy to decline. At the end of the analysis and projection modeling work, the USFWS concluded that, based on low current occupancy, high site extinction probability, low colonization probability, and a high probability of extinction in the future across the range in 50 yr time, listing the subspecies as 'threatened' under the ESA was warranted. The data analysis and projection modeling herein directly informed the decision-making process.