Aim Our specific objectives were to (1) estimate probabilities of detection of each species, after taking into account survey-specific covariates, and (2) investigate the influence of site-specific covariates on owl and nightjar abundance, integrating effects of imperfect detection. Innovation We conceived a survey protocol to estimate probabilities of detection and estimates of abundance of owls and nightjars in a large area, the Basque Country, northern Spain. Main Conclusions Our results show that detection probability was strongly influenced by playback broadcast and by observer experience. Date irregularly affected species according to their reproductive periods, and we also found that vocal activity gradually diminished proportionally to the hour after sunset. Tawny owl (Strix aluco) was the most abundant and widely distributed species. Its abundance was positively related to forest areas (mainly pine timber forests) and decreased in large urban and agricultural areas. Open space species were less common. Barn owls (Tyto alba), little owls (Athene noctua), Eurasian scops owls (Otus scops) and long-eared owls (Asio otus) avoided forest areas, but showed different responses to agriculture, grass-fields, scrub and urban areas. Finally, European nightjar (Caprimulgus europaeus) was moderately frequent, and its abundance was favored by scrub areas and, weakly, by eucalyptus patches, whereas it was negatively affected by large forest areas. We have shown that it is fundamental to consider the effects of survey-specific covariates in the methodology design and analytical development. Our results also indicate ecological adaptations and population changes in the nocturnal bird community following an increase in urbanization and in the extent of timber plantations, and also the simplification of natural habitats.

Study area

This study was carried out in the Basque Country, northern Spain (7,234 km², lying between 42º and 43º N and 1º and 3º W; Fig. 1). There are two clearly defined areas (roughly north and south). The northern area, Cantabric region, runs along the coast of the Bay of Biscay. It has an Atlantic climate and mild temperatures with a thermic oscillation of 12ºC from the coldest to the hottest months and 1,200-2,000 mm of rainfall distributed throughout the year (www.euskalmet.euskadi.net). The landscape is mountainous and densely populated, with extensive urban and industrial areas, mainly located in valley floors and on the gentler slopes. Forestry plantations (Pinus radiata and Eucalyptus spp.) have become widespread in the last 80 years, gradually replacing grazing land for extensively-reared livestock, traditional agricultural activities, as well as a few remnants of native forest. The other area, of some 2,500 km², lies to the south and is situated in a transition area to the Mediterranean climatic region. The climate is dry with a thermic oscillation of 17ºC from the coldest to the hottest months and the landscape is dominated by arable lands, vineyards, Mediterranean scrub and holm-oak woods in the sloping areas.

Survey design and survey protocol

Survey methods were based on the methodology used in the little owl (Athene noctua) census carried out in the Basque Country (Zuberogoitia, Zabala, & Martínez, 2011), in which we conducted a large scale detection/non-detection survey encompassing all the Basque Country. However, in the present case, we considered all the owl species and also European nightjars (Caprimulgus europaeus).

First, we randomly selected 65 5x5 km Universal Transverse Mercator (UTM) squares that represented all the vegetation types of the Basque Country. These UTMs were considered our Sampling Units (SU). Second, we randomly selected eight Survey Points (SP) in each SU. The SPs were chosen according to the main habitat types at each SU, considering a minimum distance of 1 km between two SPs. All the SPs were established before the beginning of the surveys and were kept without change until the end of the study. In all, we surveyed 521 SPs which will be considered the sampling sites for our hierarchical models (Fig. 1).

We considered seven survey periods in 2018 (one per month through January to July, most of the breeding cycle of the target species; Zuberogoitia, Martínez, & Alonso, 2011). In each period we surveyed between four (minimum) and eight (maximum) SPs per SU, being the average of six SPs surveyed per SU and survey period.

In each SP we developed a three-period survey protocol for owl and nightjar census: 5 min waiting for spontaneous voices, 5 min of playback broadcast voices of only one species and 5 min waiting in silence. Therefore, we spent 15 min per survey in each SP. All the surveys were developed during the first hours after dusk. Surveying eight SP in a SU required on average 180 min, and therefore only one or a maximum of two SUs were surveyed per night. Surveys were conducted on calm and dry nights and they were suspended if it started to rain or wind increased, reducing owl detectability (Lengagne & Slater, 2002; Braga & Motta-Junior, 2009; Zuberogoitia et al., 2019).

We noted every owl species and nightjars in each survey, considering those that were detected in the first 5-min spontaneous, during the 5-min broadcast period and during the last 5-min silent period. We broadcast playback voices of only one species per survey period. The species were chosen according to the annual maximum peak of response to playback broadcasts in our study area. In this sense, we broadcast voices of eagle owls (Bubo bubo) in January (Martínez & Zuberogoitia, 2003), tawny owls (Strix aluco) in February (Zuberogoitia & Martínez, 2000; Zuberogoitia et al., 2019), long-eared owls (Asio otus) in March (Martínez, Zuberogoitia, Colás, & Macía, 2002), barn owls (Tyto alba) in April (Zuberogoitia & Campos, 1998), Eurasian scops owls (Otus scops) in the last weeks of April and first weeks of May (Zuberogoitia, Martínez, & Alonso, 2011) and little owls during June and the first half of July (Zuberogoitia, Zabala, & Martínez, 2011). There were no records of boreal owls (Aegolius funereus) in the study area, although we also included survey points with broadcast voices of this species in mountain old forest areas in January (Badosa, López, Potrony, Bonada, & Gil, 2012). Similarly, although there were no records of breeding attempts of short-toed owls (Asio flammeus) in the study area, we considered the courtship main period of this species in Spain (Onrubia, 2016) and included some surveys in open areas broadcasting voices of the species in March and April. We did not use broadcast voices of nightjars.

