Divergent responses of native predators to severe wildfire and biological invasion are mediated by life history

This dataset supports an analysis of the relative and combined effects of severe wildfires and an early-stage barred owl (Strix varia) invasion on three native forest owl species in the Sierra Nevada, California, USA. Data were collected between 2018 and 2023 using a regional passive acoustic monitoring program, complemented by manipulative (lethal barred owl removals) and natural (severe wildfire) experiments. The dataset includes species-specific detection histories for flammulated owls (Psiloscops flammeolus), great horned owls (Bubo virginianus), and northern pygmy owls (Glaucidium californicum), as well as site-level covariates describing fire severity, barred owl removal treatment, topographic conditions, and survey effort. Derived variables include wildfire severity classifications and pre- and post-removal monitoring periods. The data can be reused to run Bayesian single-species "stacked" single-season occupancy models. All sensitive location information of forest owls has been omitted.

Methods
Study system

Our study area encompassed > 6,000 km2 in the northern Sierra Nevada, California, comprised primarily of publicly-managed lands (Lassen and Plumas National Forests) and some private lands. The area included a topographically complex mountain range that varied in elevation and vegetation type (though predominantly mixed conifer forest). Three significant environmental changes occurred in this region from 2018 to 2022: 1) the North Complex Fire ( 1,220 km2 burned in late summer 2020, 599 km2 [49%] at high-severity); 2) the Dixie Fire (3,740 km2 burned in late summer 2021, 2,090 km2 [56%] at high-severity); and 3) a barred owl lethal removal experiment (removals conducted across the entire study area primarily in summer 2019 described below) (Hofstadter et al. 2022). The extent and severity of each fire were both a departure from historical fire regimes and emblematic of emerging trends in the Sierra Nevada (Cova et al. 2023). 

The lethal removal experiment was initiated after documentation of rapid population growth in the project area between 2017 and 2018 (Wood et al. 2020a); most removals were conducted in 2019 (54 birds), with 8 owls removed in 2018, 14 owls removed in 2020, and two in 2021 and 2022 (Hofstadter et al. 2022). We lethally removed barred owls and hybrids using a 12- or 20-gauge shotgun following field protocols established by Diller et al. (2014, 2016), and all removals were carried out by trained and permitted personnel from Sierra Pacific Industries and the University of Wisconsin (Hofstadter et al. 2022). Of the 80 total removals, 67 were barred owls and 13 were barred owl x spotted owl hybrids. We grouped barred owls and hybrids due to their genetic similarity and small hybrid sample size, which prevented a separate analysis for barred owls and hybrids.

Passive acoustic surveys

To characterize the effects of severe fire and invasive barred owls on forest owl communities, we conducted passive acoustic surveys in the northern Sierra Nevada from early May to late July or the onset of wildfire in 2018, 2021, 2022, and 2023, bracketing the periods in which the North Complex Fire (August to September, 2020) and Dixie Fire (July to October, 2021) burned and barred owls were removed (survey data from 2019 and 2020 was excluded because lethal removals and acoustic surveys occurred simultaneously, violating the assumption of closure). Our study relied on nocturnal surveys from 2000 to 0800 hours PDT, aligning with the activity periods of most owls and capturing the crepuscular behavior of pygmy owls (Sater et al. 2006). Although pygmy owls are primarily diurnal hunters, their frequent detections during the crepuscular period justified their inclusion in our analysis (Sater et al. 2006, Deshler et al. 2019). However, because our surveys did not extend into full daylight hours, some pygmy owl detections were likely missed.

