Environmental change, shifting distributions, and habitat conservation plans: A case study of the California gnatcatcher

Abstract Many species have already experienced distributional shifts due to changing environmental conditions, and analyzing past shifts can help us to understand the influence of environmental stressors on a species as well as to analyze the effectiveness of conservation strategies. We aimed to (1) quantify regional habitat associations of the California gnatcatcher (Polioptila californica); (2) describe changes in environmental variables and gnatcatcher distributions through time; (3) identify environmental drivers associated with habitat suitability changes; and (4) relate habitat suitability changes through time to habitat conservation plans. Southern California's Western Riverside County (WRC), an approximately 4,675 km2 conservation planning area. We assessed environmental correlates of distributional shifts of the federally threatened California gnatcatcher (hereafter, gnatcatcher) using partitioned Mahalanobis D 2 niche modeling for three time periods: 1980–1997, 1998–2003, and 2004–2012, corresponding to distinct periods in habitat conservation planning. Highly suitable gnatcatcher habitat was consistently warmer and drier and occurred at a lower elevation than less suitable habitat and consistently had more CSS, less agriculture, and less chaparral. However, its relationship to development changed among periods, mainly due to the rapid change in this variable. Likewise, other aspects of highly suitable habitat changed among time periods, which became cooler and higher in elevation. The gnatcatcher lost 11.7% and 40.6% of highly suitable habitat within WRC between 1980–1997 to 1998–2003, and 1998–2003 to 2004–2012, respectively. Unprotected landscapes lost relatively more suitable habitat (−64.3%) than protected landscapes (30.5%). Over the past four decades, suitable habitat loss within WRC, especially between the second and third time periods, was associated with temperature‐related factors coupled with landscape development across coastal sage scrub habitat; however, development appears to be driving change more rapidly than climate change. Our study demonstrates the importance of providing protected lands for potential suitable habitat in future scenarios.

