Spatial and temporal predictions of whooping crane (Grus americana) habitat along the US Gulf Coast

The challenge of conserving viable habitat while simultaneously predicting how land cover may geographically shift with future climate change has put pressure on ecologists and policy‐makers to develop near‐term (several years to a decade) ecological and geospatial predictions. This is particularly relevant for endangered species as ranges adjust to track a changing climate. The whooping crane is vulnerable to these changes, as the overwintering habitat of a small population is susceptible to climate impacts. This study mapped the historical spatial transformation of crane habitat in and around the Aransas National Wildlife Refuge. A time series of ecological niche models was developed to determine the biotic and abiotic factors correlated with crane presence and track how importance has changed over time. The results from the multitemporal models were used to predict areas along the US Gulf Coast where additional unoccupied habitat may be located for crane population expansion and model how the areas may degrade or change as sea levels rise through future climate change scenarios. Findings indicate that the percentage of emergent herbaceous wetland and water are the most important variables influencing crane presence. Sea level rise analysis indicates that potential habitat throughout the Texas–Louisiana Gulf Coast will be impacted considerably by climate change. The lack of large, continuous blocks of usable land cover could limit population expansion and future recovery efforts. However, the findings can help facilitate winter range expansion to accommodate the growing population by identifying additional areas to protect that could be used by the current wild population or experimental populations.


| INTRODUCTION
The world is at a tipping point in terms of biodiversity, with more than a million species currently under threat of extinction (Barnosky et al., 2012;Steffen et al., 2018). Assessing which species are at greatest risk from habitat loss remains a conservation priority, since loss of both species and habitats can degrade ecosystem services (Mooney et al., 2009), and species abundance can reflect habitat availability. A changing climate complicates conservation efforts by shifting habitat boundaries while simultaneously pushing aspects of the Earth system potentially past the point of reversible change (Leclère et al., 2020). As an example, the melting of the Greenland and West Antarctica ice sheets from global warming could lead to sea level rise up to 3 m (Bamber et al., 2009;Lhermitte et al., 2020), which would have detrimental consequences for coastal-dwelling species, such as wintering whooping cranes.
The whooping crane (Grus americana) is particularly vulnerable to the combined impacts of habitat loss and climate change because of critical overwintering habitat requirements of the wild population. Overhunting coupled with habitat loss and degradation from humandriven disturbances led to a rapid decline of the population, and by the 1930s, the whooping crane was at the brink of extinction. The only wild population (~20 individuals in 1937) winters in a small area of coastal estuarine habitat along the US Texas Gulf coast and migrates annually to breeding grounds in northern Canada (French Jr. et al., 2019). Today, this group, known as the Aransas-Wood Buffalo population (AWBP), is estimated to be composed of more than 500 birds (95% CI = 342.6-678.0) (USFWS, 2020) and is the only remaining viable population of wild whooping cranes.
The vulnerability of whooping cranes to human activities along the Texas Gulf coast was recognized long before there was convincing evidence for climate change. The Aransas National Wildlife Refuge (ANWR) was protected in 1937 to conserve AWBP habitat (Evans & Waring, 1993) and was declared critical habitat in 1978 (USFWS, 1978). Conservation and land management efforts have contributed to the increasing crane population with a long-term average growth rate of 4.4% (n = 80; 95% CI = 1.85-6.96%), although the population number has remained apparently stable over the past few years (USFWS, 2020). These findings indicate improvements in whooping crane viability over the past 80 years, but population numbers are still well below the 1000 individuals and 250 breeding pairs required to down-list the species from endangered to threatened (CWS & USFWS, 2007).
As the population has grown, the cranes have expanded their wintering area beyond ANWR (Harrell & Bidwell, 2019;Stehn & Prieto, 2010), but expansion has been limited to areas just outside refuge boundaries, and documented sightings in the last 60 years beyond areas immediately adjacent to the refuge have been rare. More recently, pairs have been documented using areas outside the refuge, but the most concentrated area of overwintering remains within ANWR (USGS, 2020). The vast majority of the wintering population primarily occurs in areas of low-lying coastal marsh vegetation, but limited areas of this vegetation community and reduced habitat quality may be preventing more expansive population growth (CWS & USFWS, 2007;Lewis, 1995;Metzger et al., 2020). Tidal flats are also an integral component of crane habitat and are used in higher proportion than available (Golden, 2020). Thus, there is widespread interest in identifying which vegetation communities and land cover types whooping cranes are using so that unoccupied available habitat can be protected (Metzger et al., 2020). Predicting where additional preferred habitat may be located for the future is also critical to proactively foster population growth.
A fundamental step in meeting the goals of down-listing the species is to identify the biotic and abiotic conditions that characterize crane winter habitat along the Texas Gulf Coast by building on past research (Chavez-Ramirez, 1996;Chavez-Ramirez et al., 1996;Chavez-Ramirez & Slack, 1999;Davis, 2019;Metzger et al., 2020). The historical overwintering habitat of the AWBP included the entire Texas coast and considerable portions of the Louisiana coast . Understanding how land management strategies taking place within the refuge (i.e., prescribed burning) and habitat availability along the coast influence crane use is necessary to identify whether other areas exist where unoccupied habitat can be conserved or managed to allow population expansion. It is also crucial to predict how these areas may change or disappear with rising sea levels from climate change, given the precarious location of the current realized wintering niche along the coast.
The objectives of this study were threefold. Using the historical archive of whooping crane observations in ANWR, we first mapped how the historic spatial distribution of cranes has changed over time in and around ANWR as the AWBP has grown. This exploratory step highlights trends in the use of space by the cranes over time that can inform future habitat needs. Second, we developed a time series of ecological niche models to (1) determine the biotic and abiotic factors correlated with crane presence to infer habitat availability, (2) determine how land management practices within the refuge contribute to habitat selection, and (3) track how land cover and variable importance has changed over time. Lastly, we used the results from the multitemporal ecological niche models (series of models over time) to spatially predict areas along the US Gulf Coast (Texas and Louisiana) where additional unoccupied habitat may be located for crane expansion and model how these areas may degrade or change as sea levels rise using future climate change scenarios.