Passive acoustic surveys can provide detection/non-detection data required for occupancy modeling (Wood et al. 2019). Each year we surveyed the same 265 hexagonal cells, each 400 ha, similar in size to spotted owl and barred owl territories in the region (Wood et al. 2020a). We systematically selected approximately 1 of 5 hexagonal cells from the grid in the Lassen and Plumas National Forests for surveying, ensuring they were (1) non-contiguous to minimize nonindependence among sites (e.g., detecting the same owls in multiple hexagonal cells) and (2) readily accessible by road. In each hexagonal cell, we deployed two autonomous recording units (ARUs; Swift One Recorder, K. Lisa Yang Center for Conservation Bioacoustics, Cornell Lab of Ornithology, Ithaca, NY) at locations with high topographic relief and a minimum spacing of 500 m. The elevations of sites surveyed ranged from 669 to 2262 m with a mean of 1573 m. In 2018, we deployed ARUs three times at each site for approximately one week per deployment, with an average interval of 26 days (minimum 14 days) between deployments. In 2021, 2022, and 2023, we deployed ARUs continuously for five weeks at each site. Although total survey effort was lower and the deployment schedule differed in 2018, all deployments occurred during the core breeding season (May - July) and at the same spatial locations across years. We also accounted for sampling effort differences by including sampling effort as a fixed effect in our detection models (see below). We deployed ARUs in the same locations each year. Equipped with an omnidirectional microphone, the ARUs recorded sound data at a sample rate of 32 kHz, 16-bit resolution, and gain of +33 dB.

Processing acoustic survey data

To identify forest owl vocalizations, we used the BirdNET algorithm, a deep convolutional neural network designed to identify > 6,000 species by sound, most of them birds (Kahl et al. 2021). We used a custom version of BirdNET overfit to vocalizations of the species of interest in this study. Every 3 seconds, BirdNET produces a unitless numeric prediction score from 0 to 1 for each species. The prediction score reflects its relative confidence in correctly identifying the species, with higher scores representing higher confidence (Wood and Kahl 2024).

To determine the occurrence of a species-specific vocalization, we selected minimum prediction score thresholds that would result in false positive rates of < 1% for each species at the scale of a weekly secondary sampling period. To do so, we first used the program Raven Pro (K. Lisa Yang Center for Conservation Bioacoustics, Cornell Lab of Ornithology) to manually review a random subset of at least 200 hours of audio data each year for each species, with all selected hours containing at least one BirdNET prediction with a score > 0.1. When reviewing the acoustic data, we assigned a code of ‘1’ to an hour-long sample if it contained at least one confirmed true positive or ‘0’ if it did not contain at least one confirmed true positive. In each hour-long sample, we then counted the number of BirdNET predictions above a range of scores (0.1, 0.2, …, 0.9; 0.91, ..., 0.99). Using logistic regression, we estimated the probability of a false positive in a random hour of acoustic data as a function of the number of BirdNET predictions over each score, where we treated the true positive/false positive status as a binary response (i.e., false positive =1; true positive = 0) and the number of BirdNET predictions in an hour above a given score threshold as the predictor (lme4; Bates et al. 2015). We estimated false positive rates for each potential threshold at the scale of a secondary sampling period using the equation 1-(1-FP)n, where FP was the hourly false positive rate, and n was the number of hours in a typical secondary sampling period (n = 84 hours; see below). We selected thresholds that yielded an expected < 1% false positive rate for each species by year combination.

As a result of transitioning to different microphones between 2018 and 2021 and maintaining the same microphones from 2021 onward, we decided to uniquely assess audio data from the 2018 and 2021 seasons, with call rate thresholds in 2022 and 2023 matching those from 2021. To increase the sample size before and after wildfire and lethal barred owl removals, we manually validated predictions for pygmy owls and flammulated owls with a 15% and 42% false positive rate, respectively. We initially included western screech owls and saw-whet owls in our analyses, but had to exclude them due to low detection occurrences, both at barred owl lethal removal and severely burned sites.

Creating encounter histories for occupancy modeling

For flammulated and pygmy owls, we created encounter histories from the four years of passive acoustic survey data at the level of individual ARUs, as these species have relatively small home ranges and quiet calls, meaning ARUs separated by ≥ 500m could be treated as independent sampling units. For great horned owls, a larger species with a larger home range and louder calls, we created encounter histories at the level of sampling hexagons (i.e., two ARUs), as multiple ARUs within a 400-ha hexagon may record calls from the same individual. We divided the survey season into eleven one-week-long secondary sampling periods, starting on May 7th and ending on July 23rd. We assumed population closure throughout the year. If an ARU did not record during a specific secondary sampling period, we treated that period as null. However, if an ARU was active for only part of a secondary sampling period, we included the data from that period and treated survey hours as a detection covariate to account for unequal sampling (see below).

Environmental predictors