Together, these changes in geographic distributions can lead to the local or regional extinction of species (Thomas et al., 2004) or the generation of novel communities (Hobbs et al., 2006;Ohlemüller, Walker, & Wilson, 2006;Williams & Jackson, 2007). Furthermore, there are multiple environmental factors that can influence a species' distribution, and recent research has shown that a species' ability to respond to climate change may be affected by multiple natural and anthropogenic factors (e.g., Brook, Sodhi, & Bradshaw, 2008;Swab, Regan, Keith, Regan, & Ooi, 2012).
Ecological niche modeling can provide insights into consequences of environmental change on biodiversity (Barrows & Murphy-Mariscal, 2012;Thomas et al., 2004), and it is increasingly used to evaluate the responses of biodiversity to changes in climate and other environmental attributes (e.g., Dawson, Jackson, House, Prentice, & Mace, 2011;Pereira et al., 2010;Thomas et al., 2004). Results from niche models can be used to plan for priority area selection for conservation purposes, and many habitat conservation plans created to mitigate climate and environmental change are based on future geographic ranges (summarized in Peterson, 2006). However, there can be substantial uncertainty when predicting future species' distributions due to the uncertainty of future climatic conditions, potential patterns of landscape development, changing species interactions, and even species' adaptations to new conditions. Nevertheless, many species have already undergone shifts in distribution due to habitat and climatic changes, and we can use past shifts to better understand the influence of environmental factors on a species' range and to evaluate, create, and/or adapt habitat conservation plans, to the extent that past geographic shifts show us the trajectory of change. Furthermore, analyzing past shifts may provide us with the ability to evaluate the effectiveness of existing conservation protections.
However, this region continues to experience a multitude of anthropogenic disturbances, including but not limited to urbanization, air pollution, the introduction of non-native species (summarized in Lovich & Bainbridge, 1999), and climate change. We analyzed past habitat associations and distributional shifts of the federally threatened California gnatcatcher (Polioptila californica; hereafter gnatcatcher; see Figure 1 for gnatcatcher photo) during different time periods within southern California's Western Riverside County (WRC) planning area using ecological niche modeling. We chose the gnatcatcher as a case study due to an extensive observation and survey database and because it is embedded in a regional conservation plan for WRC. The gnatcatcherpreferred vegetation type is coastal sage scrub (CSS; Atwood, 1993).
CSS communities are among the most endangered habitats in the United States with estimated losses of 60%-90% since the start of the 20th century (O'Leary, 1995). Within the remaining intact CSS, many areas have been invaded by exotic grasses (Minnich & Dezzani, 1998) through the combination of nitrogen deposition from urban southern California (Allen, Padgett, Bytnerowicz, & Minnich, 1998;Padgett & Allen, 1999) and fire (Cox, Preston, Johnson, Minnich, & Allen, 2014).
As shrub cover and shrub species diversity decline, critical CSS habitat for more than 200 plant and animal species that are currently endangered, threatened, or of "special concern" is compromised (Bowler, 2000).
To better understand the influence of past environmental changes and the importance of conservation lands in fulfilling a mandate to protect this species, we had four aims: (1) to quantify regional habitat associations of the gnatcatcher; (2) to describe changes in environmental variables and gnatcatcher distributions through time; (3) to identify environmental drivers that are associated with habitat suitability changes; and (4) to relate habitat suitability gains and losses to habitat conservation plans. Because our analyses covered pre-and postconservation planning periods, we expected to see shifts in the distribution of gnatcatchers associated with overall deteriorating environmental conditions that helped trigger conservation action in the first place. As there are many known environmental stressors within WRC, we expected to see a suite of environmental variables as drivers for shifting distributions. Nevertheless, we anticipated that urbanization would be a stronger driver than climatic variables due to greater pace of change of the former compared to the latter. We also expected that acquired conservation lands would help mitigate the negative effects of environmental change.

| METHODS
Our study location was southern California's WRC (Figure 2), at the northern extent of the California gnatcatcher's geographic range ( Figure 1). Rapid urban development beginning in the 1970s, overlaid on extensive agricultural development that started in the late 1800s, led WRC to develop one of the first multiple species habitat conservation plans in the state in 2004 (planning initiated in 1997, plans established by 2004) to protect a growing number of species, including the gnatcatcher, listed as rare, threatened, endangered, or sensitive by both the state and federal governments (Preston & Rotenberry, 2007;Preston, Rotenberry, Redak, & Allen, 2008;Western Riverside County 1997 ca.gov/Data/CNDDB/Maps-and-Data). We only used records with spatial precision of <125 m. We deleted all spatially redundant records (observations within the same grid cell) for each time period independently.
We split the observational data into three time periods: (1) 1980-1997, (2) 1998-2003, and (3) 2004-2012. Dates were chosen to represent both land management and ecological changes in WRC. The first time period was essentially preconservation planning.
Native habitat, especially CSS, began rapidly disappearing due to conversion to urban and suburban development and expanding agricultural development beginning in the mid-1980s (O'Leary, 1995.

| Partitioned Mahalanobis D 2 niche modeling
Mahalanobis D 2 is a niche modeling technique based on the standardized difference between the multivariate mean for environmental variables at locations where a species is detected relative to the values for these same environmental variables at any point in the region being modeled (Clark, Dunn, & Smith, 1993;Rotenberry, Knick, & Dunn, 2002;Rotenberry, Preston, & Knick, 2006). The more similar in environmental conditions a point is to the species' mean, the smaller the D 2 and the more "suitable" the habitat at that point is assumed to be. Habitat similarity index (HSI) values are derived from D 2 values and are generally rescaled to range from 0 to 1 (Clark et al., 1993). An HSI of 1 represents environmental conditions identical to the species' mean; HSI of 0 represents conditions most dissimilar. D 2 assumes that at least some environmental variables influencing a species' distribution have been included in the model, and it performs reasonably well in identifying suitable habitat based on species presence only (Knick and Dyer 1997;Knick & Rotenberry, 1998).
Using eigenvector analysis, Mahalanobis D 2 can then be partitioned into independent, additive components that represent independent relationships between a species' distribution and environmental variables (detailed in Rotenberry et al., 2002Rotenberry et al., , 2006. A particular advantage to using partitioned D 2 is that it focuses on variables that have a relatively low variance across a species' occurrences; variables maintaining a consistent value where a species occurs (and hence with relatively low variance) are most likely to be associated with factors limiting a species' distribution, especially compared to those taking on a wider range of values (Dunn & Duncan, 2000;Rotenberry et al., 2002). The number of partitions for a model equals the number of variables included in that model. To choose the appropriate partition that best represented the gnatcatcher's distribution, we calculated the median habitat similarity index for gnatcatcher-occupied cells for each of the 12 partitions in our model, then chose the partition associated with the largest of those medians as the best partition (Rotenberry et al., 2006). We then used the selected partition to calculate HSI values for every cell in the landscape.
F I G U R E 1 (Upper) California gnatcatcher (Polioptila californica) and (lower) its geographic range. Map adapted from birdphotos.com with original data acquired from Ridgely et al. (2005). Photo courtesy of Mark A. Chappell

| Modeling the California gnatcatcher through time
We selected habitat variables based on previous extensive literature reviews and validation procedures created for WRCMSHCP that are outlined in Preston and Rotenberry (2007) and Preston et al. (2008).
There were a total of 12 environmental variables in our models (Table 1) If a proportional change decreased a habitat below zero, we assumed a zero percent cover of that habitat; if a proportional change increased a habitat over 100, we assumed a 100% cover of that habitat. We note that this calculation assumes that changes in development affect each vegetation type to the same degree. We then calculated the WRC-wide mean and standard deviation of all 12 environmental variables for each time period (Table 2).
We modeled gnatcatcher habitat within each time period using the partitioned Mahalanobis D 2 applied to the 12 environmental variables for each period. We modeled each time period separately to account for any potential changes in the variables that might limit gnatcatcher distribution. Calculations were carried out in SAS (SAS Institute 2001) using SAS code from Rotenberry et al. (2006). We then validated the models quantitatively using the evaluate() command in the {dismo} package in R (Core Team, 2012). This analyses cross-validates models with presence only or presence/absence data.
Given a vector of presences and a vector of absences (or in our case pseudoabsences generated by evaluate()), and a vector of HSI values, confusion matrices are computed (for varying thresholds), and model evaluation statistics (AUC values) are computed for each confusion matrix/threshold. Models are considered robust if AUC is ≥0.70 (Fielding & Bell, 1997).
Once we validated models, within each time period, we calculated HSIs for every cell in the landscape. For visualization, using ArcGIS (ESRI 2012) we mapped each cell as either highly suitable (HSI ≥ 0.66), moderately suitable (0.33 < HSI < 0.66), or minimally suitable (HSI ≤ 0.33), following categories defined by Preston and Rotenberry (2007) and Preston et al. (2008). We then simplified our HSI categories into habitat that was either suitable (HSI ≥ 0.66) or unsuitable (HSI < 0.66) and created suitability change maps representing each cell's suitability transition from the first to second time period and the second to third time period. Thus, we had four transition categories: (1) cells that maintained high habitat suitability (2) cells that remained unsuitable, (3) cells that lost suitability, and (4) cells that gained suitability. Finally, we calculated the number of cells that were highly suitable or unsuitable in each time period in order to better understand the amount of area that gained or lost suitability.

| Relating environmental variables to suitability changes within grid cells
For environmental analyses, we included three additional environmental variables that are known stressors within WRC: exotic plant cover, nitrogen deposition, and fire frequency (Table 1).
Nitrogen deposition values were calculated on a 4-km grid using the Community Multiscale Air Quality model (Tonnesen, Wang, Omary, & Chien, 2007 share similar values for a variable than points further apart) creates a lack of statistical independence among cells and thus invalidates most normal tests of statistical significance (e.g., Dale & Fortin, 2002). Intuitively, because some sample points are not fully independent of one another, our "effective" sample size is less than the nominal sample size (n), increasing the chance of Type I error. To compensate, we calculated an effective sample size (n e ) following Griffith (2005;eqn. 3); for each test, we performed assuming a very high level of spatial dependency (e.g., ρ = 0.8, for a simple autoregressive model of y = ρCy + ε, where y is a variable of interest, C is a matrix of spatial weights [e.g., inverse distances], and ε represents error). Although this will lead us to underestimate statistical significance for variables with low spatial autocorrelation, this approach will provide some assurance that any apparently statistically significant results we observe are indeed significant.
Sample means and standard deviations are calculated in the usual way, then n e is used instead of n in calculating standard errors and, hence, t-statistics (Dale & Fortin, 2002). We also compared environmental variable means using the method described above between habitat that gained or lost suitability through time based on protection level.

| Assessing the effectiveness of prioritized conservation lands
Conservation lands included properties within WRC owned, man-

| Modeling the California gnatcatcher through time
For the WRC study area as a whole, temperatures across time periods only varied about 0.2°C for maxima and less than 0.1°C for minima (Table 2a). On the other hand, precipitation varied substantially, declining nearly 100 mm between the first and second periods, then increasing by nearly 50 mm between the second and third. All of the land cover classes changed substantially through time: coastal sage scrub declined 17% between the first and third periods, chaparral −6%, grassland −26%, and agriculture −36%. The most dramatic change was in developed land, which increased over 900%, most of which occurred between the first and second periods (Tables 2a,b; Fig. S1).
There were 480, 528, and 554 unique grid cells where gnatcatchers were observed for the first, second, and third time periods, respectively. Cells occupied by gnatcatchers differed significantly from the landscape as a whole with respect to a number of environmental variables (Table 2b). Occupied cells were consistently warmer and drier and at lower elevations. They also contained more coastal sage scrub but less chaparral, with higher exotic cover and generally higher total nitrogen deposition. They had less than average agriculture during the first period (when agriculture was highest) and less than average development during the third period (when development was highest). Characteristics of occupied cells were also largely consistent in composition from one year to the next, differing significantly among periods only in precipitation and development, the latter only between the first and second periods (  (Table 3); however, during each time period, more suitable habitat became cooler and higher in elevation, particularly comparing the third period to the first.
Suitable maximum temperature varied little among time periods, and the relationship between mean highly suitable HSI and precipitation appeared to fluctuate with precipitation's natural variation.
In one period, highly suitable habitat was significantly more southfacing, and in two periods, significantly more west-facing (Table 3).
Highly suitable habitat consistently had substantially more coastal sage scrub, much less chaparral, less grassland and agriculture, and higher exotic cover (Table 3). In one period, it was associated with more historical fires and less deposited nitrogen, but not in the other two (Table 3). Although high HSI cells were associated with relatively  (Table 3). Overall, results from comparing high-versus low-suitability cells largely paralleled those comparing occupied to unoccupied cells.

| Changing habitat suitability through time
The gnatcatcher lost 11.7% and 40.6% of cells with highly suitable habitat within WRC between time period transitions (Table 4a,b).
Visually, we can see the reduction in highly suitable habitat through time (Figure 3a-c), and its shift toward the southeast section of the study area (Figure 3d,e). The largest numbers of grid cells regardless of their starting and middle suitability transition states were those considered unsuitable in the third time period (Fig. S2). Of the cells that were unsuitable during 1980-1997 and became suitable during 1998-2003, 76% lost their suitability during 2004-2012 (Fig. S2).
Only 1,918 grid cells (20.7%) maintained high suitability through all time periods out of the initial 9,258 grid cells that were suitable in the first time period (Fig. S2).

| Assessing the effectiveness of conservation lands
During the first transition period, the number of highly suitable cells declined 18% within protected areas, but only 6% in unprotected ones ( T A B L E 1 Environmental variables included in the niche models and subsequent analyses. The first 12 environmental variables were used in the ecological niche modeling. All environmental variables were included in additional analyses precipitation, less agriculture, and less chaparral (Table 5). Unprotected cells that gained compared to losing suitability were also dryer and had less agriculture, but had more coastal sage scrub. During the second transition period, changes in suitability both inside and outside protected areas were similarly driven by climate; compared to losses, gains were associated with warmer minima, cooler maxima, and less precipitation, and gains were associated with less development and more coastal sage scrub inside and out ( Table 5). Cells that gained suitability in protected areas had more chaparral, whereas those outside protection had less agriculture.

| DISCUSSION
Our study shows the utility of using ecological niche models for understanding how environmental variables are related to changes in T A B L E 2 (a) Environmental variable means ± standard deviations calculated across the entire WRC (N = 74,832 250 m × 250 m cells). See Table 1 for explanation of environmental variables. (b) Environmental variable means ± standard deviations calculated for cells within WRC that were occupied by California gnatcatchers. We present means and ordinary standard deviations based on n; however, we use standard errors and degrees of freedom based on n e . *p < .01, 1-sample t test compared to landscape mean (a) habitat suitability through time. Further, we provide evidence for the importance of creating habitat plans that maintain adequate protected lands. Our data described a species' changing distribution in response to environmental change and provide a background for creating and adapting conservation plans to begin to accommodate future responses to environmental change. Over the past four decades, the gnatcatcher has lost large amounts of suitable habitat within WRC especially between the second and third time periods. The change in suitable habitat was driven mainly by landscape development across CSS habitat, but also yielded some changes in climate-related associations.
Gnatcatchers occupied cells that varied less (both through time and within occupied cells) for maximum temperature than minimum temperature, and we did not find a consistent temporal trend or significant difference in maximum temperature between cells that gained or lost suitability. However, we found significant differences between mean minimum temperatures among cells that gained or lost suitability. Thus, when predicting future geographic responses to temperature changes for the gnatcatcher, it is crucial to allow flexibility (plasticity) with minimum temperature, but restrict maximum temperature for future scenarios. We note that based on the changing response curve to the relationship between habitat suitability changes and the environmental variables through time, the gnatcatcher may also be adapting to environmental change or have a high plasticity for those variables. Additional studies are needed to disentangle adaptation from plasticity in the gnatcatcher's response to environmental change.
Although the loss of suitable cells within protected areas outpaced gains, the rate of loss was manifestly less than that outside. Thus, it is clear that gnatcatcher habitat is better maintained within protected areas than outside, but only insofar as it has experienced a less rapid decline inside than out. Changes in both climate and land use/land cover were associated with changes in habitat suitability inside and outside protected areas, but no easily interpretable differences between the two were apparent. It seems therefore that processes that modify suitability are generally the same outside versus inside, but are more extensive outside protected areas.
The mean minimum temperature of suitable habitat declined slightly through time; the lower minimum temperature is likely to be a latent effect of the movement toward higher elevation. Although this is consistent with literature showing that as temperatures rise due to climate change, bird species are known move up in elevation to stay within their physiological tolerances (Tingley, Monahan, Beissinger, & Moritz, 2009), it seems more likely that movement toward higher elevation may be a by-product of displacement as lower elevational landscapes are developed with concomitant loss of CSS. Both maxima and minima temperatures remained within the known temperature constraints of the gnatcatcher (Mock 1998). Interestingly, lower temperatures at higher elevations are likely to pose a stronger distributional limit to gnatcatchers than higher temperatures at lower eleva-

tions. Range limits in otherwise suitable vegetation types in southern
California appear to be associated with average January minima of about 2.5°C (Mock 1998), only about 2°C below what we observed.
On the other hand, summer maxima in gnatcatcher-occupied areas in Baja California Sur routinely exceed 35°C (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), well above those observed in our study area.
The relationships among the variables we considered in our analyses are complex. One important feature is that things such as topography, nitrogen deposition, fire history, and agriculture/development exert their effects on gnatcatchers only indirectly, via their direct effects on vegetation type. Gnatcatchers do not seem to be edge-sensitive (Bolger, Scott, & Rotenberry, 1997), and hence, effects of agriculture and development are mainly manifest through habitat loss. Nevertheless, even indirect effects can be rapid and profound, mediated through sudden impacts of land-use change and/or fire on vegetation type, including complete replacement of one land cover type by another. A further complication from a modeling and changedetection perspective is that topography is a fixed effect (no temporal variance, at least over relevant time scales), yet is a strong driver of local and regional patterns of climate and land use. Indeed, in our analysis, increasing elevation appears to be associated with changes in habitat suitability, but this is most like an expression of the fact that most development and agriculture have occurred below 600 m, where most CSS historically occurred as well. And although patterns T A B L E 4 (a) Number of grid cells that were highly suitable based on whether cell was within protected conservation area and (b) number of grid cells that gained or lost high suitability based on whether cell was within protected conservation area  (Preston et al., 2008).
Two things are shifting through time: the distribution of values of environmental variables (e.g., precipitation, temperature, CSS, development) and the distribution of California gnatcatchers in response to some of those variables. Although the three partitioned Mahalanobis D 2 niche models that capture bird-environment relationships were fundamentally similar, they differed in some attributes between periods. We do not believe that these represent "niche shifts" or local adaptive changes in gnatcatcher habitat relationships, but rather suggest that the nature of the limits to the distribution of the species in this region may shift between periods, compounded by lags in the response of gnatcatchers to environmental change (i.e., "ghosts of habitats past"; Knick & Rotenberry, 2000).
Such "ghosts" are most likely to occur when largely sedentary, permanent resident species such as the California gnatcatcher (Atwood & Bontrager, 2001) 1980-1997 and 1998-2003 and (e) 1998-2003 and 2004-2012 dynamics and may change a species' distribution . Recent studies by Swab et al. (2012) and  found that predicted range shifts induced by climate change were only part of the threat for the future outcome of two plant species, and other environmental stressors may play an equal or larger or interactive role in the survival of a species. Habitat loss has been the greatest threat to biodiversity, at least in the near term (Brooks et al., 2002;Groom & Grubb, 2006;Hanski, 2005). This factor was evident in the impact of expanding development and conversion of CSS within southern California on the California gnatcatcher. While the gnatcatcher lost substantial highly suitable habitat throughout the region, our study found that land outside of conservation protection by the WRCMSHCP lost more suitable habitat than habitat within conservation lands over the course of our study.
These findings emphasize the importance of providing protected areas T A B L E 5 Mean difference (±standard deviation) of grid cells that gained or lost high suitability based on whether cell was within protected conservation area for the transitions between 1980-1997 to 1998-2003, and 1998-2003 to 2004-2012. We present means and ordinary standard deviations based on n; however, we use standard errors and degrees of freedom based on n e . *p < .01, 2-sample t test comparing mean gain to mean loss beyond the present distributions, as these may mitigate at least some of the negative effects of environmental changes that do not respect political boundaries. Changing climate and the development of land due to lack of protection outside the initial critical habitat may be crucial to the long-term protection of the gnatcatcher and other species.
While southern California is considered a biodiversity hotspot, the remaining CSS especially is also experiencing a multitude of anthropogenic stressors that differentially affect an already patchy habitat type. Because of this, prioritizing future conservation lands is critical to successfully manage habitat for plant and animal populations.
During our study period, WRC adopted a multiple habitat and species conservation plan in 2004 with a goal of conserving 146 sensitive plant and animal species and their habitats (Western Riverside County Regional Conservation Authority 2017). Based on the trends of the nonprotected lands losing suitability in the most recent period, our data suggest that the gnatcatcher would have lost even larger amounts of suitable habitat without the prioritization of conservation lands.
Our study emphasizes the importance of creating habitat conservation plans and prioritizing conservation lands as well as protecting project future habitats are crucial maintain populations during periods of high environmental stress.