| The Aransas National Wildlife Refuge
ANWR is located along the southeastern Gulf Coast approximately 200 km south of Houston, Texas, USA. The dominant physical features in the region are coastal bays and inlets, along with a saltmarsh vegetation community that supports abundant terrestrial and marine life and attracts many bird species. Beyond the marshlands, the terrain is relatively flat with a mosaic of oak woodlands and low sandy prairie openings. The climate is subtropical with hot, humid summers and mild winters. Much of the area has experienced substantial land use and land cover changes, particularly over the past two decades, mainly in the form of salt marsh loss and woody plant encroachment from oaks and mangroves (Armitage et al., 2015;Brown & Archer, 1999).
Today, ANWR consists of five units totaling 46,810 ha ( Figure 1). The main unit, Aransas, comprises 19,126 ha and was established in 1937. It is surrounded on three sides by Aransas Bay, part of the Intercoastal Waterway separating the mainland from barrier islands. The Matagorda Island unit is part of a chain of barrier islands extending along the coast, with the unit comprising 22,939 ha of the island. The Tatton unit (3063 ha) protects a remnant of low upland dark soil coastal prairie. The Myrtle Foester-Whitmire unit (1392 ha) was purchased in 1993 for waterfowl and shorebird protection and native prairie restoration. The Lamar unit is the smallest, comprising 396 ha of salt marsh habitat and native coastal woodlands.
The refuge is a remnant of a coastal prairie firedependent ecosystem and is composed of fire-adapted vegetation communities that have historically burned frequently due to both natural ignitions and intentional burning by Native Americans and early settlers (Hanselka, 1980;Lynch, 1941;Stambaugh et al., 2014). Since the 1980s, ANWR has taken an active approach to managing vegetation communities used by whooping cranes via prescribed burning, which attempts to maintain woody vegetation at F I G U R E 1 Location of the Aransas National Wildlife Refuge complex on the southeastern Texas Gulf Coast, USA early successional stages and vegetation height below 1.4 m to foster visibility, aid in predator detection, and promote foraging and feeding in the uplands (Armbruster, 1990;Chavez-Ramirez, 1996;Chavez-Ramirez & Wehtje, 2012;Venne & Frederick, 2013). Since cranes use unobstructed areas for roosting and feeding (Urbanek & Lewis, 2020), burning may benefit these activities. While the majority of the whooping crane diet consists of blue crabs and clams (Hunt & Slack, 1989;Pugesek et al., 2013), upland burning may provide cranes access to supplemental food when preferred items are scarce Hunt & Slack, 1989). Since 1985, ANWR has maintained a fire return interval between 1.3 and 28 years (Golden, 2020), but there have been few efforts to quantify the impact of prescribed burning on crane habitat (Chavez-Ramirez, 1996;Hunt, 1987).

