A habitat-based approach to determining the effects of drought on aridland bird communities

Aridland breeding bird communities of the United States are among the most vulnerable to drought, with many species showing significant population declines associated with decreasing precipitation and increasing temperature. Individual breeding bird species have varied responses to drought which suggests complex responses to changes in water availability. Here, we evaluated the influence of water deficit, an integrative metric of drought stress, on breeding bird communities within three distinct aridland habitat types: riparian, pinyon-juniper, and sagebrush shrubland. We used 12 years of breeding bird survey data from 11 National Parks and Monuments in the Northern Colorado Plateau Inventory and Monitoring Network (NCPN). We used a multivariate community level approach to test for the effect of annual water deficit on the bird communities in the three habitats. We found that bird communities responded to water deficit in all three habitat types, and 70% of the 30 species-habitat combinations show significant relationships between density and variation in water deficit. Our analyses revealed that the direction and magnitude of species responses to water deficit was habitat dependent. The habitat specific responses we observed suggest that the adaptive capacity of some species depends on the habitat in which they occur, a pattern only elucidated with our habitat-based approach. The direction and magnitude of the relationships between predicted densities and annual water deficit can be used to predict the relative sensitivity of species within habitat climate changes. These results provide the first attempt to determine how the indirect effect of changes in water availability might affect aridland breeding birds in distinct habitat types. Linking breeding bird density to annual water deficit may be valuable for predicting changes in species persistence and distribution in response to climate change.The code and model to produce the water deficit values are not included here because they are in the process of being published.Study Design and Habitat Descriptions We used 12 years (2005-2015 and 2017) of breeding bird survey data to estimate bird densities in three habitat types (riparian, pinyon-juniper, and sagebrush shrubland) in 11 National Parks and Monuments of the National Park Service's Northern Colorado Plateau Inventory and Monitoring Network (Figure 1). We identified pinyon-juniper and sagebrush shrubland habitat types on land-cover maps from the Southwest Regional Re-GAP Analysis Project (Lowry et al. 2005) or using vegetation-association maps (Daw et al. 2017). Riparian habitat occurs along perennial streams as narrow strips of habitat, surrounded by dry uplands, and contains multiple layers of canopy, which may experience different degrees of dryness depending on ground water levels and rooting depths. The riparian habitat is dominated by willow (Salix spp.), Tamarisk (Tamarix ramosissima), and isolated stands of Fremont cottonwood (Populus fremontii) and boxelder (Acer negundo). Pinyon-juniper habitat is dominated by two-needle pinyon (Pinus edulis) and junipers (Juniperus spp.), which varied in relative abundance. The shrub layer in pinyon-juniper habitat varies throughout the NCPN, but is often dominated by sagebrush (Artemisia spp.), mountain mahogany (Cercocarpus spp.), jointfir (Ephedra spp.), and cliffrose (Purshia spp.). Sagebrush shrubland habitat is dominated by sagebrush, primarily big sagebrush (Artemisia tridentata) and prairie sagewort (A. frigida), with rabbitbrush (Chrysothamnus spp.), greasewood (Sarcobatus spp.), and other shrub species interspersed. Riparian transects ranged in elevation from 1,283 – 1,901 m, pinyon-juniper transects ranged in elevation from 1,393 – 2,402 m, and sagebrush shrubland transects ranged in elevation from 1,666 – 2,447 m. Transects within each habitat type were established across a large latitudinal gradient to allow for habitat comparisons and to guard against variation that occurs across broad geographic gradients (Figure 1). Water Deficit Estimates Water deficit is a measure of drought stress and is the amount of additional water vegetation would use if it was available (Stephenson 1998). The use of water-year, defined as October 1 – September 30 of the following year, in the estimation of water deficit is routine in studies of vegetative responses to climate in the desert Southwest, as it helps account for precipitation legacies that influence plants during the growing season (Reichmann et al. 2013, Bunting et al. 2017, Thoma et al. 2019). We used a monthly water balance model to estimate water-year water deficit using temperature and precipitation at the center of 45 bird survey transects (see below) following the methods of Lutz et al. (2010). Precipitation was partitioned into soil moisture, the quantity of water stored in the top meter of soil at the end of each month, and runoff which included overland flow plus infiltration below the rooting zone. Runoff is the proportion of water that is not available for plant growth. Maximum storage in the top meter of soil was defined by water holding capacity obtained from soil surveys (Soil Survey Staff 2019). Potential evapotranspiration (PET, mm) was the amount of water that could be evaporated or transpired with available energy if soil moisture was unlimited. Actual evapotranspiration (AET, mm) was the estimated monthly loss of water from soil via evaporation and transpiration, limited by availability of soil moisture. Water deficit (mm) was calculated as the difference between PET and AET (Stephenson 1998). We used Daymet daily temperature and precipitation data at a 1-km grain as the climatic input to the model (Thornton et al. 2016). These data were co-located with bird point count transects and thus represent the local elevation and latitude effects on precipitation and temperature at 1-km resolution. Within each 1-km grid cell of a temperature and precipitation time series the water balance model calculated heat load due to slope and aspect obtained from 30 m digital elevation model at the center of each transect (U.S. Geological Survey, 2017). Heat load was calculated as a scaling factor used to adjust PET up on south aspects or down on north aspects, thus accounting for slope and aspect interactions (Lutz et al. 2010 after McCune and Keon 2002). Breeding Bird Community We sampled the breeding bird community at 675 unique sampling locations on 45 transects during the breeding season (May 1st through July 15th), the timing of which varied by transect location. The elevation and latitudinal position of each transect was considered during the scheduling of each field season such that all transects were surveyed within their peak breeding season and after the passage of migratory birds. Transect orientation in pinyon-juniper and sagebrush shrubland habitats were randomly determined following protocol procedures outlined in Daw et al. (2017), the riparian transects surveyed in this study consisted of narrow strips of habitat, requiring transect placement to be centered along the length of the habitat and to follow the natural orientation of the riparian zone. All riparian habitat was narrow enough that the entire width of the riparian zone was included in the point radius of each survey (see below) and that dry upland habitat was also included in surveys of riparian habitat at some locations. Each transect consisted of 15 point count locations, with points spaced 250 m apart. Observers visited one transect per day, and initiated sampling within 30 minutes before sunrise and completed the survey within five hours after sunrise. Observers walked the transect, navigated to each point using a GPS unit and completed a five-minute point count survey. We did not conduct surveys if winds exceeded a four on the Beaufort Scale (13 – 18 mph) or precipitation was more than a drizzle. During each point count, observers recoded the species and distance (m) to all individuals detected. Observers estimated the distance to each bird or cluster of birds using a laser rangefinder (Simmons LRF 600; Simmons Outdoor Products, Overland Park, KS, USA). From 2005 – 2013, field teams surveyed each transect twice annually. Starting in 2014, survey effort was reduced to one visit per year because a sufficient number of detections to estimate density had been obtained for most bird species detected and to reduce costs. The number of completed surveys varied by year, but >75% of possible surveys were completed annually. For more detailed information about survey protocols, refer to Daw et al. (2017). We first estimated transect-level detection-adjusted densities for each bird species using the Distance package (Miller et al. 2019) in program R (3.6.1, http://www.r-project.org, accessed July 10, 2019). To avoid including double-counted individuals in our analyses, we only included detections within 125 m from the point center in our analyses. The total area covered by each transect was 73.6 ha (i.e. 15, 125 m radius points per transect). Next, we used the detection-adjusted density estimates in a multivariate analysis to determine the bird community response to changes in water deficit using the manyglm function in the R package mvabund (Wang et al. 2012). The mvabund package provides a novel set of hypothesis testing tools that is a flexible and powerful framework for analyzing multivariate abundance data (Wang et al. 2012). We used the manyglm function to simultaneously fit general linear models to each species using water deficit as the common predictor variable (Wang et al. 2012). This model-based approach alleviates the mean-variance relationship problem associated with distance-based community level metrics (Warton et al. 2012). To meet the data structure requirements of mvabund, we sought to maximize the number of transect-level density estimates available using Distance by constraining our analyses to only include the 10 species within each habitat type with the greatest number of transect-level density estimates available (Table 1). Species with the greatest number of detections allow for the greatest number of transect-level density estimate, therefore, the species included in each habitat are not the species with the highest densities, rather the species with the greatest number of detections. For species-transect-year combinations where detections were too few to estimate density, we summed the unadjusted counts and divided by the average detection probability for that species-year. This was done for 6% of the density estimates to satisfy the need for a response variable for all species-transect year combinations. Because there is likely a time lag in the bird community response to water deficit driven changes in vegetation, we compared two time-lag models (a one-year lag and a two-year lag effect of annual water deficit) against a null model within each habitat and selected the model with the lowest Akaike Information Criterion (AICc, Burnham and Anderson 2002). We did not explore other lags because vegetation response to legacy effects of precipitation shortfalls revert to average condition within two years in semiarid grasslands (Thoma et al. 2016). Once we determined the most appropriate time lag within each habitat type, we then tested for community-wide responses to annual water deficit within each habitat using ANOVA on the 'manyglm' object and assessed the resulting likelihood ratio tests (LRT; Warton 2011) and resampled P-values (Wang et al. 2012). If we detected a significant community level response to water deficit within a habitat type, we then used p.uni = "adjusted" argument to determine the effect size and direction of the response for each species. Finally, we predicted the relationships between species' densities and annual water deficit using the 'predict' function on the best model within each habitat type