Advancement in timing of important life history events for birds due to climate change presents conservation and monitoring challenges. Song and other vocal activity are strong indicators of avian phenology because they correspond to territorial defense and mate attraction during the breeding season. We combined data from 2 projects using point counts and automated sound recorders to survey passerines during the breeding season in northern California, USA (553 sites, 2009–2011). We used multi-species occupancy modeling to estimate how detection probability based on vocalizations varied over the course of the breeding season. We estimated dates of peak vocal activity, which we reasoned were indicative of reproductive phenology. We demonstrated a strong unimodal relationship between vocal activity and survey date for 8 focal species for which average detectability peaked on June 21 (90% CI: June 19–June 22). Although we found no difference in peak dates of average detectability for migrants vs. residents, the variability of this estimate was lowest for Neotropical migrants compared to residents and elevational migrants. This finding supports previous research suggesting that long-distance migrants may be less flexible in adapting their phenology and more vulnerable to climate change. For an 80% power standard, we found that repeating our level of survey effort on an annual basis would allow detection of an advancement of average peak vocal activity by as small as 2.2 days over 10 yr for the focal species. We could detect smaller average shifts of 0.8 day decade−1 for this group or 1.6 days decade−1 for all passerines over 20 yr. Monitoring vocal phenology of birds through occupancy modeling of survey data from the breeding season is an efficient approach to assessing climate change impacts because species occupancies and measures of community-level diversity can be simultaneously estimated.
The evidence of climate change impacts on birds is mounting both globally and regionally (Parmesan and Yohe 2003, Tingley et al. 2012, Stephens et al. 2016). Effects on avian phenology will be particularly troubling because of the potential for mismatch between environmental conditions associated with egg laying, as well as asynchronies in the rates of climate change on winter and summer ranges used by migratory species (Stenseth and Mysterud 2002, Mayor et al. 2017). Much research has documented temporal shifts in migration arrival and egg laying, including some of the underlying mechanisms for how increasing temperature affects the reproductive success of birds (Stephenson and Bryant 2000, Crick 2004, Both et al. 2009). From a meta-analysis of these types of data, Parmesan (2007) estimated an overall spring advancement of bird phenology of 3.7 days decade −1 across the northern hemisphere, but this rate is an average, not a uniform trend. Mayor et al. (2017) provide examples of strong regional variation in phenological shifts. Therefore, regional monitoring of the phenologies of multiple species is warranted.
Vocal activity in birds is directly related to timing of important breeding behaviors, and has been shown to respond to climate change (Møller et al. 2010, Rubolini et al. 2010). For many species, vocalization rate is related to pairing status, peaks during nest building, and functions in territorial defense against rival males (Slagsvold 1977, Catchpole and Slater 2008, McClure et al. 2011). Moreover, patterns in vocal activity can be highly correlated with seasonal or daily variation in precipitation and temperature (Gordo et al. 2008, McGrann and Furnas 2016). Therefore, vocal activity can be expected to be a good indicator of how breeding phenology is responding to long-term climate trends (Rubolini et al. 2010). Strebel et al. (2014) argue that vocal activity may be a better measure of reproductive phenology than migration arrival dates, which are less directly related to nesting and ignore the contribution of resident species.
Vocalizations are also the primary means by which observers in the field identify birds using standard point counts (Ralph et al. 1995) or automated recorders (Shonfield and Bayne 2017). Both techniques have been applied to efficiently survey birds across large geographic areas (McGrann et al. 2014, Furnas and Callas 2015). The data collected using either of these methods are particularly amenable to occupancy modeling, which requires a large number of sites and temporal replication of surveys at each site during a period of population closure (MacKenzie et al. 2006). Multi-species occupancy modeling provides a means to get the most out of multi-species surveys like those provided by point counts and automated recorders, because data from numerous species can be combined in the same model to increase precision (DeWan and Zipkin 2010, Iknayan et al. 2014). Detection probability in occupancy modeling is usually considered a nuisance variable, but with careful consideration inferences can be drawn with respect to animal behavior. Indeed, in other contexts, activity levels of animals observed at survey locations is considered an important behavioral metric (Rowcliffe et al. 2014). Using occupancy modeling, McGrann and Furnas (2016) drew inference about potential climate vulnerabilities of Neotropical migrants based on temperature effects on detection probability. Of greater interest to the topic of avian phenology, Strebel et al. (2014) proposed the use of date of highest singing activity, quantified by detection probability in occupancy modeling, as a proxy for breeding phenology.
We expanded upon the recommendations of Strebel et al. (2014) to use detection probability in occupancy modeling to draw inferences about the phenological response of birds to climate change. We combined data from 2 projects using point counts and automated recorders to survey passerines during the breeding season in northern California, USA. We integrated these data in a multi-species occupancy model in which we used detection covariates on survey date and its quadratic term to directly estimate the date of peak vocal activity and its uncertainty. We calculated peak dates for 8 focal species known to have strong unimodal singing phenology. To evaluate potential differences in phenological responses among migratory guilds, we calculated average peak date for Neotropical migrants, elevational migrants, and resident species. Similarly, we computed averages for the focal species and for all passerines as a group. We also calculated statistical power to detect changes in peak dates and considered the implications of these results for long-term monitoring.
Study Area and Design
We analyzed bird survey data from an area of northern California, USA, comprising the Klamath Mountains and Southern Cascades ecoregions (39,263 km2). These areas are predominantly covered by conifer- and hardwood/conifer-dominated forests punctuated by meadows and pockets of chaparral, with alpine habitats at higher elevations (Schoenherr 1992). For this study area, we included bird surveys from the Ecoregion Biodiversity Monitoring project (EBM; Furnas and Callas 2015). Those survey sites sampled a mix of public and private lands, but although site selection was generally random, there was less sampling at higher elevations and in wilderness areas where access by roads was limited. To remedy this potential bias, we added bird surveys conducted along a wilderness trail system bisecting the study area (Figure 1, Table 1). The Pacific Crest National Scenic Trail (PCT) is a recreational hiking and equestrian trail that traverses the mountain ranges of the Pacific cordillera from Mexico to Canada. We used survey sites along a 571 km segment of the PCT in northern California. The sites were ∼10 minutes hiking distance apart separated by approximately 500–700 m (McGrann and Thorne 2014).
Summary biophysical conditions at sites where birds were surveyed by the Ecoregion Biodiversity Monitoring (EBM) project, Pacific Crest Trail (PCT) project and at random locations representative of the study area.
The EBM surveys occurred at 185 sites each surveyed in a single year, May 9–July 7, 2009–2011. The automated recorder protocol used at these sites is described in detail by Furnas and Callas (2015). Briefly, we used a lightweight digital voice recorder with factory provided stereo microphones (Olympus DS-40 and DS-61, Olympus Corporation, Center Valley, Pennsylvania, USA). We placed the digital recorder inside a protective container on the ground away from sound-attenuating obstructions (large tree bole, large rock, shrub, etc.) where possible. The devices were programmed to record for 5 min, 3 times each day, for 3 consecutive days at 30 min before local sunrise, sunrise, and 30 min after sunrise. A biologist experienced in California bird identification listened to and viewed the spectrogram of each recording using Raven Pro software (Cornell Lab of Ornithology, Ithaca, New York, USA) to identify species by song, call, or other aural cue.
The PCT surveys occurred at 368 sites surveyed May 16–July 22, 2010. As described in detail by McGrann and Furnas (2016), we followed a standard methodology for fixed-radius point counts (Ralph et al. 1995). Observers hiked and conducted point counts along an 8–11 km trail segment each morning. Observers then walked a trail segment of similar length each afternoon to reach the start of the next set of morning surveys. Observers completed 5 min counts at each site and counted all birds seen or heard on fixed-radius (50 m) circular plots. Two observers separately completed a point count on different dates. The first observer completed the first point count and deployed an automated recorder at 157 of the sites following the same protocol as for EBM. The second observer followed 3 days later, completed a second point count, and retrieved the recorder at locations where deployed.
We included elevation and canopy cover as site covariates in our occupancy modeling. Using a geographic information system (GIS), we extracted elevations at each site from the National Elevation Dataset (10 m resolution; U.S. Geological Survey, http://ned.usgs.gov). Average canopy cover within a 200 m radius surrounding each survey site was calculated using data representing vertically projected percent cover of the live tree canopy layer for 30 m grid cells (Toney et al. 2009).
Furnas and Callas (2015) found that the effective distance for detecting species using recording devices of the design we employed was ∼40–50 m, similar to the 50 m fixed distance used in our point counts. Therefore, we combined data from the 2 point counts and the 9 recordings, producing 11 temporal replicates at sites where both methods were used. Multiple visits at each site allowed us to construct a detection history for use in occupancy estimation techniques (MacKenzie et al. 2006, Iknayan et al. 2014). We limited modeling to 87 species detected at least once in our surveys and belonging to Passeriformes (see Table 2 for a subset of these species, or Supplemental Material Data S1 for the full list). This decision allowed us to confine our analysis of peak vocal activity to a taxon in which sexual selection has shaped bird song to function as one of the most important mechanisms males use to attract mates and to defend territories during the breeding season (Catchpole and Slater 2008).
Occupancy represents the proportion of a study area in which a particular species occurs or the probability that a given location is occupied by the species (MacKenzie et al. 2006). Occupancy modeling allows simultaneous estimation of detection and occupancy probabilities using temporally replicated surveys occurring over a time period for which occupancy is assumed to remain constant. We used multispecies occupancy modeling to control for expected heterogeneity in detection probability among species, sites, temporal survey replicates, and survey methods (Dorazio and Royle 2005, Zipkin et al. 2009). In particular, we combined the survey data from numerous species to fit a Bayesian model wherein differences among species were treated as hyperparameters varying as random effects about mean values for all species. To provide unbiased estimates of species richness, our model included data augmentation to include the effects of rare and cryptic species that were present but never observed (Iknayan et al. 2014). This approach is analogous to the rarefaction of a species accumulation curve to estimate total richness (Colwell and Coddington 1994).
We modeled whether a species i drawn from the augmented dataset occurred in our avian community: wi ∼ Bernoulli (Ω i); whether this species occurred at site j: zi,j ∼ Bernoulli (wi*ψi,j); and whether it was detected in survey k: yi,j,k ∼ Bernoulli (zi,j*pi,j,k). Covariates on occupancy (ψ) and detection probability (p) were modeled via logistic regression. Hyperparameters were placed on all covariates. Data augmentation consisted of 60 additional, unnamed species added to the detection history (y) used in the models. These species were added to the detection history as rows of zeros for each survey site.
We fit a multi-species occupancy model including covariates based on previous evidence from the literature about the important factors likely to affect occupancy and detection probabilities of passerine birds in northern California forests. Based on the findings of Furnas and Callas (2015), the detection component included an intercept representing either one of the 3 recorder survey times (30 min before sunrise, sunrise, 30 min after sunrise) from both projects or either of the 2 point count observers from the PCT for a total of 8 mutually exclusive intercept terms. We also included tree canopy cover to control for variation in the ability of our surveys to detect bird vocalizations in areas of greater vegetation density (Pacifici et al. 2008, Yip et al. 2017).
After controlling for the factors listed above, we included survey date and its quadratic term in the detection component of the model. This pair of covariates allowed us to evaluate our primary research question of how differences in the phenology of vocal activity varied by species. Inclusion of a quadratic term facilitated our investigation of whether detectability for some species was characterized by a unimodal pattern that peaked on some date during the course of the breeding season. Whether with song or other calls, most of the birds we modeled were more active vocally during the breeding season. We assumed that detection probably in our data was a direct measure of vocal activity at each site because the vast majority of detections were based on the aural identification of birds. All of the detections based on recordings were aural, while for the point counts in the field, observers indicated whether each detection was aural only, aural and visual, or visual only. After combining the automated recorder and point counts data (n = 12,083 detections), 98.2% of the detections were aural, 1.3% were both aural and visual, and 0.5% were visual only.
To address spatial differences in species occurrences related to biophysical variables, we included elevation, tree canopy cover, and their quadratic terms in the occupancy component of the model. The inclusion of hyperparameters and quadratic terms in the model allowed the shapes of these associations (positive or negative slope vs. unimodal) to vary by species. We expected differences in elevation and canopy cover to reflect variation in habitat and to be associated with differences in bird species composition (Beedy 1981, Siegel and DeSante 2003, McGrann et al. 2014, Furnas and Callas 2015). All continuous covariates on both detection probability and occupancy (elevation, canopy cover, date) were standardized.
Posterior distributions of derived quantities can be readily computed in Bayesian models. Excluding extra species from data augmentation, we computed dates of maximum detection probability for species for which detection probability displayed a unimodal pattern. We first identified those species with a quadratic term on survey date with a 90% credible interval that did not overlap zero. We further limited species to a group of focal species displaying a strong unimodal pattern for which the standard error on our estimate of peak date was <4 days. We limited the range of possible peak dates evaluated to the months of May, June, and July to match the temporal span of our surveys. To investigate community-level patterns, we computed average peak dates for all passerines, the group of focal species displaying a strong unimodal pattern, and migratory guilds (i.e. resident, elevational, Neotropical). We used Birds of North America species accounts ( http://bna.birds.cornell.edu) to categorize species as residents if they are sedentary and remain in the same climate year-round, elevational migrants if they winter at low-elevation climates and breed in high-elevation climates in California, and Neotropical migrants if they winter in the Neotropics (Mexico, Central America, and South America).
We also calculated average occupancies of species across the study area. We used a model-based rather than design-based approach to inference because we did not assume the combination of the EBM and PCT sites represented a random sample from the study area (Gregoire 1998). Specifically, we predicted occupancies at 1,000 random locations throughout the study area using GIS. We extracted tree canopy cover and elevation covariates values for each of the “virtual” locations using GIS, calculated predicted occupancies based on covariate associations from our modeling, and took the average for each species. Last, we computed average site-level species richness and total area-level species richness for the study area (alpha and gamma diversity; Magurran 2004). Using a similar procedure to that used for peak dates, we estimated the elevation and tree canopy cover at which site-level species richness was maximized.
All models were solved through a Markov Chain Monte Carlo (MCMC) algorithm (Link et al. 2002) implemented in JAGS (4.2.0; Plummer 2003) accessed via R statistical software (3.4.1; R Core Team 2013) with the jagsUI package (Kellner 2015). Uninformative priors were assumed for all parameters. Five independent chains each of 5,000 samples were run with a burn-in period of 2,000 and a thinning rate of 5. Effective mixing of these chains was assessed by means of the Gelman-Rubin convergence statistic (R̂ < 1.1; Gelman et al. 2004). The survey data and the R code we used showing additional technical details on our model are provided as an online supplement (see Supplemental Material Data S1).
We used the model results to assess statistical power to observe changes in dates of peak vocal activity over monitoring timeframes of 10 and 20 yr. We evaluated power to observe shifts in peak dates for each of the focal species. We also evaluated power to detect shifts in the group average peak dates for the focal species, each migratory guild, and all passerines as a group. We used simulation to assess power to detect temporal trends (Purcell et al. 2005, Nielsen et al. 2009, Meyer et al. 2010). We simulated a time series where the true starting peak date was what we estimated from occupancy modeling. We then modeled a state trend as a steady annual rate of advancement (r) from the starting date (D1) such that Dt = [1–r(t–1)] D1 for year t.
We stochastically generated the monitored peak date values about the time series of declining true values by assuming a normal sampling distribution, Dobserved ∼ Norm(Dtrue, σ2 D true). We used the standard errors from our models and adjusted them downward as Dtrue decreased such that the coefficient of variation was held constant. We ran 10,000 Monte Carlo simulations (Metropolis and Ulam 1949) for each species or group and tested for a simple linear trend via ordinary least squares regression. We calculated power as the proportion of simulations where we could reject (P < 0.1) the null hypothesis of a zero or positive slope for the trend line. Through this method of simulations, we determined the minimum rate of phenological advancement that could be detected for each species or group that exceeded an 80% power standard. Our choice of this standard (α: Type I error = 0.1, β: Type II error = 0.2) was based on the recommendations of Bart et al. (2004) for avian monitoring programs. For consistency, we report 90% credible intervals for all parameter estimates from occupancy modeling.
Average elevation and tree canopy cover at EBM sites was similar to averages for the entire study area, but average elevation at PCT sites was 350 m higher than for the study area (Table 1). Average detection probability was highest at sunrise for the automated recorders (Table 3). Detection probabilities varied considerably by species, survey method, and personnel, but automated recorders provided ∼27% higher detection probabilities on average than point counts (Figure 2). Average detection probability (per 5 min survey) was also slightly higher for Neotropical migrants (0.24, 90% CI: 0.23–0.25) than for either elevational migrants (0.21, 90% CI: 0.20–0.22) or residents (0.21, 90% CI: 0.19–0.22). On average for all modeled species, there was a strong unimodal relationship between survey date and detection probability (i.e. quadratic term is credible; Table 3, Figure 3).
Community-level parameter estimates from multi-species occupancy modeling of automated recorder (AR) and point count (PC) avian surveys from the Ecoregion Biodiversity Monitoring (EBM) and Pacific Crest Trail (PCT) projects in northern California, USA, 2009–2011.
There were also unimodal relationships with occupancy for elevation and tree canopy cover (Table 3). Holding canopy cover to its mean value, species richness of passerines peaked at 1,206 m elevation (90% CI: 1,050–1,350). Holding elevation to its mean value, species richness of passerines peaked at 33% canopy cover (90% CI: 28–37). We estimated 34 passerine species with occupancies ≥0.1 (Figure 4). Average site-level species richness of passerines was 10.3 (90% CI: 10.2–10.5). Total, area-level species richness of passerines was 110 (90% CI: 92–132).
Based on 90% credible intervals for the quadratic effect of survey date on detection probability, we found evidence of peak dates of vocal activity in 33 species (see Supplemental Material Data S1). The pattern was most discernable in 8 focal species for which the standard error in peak date was <4 days (Figure 3). These species were Hutton's Vireo (Vireo huttoni), Hermit Thrush (Catharus guttatus), Dark-eyed Junco (Junco hyemalis), Nashville Warbler (Oreothlypis ruficapilla), MacGillivray's Warbler (Geothlypis tolmiei), Yellow Warbler (Setophaga petechia), Western Tanager (Piranga ludoviciana), and Black-headed Grosbeak (Pheucticus melanocephalus). The group mean for these focal species (June 21, 90% CI: June 19–June 22) was later than the average for all species (June 15, 90% CI: June 12–June 17). Although we found no difference in peak dates among migratory guilds, there was less variance in our estimates of peak date of Neotropical migrants compared to elevational migrants and residents (Figure 3).
For a program of repeating the surveys we reported on an annual basis over 10 yr, power analysis showed that only Western Tanager could be monitored at 80% power to observe a shift in peak vocal activity <5.0 days decade−1. However, the range of shifts that could be detected over 20 yr was 1.5–3.8 days decade−1 for each of the focal species. For the focal species, Neotropical, and all-passerines groups, it would be possible to detect average advancements of ∼2–5 days decade−1 over 10 yr (Table 4). If annual surveys were sustained for 20 yr, a shift of ∼1–2 days decade−1 could be detected for these same groups. Our power analysis suggests it would not be possible to monitor either residents or elevational migrants for average shifts <5.0 days decade−1 for timeframes ≤20 yr.
Power analysis for using bird surveys to detect advancement in phenology of avian vocal activity due to climate change.
Using multi-species occupancy modeling, we were able to identify reasonably precise dates of peak vocal activity in 8 focal passerine species (SE < 4 days). We treated peak in vocal activity as a signal of breeding phenology. Furthermore, using power analysis we demonstrated feasibility to monitor group average advances in peak vocal activity of 0.8–4.4 days decade−1 over timeframes of 10–20 yr for the focal species, Neotropical migrants, and all passerines as a group. We also showed that each of the focal species could be monitored individually (≤ 3.8 days decade−1). These levels of precision are comparable to effect sizes already documented for advancement of migration or breeding phenology of several avian species from the same region. Brown et al. (1999) found that nesting and egg laying of Mexican Jay (Aphelocoma wollweberi) advanced an average of 10 days from 1971 to 1998 (3.7 days decade−1). Dunn and Winkler (1999) showed that the egg-laying dates of Tree Swallow (Tachycineta bicolor) throughout North America advanced by 9 days during 1959–1991 (2.8 days decade−1). Inouye et al. (2000) demonstrated dates of first sighting of American Robin (Turdus migratorius) in Colorado moved forward 14 days during 1981–1999 (7.8 days decade−1). These results are also comparable to trends in advancement of avian vocal activity in northern Germany; Rubolini et al. (2010) found that 56 species of short- and long-distance migrants advanced their first singing date by 3.2 days decade−1 on average over an approximately 30-year period. At the global scale, Parmesan (2007) conducted a meta-analysis that found an average phenological advancement of 3.7 days decade−1 for 41 bird species throughout the northern hemisphere.
We found that the focal species characterized by strong unimodal patterns in vocal activity tended to peak later (6 days on average) than other species (Figure 3). One possible explanation is that birds displaying tightly constrained reproductive phenologies may be less able to adapt to changing climatic conditions by advancing their phenologies. For example, Rubolini et al. (2010) found that Afrotropical migrants advanced their date of first song less than resident species. Strebel et al. (2014) used citizen science observations collected in Switzerland to show peak detectability dates in March–April for residents versus April–May for long-distance migrants. Although we were not able to find evidence of different peak dates among migratory guilds, we did find less variance in peak dates among Neotropical migrants compared to elevational migrants and residents. Tropical migrants that breed in temperate regions may have evolved limited phenotypic and behavioral plasticity to changes on their temperate breeding grounds. Compared to resident populations, migrants, having evolved behavioral adaptations that are suitable to tropical and temperate habitats, may have limited ability to adjust behaviors to rapidly changing temperatures on their temperate breeding grounds (Pulido and Widmer 2005). Furthermore, climate change–induced mismatches between environmental conditions on winter and summer ranges may further constrain the ability of migrants to adjust their breeding phenology (Mayor et al. 2017). This inflexibility among migrants is consistent with our observation of reduced variance in peak dates and higher overall detectability of Neotropical migrants compared to other species. For example, when long-distance migrants arrive on summer ranges, males may have less time to establish territories and find mates. Therefore Neotropical migrants sing on hotter days when residents do not, even though this activity is metabolically expensive (McGrann and Furnas 2016).
Resource partitioning is another possible reason why some species display a more tightly constrained reproductive phenology than others, especially among sympatric species on summer range habitats where foraging resources are seasonal (Fretwell 1972). Our results for wood warblers (i.e. Parulidae) support this reasoning; they were the most represented taxon among the 8 focal species and 5 of 6 warblers with occupancies >0.1 (Figure 4) showed some evidence of a unimodal pattern in vocal activity (e.g., Survey date2 in Table 3). Warblers are more likely than other avian taxa to occur in sympatry, which may be linked to greater resource partitioning, especially in foraging behavior and the timing of nesting (MacArthur 1958, Lovette and Bermingham 1999).
Improved monitoring of avian diversity is required to inform effective conservation planning responding to climate change and other stressors (Manly et al. 2005, DeWan and Zipkin 2010, Haughland et al. 2010). Estimating peak dates of vocal activity is a more direct measure of reproductive phenology than migration arrival dates, which excludes resident species (Slagsvold 1977, Strebel et al. 2014). Multi-species occupancy modeling is an efficient tool for combining data from multiple species to simultaneously estimate peak detection dates as well as a variety of other parameters of conservation interest (Iknayan et al. 2014). While watching for advances in phenology, many of the same species can be monitored for coincident declines in occupancy (∼2.5% per year; Furnas and Callas 2015).
To our knowledge, this study provides the first comprehensive assessment of baseline occupancies of passerines over a ∼40,000 km2 region of northern California (Figure 4). We did so by integrating data from 2 large surveys with different sampling coverages, and by using regression to extrapolate average occupancies over the region. Occupancy is likely a good proxy for abundance for point count and automated recorder surveys of breeding birds (Furnas and Callas 2015), and detection/non-detection data provided by these methods could be used to directly model abundance (Royle and Nichols 2003, Royle 2004). With better information about the effective area of surveys, those estimates of abundance could be transformed into estimates of density and population size (Darras et al 2016, Sollmann et al 2016).
The elevation and canopy cover conditions at which we found the highest species richness values could also be monitored for changes over time. We note that the canopy cover optimum was below the average for the study area. This suggests that some mix of climate change and fire suppression may already have reduced avian diversity and populations through the creation of denser forest conditions (Stephens et al. 2009, McIntyre et al. 2015). However, the current optimums for both canopy cover and elevation point toward possible adaption strategies including the movement of avian communities upslope, especially under management designed to promote structurally diverse forests characterized by a greater variety of canopy cover conditions.
Monitoring of phenology and population statuses would benefit from larger data sets with greater sample sizes over larger areas and longer timeframes. We demonstrated the combination of survey data from 2 projects to improve inferences across a large study region. The EBM data provided inference across public and private ownership and over a larger geographical area, whereas the PCT data accessed higher elevations and wilderness areas not well sampled under EBM. We showed how similar methods and protocols made combining data easier and how we adjusted our modeling structure to address differences that remained. In the near future, large amounts of citizen science data collected with the aid of smartphone applications and other types of personal technology could augment data collected by professional scientists (Dickinson et al. 2012, Bonney et al. 2014). Integrated modeling methods can be used to efficiently combine these data to improve accuracy and precision while adjusting for potentially large biases in sampling effort (Handel and Sauer 2017, Pacifici et al. 2017). Furthermore, by controlling for vocal phenology in the detection component of occupancy models, researchers can address concerns about potential confounding of population trends being monitored (Strebel et al. 2014, Handel and Sauer 2017).
Although we identified an overall, community-level unimodal seasonal pattern in vocal activity (e.g., Survey date2 in Table 3), we were only able to identify precise dates of peak vocal activity in 8 focal species. We suspect that day-to-day variation in weather and other environmental factors potentially affecting singing behavior of birds added statistical noise to our data, but that a larger sample size would likely identify peak dates for a larger set of individual species. As the availability of more data and the ability to model them increases in the future, we believe our current power analysis conclusions are conservative. The consequence of this is that we will likely be able to monitor more species over smaller areas, which will be of critical importance because the effects of climate change are expected to be regionally heterogeneous and disruptive (Stralberg et al. 2009, Tingley et al. 2012, McGrann et al. 2014, Mayor et al. 2017). For example, we expect that peak dates of vocal activity are more complex than the estimates we provided for our study area. They likely vary geographically in association with temperature, latitude, and elevation, nuances that could be incorporated into modeling that includes more data from a larger area. One caveat is that we did not use biophysical covariates and model-based inference to adjust our estimates of peak vocal activity for the study area as we did for occupancy and species richness. However, we recommend this approach be taken as more data on phenology become available.
Our use of automated recorders increased site-level detection probabilities of our bird surveys, which allowed us to get more precise estimates of the dates of peak vocal activity. This was a direct consequence of the numerous survey replicates made at each site. We could have achieved a similar increase in precision through additional point counts at each site, but it would have been challenging for us to access each of our survey sites on 3 consecutive mornings. Others have noted the advantages of automated recorders, including limited availability of expert point count observers, synchronization of survey times among sites, and a verifiable record of bird detections (Hobson et al. 2002, Rempel et al. 2005, Shonfield and Bayne 2017). In our case, we used both methods, which further increased precision by allowing us to sample over a greater diversity of times of the morning. Our results show that detection probability varied by species, time of morning, survey method, and survey personnel, details that will be of value to others designing programs to monitor avian phenology and population status (Figure 2).
The use of occupancy modeling of bird detections to monitor avian reproductive phenology has been demonstrated by us in California and others in Switzerland (Strebel et al. 2014). The impacts of climate change on birds are expected to be heterogeneous and disruptive, and migrants may be particularly vulnerable (McGrann and Furnas 2016, Mayor et al. 2017). Therefore, coordinated monitoring of phenology and population status across regions is needed to inform effective conservation. Although integrative modeling approaches can be used to incorporate a greater range of data, greater standardization of survey methods will also be important. Automated recorder survey methods and sound recordings from other sources may help improve standardization.
We thank A. Engilis who interpreted the EBM recordings to identify species detected. We also thank R. Landers, who was one of the point count observers, and A. McGrann, who provided logistical support in the field for the PCT project.
Funding statement: The California Department of Fish and Wildlife and the US Fish and Wildlife Service jointly funded this research through State Wildlife Grants F08AF00126 and F12AF00829.
Ethics statement: We followed passive, noninvasive survey methods (e.g., point counts and automated recorders).
Author contributions: B.J.F. and M.C.M conceived the research, designed surveys, supervised the research, collected data, and wrote the paper. M.C.M. interpreted recordings from the PCT project to identity species detected and was one of the point count observers. B.J.F analyzed the data and orchestrated funding.