| DATA AND ANALYSES
Exploratory spatial data analysis (ESDA) is key for generating a baseline understanding of which areas of ANWR cranes have been using and how their space use has evolved over time. We address this gap through two spatial statistical analyses (mean center and kernel density) interpreted in the context of land conservation efforts. Next, we develop a time series of ecological niche models using MaxEnt (Phillips et al., 2006) to uncover the environmental indicators most related to crane habitat use and chart how those variables have shifted in importance over time. Using the time series model results, we forecasted areas along the Gulf Coast where potentially unoccupied habitat might be located for cranes to expand their winter range to accommodate the increased population size required for down-listing the species. Crane survey methods (described below) changed in 2012, and prescribed burn records are only available back into the 1980s, so our analysis spans the 20-year period 1990 through 2009.

| Crane surveys
From 1950 to 2010, the United States Fish and Wildlife Service (USFWS) conducted aerial surveys during the overwintering period from October to May to monitor crane population size and document their locations. Surveys were conducted biweekly (as permitting) by a two-man crew flying over the study area (Stehn & Taylor, 2008). Flights were conducted at a low enough altitude (61 m a.s.l.) and a slow enough average speed (167 km/h) to permit identification of cranes by leg bands and size (Stehn & Taylor, 2008). Observers marked crane locations on a physical map and recorded attributes such as the number of cranes/juveniles present. Flight numbers varied by season (Table 1), with a median of 21 flights conducted annually. The aerial survey methodology has been scrutinized for the lack of standardized sampling protocols and other effects that result in imperfect detectability and unquantifiable precision such as observer fatigue, blind spots on the plane, weather limitations, and observer confusion with other white birds that overwinter at Aransas (Strobel et al., 2012;Strobel & Butler, 2014). However, it is the only consistently captured, multi-decadal, historical survey of the crane population in ANWR, so despite the drawbacks, there is inherent value in the dataset as it offers a rich archive for describing and modeling the general spatial distribution of the population.
Approximately 21,000 locational observations of crane groups, including 48,000 cranes, were made during the 384 surveys that occurred during the 20-year study period. Given that the current population size is around 500 birds, the same individuals were observed multiple times through the repeated counts. Since the focus here is on analyzing geographic distribution, we randomly selected 500 use locations in each survey year for analysis to ensure an equal number of samples were used for the spatial distribution analyses described below.

| Spatial distribution analysis
Mean center analysis quantifies central tendency and is the geographic (x, y) center of a set of features. When computed over time, it can provide a spatially relevant statistical measure of population movement. For the AWBP, which returns to the same spot year after year, subtle but sustained shifts in the geographic distribution can provide insight into where new territories are being established, which is a key component of future conservation efforts. Since each location observation can include multiple individuals, we computed weighted mean center as: where x i and y i are the coordinates for observation i, w i is the weight for observation i (the number of birds tallied for that observation), and n is the total number of observations. We used the subsample of 500 observations to reduce outlier effects.
Kernel density estimation is a nonparametric method for estimating the probability density function of a variable and smoothing data trends in order to make inferences about the population from a sample (Wiegand & Moloney, 2014).
We used kernel density to map "hotspots" from the sampled crane presence locations. These density surfaces are particularly useful when working with ecological observations because they can provide an estimate of the probability of finding the object of interest in an area. When interpreted in conjunction with the mean center, they can indicate how the distribution is spreading over time. Kernel density is estimated as (Silverman, 1986) where K is the probability of occurrence, X i is an independent sample drawn from an unknown density (crane locations), n is the number of independent points, and h is bandwidth. The bandwidth is a key consideration because it exhibits a strong influence on estimates. We determined kernel bandwidth based on typical territory sizes (Stehn & Prieto, 2010), and using the upper range of average territory size (314 ha), we computed the approximate territory radius as 1 km. We held this bandwidth constant across all years of analysis so surfaces could be compared. Although there has been a shift to slightly smaller territories over time, the average has not changed significantly (Stehn & Prieto, 2010); therefore, the bandwidth encompassed the range of changes in territory sizes over time.

