Spatiotemporal patterns in Golden-cheeked Warbler breeding habitat quantity and suitability

,


INTRODUCTION
The Golden-cheeked Warbler (GCWA; Setophaga chrysoparia) is a Neotropical migratory songbird that is listed as endangered under the U.S. Endangered Species Act (ESA) and is heavily impacted by forest conversion. The species is a habitat specialist, breeding exclusively in the Edwards Plateau of central Texas, commonly referred to as the Texas Hill Country, and preferring mature, mixed ashe juniper (Juniperus ashei) and oak (Quercus) woodlands (Pulich 1976, Ladd and Gass 1999, Long et al. 2016). These tree species provide critical material for nest-building and shelter for main GCWA food sources , Marshall et al. 2013. A wealth of literature assessing the influence of habitat factors on measures of GCWA survival (e.g., presence, reproductive success, nest survival, and density) finds that forest composition, age, and patch size are important to species success (Pulich 1976, Shaw and Atkinson 1990, USFWS 1992, Jetté et al. 1998, Magness et al. 2006, Diamond 2007, Klassen et al. 2012, Stewart et al. 2014, Reidy et al. 2017, Colón et al. 2019, Trumbo et al. 2021. Their specialized preference for this already rangerestricted forest type, as well as the high rates of habitat loss from urban developments and transportation infrastructure, led to the listing of GCWA as endangered by the U.S. Fish andWildlife Service in 1990 (USFWS 1990). Since listing, the threat of habitat loss has persisted with an estimated 29% reduction in total GCWA breeding habitat between 2000 and 2010, and with most remaining, intact habitat occurring in the southern portion of their breeding range (Duarte et al. 2013). As of 2020, a very small portion (1.5%) of Texas lands are protected or managed in a way that is consistent with biodiversity conservation, designated as GAP code 1 or 2 from the Protected Areas Database of the U.S. (USGS GAP 2020, Dreiss andMalcom 2022).
Additionally, human impact is on the rise in Texas landscapes and may compromise habitat quality. In recent years, Texas had the largest increases in population of any state in the U.S. (U.S. Census Bureau 2020). In the 30-year period that was studied, the state's population grew 73.9% and growth is projected to continue at a steady rate to 2050, with 88.3% increase in the next three decades (Texas Demographic Center 2019). There are near 40 urban census areas that overlap the GCWA breeding range. Of those, Austin and San Antonio urban areas have the largest populations and the highest population density (U.S. Census Bureau 2020). They are also undergoing significant population growth, with 22% and 8% increase in Austin and San Antonio cities, respectively (U.S. Census Bureau 2020), and sprawl events. The two urban areas are located along the southeastern border of the GCWA range (Fig. 1). Population growth for this area of the state is projected to continue at one of the highest rates in the country (Texas Demographic Center 2019). As such, the threat of further habitat loss and degradation remains, and continued monitoring and evaluation of landscapes and populations is necessary to understand the progress of GCWA recovery (Eichenwald et al. 2020). Regular reviews of the best available science and commercial information are a mandate of the ESA, section 4(c)(2), and are conducted once every five years to examine population trends, threats to recovery, and accuracy of the listing. Science that demonstrates a range-wide understanding of available breeding habitat conditions and distribution is critical to engaging federal action for proper protections for species recovery (La Sorte et al. 2015). Though coarse, this information may be used to estimate population size and viability and, in the case of habitat conversion, can also serve to evaluate the status and trends of major threats to species recovery (McGowan et al. 2017). For GCWA, analyses are generally focused on habitat quantity, i.e., the amount of habitat, but additional nuances of habitat suitability, i.e., the potential capacity of a given habitat to support a species (Peak and Thompson 2013, Peak and Thompson 2014, Reidy et al. 2016, Reidy et al. 2017, Reidy et al. 2018) may result in a more refined understanding of population dynamics and can help managers prioritize habitats for conservation and restoration. Although reporting on changes in habitat quantity can aid understanding of absolute habitat loss from conversion, it may not be reflective of impacts from changes in habitat configuration or fragmentation over range-wide scales. Most important, regular updates to species habitat suitability, quantity, location, and use are necessary to understand longer-term temporal trends for better informed conservation efforts and for consideration in regular federal assessments revisiting species listing and other decisions. As such, up-to-date spatiotemporal trends in both the suitability and quantity of species habitat at range-wide scales are important for informing conservation management and supporting continued protections for threatened and endangered species.
The last five-year review for GCWA was in 2014 (USFWS 2014). Since then, there have been legal challenges regarding species recovery and listing status (General Land Office of the State of Texas v. United States Department of the Interior et al. 2020, USFWS 2021. Currently available science on GCWA breeding habitat quantity assesses temporal trends on short timeframes (no more than a decade) and is already outdated (Duarte et al. 2013). Additionally, available analyses of habitat suitability are restricted to fractions of the breeding range (Loomis Austin 2008, Heger andHayes 2013). There is a need for an amended breeding habitat assessment to reflect recent landscape changes and data availability as well as to inform upcoming species status assessments and the next steps in conservation planning.
The objective of this analysis was to conduct a comprehensive spatiotemporal model of GCWA breeding habitat distribution. As part of this study, we used geospatial data spanning 1985 to 2018 in order to (1) update spatiotemporal patterns in GCWA habitat quantity; and (2) analyze patterns in suitability over GCWA breeding range. We assessed both changes in evergreen and mixed forests within the breeding range, and juniper-oak forest more specifically.

METHODS
This study focused on spatiotemporal changes to habitat for the entirety of the GCWA breeding range as defined by the U.S. Geological Survey's Gap Analysis Project (Fig. 1, Table 1). Breeding and nesting activities are confined to central Texas, where habitat varies in tree density and composition. Juniper-oak woodland is more commonly found in the southern and eastern portions of the breeding range (Duarte et al. 2013). Nesting habitat is generally defined by the tree-species composition. GCWA nest in habitat made up of mature ashe juniper and a combination of other species, such as live oak (Quercus fusiformis), bastard oak (Quercus sinuata), Texas oak (Quercus buckleyi), post oak (Quercus stellata), blackjack oak (Quercus marilandica), lacey oak (Quercus laceyi), sugarberry (Celtis laevigata), Texas ash (Fraxinus albicans), cedar elm (Ulmus crassifolia), black cherry (Prunus serotina var. eximia), pecan (Carya illinoinensis), and little walnut (Juglans microcarpa; USFWS 1992). GCWA occupancy and density have been shown to be positively related to total woodland, especially mixed juniper-oak woodland, canopy height, and patch size (Wahl et al. 1990, Magness et al. 2006, Peak and Thompson 2013, Reidy et al. 2016. GCWA also showed a negative response to increasing urban cover in the landscape (Reidy et al. 2016), which is especially concerning given the rapid population growth near some of the largest tracts of breeding woodland.

Data acquisition
Given the primary goal of assessing change over time, we considered publicly available datasets on land cover and ecological systems with spatial information for more than a single time step. We recognize there are other available datasets that may be effective for other analyses and comparisons in GCWA habitat, e.g., Texas Parks and Wildlife's Ecological Mapping Systems of Texas. Analyses herein focus on locations inside the GCWA breeding range as defined by the U.S. Geological Survey's Gap  Sohl et al. 2018Sohl et al. 1985Sohl et al. used (1938Sohl et al. -1992 250m Habitat suitability National Land Cover Database Multi-Resolution Land Characteristics Consortium 1992Consortium , 2001Consortium , 2006Consortium , 2011Consortium , 2016 30m

Spatial analyses Habitat quantity
We used Google Earth Engine implementation of the LandTrendr algorithm (Kennedy et al. 2018) to identify land-use changes within the GCWA breeding range for every year between 1985 and 2018 from Landsat imagery via breakpoints in temporal trends of the Normalized Difference Vegetation Index (NDVI; see Eichenwald et al. 2020). The algorithm was used to identify areas where vegetation changed quickly over a short period of time, including from natural or prescribed burns. Additional datasets were used to focus on specific habitat types. First, we used this data to calculate quantity of all evergreen and mixed (i.e., neither deciduous nor evergreen species are more than 75% of dominant tree cover) forest areas within the GCWA breeding range for every year between 1985 and 2018. These calculations were compared with historic forest disturbance data from the National Land Cover Dataset (NLCD; Dewitz 2019) provided for every two to three years from 1986 to 2019. It was assumed that all pixels labeled as either "never experiencing a disturbance" or "experiencing a disturbance after 1985" were forested in 1985.
Second, we used LANDFIRE existing vegetation type data to calculate changes specific to forested habitats with mixed ashe juniper and oak overstories: Edwards Plateau mesic canyon, Edwards Plateau riparian forest and woodland, warm temperate urban mixed and evergreen forest (NatureServe 2018). LANDFIRE's existing vegetation type data are available for 2001, 2008, 2010, 2012, 2014, and 2016. The vegetation datasets were developed using new plots, processing methods, and classification models starting in 2016 (see landfire.gov/lf_remap.php). More specifically, the classification was expanded to disaggregate some vegetation types. To account for changes to vegetation categories, we used a crosswalk between legends to find the equivalent vegetation system value, a four-digit code given to each ecological system (see landfire.gov/remapevt_assessment.php). For data on mixed/evergreen (from NLCD) and juniper-oak forests (from LANDFIRE) separately, we calculated the area of disturbed and undisturbed habitat throughout the entire breeding range, in urban and non-urban portions of the range, and in the metropolitan areas of Austin and San Antonio separately, as defined by U.S. Census urban-area boundaries for every year the landscape data were available. U.S. Census boundaries were kept fixed to the most recent year (2018) to allow for capturing changes made to all lands that ultimately fall into the fullest urban extent.

Habitat suitability
We applied a habitat suitability framework to determine location and acreage of suitable GCWA habitat for the entire breeding range in 1985 and 2016, the most current year of the National Land Cover Dataset (Heger andHayes 2013, Dewitz 2019). The framework for scoring habitat suitability is based on a large body of literature citing forest composition, landscape fragmentation, and edge effects as related to GCWA presence, survival, and breeding success (e.g., Pulich 1976, USFWS 1992, Jetté et al. 1998, Magness et al. 2006, Peak 2007, Long et al. 2016, Reidy et al. 2017, Reidy et al. 2018, Colón et al. 2019. Based on the scoring, intact mixed or evergreen forests part of a larger forested landscape should receive the highest score (most suitable). The following scores were combined such that each location (i.e., 30 m pixel) in the final output had a score in the range of 0 (not suitable) to 4 (most suitable). The suitability analysis was done for mixed/evergreen forests in 1985 and 2016 (from NLCD) and for juniper-oak forests in 2001 and 2016 (from LANDFIRE), separately. The general process for scoring was based on three main characteristics (Fig. A1.1, Appendix 1).
We reclassified land-cover types to two general groups. Pixels assigned a mixed or evergreen forest type and deciduous forest pixels within 100 m of mixed/evergreen forest were included in the analysis. All other land-cover types were considered not suitable for breeding and were assigned a score of 0. For analysis on juniper-oak forests, aforementioned existing vegetation types from LANDFIRE dataset were analyzed, and all others were assigned a score of 0.
We took landscape context into consideration by using neighborhood statistics to determine the percent of forest cover in a 218 m radius of each mixed/evergreen forest pixel. Generally, patches below this size are not likely to successfully support breeding (Butcher et al. 2010). The analysis assigned each pixel a value indicating the proportion of the surrounding neighborhood that is classified as mixed/evergreen (when using NLCD) or juniper-oak forest (for LANDFIRE). Pixel values were reassigned a score between 0 and 4: areas that were 80-100% forested received the highest score (4) and areas 0-20% forested the lowest (0). To account for potential edge effects, pixel scores were docked 1 point if they were located within 90 m of the edge of the forest habitat.
This model was previously tested for Travis County and validated by comparing model predictions of habitat suitability to that of GCWA sightings in Balcones Canyonlands National Wildlife Refuge (Heger and Hayes 2013). We used this framework to assess the suitability of breeding-range habitat in 1985 and 2016 using land cover data (Sohl et al. 2018, Dewitz 2019. However, the historical (1985) land use and land cover data for the contiguous U.S. has a coarser resolution than current models (Sohl et al. 2018; Table 1). To account for differences in data resolution, we resampled and masked the land cover data from 1985 by historic forest disturbance data (Dewitz 2019), assuming that all pixels labeled as either "never experiencing a disturbance" or "experiencing a disturbance after 1985" were forested in 1985.
For each dataset, we calculated the amount of habitat by score throughout the entire breeding range. Descriptive statistics were also generated to compare results in urban areas and specifically for Austin and San Antonio metropolitan areas.
We assessed changes in patch size for more suitable breeding habitat areas. This was done by grouping pixels by score and proximity. For each pixel with a score of 3 or 4, neighboring pixels were identified and those with a similar score (3 or 4) that were connected by a common edge or corner were included in the patch. Area of land covered by each patch was calculated by counting up the number of pixels in each patch and converting to hectares. Average patch size was calculated for mixed/evergreen forests in 1985 and 2016 (NLCD) and for juniper-oak forests in 2001 and 2016 (LANDFIRE). Differences between years were analyzed using an unpaired sample t-test where unequal variances were assumed.

Protected areas overlay
We combined the results from the habitat quantity and suitability analyses with the protected areas database of the U.S. (PADUS 2.0) to determine the amount and quality of habitat that is currently protected. U.S. Geological Survey's Gap Analysis Program (USGS GAP) codes are specific to the management intent to conserve biodiversity (USGS GAP 2020). GAP 1 and 2 areas were managed in ways typically consistent with conservation and were considered protected in this context.

Mixed/evergreen forest cover
Prior to 1985 there were over 1.83 million forested hectares within the GCWA breeding range, but by 2016, the forested area had decreased by 13% to under 1.6 million ha (A, Fig. 2). Forest conversion was more extreme in parts of the range in metropolitan areas, with 24% forest loss in the Austin metropolitan area and 32% loss in the San Antonio metropolitan area. The steepest rates of forest loss occurred in 2008 and 2011 in both urban and nonurban areas. Habitat suitability scores also declined during this period throughout the breeding range (A, Fig. 3). In the 1980s over 10% of the forested habitat within GCWA breeding range was intact mixed or evergreen core forests (score 3 or 4). In 2016 the same habitat made up 5% of the breeding range, indicating a 45% decrease in more suitable breeding habitat (A, Fig. 3). Remaining suitable habitat (score 3 or 4) was more fragmented in 2016 (14.6 ± 0.7 ha), with significantly smaller patch sizes than in 1985 (27.9 ± 2.0 ha) for mixed/evergreen forest ( Fig. 5; t = 1.96, p < 0.001).
Generally, more suitable habitat (score 3 or 4) was concentrated along the southeastern extent of the breeding range and some forested areas northwest of Austin. In addition, there was more suitable habitat in the northern portion of the breeding range than in 1985 (A, Fig. 4). Forests within protected areas (i.e., GAP codes 1 and 2) of the breeding range experienced fewer sharp declines (27% to 19%) in more suitable habitat (score 3 or 4) between 1985 and 2016 (27% decrease; Fig. 4).

Juniper-oak forest cover
Declines were larger for juniper-oak forests than for all mixed/ evergreen forests regardless of whether the forests were non-urban or urban areas, particularly between 2008 and 2012 (B, Fig. 2). Juniper/oak forests were comprised of 33% intact core habitat in 2001, but by 2016, only ~14% was intact core, suggesting a drastic decline in highly suitable habitat (B, Fig. 4). Remaining suitable habitat (score 3 or 4) was more fragmented in 2016 (13.0 ± 1.0 ha), with significantly smaller patch sizes than in 2001 (18.1 ± 1.2 ha) for juniper-oak forest (t = 3.19, p = 0.001). Juniper-oak forests within protected areas (i.e., GAP codes 1 and 2) of the breeding range experienced fewer sharp declines (69% to 64%) in more suitable habitat (score 3 or 4) between 2001 and 2016 (B, Fig. 4). However, protected habitats currently represent only 10% of most suitable (score 4) juniper-oak forests within the breeding range.

DISCUSSION
To inform conservation decision-making, we quantified the absolute change in mixed/evergreen forest cover and, more specifically, juniper-oak forests within GCWA breeding range and changes in breeding habitat suitability over a 30-year period. Overall, we estimated a 13% loss in mixed/evergreen forest habitat with high spatial heterogeneity in landscape conversion closely tied to human developments. This value is lower than relatively recent estimates and may be reflective of changing forest dynamics that occurred over a longer study period , Duarte et al. 2013. Declines may in fact be larger for forested habitats specific to GCWA ecology given a 19% loss of juniper-oak forests since 2001. Additionally, drastic declines in habitat suitability suggest that 13% is an underestimation of effective habitat loss, as degradation may impact the ecological viability of remaining forests for GCWA nesting. The amount of intact core forest habitat fell 45% in the 30-year period (42% for suitable juniperoak habitat since 2001) leaving more suitable breeding sites concentrated along the southeastern extent of the breeding range and in protected areas.
In recent years, Texas has had the largest increases in population of any state in the U.S. (U.S. Census Bureau 2020). Within GCWA breeding range, at least four counties are projected to see population increases of over 100% by 2050, all of which coincide with areas of more suitable habitat in the southeastern parts of the range: Williamson, Hayes, Comal, and Kendall counties. Increased development pressures in the Texas Hill Country will likely continue to drive the trends of GCWA habitat disturbance and degradation. Our data indicate that suitable forests have undergone fragmentation resulting in smaller average habitat patches. Similar trends have been reported more specifically for ashe juniper distributions across the state due to an increase in pastureland and development (Diamond 1997). Increased fragmentation results in smaller patch sizes and more forest edge, the latter of which has been linked to decreased nest survival for GCWA (Peak 2007, Reidy et al. 2009). Forest edges can leave GCWA nests less concealed and more vulnerable to nest predation (Reidy et al. 2008). Additionally, fragmentation of breeding habitat may represent barriers to dispersal of birds and important genetic material (Lindsay et al. 2008). Hence, there is already evidence of notable genetic differentiation among populations of GCWA, having important implications for management of species like GCWA that are relatively vagile but highly specialized Grey color gradient indicates the habitat suitability score (0-4) of the first time step, with dark grey representing areas of higher suitability (score 3 or 4). Red represents pixels that were more suitable breeding habitat (score 3 or 4) in the first time step, but became less suitable by 2016 (score 0-2). Mixed/evergreen forest analyses use National Land Cover Database and juniper-oak forest analyses use LANDFIRE existing vegetation type.
in their habitat preferences. Restoration and protection of large and connected patches may be the best option for conserving or recovering such species (Young and Clarke 2000).
We found that only 10% of the species' most-suitable habitat is in protected areas, creating both challenges and opportunities. These lands, because they are managed in ways consistent with biodiversity conservation, generally represent more suitable habitats with fewer human disturbances (Dreiss and Malcom 2022): where 27% of the juniper-oak forests throughout the entire breeding range had highly suitable scores in 2016, 64% of the forest in protected areas received the same scores. Our findings indicate that protected areas within GCWA breeding range also exhibited declines in suitability, but degradation was buffered relative to the overall range. As human populations grow and landscape conversion continues, protected areas are expected to grow in importance. A preliminary analysis of GCWA occupancy models from Morrison et al. (2010) also demonstrates the importance of protected areas to GCWA success: areas with at least 70% probability of occupancy make up 13% of the breeding range, but 62% of protected areas. Public protected areas can play a central role in habitat conservation efforts because they are more amenable to the application of broad-scale management strategies that more closely align with species conservation. However, the extent to which public protected areas can benefit migratory bird populations depends on how well protected areas are represented within the breeding range (La Sorte et al. 2015). Currently, areas managed for conservation (GAP status 1 and 2) represent 3.23% of the breeding range. Lands with more intermediate mandates (GAP status 3) provide a higher degree of flexibility for the implementation of management recommendations more closely aligned with maintaining biodiversity. However, these lands are also limited in the Edwards Plateau in central Texas (1.71%). Expansion of protections to key habitats would require that resources be spent by conservation agencies on land acquisition, conservation easements, or providing incentives for conservation on private lands. Our findings demonstrate a need for strengthening current conservation measures and improving private lands protections within GCWA habitat to ensure higher availability of more suitable habitat, which should aid in species recovery. Newer proposals, such as Executive Order 14008 (2021), to protect at least 30% of U.S. lands and waters by 2030 to address the biodiversity and climate crises may provide additional opportunities for land designations and conservation efforts for imperiled species like GCWA. Whereas most GCWA habitat conservation dollars have been spent conserving GCWA breeding habitat on the outskirts of the cities of Austin and San Antonio, our findings support previous work demonstrating higher rates of habitat conversion near metropolitan areas (Duarte et al. 2013). Given the scarcity of public lands, the distribution of intact forest habitat, and the relatively high habitat loss and degradation in and around metropolitan areas, future GCWA habitat conservation efforts should be more focused on supporting current protected areas and expanding protections to suitable habitats in the Balcones Canyonlands and Fort Hood areas and regions west of San Antonio and Fort Worth. It will be particularly important to work collaboratively with land owners to strengthen protections on private lands (e.g., through conservation easements). Additionally, projected species distribution models reflecting climate change impacts on tree species indicate that the Texas Hill Country will continue to be a stronghold for ashe juniper, with potential for population stabilization and maybe even increase/spread to the northeast (McKenney et al. 2007). This suggests that efforts to conserve or restore quality GCWA habitat will have long-term benefits.
We recognize the limitations of some of our analyses, which equate all available mixed or evergreen forests within the breeding range. However, our complementary analyses specific to juniper-oak forests, even though available for only half of the time period, tell a similar story. Further advances in available technology (e.g., LiDAR) may clarify spectral or structural differences in forest composition that could improve spatial analysis of GCWA suitable habitat. Regardless, overall classification accuracies of the habitat loss dataset followed methods that average a mean absolute error of less than 3% (Kennedy et al. 2018). Additionally, datasets used for assessing habitat suitability, though they represent the most current version available, are already out of date. Collectively, this indicates that we have likely underestimated the rate of habitat loss. Finally, although widely used, there are key limitations to using a habitat suitability index. Assignment of qualitative scores to habitat quality may be unreliable depending on the relevance of selected habitat variables, quality of the geospatial data, and reliability of the wildlife-habitat relationships (Cole andSmith 1983, Roloff andKernohan 1999). Ultimately, our suitability index serves the intent of quantifying the value of habitats to aid management considerations and alternatives in species-specific conservation and restoration. However, we emphasize the importance of validation to determine the reliability and utility of habitat suitability models. The model used in this analysis was originally built using a subset of the GCWA range and validated using GCWA sightings in Balcones Canyonlands National Wildlife Refuge (Heger and Hayes 2013). However, future work validating the predictions of habitat quality across the entire breeding range is needed to determine the model's reliability. In the context of federal species listing and review, landscape-change analyses, habitat identification and classification, and the characterization of trends over time must be considered. The metrics used in this study are meant to facilitate such consideration: they can be applied to multiple scales and interpreted by non-GIS audiences, helping diverse stakeholder groups to engage in the conservation decision-making process. Although current limitations in data, technology, and metrics may influence interpretation of landscapes, the intent is for future research to continue to improve upon this methodology and on our understanding of the changing habitat.
Human landscape modification is likely to continue in the Texas Hill Country, but conservation and land management actions can be taken to minimize further habitat loss and degradation in GCWA breeding range. This information will assist researchers and managers in prioritizing range-wide breeding habitat conservation efforts and highlights the significant role land management for conservation biodiversity plays on the landscape. There remains a need to grow the network of protected areas for GCWA restoration. Further, continued, regular spatiotemporal assessments of habitat quantity and suitability are necessary to assess changes to species potential for persistence and extrapolate population viability, given these dynamics.
Responses to this article can be read online at: https://www.ace-eco.org/issues/responses.php/2245