Monitoring of bird populations
We monitored nest box populations of great tits during the breeding season in 2013–2020 in five urban parks. The parks were located within the central area of the city of Malmö (55°35’24”N 12°59’19”E; Fig. 1a), which is the 3rd largest city in Sweden with 350 000 inhabitants (SCB 2021). The parks ranged between 3 to 45 ha in size and housed a total of 400 nest boxes, spread evenly throughout them (Fig. 1b). The parks were characterized by a mixture of tree species, amenity grass, ponds and urban infrastructure such as paved roads, paths, lamp-posts and buildings. While great tits account for most breeding attempts in the nest box population, a significant number of blue tits also utilize the boxes, whose breeding season overlaps with that of great tits. A few observations of other passerine birds have also been recorded (see Appendix S2).
During the breeding season of great tits, which locally occurs April-June, the nest boxes were checked weekly. When eggs were found, the laying date was estimated by back-calculating from the number of eggs, as great tits typically lay one egg per day (Perrins 1970). When incubation was confirmed, the date of incubation start was calculated in a similar way, using the clutch size to estimate when incubation likely began. Nests were then left undisturbed for 12 days from the estimated incubation start and subsequently checked daily for hatching. Once nestlings reached 14 days from hatching, they were weighed and ringed. In 2016, not all nest boxes were monitored due to logistical constraints, and we therefore excluded the data from this year.
Territory mapping
To investigate which aspects of the local tree composition affected the breeding success of urban great tits, we used the georeferenced tree database of Malmö, covering publicly managed land, as of August 2019 provided by the city of Malmö (pers. comm. Tim Delshammar). This database contains spatially explicit information on individual tree species, age, height, and crown radius. Using R packages rgeos and sp (Pebesma and Bivand 2005; Bivand et al. 2013; Bivand and Rundel 2020), we created a GIS polygon layer of the spatial coverage of all tree crowns by creating individual buffers of the crown radius for each tree based around the geographical position of the tree. Using a 35 m radius around each nest box, all tree canopies reaching within the radius were identified as part of the local tree composition. The 35 m cut-off point was chosen a priori as tits have been shown to hold territories within this area and tend, on average, to travel this distance to forage (Krebs 1971; Stauss et al. 2005; Jarrett et al. 2020).
We characterized urban territories based on a number of hypothesis-driven factors of the biotic environment: the number of common oak trees (Quercus robur), European beech trees (Fagus sylvatica), birch trees (Betula pendula and B. pubescens), tree diversity and number of non-native tree individuals. We included oak trees as they have been found to be a primary source of caterpillars and modulate the breeding phenology of great tits in forest populations (van Noordwijk et al. 1995; Visser et al. 2006). We included beech trees as this was by far the most common species within the parks (Jensen et al. 2021). Birch species were also common and are known to support relatively high caterpillar productivity (Eeva et al. 2000).
An inverse Simpson diversity index was included as we hypothesized that tree diversity could increase the number of potential different food sources within the territory (Ehrlich and Raven 1964). The abundance of non-native tree species was included as they host significantly fewer invertebrates compared with native trees (Jensen et al. 2021) and have been found to negatively affect avian reproduction (Narango et al. 2018). Non-native trees were defined as species introduced after the 13th century (Essl et al. 2018). While the tree composition could potentially change over the 8-year period the study was conducted, the average age of the trees located within the parks is over 80 years, and it is therefore unlikely that any major changes in the territories had occurred within the time-window of our study.
Statistical analysis
All statistical analyses were conducted using R version 4.0.3 (R Core Team 2020). We constructed Generalized Linear Mixed Models (GLMMs), using the glmmTMB function from the glmmTMB package (Brooks et al. 2017), to analyze the correlation between local tree composition and breeding attempts, breeding onset, offspring survival and weight. We selected models through backward elimination whilst respecting marginality between interactive terms (see below; (Zuur et al. 2007, 2009), where all models had the same initial structure where possible. We inspected model residuals manually, in addition to using the function testResiduals of the DHARMa package (Hartig, 2020), to verify model assumptions (variance homogeneity and acceptable normality of model residuals).
To investigate which aspects of the tree composition influenced the likelihood of a breeding attempt, data from all boxes from the 7 years was used, totaling in a sample size of 2 800 cases. Note that we considered breeding attempts from both blue and great tits for this specific measure only, since not all breeding attempts could be confidently assigned to a species (e.g. if an observer could not identify nest ownership from eggs alone). Moreover, we expected interspecific competition for nest boxes, given the overlapping breeding seasons. A nest box was hence considered occupied when at least one blue or great tit egg was observed for the given year. We constructed an initial GLMM, with a logistic binomial distribution, with breeding attempt as the response variable and number of oaks, beech and birch trees, number of non-native trees, tree diversity, and year as fixed effects. The two-way interactions between year and all other main effects were also included. Neither lay date nor its interactions were included in the analysis on breeding attempts, as lay date could not be determined in all cases for this subset of data. The final GLMM, selected through backward elimination, included year, number of European beech trees and non-native tree individuals. Nest box was nested within park as a random structure, to account for the repeated sampling across years and potential differences between parks.
The response variable in the model for effects of local tree composition on breeding onset was the mean-centered laying date (day of first egg) per year for all great tit clutches that reached incubation, resulting in a sample size of 464 clutches. We only included eggs from completed clutches, as the species can be hard to determine up until this point, when the identity of the female is confirmed. The initial model had the same structure as above: number of oaks, beech and birch trees, number of non-native tree individuals, tree diversity, and year, together with the two-way interactions between year and all other fixed effects, were included. The final GLMM consisted of the number of common oak and European beech trees as fixed effects and nest box nested within park as a random factor. A Gaussian distribution was specified for this model. Inspection of model residuals revealed that variance was not homogenous between years in the model and the dispformula argument of the glmmTMB function was used to correct for the heteroscedastic variance (Brooks et al. 2017).
Great tit offspring survival, defined as the probability of an egg reaching day 14 as a nestling, was studied for all eggs of completed clutches where incubation onset was confirmed. In total, 3 386 eggs from 464 clutches were included in the analysis. The initial model included number of oaks, beech and birch trees, number of non-native tree individuals, tree diversity, year and lay date (mean centered per year) as fixed effects. The two-way interactions between year and lay date and all other factors were also included as fixed effects. The final GLMM included year, number of oak trees, tree diversity and lay date, together with the interactions between lay date and year, oak trees and lay date, and tree diversity and lay date. A random factor of clutch ID nested within park was included. A logistic binomial distribution was specified for the model.
To investigate the effect of local tree composition on nestling weight, the individual weight at day 14 of 1 270 great tit nestlings, belonging to the 233 broods that reached this point in the 7 years, were included. The initial model had the same structure as above: number of oaks, beech and birch trees, number of non-native tree individuals, tree diversity, year and lay date (mean-centered per year) together with the two-way interactions between year and lay date and all other main effects. The final model, selected through backward elimination, included year, number of beech trees and non-native trees and lay date, together with the interactions between lay date and year, as well as beech trees and year, as fixed effects. A random factor of clutch ID nested within park was also included. A Gaussian distribution was used for the model. Inspection of model residuals revealed that variance was not homogenous between years in the model, and as above, the dispformula argument of the glmmTMB function was used to correct for the heteroscedastic variance (Brooks et al. 2017).
Since Park 5 had a significant proportion of the local non-native tree abundance located near the edge of the park, close to a road, we excluded the nestboxes along the park edge (13 nestboxes east of the main waterbody; Fig. 1b) from the models where non-native trees were included as a fixed effect (breeding occupancy and nestling weight). However, this potential spatial confounder did not alter the results qualitatively when excluded, and the full data set was therefore used. Also note that there were no significant difference in age or size between native and non-native trees in the study system (Jensen et al. 2021). Clutch size was not included in models for nestling weight or offspring survival; all significant results were however retained for both final models when we added clutch size to control this did not influence the results. To control for co-linearity between factors, variance inflation factors (VIFs) were calculated for all initial models, excluding interactions, using the check_collinearity function of the performance package (Lüdecke et al. 2020). No VIFs > 2 were found (Appendix S3), indicating negligible co-linearity between fixed effects (Zuur et al. 2007). Models were compared with the anova function in the backward elimination process (Appendix S4). P-values were obtained with function Anova (car package) on final models using Type III Wald chi-square tests. R2-values were calculated through the r2 function of the performance package.