| Ecological niche modeling
Ecological niche models (a.k.a. species distribution models [SDMs]), are used to model habitat availability and predict species occurrence (Elith et al., 2011;Elith & Leathwick, 2009;Franklin, 2010). MaxEnt (Phillips et al., 2006) is one of the most popular and highest performing software packages for presence-only modeling (Elith et al., 2011;Merow et al., 2013;Phillips & Dudík, 2008). MaxEnt seeks to maximize the entropy (i.e., uniformity) of the geographic distribution of a species constrained by environmental covariates, and, when conceptualized as an inhomogeneous Poisson process, the cumulative log-log output can be interpreted as the probability of presence (Franklin, 2010;Phillips et al., 2017). The model parameterization is specified

| Environmental covariates
A suite of seven predictor variables (proportion emergent herbaceous wetlands, proportion water, distance from protected areas, normalized differenced vegetation index, distance to woody wetlands, cumulative burn index, and distance to developed areas) was constructed for each of the 20 years of analysis (1990-2009) based on observational and empirical evidence of crane habitat preferences from the literature. Variables were resampled when necessary to 30 m spatial resolution for modeling purposes. For all variables, water areas greater than 2 m deep were masked using bathymetry data obtained from the Texas Natural Resources Information System, since cranes prefer shallow water, with average depths between 10 and 92 cm (Austin & Richert, 2005;Kang & King, 2014). The bathymetry dataset represents a single snapshot from 2013, so we made sure this mask did not preclude any observations from the model as it is possible certain depths at the time of observation were different due to tide and sea level changes.

| Normalized difference vegetation index
The normalized difference vegetation index (NDVI) is a measure of vegetation greenness that provides information on density and vegetation type using red and near infrared (NIR) spectral bands. Since cranes prefer shallow water and wetlands for foraging and roosting (Austin & Richert, 2005) and avoid tall, dense vegetation (Armbruster, 1990), we expected crane presence to be associated with low to moderate levels of NDVI but not high values that would signal denser vegetation (Jensen, 2005). We derived NDVI for each survey year from Landsat 4-5 TM and 7 ETM+ using Google Earth Engine. We computed NDVI using the median reflectance values for all pixels from cloud-free images captured during the overwintering period (October 15-May 15) of each year.

| Land cover variables
Cranes generally prefer areas distant from forested wetlands and human developments (Austin & Richert, 2005;Pearse et al., 2017). Using the National Land Cover Database (NLCD), which spans seven comparative epochs (2001,2004,2006,2008,2011,2013, and 2016), we isolated developed areas, woody wetlands, emergent herbaceous wetlands, and water. We did not include agricultural lands as they are not essential for this population of overwintering cranes, although they may be beneficial in years of low crab availability (Armbruster, 1990;Golden, 2020). For each epoch, we computed the proportion of water and emergent herbaceous wetlands within a 1 km radius from each cell using a spatial convolution filter. We computed the distance from each cell in the study area to the closest developed land cover and woody wetlands using Euclidean distance.

| Protected areas
Evidence suggests that cranes may select protected regions as they have been observed in protected areas disjointed from ANWR without using areas in between (Golden, 2020). Using the USGS database of marine and terrestrial protected areas, we developed a variable representing the distance from every cell in the study area to the closest protected area. The protected areas database includes approximately 95% of federal lands and waters and 60% of state, regional and local parks, and other preserved lands. We removed areas under 2.02 ha to eliminate sliver polygons and smaller parks, considering the species' large wintering territory size. In cases where data were missing, we researched dates when the protected area was established and created separate distance rasters for each year of analysis.

