The city of Xalapa is located in central Veracruz, Mexico (19° 31′–19° 36′ North, 96° 55′–96° 59′ West, elevation 1120–1720 m; rainfall 1100–1600 mm/year). The urbanized area of this city is ~64 km2 and has 488,531 inhabitants (INEGI 2020). Botanists have reported 140 tree and shrub species in the streetscape of Xalapa, with exotics representing more than half of the recorded species, and native species distributed unevenly throughout the city. Among the most conspicuous tree species are weeping fig (Ficus benjamina), Chinese hibiscus (Hibiscus rosa-sinensis), paper flower (Bougainvillea glabra), rhododendron (Rhododendron sp.), Mexican cypress (Cupressus lusitanica), and little-leaf boxwood (Buxus microphylla). Despite its exponential growth and substantial vegetation loss, Xalapa is home to more than 340 bird species that includes 31 raptor species (González-García et al. 2014, 2016), and is considered as a natural laboratory for ecological studies by urban ecologists (MacGregor-Fors et al. 2015). Recent studies show that Xalapa’s woodland covers around 40% of its territory (Falfán et al. 2018). Land uses in Xalapa include cloud forest remnants, coffee, sugarcane and corn fields, grasslands, wetlands, water bodies, and urban areas (Benítez Badillo 2011; Lemoine-Rodríguez et al. 2019).
Urbanization gradient and sampling sites
To characterize the urbanization gradient, the city's vegetation was categorized based on the Normalized Difference Vegetation Index (NDVI) from a 2019 Sentinel 2A satellite image (ESA 2019). The NDVI allows for the generation of an image showing the greenness (vegetation reflectance indicative of relative plant biomass). This index takes advantage of the contrasting characteristics of two bands of a multispectral image: the chlorophyll pigment absorptions reflect in the red band and the high plant reflectivity in the near-infrared band (Nageswara Rao et al. 2005). This process was performed in ArcMap 10 (ESRI 2012) using the NDVI tool from the Image Analysis module.
First, it was necessary to join the spectral bands of the Sentinel satellite image. For this purpose, we used bands 2, 3, 4, and 8 with a spatial resolution of 10 m. Second, we cropped the resulting image over the urban area of Xalapa at a scale of 1:68,000 to exclude areas that were not of interest for the analysis. We generated a polygon at a scale of 1:80,000 to make the NDVI corresponding to the urban boundary of Xalapa that was modified from the polygon of Lemoine-Rodríguez et al. (2019) and INEGI (2010).
The highest percentage of urbanization obtained from the NDVI was 77% since this city’s gray, impervious surface has interspersed vegetation on city parks, ravines, ridges, and sidewalks. For this reason, the NDVI does not reach 100%. From this urbanization continuous value, we generated 10 categories of the gradient (from 1–10), where 1 is 77% urban cover, and 10 is <7.7% urban cover. It is worth mentioning that we decided to include in the 10th category two preserved areas completely away from the urbanization of the city to obtain the highest possible NDVI greener value. Subsequently, we generated 2 random sampling sites in ArcMap for each urban gradient category that were at least 2 km apart (20 sampling sites in total; Figure 1). For each sampling site, we delimited a 500 m radius around the center using ArcMap 10, in order to characterize the surface composition (in ha) of the urbanized to green area of each site.
A pilot test was carried out during May 2019 in the selected sampling sites in order to evaluate the logistics of the project and preliminarily detect which raptors were present. We visited each site two times during three field seasons: breeding (from June–July 2019), fall migration (from August–November 2019), and winter (from December 2019–February 2020). We surveyed each site twice a day, from dawn–11:00 and from 18:00–23:00 h (Bird et al. 2018).
For diurnal and nocturnal raptor censuses, we conducted systematic sampling at each of the selected 500 m radius circular plots (Bildstein and Bird 2007). For species identification, we used (Clark and Schmitt 2017), cellphone applications such as iNaturalista (CONABIO 2019) and Merlin (The Cornell Lab of Ornithology 2020), 10×42mm Celestron Nature Dx binoculars, and a Nikon D7500 camera with a 200–500 mm lens. Whenever possible, we took photographs to document our records.
For nocturnal raptors, whose study is difficult because of their nocturnal and elusive habits that make direct binocular observations impossible (Aguilar et al. 2001), we used audio playback recordings. We obtained recordings from the xeno-canto.org repository (Hardy and Morrison 2000; Enríquez and Rangel-Salazar 2001) and the Merlin application for cellphones. When possible, we used flashlights to allow the observation and identification of the species.
At each observation site we recorded the site number, urban gradient category of the site, survey number, day of the survey (hereafter “day”; recorded as an ordinal number), time of the day (diurnal or nocturnal), site coordinates, altitude, raptor season (breeding, fall migration, and winter), estimated vegetation proportion (native and exotic), presence or absence of water bodies, pedestrians, pets, raptor species, type of detection (visual or auditive), detection/non-detection of the species, behavior (hunting, perched, nesting, courtship, flying, and gliding), and perch type used (tree or building; Berthiaume et al. 2008; Conway et al. 2009; Dowling et al. 2012). To document pedestrians and pets at each sampling site, we assigned a subjective category from 0–4 (e.g., pedestrians: where 0 was the absence of the variable, 2 will be moderate, and 4 intense activity; Hogg and Nilon 2015).
A widely used sampling design to determine the factors or habitat characteristics that determine the presence of species in urban areas is to repeatedly survey a selection of sampling units and record detections and no detections (Rota et al. 2016). Failure to detect a species may be due to the species being truly absent from the sampling unit or simply remaining undetected MacKenzie et al. 2002). Binomial zero-inflation models (hereafter referred to as occupancy models; MacKenzie et al. 2017) are useful in modeling the factors that influence the probabilities of species occupancy in the face of imperfect detection (Melles et al. 2003; Kroll et al. 2007; Henneman and Andersen 2009; Hansen et al. 2011; Hogg and Nilon 2015; Rota et al. 2016).
We assigned each species to an abundance category that expresses the probability of detection of the species, according to Pugnali and Chamorro (2006). We also calculated relative abundance (RelAb) with the number of observed and heard individuals of each species divided by the number of observed individuals of all species at each sampling site (RelAb = (#ind_sp)/ (#ind_total per site). The categories were accidental: 0–1%, rare: 1–9%, scarce: 10–30%, common: 31–64%, frequent: 65–89% and abundant: 90–100%.). This result, multiplied by 100 (RelAb×100) corresponds to the percentage that each species contributes to the community (Smith and Smith 2001).
We used the Bray Curtis distance measurement (QS) to obtain the dissimilarity values between sampling sites and urbanization gradient. When used in presence/absence data, the Bray Curtis index is known as Sorensen's measured distance. QS is the similarity ratio and varies from 0–1. This expression easily extends to abundance rather than the presence/absence of species. The value 0 indicates full dissimilarity between the areas sampled because they do not share species, and the value 1 indicates total similarity; intermediate values are considered as <0.25 Low, 0.26–0.5 Moderate, 0.51–0.75 High, and 0.76–1 Total similarity (Bray and Curtis 1957). Raptor diversity along the urban gradient was calculated using Shannon Exponential Index (eH´), which gives all species an equal weight to their relative abundance, without favoring or penalizing either rare or frequent species: q = 1. This index is expressed as a positive number, which in most ecosystems vary between 0.5–5, although its normal value is between 2–3; values <2 are considered low, and >3 are high (Jost and González-Oreja 2012).
Modeling habitat associations
We used sampling data from each site to construct detection probability and occupancy probability models in a two-stage theoretical information modeling approach (Hansen et al. 2011; MacKenzie et al. 2017). In the first stage, we modeled the detection probability for each sampling site with each of the variables recorded per site. In the second stage, the occupancy probability was modeled by species at each sampling site by incorporating habitat characteristics to find the best fit for each species. We analyzed the candidate models with the Unmarked library in R (Fiske and Chandler 2011) and the best models were selected using the Akaike Information Criterion with correction for small samples (AICc). We also calculated Akaike's weight for each model (wi).
Following that, we calculated the over-dispersion parameter (ĉ) over 1000 parametric bootstraps of the R-occupancy data. This procedure generates 1000 random detection data sets and adjusts the global model for each of them by producing 1000 statistical Chi-square (χ2) tests. The P-value of this procedure refers to the probability that each statistical test is greater than or equal to the chi-square value calculated by the global model adjusted to the currently observed dataset (χ2obs). The dispersion parameter (ĉ) is the average of all statistical Bootstrap tests divided by the observed statistical tests (χ2obs / χ2). A value of ĉ near 1 was considered to be a well-fitting model for the data (Hogg and Nilon 2015). A value ĉ greater than 1 indicated that the data were overestimated, or that there was greater variation in the observed data than the overall model fitted to the species' occupation. We modeled the probability of detection with the variables that were considered that could influence the presence of these birds at each sampling site: day, people, pets, and proportion of native/exotic vegetation per site.
The second stage involved the creation of a set of candidate models to test the effects of measured habitat variables on observed variation in occupation (ψ) between sampling sites, given the modeled detection probabilities (Kroll et al. 2007; Henneman and Andersen 2009; Hansen et al. 2011). The habitat variables used were altitude, NDVI, presence/absence of water bodies, and composition of the 500 m radius sampling sites (urban/vegetation area in hectares).