To calculate predictor variables, we created circular buffers around all sites that represent the approximate size of each species’ home range (Giese and Forsman 2003, Bennett and Bloom 2005, Yanco and Linkhart 2018). For flammulated and pygmy owls, we buffered individual ARUs at 50 ha and 300 ha, respectively. To represent larger home ranges in great horned owls, we delineated buffers of 700 ha around the center of the sampled 400-ha hexagons. Next, we utilized severely burned area (with at least 75% overstory mortality) as estimated by the Monitoring Trends in Burn Severity database (www.mtbs.gov) to generate two wildfire predictors: 1) a binary response of either burned (1) or unburned (0) indicating whether any severe fire occurred within a home range buffer and 2) a continuous value ranging from 0 to 1 representing the proportion of the area that burned at high-severity during the Dixie Fire and North Complex Fire, providing a measure of the local extent of severe fire.

We created a binary predictor to represent the lethal removal of barred owls based on the 57 lethal removal locations of barred owls and barred owl x spotted owl hybrids from 2018 through 2022 that overlapped with our acoustic study area (Hofstadter et al. 2022). We first delineated a circular buffer of 2,004 ha around each barred owl removal location, the approximate home range size of this species in the region (Wood et al. 2020a). Then, for flammulated and pygmy owls, we classified a site (i.e., ARU) as ‘1’ if their species-specific home range buffer overlapped with a barred owl home range buffer, and as ‘0’ if no such overlap occurred. For great horned owls, we designated sites as '1' if a great horned owl's home range buffer centered around a hexagon overlapped a barred owl's home range buffer. We coded sites as '0' if the great horned owl’s home range buffer did not overlap a barred owl's home range buffer.

Modeling site occupancy probability

We developed Bayesian single-species “stacked” single-season occupancy models to assess the effect of barred owls and wildfire on site occupancy rates for the three forest owl species and controlled for imperfect detection (MacKenzie et al. 2018, Jones et al. 2020). We stacked the four-year encounter histories (2018, 2021, 2022, 2023) because we were more interested in regional owl occupancy patterns rather than turnover rates. We considered a year and site random effect on occupancy (see below) because the stacked model assumes independence among sites and therefore may underestimate error in model coefficients. We also chose this data format to increase the effective sample size, thereby enabling us to fit a larger set of occupancy predictors for less prevalent species and avoid difficulties associated with poor model fit. 

We first modeled the potential effects of five predictors as fixed effects on the probability of detecting each species. We considered ordinal date as a predictor to accommodate potential changes in vocalization patterns over the breeding season (Crozier et al. 2003). We also considered the year as a predictor to assess interannual variation in detectability and account for different microphones in 2018. We included terrain ruggedness as a predictor because owl vocalizations may travel a shorter distance in areas of complex topography. We calculated terrain ruggedness by averaging the absolute differences between the elevation of a 10 m raster grid cell and the values of its eight surrounding cells within a 390 m buffer around each ARU location, where elevation data were obtained from a digital elevation model (Hijmans et al. 2015). We treated the number of hours ARUs recorded per secondary sampling period as a predictor, given that units were deployed for variable durations. We used the barred owl detection/non-detection history obtained from our passive acoustic surveys as a predictor to account for potential changes in vocalization behavior in response to barred owl presence (Wood et al. 2020b, Rugg et al. 2023). We coded hexagonal cells as ‘1’ if a barred owl was detected at least once in a given year at either ARU or ‘0’ if no detection occurred during that year. In 2018, barred owls were acoustically detected at 31 hexagonal cells, 18 of which overlapped with follow-up barred owl lethal removal locations. We did not include background noise as a detection covariate to maintain a minimal set of covariates and because BirdNET is robust to high ambient noise (Kahl et al. 2021). No pairs of predictors were highly correlated (all Pearson’s correlation coefficients < 0.6). We incorporated all five predictors into the occupancy models for each species to minimize the potential number of models (see below).

Next, we built six individual occupancy models to assess whether barred owl removals, severe wildfires, or interactions between the two predictors affected owl occupancy rates. We treated unburned and barred owl non-removal sites as control sites for each owl species and employed an unbalanced Before-After-Control-Impact design. 