| Visibility based on prescribed burning
To examine crane response to prescribed burning at ANWR, we built a digital atlas of burn activity on the refuge based on the historical record and a satellite-based burn index (Rogan & Yool, 2001). We then developed a spatially continuous variable from that atlas, which summarizes the spatial and temporal density of burns.
To build the digital atlas, we obtained a burn database from the USFWS that included the management unit where the burn occurred (but not precise coordinates) and the burn date. We mapped burn footprints using the differenced normalized burn ratio (dNBR), (ΔNBR = Prefire NBR À Postfire NBR), with Landsat images (pre-and post-fire images for every expected burn). Pre-fire images could extend up to 3 months prior to the burn, while post-fire images were limited to 1 month afterward since vegetation regrowth can be rapid, preventing accurate identification of burns beyond that time. The normalized burn ratio (NBR) index incorporates near infrared (NIR) and shortwave infrared (SWIR) reflectance. Healthy vegetation strongly reflects NIR light, but fire destroys plant cell structure, reducing or eliminating chlorophyll production and decreasing the reflectance of these wavelengths (Jensen, 2005). Similarly, fire-altered landscapes absorb less SWIR light due to the decrease of water content in the vegetation (Rogan & Yool, 2001). Together, these attributes make NIR and SWIR bands well-suited for mapping burns.
ΔNBR differences the NBR from the pre-and post-fire images and ranges [À2, 2], with values from 0.1 to~1.35 indicative of burned areas (Key & Benson, 2005). Using these guidelines alongside visual validation, pixels with ΔNBR values from 0.15 to 1.50 were classified as burned. This range was determined to be suitable for delineating the extent of burned areas while also mitigating the effects of clouds. Contiguous burned pixels were vectorized and assigned the burn date. Processing was completed using ENVI IDL programming software, and accuracy was verified by ANWR personnel with first-hand knowledge of the burns (Harrell, personal communication). Burns were accumulated from the 5 years preceding and including the survey year using an inverse weighting system where burns conducted within the survey year were weighted 5, burns from the prior year weighted 4, etc. This allowed burns during the survey year to contribute most heavily to the model while also allowing for lagging vegetation effects from prior years' burns as woody vegetation takes several years to reestablish (Ruthven III et al., 2003). A kernel density estimator was applied to the aggregated, weighted burn polygons to estimate the probability density function of burning effects. The kernel bandwidth was 1 km as above. The result was a continuous surface of probabilistic effects of burn activities on ANWR.
It should be noted that only a single burn inventory was conducted for each year, which means any crane observed during that year will be related to the burn regardless of when the burn occurred. This decision only impacts a limited number of observations during the survey year and does not affect the preceding 4 years of burns. Therefore, the impact on the magnitude of the visual acuity variable is marginal.

| MaxEnt parameterization
Models for each year were run allowing for a maximum of 500 iterations with a convergence threshold of 0.00001, and 10,000 background points selected at random from the environmental layers matching the year of crane observations. From the 500 sampled observations, 25% were set aside for statistical validation. Model performance was evaluated using area under the receiver operating characteristic curve (AUC). AUC ranges [0, 1], values closer to one represent a more accurate model that can be interpreted as being equal to the probability that the model will rank a randomly selected presence location as being more suitable than a randomly selected background location (Merow et al., 2013). AUC is widely used as a measure of SDM performance and is preferable to other methods of classification performance, as it does not require arbitrary thresholding of results (Merow et al., 2013). Jackknife tests for variable importance were conducted and species response curves, which plot the predicted relative occurrence rate against the environmental predictor values and are an important tool for evaluating the biological soundness of the model (Merow et al., 2013), were evaluated to contextualize variable importance.
F I G U R E 2 Geographical mean center of whooping crane observations from 1990 to 2009. Aransas National Wildlife Refuge, Texas, USA 3.10 | Spatial and temporal forecasting Using MaxEnt variable importance values, temporal trends were identified for each variable and projected to 2020. A MaxEnt model matching "present" variable importance using the 2020 data was applied to the entire Texas and Louisiana Gulf coast (TX-LA coast) extending inland 75 km to capture bays and inlets. This larger area was selected since the former range of this population stretched across both states, and whooping cranes have been reintroduced to Louisiana (Pickens et al., 2017).

