DOI: 10.5061/dryad.pzgmsbckh
title: A habitat-based approach to determining the effects of drought on aridland bird communities

description: 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.
description: <b>Study Design and Habitat
 Descriptions</b> 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 (<i>Salix
 </i>spp.), Tamarisk (<i>Tamarix
 <em>ramosissima</em></i>),<i> </i>and
 isolated stands of Fremont cottonwood (<i>Populus
 fremontii</i>) and boxelder (<i>Acer negundo</i>).
 Pinyon-juniper habitat is dominated by two-needle pinyon (<i>Pinus
 edulis</i>) and junipers (<i>Juniperus</i> spp.), which
 varied in relative abundance. The shrub layer in pinyon-juniper habitat
 varies throughout the NCPN, but is often dominated by sagebrush
 (<i>Artemisia </i>spp.), mountain mahogany
 (<i>Cercocarpus </i>spp.), jointfir (<i>Ephedra
 </i>spp<i>.</i>), and cliffrose (<i>Purshia
 </i>spp.). Sagebrush shrubland habitat is dominated by sagebrush,
 primarily big sagebrush (<i>Artemisia tridentata</i>) and
 prairie sagewort (<i>A. frigida</i>), with rabbitbrush
 (<i>Chrysothamnus </i>spp<i>.</i>), greasewood
 (<i>Sarcobatus </i>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). <b>Water Deficit
 Estimates</b> 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).   <b>Breeding Bird Community
 </b> We sampled the breeding bird community at
 675 unique sampling locations on 45 transects during the breeding season
 (May 1<sup>st</sup> through July 15<sup>th</sup>),
 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 &gt;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
 <i>Distance</i> 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 <i>manyglm</i> function in the R package
 <i>mvabund</i> (Wang et al. 2012). The
 <i>mvabund</i> 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
 <i>manyglm</i> 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 <i>mvabund, </i>we sought to maximize the number of
 transect-level density estimates available using <i>Distance
 </i>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 <i>P</i>-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.
description: The code and model to produce the water deficit values are not
 included here because they are in the process of being
 published.