DOI: 10.5061/dryad.hqbzkh1g6
title: Year-round monitoring at a Pacific coastal campus reveals similar winter and spring collision mortality and high vulnerability of the Varied Thrush

description: Bird-window collisions are a leading cause of direct anthropogenic avian
 mortality, yet our state of knowledge regarding this threat relies heavily
 on eastern North American studies. Seasonal patterns of collision
 mortality may differ along the Pacific coast, and western North American
 species remain understudied. We therefore surveyed a stratified random
 sample of 8 buildings for collisions at the University of British
 Columbia, Vancouver, Canada over 45-day periods during 2 winters, 1
 spring, 1 summer and 1 fall season between January 22, 2015 and March 15,
 2017. After accounting for the rate of scavenging and efficiency of
 observers in finding carcasses, we estimated that 360 collision fatalities
 (95% C.I.: 281 to 486) occurred over 225 days of monitoring. Collision
 mortality was highest in fall, but in contrast to most published research,
 collision mortality was intermediate in both winter and spring, and was
 lowest in summer. In winter 2017, we performed point count surveys to
 assess whether individual species are disproportionately vulnerable to
 collisions when accounting for population size, and found that the Varied
 Thrush (Ixoreus naevius) was 76.9 times more likely to collide with
 buildings, relative to average species vulnerability in winter. To our
 knowledge, this is the first study to report the Varied Thrush as a
 species that is disproportionately vulnerable to collisions. Further
 studies are needed to assess the vulnerability of Western North American
 species and subspecies and to determine whether similar patterns of
 seasonal collision mortality are found elsewhere.