| Climate-related variables for spatial forecasting
Sea level rise data for three scenarios (0.3, 0.6, and 0.9 m) were obtained from the National Oceanic and Atmospheric Administration Office for Coastal Management for TX-LA coasts (https://coast.noaa.gov/ slrdata/). Data were mosaicked and re-projected to match the datasets used in modeling, resampled to 30 meters using bilinear resampling, and reclassified to isolate pixels that would be underwater given a 0.3,

| Results from ESDA
We observed a small, but clear and consistent southward shift in crane distribution starting in 2001 through 2009 (Figure 2). Part of Matagorda Island was protected in 1983, but the entire island was not conserved until 1994. It is during the period after complete conservation that we observe clear movement of the population toward the island; but while Matagorda Island extends north, the population center shifts south almost 4 km over a 10-year period (Figure 2). The area directly south of Matagorda Island is San Jose Island (Figure 1), which is not protected as part of ANWR but is privately owned and managed primarily for wildlife.
We observed two key findings through the kernel density estimate results. First, the area around ANWR used by cranes remained generally stable over two decades, although observations expanded approximately 15 km north on Matagorda Island during the study period ( Figure 3). More subtly, results showed increasing densities of observations on Matagorda Island, with hotspots growing south of the protected areas (Figure 3b-e). While there has been some movement to the north, most range expansion appears to be trending south. Additionally, over time, the density histograms ( Figure 3) shifted from heavy tailed distributions indicating greater areas of less density in 1990 to a more equal distribution across the density classes in 2009. This shift can be seen in the density maps as the diffusion of the pockets of very high densities along the Aransas mainland and Matagorda Island and the growth of the hotspots to the south.

| Ecological niche modeling results
Model fits were consistently high across the 20 years with a mean AUC of 0.934 and a standard deviation of 0.013.   The lowest AUC, occurring in 2000, was 0.9136, which is still very good in terms of model accuracy. For all 20 models, regularized gains for all variables were lower when that variable was omitted from the model compared to the model with all variables, indicating that the set of variables selected for modeling all held key information for predicting crane presence. The relative contributions of each variable to the overall model show a clear separation between the proportion of emergent herbaceous wetlands (hereafter, "proportion EHW") and the proportion of water and the rest of the variables. These two variables were consistently the most important across the 20-year time span and their importance generally increased over the same period. While the protected area variable did show similar values to the proportion of water during the early model years, with crane occurrences more likely in areas within or close to protected areas, over time the importance of this variable decreased while the importance of water and proportion EHW increased. Jackknife tests support these findings with proportion of water having the highest regularized gain when it was the only variable considered as well as the lowest gain when it was omitted from the models, indicating it holds key information not found in the other variables. Species response curves indicated that cranes are selecting areas where EHW comprises about 20-80% of the area, and water comprises less than 80%. When the proportion shifted to greater than 80% water, crane presence declined. Presences were also highest when woody wetlands were at least 10,000 m away. These variable responses were consistent over the 20-year period.
The remaining variables had considerably lower variable importance scores that remained stable or decreased over time (Figure 4). The cumulative burn index, the NDVI, distance to developed areas, and distance to woody wetlands all had variable importances less than 10% for most of the study years. Modeling results for a single year (1999) are visualized as an example ( Figure 5) and show areas of high probability of crane presence along the coastal edges of Aransas and Matagorda Island, with other suitable areas along the inner bays, inlets, and lakes. Areas of highest probability (>0.95) were limited to Aransas, Tatton, a small area on the western side of F I G U R E 6 Habitat suitability along the Texas and Louisiana Gulf Coast from a spatial projection of the MaxEnt model using 2020 data. Existing national wildlife refuges are labeled Matagorda Island, and a small area on the western coast of San Jose Island.

| Results of spatial forecasting and sea level rise analysis
Values projected from the variable importance trends for 2020 most closely matched the 1999 model (Table 2), which was then used to extrapolate the model throughout the Gulf Coast. The extrapolated model identified more than 1.7 million ha of available habitat (>0.25 probability of presence) along the TX-LA coasts and approximately 45,868 ha of highly suitable habitat (>0.95 probability of presence). The largest areas of habitat were located along the Louisiana coast, particularly in the Chenier Plains and Mississippi Delta regions (Figure 6, insets).  Categories are based on the habitat suitability index (probability of crane presence).
The sea level rise analysis indicated that available habitat throughout the Texas-Louisiana Gulf Coast will be impacted considerably by climate change (Table 3). Approximately 79% of the lowest habitat suitability class (0.25-0.5 probability of presence) will be impacted by a 0.3-m rise, and 89% of this class will be impacted by a 0.9-m rise. For the highest suitability class (>0.95 probability of presence), 97% of the area will be impacted by the conservative 0.3-m scenario, and 100% of this class will be impacted by the 0.6-m scenario. Specifically, the areas of the Louisiana Coast around the Chenier Plains and the Mississippi Delta regions will lose most of their currently suitable whooping crane habitat under the most conservative scenario (Figure 7).

| DISCUSSION
Protection under the US Endangered Species Act has allowed the AWBP to recover from a low of fewer than 20 birds to a recent count of more than 500 individuals, but another doubling of the population is needed to reach the numbers needed for down-listing from endangered to threatened (CWS & USFWS, 2007). Lack of available habitat may prevent the population from meeting this milestone in the future (Lewis, 1995;Metzger et al., 2020;Stehn & Prieto, 2010). The results of this study highlight several habitat-related variables related to whooping cranes overwintering habitat and indicate where potential future unoccupied habitat may exist for the whooping crane.
The exploratory spatial data analyses identified relative stability in the spatial extent of the crane distribution over time, with clustering occurring in and around ANWR. An increase in population size likely contributed to expansion north and south of the refuge. Territoriality may provide a partial explanation for the clustering observed; as new pairs form, subadults establish territory as close as they can to their parents (Stehn & Prieto, 2010), so maintaining preferred habitat immediately surrounding the refuge is imperative for future AWBP viability. With over 200 breeding pairs in the latest count (Harrell & Bidwell, 2019;USFWS, 2020) and a current average territory size of 172 hectares, additional areas will be needed (Gil-Weir et al., 2012) to support continued expansion. Beyond San Jose Island, there is not contiguous unoccupied available habitat to the south (Figure 6), suggesting that the continued expansion of the birds in this direction may not be sustainable in the long term.
While we did not address habitat fragmentation directly in our study, our findings suggest it could limit population expansion and future recovery efforts. Large areas of potentially available habitat are located north of ANWR in the Galveston, Texas area ( Figure 6) near three existing National Wildlife Refuges (NWR), the closest which is approximately 65 miles away. While this distance is not prohibitive for a crane to traverse, prior research suggests that rather than exploring new territory north and south of the refuge such as these areas near Galveston, the cranes may instead be reducing their territory sizes (Harrell & Bidwell, 2019;Stehn & Prieto, 2010) to remain close to ANWR. While the available habitat areas around Brazoria, San Bernard, and Big Boggy NWRs are used by whooping cranes, they are separated from current ANWR habitat (USGS, 2020), indicating unlikely use as territories in the immediate future, although they remain important areas of crane use.
Focusing on inland wetland conservation will likely become necessary for whooping crane overwintering survival as coastal habitat is lost. If sea level rise remains conservative, habitat that will not be submerged under the 0.9-m scenario can help provide habitat and grow the population of non-migratory whooping cranes. If the 0.9-m estimate proves to be too conservative, management, and conservation efforts will likely need to shift the focus to protecting and studying inland wetlands. The habitat types used in breeding areas and along migratory routes include inland wetlands (Armbruster, 1990;Timoney, 1999), which may provide insight to what type (size, characteristics, etc.) of overwintering wetlands can be protected that provide for both habitat requirements and habitat preferences of the cranes. Cranes have recently been observed overwintering in locations distant from ANWR (Wright et al., 2014), and a non-migratory population in Louisiana has been documented using inland wetlands associated with agricultural land (Pickens et al., 2017). However, the non-migratory population in particular uses agricultural wetlands (e.g., rice and crawfish aquaculture) (Pickens et al., 2017), while the agricultural areas around ANWR are primarily corn and cotton (USDA, 2021). Developed areas may also become an issue for whooping crane conservation in the future as human population densities in coastal regions are projected to increase (Neumann et al., 2015), leading to competition for available land space, directly impacting available habitat required by overwintering whooping cranes.
Whooping cranes spend most of their time at their overwintering habitat foraging (LaFever, 2006), and their preference for crabs (Hunt & Slack, 1989;Pugesek et al., 2013) likely influences habitat selection. Preference for habitat with emergent herbaceous wetland between 20%-80% coverage and water not exceeding 80% suggests that a juxtaposition of these habitat types is important. The amount of wetland habitat is directly related to the availability of crabs, and past research has found that crab size was lower in inner marsh areas compared to the connected ponds of the outer marsh area around the refuge (Hoeinghaus & Davis, 2007), where the proportion of water is greater. Our results support the utilization of these areas, with some of the highest suitable habitat and crane densities in these outer marsh areas ( Figure 5). Additional research is needed to determine what factors may be related to the configuration or interspersion of the emergent wetlands and water.
ANWR is just one of nearly 20 NWRs located along the US TX-LA coast, several of which fall within areas identified as available whooping crane habitat in this study ( Figure 6). Since the mission of the NWR system is to create a "network of lands and waters dedicated to conserving America's rich fish and wildlife heritage" (USFW, 2022), an opportunity exists to further assess how other refuges along the coast might be leveraged for whooping crane conservation. Of particular importance are Brazoria and San Bernard near Galveston, Texas, mentioned earlier. However, while more recent GPS data (USGS, 2020) has indicated that the birds have established territories in a "leapfrog" fashion, they tend to slowly expand from the core area. Establishing additional protected habitat that is contiguous with the existing preserved areas within ANWR will be needed to support continued growth. The relevance of protected areas for whooping crane conservation should not be undervalued. Our results indicated an uptick in the importance of this variable in the later years, demonstrating that protected areas are still very much needed to support conservation and delisting goals for threatened and endangered species.
Extending our species distribution model across the larger US TX-LA coast identified more substantial areas of available habitat in Louisiana, which may be notable for the future of experimental and non-migratory populations (CWS & USFWS, 2007;Pickens et al., 2017;Thompson, 2018). The areas of note are located along the Chenier Plain, where a non-migratory population of whooping cranes currently exists (Pickens et al., 2017), and in the Mississippi Delta region. While this study predicted larger swaths of available habitat in these areas, fragmentation of these regions from ANWR may limit actual availability. Predicted sea level rise, which is projected to exceed the conservative 0.3 m value used here (Bamber et al., 2009;Lhermitte et al., 2020), will also likely significantly reduce the currently available habitat and increase distances between habitat areas. There is a need to focus on restoring wetlands further inland to offset the future loss of coastal wetlands, as well as to evaluate the condition of inland wetlands that already exist. However, our analysis did not factor in land cover transitions resulting from sea level rise (e.g., an increase in wetlands), which could create additional habitat for cranes.

| Limitations
In addition to the whooping crane survey methods reported above, there are several limitations to note. We only had access to reliable land cover types, without specific data on food sources across the study areas. We can infer that certain variables (water, wetlands, burning) are likely related to food resources that drive where cranes were observed, but we did not test this directly. In the future, it would be useful to include a direct measure of crabs into models. We also did not include a measure of agriculture, which has been found to be important in other whooping crane studies (Pickens et al., 2017), but it has not been found to be as highly relevant for the overwintering population (Golden, 2020), except in years of extreme draught or low food availability . Future research should investigate the importance of inland agricultural areas for habitat provision of the AWBP cranes, especially as sea level rise causes a contraction of existing preferred areas. The use of future climate change scenarios should also be considered to better understand how future habitat availability may change. The Gulf coastal region is hydrologically dynamic, and static observations without timestamps make it difficult to determine the state of the bay during the observation. Changing conditions that can be incorporated into studies that employ GPS tracking include salinity, water depth, and other tidal conditions. Some studies have found salinity related to whooping crane behavior, resource availability, and habitat quality (Tiegs, 2017); however, there is not a long-term record of salinity and we could not incorporate it into our analyses. Additionally, the study did not consider potential land cover transitions that might result from sea level rise, such as the conversion of grasslands to wetlands. In short, there is a large amount of uncertainty that surrounds regional habitat availability under climate change. It is possible that the models used here may underestimate future habitat availability, particularly as we did not assess possible habitat expansion.
A facet of whooping crane landscape ecology that was not assessed in this study but may be important to crane habitat selection is the spatial configuration of water and emergent herbaceous wetland land cover classes. Configuration involves the spatial arrangement and spatial distribution of land cover classes (Frazier, 2019), and it impacts factors such as the size and shape of the patches of land cover patches as well as the amount of edge interface between the wetlands and water. The interposition and juxtaposition of these two land cover classes in the highly suitable areas may provide additional insights into the habitat characteristics cranes are selecting for along the Texas coast.

ACKNOWLEDGMENTS
This work was supported by the Bollenbach Endowment and the Groendyke Endowment at Oklahoma State University. Amy Frazier's participation in the project was supported by NSF Grant 1934759. The authors thank the United States Fish and Wildlife Service for providing the whooping crane location data and prescribed burn data.