Provisioning latency increases closer to roads and is associated with species-specific reproductive success in two urban adapters

Most research on urban avian ecology has focused on population- and community-level phenomena, whereas fewer studies have examined how urbanization affects individual behavioral responses to a sudden and novel stimulus, and how those translate to fitness. We measured between-individual variation in provisioning latency in two urban adapters - great tits and blue tits - in response to an infrared camera installed in the nestbox, encountered when offspring in the nest were at the peak of food demand (9–10-days old). For each nestbox, we quantified urbanization as intensity in human activity, distance to road and proportion of impervious surface area. In both species, provisioning latency increased significantly closer to roads. Moreover, increased provisioning latency when exposed to a novel object was associated with higher reproductive success in great tits whose nestboxes were surrounded by high amounts of impervious surface. In contrast, increased provisioning latency was consistently associated with lower reproductive success in blue tits. Our results suggest that provisioning latency changes in relation to the environment surrounding the nest, and may be context- and species-specific when exposed to a novel stimulus, such as a novel object in the nest. To better understand the role of initial behavioral responses towards novelty across an individual's lifetime and, ultimately, its impact on fitness in the urban mosaic, further research explicitly testing different behavioral responses across the entire breeding cycle in wild model systems is needed.The file entitled "Repeatability_data" contains data used to test repeatability in provisioning latency. (Statistical analyses, point 1) The file entitled "Data_Provisioning" contains all data used for the analyses described in 2.6 (Statistical analyses, points 2,3 and 4). The file README contains detailed info on each variable reported in the datasets mentioned above. Missing data are recorded as "NA".Funding provided by: Narodowa Agencja Wymiany AkademickiejCrossref Funder Registry ID: http://dx.doi.org/10.13039/501100014434Award Number: PPN/IWA/2019/1/00069Funding provided by: Narodowe Centrum NaukiCrossref Funder Registry ID: http://dx.doi.org/10.13039/501100004281Award Number: 2017/25/N/NZ8/02852Funding provided by: Narodowe Centrum NaukiCrossref Funder Registry ID: http://dx.doi.org/10.13039/501100004281Award Number: 2014/14/E/NZ8/00386Funding provided by: Narodowe Centrum NaukiCrossref Funder Registry ID: http://dx.doi.org/10.13039/501100004281Award Number: 2016/21/B/NZ8/030822.2 Breeding data In both years, nestboxes were inspected at weekly intervals from mid-March until the end of June to determine first-egg date and clutch size. Hatching date (recorded as day 1) was determined by checking the nestboxes one day before the expected hatching (c. 12 days after the last egg of the clutch was laid), and every two days after that, until hatching occurred. When nestlings were between 9 and 10 days of age, infrared-cameras were installed inside the nestboxes as part of a project investigating avian diet across the urban mosaic (Corsini et al., in prep). Parents were caught, measured and ringed when nestlings were ca. 11-13 days old, and aged on plumage characteristics as first-year breeders (First year) or older individuals (Older). At day 15, all nestlings were measured and ringed. Ca. 25 days after hatching, nestboxes were checked to record the number of nestlings that died before fledging, and to calculate fledging success as the proportion of eggs producing fledglings. 2.3 Measuring individual neophobia Neophobia was measured as an individual's latency time to first entrance (FE) following its first approach (FA) of the nestbox after installing the infrared camera (see Figure 2). To account for the influence of disrupting factors in the nestbox surroundings, we also recorded the number of times an individual flew off before its first entrance (Number of flights). The infrared-cameras (model SONY HK100427) were installed in the morning (between 6:26 and 9:33 AM) when nestlings were 9 or 10 days old. The time (in minutes) between each bird's first approach (FA, Figure 2, Figure S2) towards the camera (i.e., standing or flying next to the nestbox entrance, facing the camera) and its first entrance (FE, Figure 2, Figure S2), was used as a proxy of neophobia for comparison with previous studies (see 37,50). Because of clear plumage differences between male and female great tit, we could record first entrances separately for each parent. As this was not feasible for blue tits – infrared cameras recorded black and white videos – we noted latency only for the first bird entering the nestbox. Total number of videorecorded nestboxes by species, study site and year is reported in Table S1. 2.4 Characterizing urbanization We characterized the amount of urbanization surrounding each nestbox by measuring human or canine presence within 15m, distance to the closest road and percent area covered by impervious surface (i.e., % ISA). Briefly, human presence was obtained from repeated surveys on the ground aimed at detecting all humans and dogs within a 15m radius around each nestbox during 30 seconds long counts. This value was then standardized by dividing the total amount of humans and dogs detected by the number of counts performed within each study site (detailed in 48). The distance to the closest road (i.e., for vehicular use) was recorded in meters using the distance matrix tool in QGIS (see methods in 49). Percentage of ISA was derived from opensource remote sensing imagery data (Copernicus Land Monitoring, https://land.copernicus.eu/sitemap); specifically, we used the "imperviousness" layer - which included all built up areas and was quantified using a 20m-pixel resolution raster file - averaged in a 100m radius around each nestbox in QGIS using the Zonal statistics tool, as detailed in 51. A Principal Component Analysis on these three proxies of urbanization after their standardization, revealed that PC1 explained 60% of the variance. PC1 comprehended the variables Distance to Road and ISA as main contributors (Figure S3a and b). 2.5 Weather data Weather data over the two-year period were provided by the Polish Institute of Meteorology and Water Management (IMGW – PIB). Daily temperature (°C) and rainfall (mm) were derived from two climatic stations: Warsaw Okecie and Legionowo (here used as references for study sites located within and outside the city, respectively). For each nest, we used the average-daily temperature and the rainfall on the date of recording. 2.6 Statistical analyses All statistics were computed in R (v. 3.6.2, 52), R-packages and functions are detailed below. Figures and plots were generated through the R-package "ggplot2" (v.3.2.1 53) and the opensource software Inkscape (v. 1.0.2 54). In the full dataset, two great tit nests (4 birds) and one blue tit nest (1 bird) were identified by ring numbers as second broods, thus, excluded from all the analyses, as birds were already exposed to the camera from the first breeding attempt that occurred in the same year. As in great tits, the second bird of the pair was often influenced by the partner while approaching the camera (i.e., the second bird entered the nestbox faster after the first one did, MC and PL personal obs., Figure S4), the behavioral dataset of both species was restricted to birds that accessed the nestbox first. Statistical analyses were performed as follows: 1) We used the Test-Retest approach to infer whether neophobia was consistent within individuals measured in two distinct breeding years. Out of all birds measured, 9 banded great tits were caught in 2018 and 2019 (N = 18 observations). 2) We used a series of One-Way ANOVAs to test whether neophobia (in min) differed between species, age-classes and sexes. We subsequently tested species-specific differences in terms of their exposure to human presence, distance to roads and % of ISA in the nestbox surroundings. We further tested whether neophobia was associated with the number of times a bird flew away from the nestbox before the first entrance. 3) To test whether neophobia is mediated by urbanization, we built two distinct Linear Mixed – Effects Models (LMMs, "lmer" function in the R-package "lme4", v. 1.1-21) where both species were analyzed jointly, and were characterised by the same model structure except for the proxy of urbanization, which was fitted as explanatory variable as either Distance to Road or ISA. In each model, Neophobia was log-transformed to meet models' assumption and fitted as the response variable, while the variables: Distance to Road (or ISA), Human presence, rainfall in mm (Rain), Temperature and Date of Recording (1st of April coded as 1), were fitted as continuous predictors, and the variables Species and Year as two-levels factors. Site was set as random effect to control for non-independence of nestboxes sampled within the same study area. Because the variables Temperature and Date of recording were positively correlated (rPearson (113) = 0.67, p < 0.001), Temperature was removed from the global models. We used the vif function in the R-package "car" (v. 3.0-11) 55 to check possible multicollinearity issues between continuous predictors. Despite the fact that Variance Inflation Factors (VIF) were below 2 for Distance to Road and ISA, we detected a strong negative correlation between these two terms (r = -0.6; see Pearson's correlation reported in Table S2). Consequently, these two proxies of urbanization were always analyzed separately; in contrast, because Human presence was only weakly correlated to ISA and Distance to Road (r = 0.4 and r = -0.2, respectively, see Figure S3a and Table S2), it was fitted as continuous predictor in each global model. Models comparing the separate effects of Distance to Road and ISA on neophobia were also used because roads may have a major impact on wildlife (e.g., via habitat loss, but also due to the increased levels of sound, light and chemical pollution associated to them, which may extend far away from their location 56). Thus, analyses and results relative to each proxy are reported and discussed in parallel. We performed backward model selection in each global model (via the step function in the R-package lmerTest v. 3.1-3 57). We used the DHARMA package to check the global models fit and assumptions 58 (global model diagnostics shown in Figure S5). 4) To test whether neophobia in parents was associated with reproductive success and whether this relationship was mediated by urbanization (here, again, separately tested as Distance to Road or ISA), analyses were run for the two species separately, to account for species-specific differences in life-history. We ran Generalized Linear Mixed-Effects Models (GLMMs, using the glmer function in lme4), where the number of fledglings weighed by clutch size (i.e., here referred to as Fledging Success) was fitted as response variable using a binomial distribution. In great tits, one nest was excluded from the analyses (due to missing information on clutch size, N = 64). The interaction-term between neophobia and the proxy of urbanization (as ISA or Distance to Road, as explained in point 2), along with the continuous variables Date of Recording, Rain and Year, were fitted as predictors in each model. The factor Site was fitted as random effect to control for non-independence of avian fitness in nests sampled within the same study site. With the only exception of the great tit global model testing whether neophobia and fledging success were mediated by Distance to Road (where the inter group site-variance was zero, thereby necessitating the removal of random effects in the model), all 3 remaining global models included Site as random factor (Table 1, Table S7). We checked for multicollinearity issues: VIF- scores were calculated in all models and never exceeded 2. Each global model was used to generate a set of models with all possible combinations of the predictors (R-package MuMIn v.1.43.15, see 59). The best models were graded according to the Akaike's information criterion 60 (AICc). Model-averaged coefficients for a subset of models (∆AICc < 2) were extracted. As some Akaike weights of best models were below 0.9, and therefore high model selection uncertainty existed, full-model averaging was used 61. We calculated upper and lower bounds of the 95% confidence intervals (CI) for each parameter. We applied Z-score scaling to each continuous predictor for clarity of parameters estimates