Spatiotemporal variation in range-wide Golden-cheeked Warbler breeding habitat

Habitat availability ultimately limits the distribution and abundance of wildlife species. Consequently, it is paramount to identify where wildlife habitat is and understand how it changes over time in order to implement large-scale wildlife conservation plans. Yet, no work has quantified the degree of change in range-wide breeding habitat for the Golden-cheeked Warbler (Setophaga chrysoparia), despite the species being listed as endangered by the U.S. Federal Government. Thus, using available GIS data and Landsat imagery we quantified range-wide warbler breeding habitat change from 1999–2001 to 2010–2011. We detected a 29% reduction in total warbler breeding habitat and found that warbler breeding habitat was removed and became more fragmented at uneven rates across the warbler’s breeding range during this time period. This information will assist researchers and managers in prioritizing breeding habitat conservation efforts for the species and provides a foundation for more realistic carrying capacity scenarios when modeling Golden-cheeked Warbler populations over time. Additionally, this study highlights the need for future work centered on quantifying Golden-cheeked Warbler movement rates and distances in order to assess the degree of connectivity between increasingly fragmented habitat patches.


INTRODUCTION
It is well accepted that habitat loss and fragmentation caused by anthropogenic modifications to the landscape are primary drivers of wildlife population declines. Although wildlife management actions are usually implemented at the site (or local) level, recovery plans for species of concern are often structured at the range-wide (or regional) level. Thus, monitoring large-scale habitat availability is paramount in wildlife conservation efforts. This information allows researchers and managers to identify factors that may limit the distribution of the species, to create population models based on realistic carrying capacity scenarios, and to prioritize areas for habitat conservation efforts based on habitat availability, habitat patch metrics, as well as previous rates of habitat loss and fragmentation.
The Golden-cheeked Warbler (Setophaga chrysoparia; hereafter, warbler) is a migratory songbird that breeds exclusively in the mature (generally 4.6 m tall and 12.7 cm in diameter at breast height) juniper-oak/deciduous woodlands of the Texas Hill Country and Balcones Escarpment (Pulich 1976, Kroll 1980, Ladd and Gass 1999. In 1990, the species endemism coupled with perceived high rates of habitat loss across the landscape resulted in the species' emergency listing as endangered by the U.S. Federal Government (United States Fish and Wildlife Service [USFWS] 1990). Despite the species' federal listing status, most of the warbler's breeding range is comprised of privately owned property. Thus, nearly all breeding habitat can be subject to removal for multiple anthropogenic reasons (i.e., direct human removal due to perceived negative effects of Ashe juniper [Juniperus ashei] on livestock herbivory and groundwater availability, urban sprawl, land-use change, and prolonged drought).
In recent years, a suite of warbler habitat suitability models (hereafter referred to as warbler habitat maps) have been developed for the entirety of its breeding range (reviewed in ). Yet, the scope of these warbler habitat maps was to delineate a static description of the locations and metrics of habitat patches in central Texas. Groce et al. (2010) quantified woodland loss/gain across the warbler's breeding range from 1992 to 2001 using the National Land Cover Database (NLCD). However, not all woodland areas in central Texas should be considered warbler breeding habitat. Further, information concerning more recent habitat dynamics across the entire warbler breeding range is lacking. An up-to-date understanding of the dynamics of warbler breeding habitat patch metrics is particularly important since increased forest edge, decreased habitat patch size, and reduced connectivity of habitat patches have been linked to deleterious effects on warbler reproductive success (Peak 2007, Reidy et al. 2009, increased genetic differentiation among warbler populations (Lindsay et al. 2008), a decline in genetic diversity over time (Athrey et al. 2011), and reduced warbler occupancy of habitat patches (Collier et al. 2010, Warren et al. 2013. Notably, habitat patches with lower occupancy probabilities have been reported to support lower warbler abundances   Zitske et al. 2011), and therefore, have the potential to negatively affect Golden-cheeked Warbler survival.
The objective of this study was to model rangewide warbler breeding habitat dynamics from 2000 to 2010 using available GIS data and Landsat imagery. In order to manage the recovery process for warblers, the species' recovery plan delineated eight recovery units that were based on geology, vegetation, and watershed boundaries (USFWS 1992). Therefore, we quantified patch metrics in both time steps for each warbler recovery unit in order to assist researchers and managers in prioritizing breeding habitat conservation efforts for the species, which has recently been highlighted as a research need for warbler conservation (Hatfield et al. 2012).

Study area
This study pertains to the entirety of the warbler's breeding range. As stated above, the warbler's breeding range is confined to central Texas, USA and has been divided into eight recovery units to manage the recovery process for the species (Fig. 1). Warbler breeding habitat varies in density across the species' breeding range with a sparseness of habitat occurring in the far northern and far western portions of the range. Tree species that are widespread across the Texas Hill Country and Balcones Escarpment include Ashe juniper, live oak (Quercus fusiformis), Texas oak (Q. buckleyi ), post oak (Q. stellata), blackjack oak (Q. marilandica), Lacey oak (Q. glaucoides), shin oak (Q. sinuata), sugarberry (Celtis laevigata), Texas ash (Fraxinus texensis), cedar elm (Ulmus crassifolia), black cherry (Prunus serotina), pecan (Carya illinoinensis), and black walnut (Juglans nigra).

Data acquisition
Landsat 5 Thematic Mapper (TM) images at a 30-m spatial resolution were downloaded from the USGS Earth Explorer website (earthexplorer. usgs.gov/). Images were selected if they had ,10% cloud cover and were acquired in two seasons, winter (December-January) and summer (late April-July). Using images during both of these periods was desirable because it allowed for identification of areas occupied by evergreen/ mixed vegetation vs. deciduous vegetation since both types of vegetation will be green in the summer and deciduous vegetation will be either senesced or defoliated in winter images. Eighteen images were used per time step (we used two seasons worth of imagery per time step and the total number of images needed to cover the warbler's breeding range is nine, see Morrison et al. 2010: Fig. 1-5). For the first time step, images used for analysis were obtained from 1999-2001 and for the second time step, images were obtained from 2010-2011. This was necessary since not all Landsat 5 scenes were acquired on days with ,10% cloud cover within the restricted seasons in 2000 and 2010. Refer to Table 1 for a summary of image acquisition dates for each time step of the analysis.

Spatial analyses
Images were pre-processed in ERDAS Imagine 2011. Image pixel values were converted from brightness values to at-sensor reflectance in the model builder within Imagine using the radiometric coefficients and conversion algorithms presented in Chander et al. (2009). Multi-date image composites were then created by stacking  23-Jan-00 15-Jul-99 04-Dec-10 24-Apr-10 29_39 13-Dec-99 24-May-01 09-Jan-10 05-Jun-11 Note: Images are at a 30-m resolution, selected if they had ,10% cloud cover, and were acquired in two seasons, winter (December-January) and summer (April-July).
v www.esajournals.org the reflectance images from different seasons that have the same row/path and were within the same time step.
To mitigate potential classification accuracy issues due to temporal differences attributed to cloud cover (e.g., adjacent scenes acquired during the same season, but different years), we classified the composite path/row images prior to creating a mosaic of the thematic output for each time step. Land cover was classified for each of the 18 multi-date image composites using an unsupervised classification with the ISODATA clustering algorithm into 20 statistically different clusters, with a maximum of 75 iterations and convergence threshold of 0.95 for both time steps. When a cluster was determined to be indistinguishable between two or more classes, clusterbusting was used (i.e., clusters that could not be identified were extracted and a separate unsupervised classification was performed). Clusters were then combined into three thematic land cover classes, evergreen/mixed woodlands, deciduous woodlands, or other (comprised of urban, water, barren, etc.). A pixel was considered woodland if at least half of the pixel was covered by woody vegetation, which was assessed by visual observation. Classified images were then mosaicked for each time step.
Two independent accuracy assessments (hereafter referred to as the accuracy assessment performed by Interpreters (a) and (b)) were conducted for each mosaic by treating NAIP color infrared imagery (1-m resolution; 2004 and 2010) as the reference source. We used 2004 NAIP imagery to assess the accuracy of the 1999-2000 mosaic because 2004 is the closest (in time) available NAIP imagery that covers the entire warbler breeding range. Using the multinomial distribution, we determined that 106 and 102 pixels were needed to assess classification accuracy with an expected accuracy of 90% and precision of 10%, respectively, during both time steps (Tortora 1978, Congalton andGreen 1999). Measures of classification accuracy were overall classification accuracy, overall Kappa coefficient of agreement, and user's and producer's accuracy. User's accuracy is the probability that a pixel in the classified image represents the same land cover class on the ground, and is used as a measure of the reliability of the classification. Producer's accuracy indicates the probability that a reference pixel (i.e., the corresponding group of pixels in the NAIP imagery) was correctly classified. Thus, producer's accuracy provides information on how well a certain class/category has been classified.
Thematic land cover classes were then separated into two categories: ''habitat'' and ''other''. We classified habitat as any pixel that was originally classified as evergreen/mixed woodland as well as any pixel that was originally classified as deciduous woodland that was within 90 m of an evergreen/mixed woodland pixel. All remaining pixels were considered ''other''. We chose this warbler breeding habitat classification scheme because it has been shown to do well when compared to warbler detection/ non-detection surveys (N. Heger, personal communication). Habitat patches were then identified using an eight neighbor rule. All patches ,15 ha were not considered breeding habitat because it has been found that warblers need approximately 15-20 ha to successfully fledge young .
In post-classification change detection analyses, misclassification of pixels can be misleading if errors in classifications are confused with landscape changes in habitat. This may lead to biased estimates of habitat gains and losses. However, methods to account for these biases require knowledge of misclassifications of the change (Veran et al. 2012). We cannot determine with certainty that change in warbler habitat at the pixel level is actual change or misclassification since we are using a modeling approach to delineate warbler habitat from classified woodland maps. Therefore, in order to determine if the change we show here is driven by actual change, and not primarily a function of errors in classification, we ran two independent accuracy assessments (hereafter referred to the accuracy assessment performed by Interpreters (a) and (b)) on a change image between the two classified woodland map mosaics. Once again, we used the multinomial distribution with an expected accuracy of 90% (and a precision of 10%) and found it was necessary to check 158 pixels within the woodlands change image. We used the identical measures of classification accuracy as before. If change between evergreen/mixed woodlands, deciduous woodlands, and other between the two time steps has v www.esajournals.org reasonable accuracy then we assumed change in warbler habitat that were based off these woodland maps is actual change in warbler habitat, not compounded classification errors.
A post-classification change detection between the warbler habitat maps for each time step was conducted in Imagine to quantify how habitat changed across the landscape between the two time steps. Areas identified as ''converted to habitat'' were not considered habitat since these images were roughly 10 years apart and it likely takes several decades for woodlands to mature and be considered viable warbler breeding habitat (Reemts and Hansen 2008). Once these pixels were removed from the second time step's habitat map, patch sizes were reassessed to be sure all patches were 15 ha.
Habitat patch metrics were then quantified for each recovery unit in both time steps using FRAGSTATS v4 (McGarigal et al. 2012). Metrics were selected to account for class-specific area, shape characteristics, and aggregation and included mean habitat patch area, total habitat area, number of habitat patches, mean habitat patch shape index, and patch aggregation index. Mean habitat patch size, number of habitat patches, and total habitat area are the base metrics for habitat change that we considered. Mean patch area and total patch area provide information about the composition of the landscape and can be used to assess changes in habitat proportion over time. The number of habitat patches was selected to represent potential subdivision or aggregation of patches over time. An increase in the number of habitat patches combined with a decrease in mean habitat patch size indicates habitat patches are becoming fragmented. Whereas, a decrease in mean habitat patch size with a lack of change in the total number of habitat patches indicates some habitat patches are becoming fragmented, while other habitat patches are being entirely depleted. A decrease in total habitat area indicates habitat loss. Mean shape index is similar to perimeter-area ratio in that it is a measure of patch shape complexity, however, the mean shape index corrects for the size of the patch. The mean shape index is one when all habitat patches are circular and increases as habitat patches become more irregular. Notably, as a habitat patch becomes more irregular the core area of that patch decreases. Aggregation index characterizes the frequency with which like classes are adjacent and provides a spatial component to interpreting landscape fragmentation. Values are low (i.e., 0) when the frequency of like classes is low (indicating habitat fragmentation) and high (i.e., 100) when like classes are maximally aggregated. All habitat patch metrics were quantified using an eight neighbor rule.

RESULTS
Accuracy assessment results for the two land cover classifications are summarized in Table 2. The overall classification accuracy for the first time step was 93.4% or 83.2%, depending on the interpreter, (a) or (b). Producer and user accuracies for the evergreen/mixed woodland class for this time step were 90.3% and 93.3% and 70.3% and 86.7% for Interpreters (a) and (b), respectively. Kappa coefficients of agreement were 0.87 and 0.67. Accuracy assessments for the second time step demonstrated similar results, with overall classification accuracies of 93.1% and 86.3%. Producer and user accuracies for the evergreen/mixed woodland class for this time step were 84.6% and 95.7% for Interpreter (a) and 70.0% and 91.3% for Interpreter (b). Kappa coefficients of agreement for the second time step were 0.86 and 0.74.
Accuracy assessment results for the change detection between the two time steps are reported in Table 3. Overall classification accuracies were 80.4% and 70.1% for Interpreters (a) and (b) and Kappa coefficients of agreement were lower than for the individually classified maps in each time step (0.68 and 0.58). For Interpreter (a), producer and user accuracies for evergreen/mixed woodland to evergreen/mixed woodland were quite high; 70.0% and 100.0% respectively. However, these same percentages for Interpreter (b) were lower, with only 65.7% for producer accuracy and 82.1% for user accuracy.
As expected, the results indicate warbler habitat has changed (Fig. 2). From 1999-2001 to 2010-2011, the total amount of delineated warbler habitat decreased from 2,219,168 ha to 1,578,281 ha (29% loss). As mentioned above, we quantified habitat patch metrics for each v www.esajournals.org warbler recovery unit in each time step (Fig. 3). It is worth noting that changes in habitat patch metrics are not directly comparable between recovery units because the units are different sizes and differ in overall shape. Therefore, one should not expect each unit to have the identical habitat patch metrics. Instead, we should exam-ine how habitat patch metrics changed over time within a recovery unit.

DISCUSSION
We were able to quantify the change in rangewide warbler breeding habitat and describe Notes: Interpreter (a) accuracy assessments were performed by A. Duarte. Interpreter (b) accuracy assessments were performed by J. Jensen. Points were generated using random sampling and were assessed using visual interpretation. Ellipses indicate classes that were not represented in the assessments because of the low number of total pixels in those particular classes. Notes: Interpreter (a) accuracy assessments were performed by A. Duarte. Interpreter (b) accuracy assessments were performed by J. Jensen. Points were generated using random sampling and were assessed using visual interpretation.
v www.esajournals.org pertinent warbler habitat patch metrics for each of the species' recovery units over a roughly 10 year period. There was an estimated 29% loss in range-wide warbler breeding habitat from 1999-2001 to 2010-2011. Indeed, such a large reduction in breeding habitat is a function of the additive influence of direct warbler breeding habitat loss, our minimum habitat patch size criterion, and that we did not consider gained warbler breeding habitat as habitat. We chose these criteria because we set out to examine the dynamics of habitat patches that were large enough to support warbler reproduction, and it likely takes several decades for woodlands to establish and become viable warbler breeding habitat. In all recovery units, we found the total habitat area, mean habitat patch size, and habitat patch aggregation index decreased from the first time step to the second time step. Percent loss in habitat for recovery units one through eight was 30,26,27,38,23,36,36, and 17, respectively. Reduction in mean habitat patch size was more pronounced in recovery units five, six, and eight. Aggregation index had a greater decrease in recovery units five and six. Mean habitat patch shape index was variable, with the largest change occurring in recovery units five and seven. The total number of habitat patches showed little change in most recovery units (exceptions were units five, six, and eight). Collectively, these patterns suggest that warbler reproductive success and overall carrying capacity has been dramatically reduced v www.esajournals.org during this time period. Moreover, the patterns identified here raise concern regarding habitat connectivity. Alldredge et al. (2004) modeled warbler populations between the Fort Hood Military Reservation and the city of Austin and found movement between habitat patches was needed in order for populations to persist. Furthermore, reduced connectivity of habitat patches may lead to increased genetic differentiation among warbler populations and a decline in genetic diversity over time (Lindsay et al. 2008, Athrey et al. 2011). Yet, movement rates and distances have not been quantified for this species. Interestingly, a majority of warbler habitat conservation dollars have been spent conserving warbler breeding habitat on the outskirts of the city of Austin, Travis County (recovery unit five). Yet, our findings suggest warbler habitat is being removed and becoming more fragmented at equal, if not greater, rates (than in recovery unit five) on the outskirts of the city of San Antonio, Bexar County (recovery unit six) and throughout recovery unit eight. Given the scarcity of nonprivately owned warbler breeding habitat in recovery unit eight (Hatfield et al. 2012), the greater number of warbler territories that can be supported per hectare in the southern portion of the warbler's breeding range , and the relatively high amount of habitat loss and fragmentation occurring in recovery unit eight, future warbler habitat conservation efforts focused in unit eight might be a more effective range-wide breeding habitat conservation strategy in order to preserve large breeding habitat patches that support a greater number of warblers while they are still available.
Notably, our image for the second time step classifies woodlands in central Texas, and by extension potential warbler breeding habitat, with high accuracy. Thus, it can serve as yet another warbler breeding habitat map to be used in large-scale warbler conservation planning. Each of the more recent warbler breeding habitat maps (i.e., the one produced here, Diamond et al. 2010, andCollier et al. 2012) uses a classification scheme that is considered acceptable to delineate warbler breeding habitat. Yet, each warbler habitat map has its own advantages, disadvantages, and inherent assumptions. Thus, careful consideration should be conducted to determine which warbler habitat map better meets the requirements of a specific individual conservation objective in order to decide which warbler habitat map to use. Or perhaps all warbler habitat maps could be used in a multi-model inference framework (e.g., ensemble models weighted by differences in imagery age, target metric [occupancy, reproductive success, etc.], and/or other planning value) to identify breeding habitat patches of high conservation value. However, it should be noted that each of these warbler habitat maps (including the one presented here) is based off imagery that can be considered somewhat outdated. The Diamond et al. (2010) warbler habitat map is derived from the vegetation classifications done as part of the Texas Ecological Systems Classification Phase 1. The images used in this classification were acquired from 2005 to 2007. The Collier et al. (2012) warbler habitat map is based off 2007-2008 Landsat TM images. Although we recognize there will always be time lags between satellite image acquisition dates and image classifications, especially at the scale of the warbler's entire breeding range, all warbler breeding habitat maps are based on imagery taken prior to a significant, landscape altering event. In the summer of 2011, central Texas experienced severe drought conditions that resulted in die offs of many oak and juniper trees (Kukowski et al. 2013). Thus, current warbler habitat maps may have large errors of commission (i.e., false positives), and should be updated, when feasible, for more reliable large-scale warbler conservation planning to take place.
Of particular importance with regard to these results are the differences in classification accuracy and sources of uncertainty that are introduced to the analysis at different stages. For example, though every effort was made to download temporally consistent satellite image tiles, minor differences in green-up and senescence attributed to temperature or timing of precipitation may have resulted in mislabeled clusters for the individual time steps. Additionally, the use of NAIP imagery as the reference dataset, presents particular challenges with regard to visual interpretation and its use as a reference source. For instance, NAIP imagery for a region is frequently acquired over the course of days to weeks and during inconsistent satellitesun geometries, whereby shadows and shading can be introduced to the final image product such that overall species composition looks more dense or spectrally similar to different land cover classes (i.e., Ashe juniper and broadleafed evergreen species such as live oak).
Another consideration with regard to land cover classification is the differentiation between forest/woodland and shrubland. Ashe juniper shrublands may appear spectrally similar particularly to the evergreen/mix woodland class, though the difference is based on height, an attribute not available with Landsat imagery. Thus, individual class accuracies may have been v www.esajournals.org influenced by interpreter-specific identification of woodland and shrubland and may have led to artificial gains and losses of warbler habitat. Using LiDAR technology to remotely measure vegetation height could help distinguish woodlands from shrublands and has been shown to do well in predicting warbler occurrence at the site level (Farrell et al. 2013). Still, these types of data are not available for historic imagery spanning the warbler's entire breeding range. Further, optimal LiDAR based vegetation height metrics that accurately represent forest canopy characteristics at a 30 m resolution have not been determined (Jensen et al. 2013). Therefore, LiDAR based vegetation metrics were not included in this classification scheme.
Regardless, overall classification accuracies of 93.4% (or 83.2%) and 93.1% (or 86.3%) are high, and the overall Kappa coefficient of agreements (0.87 [or 0.67] and 0.86 [or 0.74]) indicated we had moderate-to-strong agreement between our classified maps and the reference dataset. Further, we had decent classification accuracy for the woodland change image, particularly given the geographic extent of the image. Collectively, this indicates we were able to quantify actual change in warbler breeding habitat over time.
It should be noted, however, that in the context of landscape change analyses, land cover classification, habitat identification, and the selection of metrics used to characterize change over time must be considered. Our selection of landscape metrics was guided by the need to quantify changes in warbler breeding habitat area, patch shape, as well as potential fragmentation within each of the species' recovery units. Further, we selected metrics that we considered would be easily interpretable to a wide audience of resource managers and policy makers, regardless of biogeographic, ecological, or geospatial expertise and training. Nonetheless, selection of metrics may in some cases influence how landscape change is quantified and interpreted, particularly when the analysis is focused on characterizing or understanding underlying landscape processes and species modeling.
In summary, we found that there was a large reduction in range-wide warbler breeding habitat and that warbler breeding habitat was removed and became more fragmented at uneven rates across the warbler's breeding range. This infor-mation will assist researchers and managers in prioritizing breeding habitat conservation efforts across recovery units, and highlights the significance in estimating warbler movement rates and distances in order to assess the degree of connectivity between increasingly fragmented habitat patches. Further, it would be useful to incorporate these rates of habitat change into warbler population models in order to extrapolate how long the species can be sustained given current habitat dynamics.