description: <strong>Study Area</strong> We
 conducted our research at the 420-ha Vancouver campus of the University of
 British Columbia (hereafter referred to as UBC), along the Pacific coast
 of southern British Columbia (BC), Canada (49.2606° N, 49.2606° W). The
 campus is situated within the Fraser Lowlands (elevation 0 - 152 m)
 bordered by the Strait of Georgia and the Coast Mountains (Figure 1). UBC
 is located on a peninsula, separated from the City of Vancouver by Pacific
 Sprit Park, a 874-ha park comprised mainly of 70 - 130 year-old
 second-growth forest (Metro Vancouver 2019). Tree canopy at UBC is
 predominantly coniferous; dominated by western red cedar (<em>Thuja
 plicata</em>), Douglas-fir (<em>Pseudotsuga
 menziesii</em>), western hemlock (<em>Tsuga
 heterophylla</em>), and grand fir (<em>Abies
 grandis</em>). ~58% of campus is in
 park-like open condition, with mowed grass, a number of coniferous and
 evergreen broad-leaf plantings, and scattered deciduous trees; primarily
 red oak (<em>Quercus ruba</em>) and maple species
 (<em>Acer</em> spp.) (Dyck 2016). The campus and surrounding
 forest lies within the Coastal Western Hemlock (CWH) Biogeoclimatic Zone
 (DataBC 2018) characterized by warm dry summers and mild wet winters, with
 a climax forest dominated by western hemlock (<em>Tsuga
 heterophylla</em>), as well as Douglas-fir
 (<em>Pseudotsuga menziesii</em>) and western red cedar
 (<em>Thuja plicata</em>) (Pojar et al. 1987, 1991). Broadly,
 the region is part of the temperate rainforest biome, also known as the
 Northern Pacific Rainforest Bird Conservation Region, composed of similar
 climate, habitat types and bird species, extending from the Gulf of
 Alaska, through to Northern California (Alaback 1996, North American Bird
 Conservation Initiative 2000, Rich et al. 2004).
 <strong>Building Selection</strong> We chose 8 buildings for
 bird-window collision monitoring at UBC, using a randomized design,
 stratified by building height and extent of surrounding vegetation (Figure
 2; see Hager and Cosentino 2014, Hager et al. 2017). We used a map layer
 of 229 UBC building footprints and size-related attributes of all
 buildings &gt;1 story, excluding parking garages (University of
 British Columbia 2015), to divide buildings into two categories: “short”
 buildings (2-4 stories) and “tall” buildings (5-12 stories). We used
 ArcMap 10.3.1 to generate a 50-m buffer around buildings and digitized all
 vegetation using the World Imagery base map (ESRI 2015) to calculate %
 vegetation within 50 m of each building. Within each building height
 class, we then randomly selected 2 buildings with &gt;50% vegetation
 cover within 50 m of the building and 2 buildings with &lt;50%
 vegetation cover within 50 m, as defined in Hager and Cosentino (2014).
 Excluding 1-story buildings eliminated most maintenance sheds and utility
 buildings with few windows. We further specified that study buildings must
 be a minimum of 100m apart, which resulted in one building being removed
 from the initial selection. We excluded a second building from our initial
 selection due to construction fencing which prevented access to the
 majority of the building perimeter.
 <strong>Collision Surveys</strong> We followed protocols outlined in
 Hager and Cosentino (2014) and Hager et al. (2017) to conduct a total of
 155 bird-window collision surveys between January 22, 2015 and March 15,
 2017 at our 8 study buildings. Surveys were conducted over 5 45-day
 periods consisting of 2 winter seasons (January 22, 2015 - March 7, 2015
 and January 23, 2017 – March 15, 2017), 1 fall season (September 9, 2015 -
 October 23, 2015), 1 spring season (April 15, 2016 - May 29, 2016 and 1
 summer season (June 16, 2016 - July 30, 2016). We aligned our spring and
 fall collision monitoring with the calendar dates when a local banding
 station records the highest numbers of migrants (WildResearch 2015). We
 performed a clean-up survey the day before each monitoring period to
 remove all evidence of bird-window collisions that accumulated prior to
 survey periods, including carcasses and feathers. In winter 2015 we
 surveyed every second day for 45 days. For the second winter season and
 all other seasons we conducted daily surveys for 21 days followed by
 surveys every second day for an additional 24 days; for a total of 45
 days. We stopped surveys for 6 days during the second winter season due to
 a heavy snowfall that prevented us from reliably finding collision
 evidence. Surveys resumed after a clean-up survey to remove any evidence
 that accumulated while surveys were halted, and the survey period was
 extended accordingly. We conducted collision surveys in
 mid-afternoon as recommended by Hager and Cosentino (2014), since window
 collisions often peak in the morning and continue at a lower rate later in
 the day (Klem 1989, Kahle et al. 2016). Two surveyors worked concurrently,
 walking around the perimeter of buildings in opposite directions,
 searching for carcasses within 2 m of the building façades. We defined a
 façade as the exterior face of a building, generally facing the same
 direction. In some cases, the complexity of a building footprint required
 us to consider several portions of a building footprint as one façade,
 e.g., several sides of an alcove, if we could not reliably discern which
 portion of the façade that a bird may have hit. One surveyor, ANP,
 remained constant throughout all 5 collision monitoring periods, while the
 second surveyor(s) differed seasonally. Evidence of a collision included
 whole carcasses and partial carcasses (e.g. wing, feet, bill and/or
 feather piles) (Hager and Cosentino 2014, Hager et al. 2017). We
 conservatively specified that feather piles must include a minimum of 10
 tail, wing or body feathers confined to 50 cm diameter circular area
 (Ponce et al. 2010). Only two stunned/injured birds were found during
 surveys, and both died shortly after discovery. We removed all carcasses
 immediately following each survey to ensure double counting did not occur,
 and identified all intact carcasses to species. We also collected partial
 carcasses, including all feathers, and where possible these were later
 identified to species-, genus- or family-level, using online resources
 (U.S. Fish and Wildlife Service Forensics Lab 2018), hard copy guides
 (Scott and McFarland 2010), or by comparing to specimens from the UBC
 Beaty Biodiversity Museum. <strong>Correcting For
 Biases Due To Scavenging And Missed Carcasses By
 Observers</strong> Collision frequency estimates
 can be negatively biased due to carcass removal by humans and scavengers
 prior to their observation by surveyors, and imperfect detection of
 carcasses by surveyors (Klem et al. 2004, Bayne et al. 2012, Riding and
 Loss 2018). Therefore, it is important to design trials to estimate bias
 parameters when comparing building mortalities among studies, or among
 seasons and regions of interest (Loss et al. 2015). To account for carcass
 removal and imperfect detection, we conducted carcass persistence and
 searcher efficiency trials, then incorporated the results into statistical
 models using the Generalized Mortality Estimator in R package
 <em>GenEst</em> (Dalthorp et al. 2018a, b) to correct our
 collision mortality estimates (details below). We tested both
 intercept-only and season-specific models, but had insufficient data for
 both carcass persistence and searcher efficiency trials to allow us to
 apply bias corrections at the species-level, by building perimeter
 substrate or by other potential covariates. <strong>Carcass persistence
 (CP). </strong>We conducted 103 carcass persistence (CP) trials to
 quantify bias due to scavenging of carcasses. Trial periods were chosen
 based on field worker availability and logistics of conducting frequent
 carcass checks; February 16 - March 7 in winter 2017 (n = 33), April 20 -
 April 28 in spring 2016 (n = 28), July 12 - July 30 in summer 2016 (n =
 28), and September 10 - October 16 in fall 2015 (n = 14). We used
 carcasses from a range of previously window-killed birds, ranging in size
 from 4g (Anna’s Hummingbird, <em>Calypte anna</em>) to 128g
 (Steller’s Jay, <em>Cyanocitta stelleri</em>). Carcasses were
 stored in deep freeze and held under our Canadian Wildlife Service Salvage
 Permit, as outlined in the Ethics Statement. We placed thawed carcasses at
 a range of times throughout the day at a randomly selected study building
 and at a randomly-selected façade. The location along each façade was also
 chosen randomly, but was no closer than 2m from a corner or edge of an
 adjacent façade. Carcasses were placed on a variety of substrates,
 depending on the randomized location chosen, including concrete, asphalt,
 river rock, bark mulch, grass, or under bushes and ground cover. Collision
 monitors were alerted to the carcass persistence study and to the location
 of CP trial carcasses. However, as a precaution, data were checked
 carefully following all collision surveys to ensure that carcasses from CP
 trials were not inadvertently included in the collision mortality dataset.
 We checked carcasses every 5 – 6 hours after initial placement on the
 first day until early evening, in early morning and mid-day on the second
 day, and at mid-day on subsequent days until carcasses were no longer
 present (presumably removed by scavengers or campus maintenance staff), or
 until the end of the study period, whichever arose first. Carcass removal
 time was the interval between placement of a carcass at a façade until it
 was removed by a “scavenger” and no longer detectable by observers, as
 defined by Riding and Loss (2018). We used
 bias-correction functions within the R package <em>GenEst</em>
 to fit a carcass persistence (or parametric survival) model to estimate
 the amount of time a carcass would persist, given the conditions under
 which it arrived (Dalthorp et al. 2018a, b). We tested whether a
 persistence model that explicitly varies depending on the season, or one
 that depends only on the intercept, best described the detection
 probability of a carcass over time. We fit one-parameter exponential and
 two-parameter (location and scale) Weibull models, and tested whether
 persistence parameters varied categorically by season. The location
 parameter in the exponential and Weibull models allowed covariates to
 shift the mean of the probability density function, whereas the additional
 scale parameter in the Weibull model allowed covariates to affect the
 degree of spread of the distribution. Since the sample size was small
 (n/K &lt; 40) we
 chose the model with the lowest 2<sup>nd</sup>-order AIC
 (AICc) goodness-of-fit statistic across models. <strong>Searcher efficiency
 (SE). </strong>We estimated searcher efficiency by having a
 non-surveyor place a total of 29 thawed bird carcasses during the winter
 2015 (n = 1), winter 2017 (n = 20), spring 2016 (n = 3), summer 2016 (n =
 3), and fall 2015 (n = 2) study periods. Carcasses ranging in size from 4
 g (Anna’s Hummingbird) to 160 g (Northern Flicker, <em>Colaptes
 auratus</em>) were placed at randomly-selected buildings, façades,
 and façade locations as described above for CP trials. Carcasses were
 removed from the trial if they were detected by at least one of the
 surveyors. Carcasses that were not detected on the first day were checked
 prior to each subsequent survey until they were located by observers,
 removed by a scavenger, or until the end of the collision monitoring
 period, whichever arose first. We used the R package
 <em>GenEst</em> to model searcher efficiency as a function of
 the probability that carcasses are located by observers on successive
 searches (Dalthorp et al. 2018a, b). The probability of finding a carcass
 on the first search after carcass arrival was denoted by the
 parameter, <em>p</em>, and <em>k</em> was the
 change in probability of finding a given carcass on the second and
 subsequent searches until the end of the collision monitoring period
 (Dalthorp et al. 2018a). We fit two models to estimate the observer bias
 parameters, including a model that depended only on the respective
 intercepts, and another model that allowed searcher efficiency to vary by
 season. We selected the best model by comparing the goodness of fit of
 models using AICc.   
 <strong>Bias-corrected mortality estimation:</strong><strong> a two-stage parametric bootstrap approach.</strong> We used the function estM in the R package <em>GenEst</em> to produce a bias-corrected estimate of mortality (Dalthorp et al. 2018b). EstM estimated each carcass’ contribution to the total mortality in each search interval by accounting for the carcasses missed due to scavenging/imperfect carcass persistence (CP), and imperfect detection/searcher efficiency (SE) (Dalthorp et al. 2018a, b). We assumed that most birds that collided with buildings fell within 2 m of building façades (Hager et al. 2013, Hager and Consentino 2014, Riding et al. 2020), and therefore set the parameter density weighted proportion (dwp) or the proportion of carcasses expected to arrive in the area searched by observers (Simonis et al. 2018), to 1. The model uncertainty from the two bias correction models (carcass persistence and searcher efficiency models) was propagated into estimates of total mortality using parametric bootstrap methods (Efron and Tibshirani 1986). Specifically, function estM in the R package <em>GenEst</em> drew 1000 sets of estimated carcass counts from the asymptotic distributions of the maximum likelihood estimators of parameters from the two bias-correction models. The result was a matrix of carcass contributions where each row corresponded to a carcass, and each column corresponded to a simulated realisation of that carcass’ estimated contribution to total mortality. These simulated outputs resulted in bias-adjusted mortality estimates that included corrections for differences in search interval (time-increments between searches on the same building façade), corrections for seasonal difference in carcass persistence, and a correction for imperfect searcher efficiency. <strong>Assessment of Relative Vulnerability to Collisions in Winter 2017</strong> After 4 seasons of collision monitoring, we noted that some species were particularly prevalent in the collision mortality data. However, assessment of relative vulnerability to collisions requires an estimate of population size of species within the study area, as well as species-specific mortality counts. Therefore, in winter 2017, we conducted avian point count surveys during our collision monitoring period, and corrected these counts for imperfect detection. We used raw collision mortality counts to calculate relative vulnerability of species to window collision mortality on UBC campus as we did not have sufficient data to generate species-specific bias-corrected mortality counts as described above. To determine how population size influenced the likelihood of collision by species, we applied a two-step approach: 1) we estimated the change in total number of birds by species over time in response to population-level processes of mortality, immigration and emigration, and then; 2). we regressed mortality estimates on the corrected population counts to determine a species-specific vulnerability rating. <strong>Winter 2017 point count surveys.</strong> We conducted 5-minute point counts at a randomly chosen façade of each building approximately weekly between January 25, 2017 and March 12, 2017, for a total of 72 point counts over 9 survey days throughout the winter 2017 collision monitoring period. We minimized factors that could contribute to detection error by: 1) recording all birds seen or heard by call notes within a small area to maximize the probability of detecting all individuals (i.e., a 50-m radius semi-circle, from the edge of building façades), and excluding birds that were only seen flying over; 2) conducting surveys during the morning hours when human foot traffic was low, and only in weather conditions favorable to maximum detection of all individuals (i.e., with no precipitation and winds below 13 km/hr), and; 3) using the same observer (KLD) to conduct all point count surveys to minimize variation in observer bias. <strong>Population size estimates. </strong>We calculated the population size for each species using the Dail–Madsen model (Dail and Madsen 2011), which is a generalized form of the Royle (2004) binomial N-mixture model developed for open populations (i.e., the closure assumption is relaxed, allowing abundance to change over time). We chose this model for our winter point count data because it 1) assumes that abundance patterns are determined by an initial territory establishment process followed by gains and losses resulting from mortality, movements into the population (recruitment) and movements out of the population, and; 2) accounts for imperfect detection probability. For all species, we assumed that individuals could move across all sites (building locations) and all days within the winter 2017 collision monitoring and point count period. Thus, the number of gains across days was directly proportional to the number of survivors (individuals occupying winter territories), so we applied an auto-regressive model to recruitment, <em>γ</em>*<em>N<sub>s</sub></em>[<em>i</em>,<em>t</em>-1]; where <em>γ</em> = recruitment, N = the number of individuals per site, i, for each species, s, on survey day, t. As a result of the relatively short detection distance, we used a starting value of 1 for probability of detection. We applied a Poisson distribution to all species in which the total number of individuals detected across all 9 surveys, N<sub>stot</sub>, was greater than 10. For species with <em>N</em><sub>stot</sub> &lt; 10, we applied either a negative-binomial or zero-inflated Poisson distribution to a simple intercept-only model, and restricted immigration to zero. We then used the parameter estimates from this simple model as starting values for the more complex, open count model that allowed immigration across survey locations and days. For all models, we set the upper bound for discrete integration, <em>K</em>, at the greater value of <em>N</em><sub>stot</sub> + 1, or 10. We summed across all sites the resulting site-specific population size estimates corrected for detection probability and movement into the study area, divided by the number of survey-days (9), and multiplied by number of buildings surveyed to obtain population estimates for each species, <em>N</em><sub>corr</sub>. We calculated the root mean squared error on the residuals of each model to estimate residual standard error for each species’ estimate of N<sub>corr</sub>. We estimated all parameters and obtained corrected population size estimates using the function, ‘pcountOpen’ and back-transformed the estimates for detection probability from each model using the function, ‘plogis’ in the R package, <em>unmarked</em> (Dail and Madsen 2011, Fiske and Chandler 2011, Hostetler and Chandler 2015, R Development Core Team 2018). <strong>Calculation of relative vulnerability to collisions by species.</strong> Next we used our raw mortality counts per species for winter 2017 and species population estimates at our 8 study buildings (N<sub>corr</sub>) to assess relative vulnerability of species to collisions during the winter 2017 season. We fit a simple linear regression to log<sub>10 </sub>(X + 1) of the raw winter 2017 mortality counts on log10 of the population size estimates, N<sub>corr</sub> (Arnold and Zink 2011, Loss et al. 2014). The use of the (X + 1) transformation on mortality counts allowed us to account for zero mortality for species that were observed in point counts but not in the mortality counts. We set regression coefficients for population size, N<sub>corr</sub> to 1.0, calculated residuals, then raised 10 to the power of the absolute value of the residuals, for all species observed in our point count surveys. A slope of 1.0 assumes that mortality due to collisions for a given species is proportionate to their population size, as described by Arnold and Zink (2011). Vulnerability values represent the factor by which a given species is more (+) or less (–) vulnerable to collisions relative to the expected value for a species with average collision risk (residual = 0).