Playback records were built according to our own experience (e.g., Zuberogoitia & Campos, 1998). We included a mix of territorial and mating voices of males and females of each species in tracks of 5 min (just the time needed for each survey). Voices were downloaded from xeno-canto (https://www.xeno-canto.org/), choosing only those clear and adequate for each case, recorded on other parts of the species range, and according to our previous experience. The volume of the playback broadcast was enough to be heard by an observer at 300 m but not as much as to produce distortion noises at close distances.

Variables for the analyses

The number of individuals (abundance, from 0 to a maximum of 9 individuals) of each species detected per SP was considered the response variable, a zero-inflated Poisson random variable (Table 1). We recorded two types of predictive variables, those that could affect species detectability (survey covariates) and those that could affect our ecological response variables, i.e., species abundance (site covariates).

We chose five survey covariates that could affect detectability. We noted whether individuals responded to conspecific broadcast voices or produced spontaneous calls (BROAD). We also considered the experience surveying owls of each observer (EXPER), being “0” for those researchers who had no experience surveying owls and “1” for those researchers that had developed owl surveys (see e.g., Zabala et al., 2006; Zuberogoitia, Zabala & Martínez, 2011). We considered linear and quadratic effects of the Julian days (1 January = 1, 31 December = 365; DATE, DATE ^˄2) to control for seasonal effects. Survey hour (HOUR) was measured counting each hour from the sunset in the study area (https://tierra.tutiempo.net/calendario/calendario-solar-de-euskadi-sp019021.htm) to the survey moment. Surveys should have been conducted on dry and calm nights, although it was impossible to effectively control some parameters as wind (WIND) and temperature (T), and therefore we obtained detailed information about the wind speed (km/h) and T (ºC) of the 15 min survey period at each SP. We obtained these data from the nearest meteorological stations (n = 27) to each SP (http://www.euskalmet.euskadi.eus).

We also selected 11 site-specific covariates that could affect species abundance. The regional climatic situation (REG) of the SP, being Cantabric region (1), Subcantabric region (2) and Mediterranean region (3) (for more details see Zuberogoitia, Zabala, & Martínez, 2011) and the altitude (ALT) m a.s.l. of the SP. Moreover, we established a 565 m radius obtaining 1 km² buffer area for each SP. The 521 buffer areas were overlapped with vegetation and urban digital maps (www.geoeuskadi.eus) to obtain the percentage of vegetation types (FIELD- grass-fields; AGR- agriculture area; SCRUB- scrub and heather areas; DEC- deciduous forest; HOLM- Holm-oak forest; PINE- pine plantations; EUC- Eucalyptus plantations; FOR- forest surface, including timber plantations; URB- Urban area). For open areas species, barn owl and little owl, we joined forest types (DEC, HOLM, PINE and EUC) in a unique covariate (FOR), whereas for the rest of species we considered the forest types in the analysis.

Data analysis

We developed binomial N-mixture occupancy models, in which we estimated abundance and detectability as a function of site-specific and survey-specific covariates using the log link function (Fiske & Chandler, 2011). Our sampling design considered 521 sites (SPs) in which we recorded the numbers of individuals for each species, considering from 4 to 8 temporal replicates for sites. Observations were generated by a combination of (1) a state process determining abundance (i.e., counts) at each site and (2) a detection process that yields observations conditional on the state process (MacKenzie et al., 2002; Royle, 2004; Kéry, Royle, & Schmid, 2005).

Given the large number of potential candidate models to evaluate abundance and detection probabilities, model fitting was conducted following a two-phase approach. First, we performed model selection for detectability models in a hierarchical process considering all the possible combinations of the survey-specific covariates, including the null model, and keeping the abundance component constant; second, we performed abundance model selection in a hierarchical process considering all the possible combinations of the site-specific covariates, keeping the component for detection probability constant (Comte & Grenouillet, 2013). For abundance models, we rescaled the variables by subtracting the mean and dividing by the standard deviation (Hedlin & Franke, 2017) and we considered the quadratic function of the urban area (URB^˄2). Finally, we performed model selection with the combination of the best detectability and abundance models. Models were ranked using the difference in Akaike’s information criterion (Δ AIC) between each model and the best model supported (i.e., the model with the smallest AIC; Burnham & Anderson, 2002). We used the pcount function from the unmarked package (Fiske & Chandler, 2011), considering latent abundance as a Poisson distribution, and detectability as binomial distribution, as well as the identifiability problems described by Kéry (2018). We generated abundance maps of the species in R using the package “raster” (https://CRAN.R-project.org/package=raster).

We included the raw database used for the analysis for the six species considered (first sheet) and also the original database (in Spanish) with other interesting data of all the owl species and European nightjar. Locations of all the species were also included in www.ornitho.eus