To test for severe fire effects independent of time since burn, we created the Fire Time-Constant Model: 

                    logit(Ψi) = β0 + β1burn + β2allpostburn + μ1[1|year] + μ2[1|site] 

where burn was an indicator predictor for burned sites (i.e., burned = 1; unburned = 0), allpostburn was a continuous predictor representing the proportion of the area that burned at high-severity at any point post-fire, year was a random effect indicator for survey year, and site was a random effect for individual sites (i.e., ARUs for flammulated and pygmy owls and hexagonal cells for great horned owls). 

To test for severe fire effects while accounting for different years post-burn, we created the Fire Time-Dependent Model: 

                     logit(Ψi) = β0 + β1burn + β2burn1 + β3burn2 + β4burn3 + μ1[1|year] + μ2[1|site]

where burn was the same as in the Fire Time-Constant Model and burn1, burn2, and burn3 were continuous predictors representing the proportion of the area that burned at high-severity relative to the number of years post-fire (i.e., all pre-fire sites and unburned sites post-fire = 0).

To test for barred owl effects independent of time since removal, we developed the Barred Owl Time-Constant Model: 

                        logit(Ψi) = β0 + β1site-type + β2 allpostlethal + μ1[1|year] + μ2[1|site]

where site-type was an indicator predictor for barred owl removal sites (i.e., removals = 1; no removals = 0 based on the criteria described above), allpostlethal was an indicator predictor for removal sites at any point post-removal, year was a random effect indicator for survey year, and site was a random effect for individual sites (i.e., ARUs for flammulated and pygmy owls and hexagonal cells for great horned owls). 

To test for barred owl effects while accounting for years since removal, we built the Barred Owl Time-Dependent Model: 

                        logit(Ψi) = β0 + β1site-type + β2 lethal1 + β3lethal2 + β4lethal3^+^ + μ1[1|year] + μ2[1|site]

where site-type was the same as in Barred Owl Time-Constant Model and lethal1, lethal2, and lethal3+ were indicator predictors for removal sites one, two, and three or more years post-removal (i.e., removal sites post-removal = 1; all pre-removal sites and non-removal sites post-removal = 0). We evaluated for the effects of barred owl lethal removal parameters by comparing native owl occupancy at removal sites to occupancy at those same sites in the years following removal. Since lethal removals happened in different years, a direct before-and-after comparison between control and impact sites was not possible. To address this, we used a contrast analysis to measure overlap between posterior distributions at removal sites (i.e., site-typelethal) and post-removal sites (i.e., lethal1lethal, lethal2lethal, lethal3+lethal, and allpostlethallethal) (Pastore and Calcagnì 2019). We considered differences between posterior distributions to be biologically meaningful if their overlap was less than 15%.

To test for interactions between fire and barred owl removals independent of time since burn, we created the Interaction Time-Constant Model:

                      logit(Ψi) = β0 + β1site-type + β2burn + β3allpostlethal*burnYrs + μ1[1|year] + μ2[1|site]

where site-type and burn were indicator predictors for barred owl removal and burned sites (i.e., removal or burned site = 1; no removals or unburned site = 0 based on the criteria described above) and allpostlethal*burnYrs indicated sites where barred owl removals preceded wildfire, measured for all years post-fire (i.e., burned removal site post-burn = 1; all other sites = 0).

Finally, to test for interactions between fire and barred owl removals while accounting for time since burn, we developed the Interaction Time-Dependent Model:

                     logit(Ψi) = β0 + β1site-type + β2burn + β3allpostlethalburn1 + β4allpostlethalburn2 + β5allpostlethal*burn3

                                    + μ1[1|year] + μ2[1|site]

where site-type and burn were the same as in Interaction Time-Constant Model and allpostlethalburn1, allpostlethalburn2, and allpostlethalburn3 were indicator predictors for sites with interactions between barred owl removals and severe fire depending on the years since the fire occurred (i.e., allpostlethalburn3 = 1 indicates a site where a barred owl removal was followed by high-severity fire in the third year after the fire, with the covariate representing a binary indicator rather than a continuous fire severity measure). We performed occupancy analyses using the R package “ubms” (Kellner et al. 2022) and used its predict function to estimate and visualize occupancy probabilities.