Are behaviour and stress-related phenotypes in urban birds adaptive?

Urbanisation is a world-wide phenomenon converting natural habitats into new artificial ones. Environmental conditions associated with urbanisation represent great challenges for wildlife. Behaviour and stress tolerance are considered of major importance in the adaptation to novel urban habitats and numerous studies already reported behavioural and stress response phenotypes associated with urbanisation, often suggesting they represented adaptations, while rarely demonstrating it. The main goal of this study was to test the adaptive nature of urban shifts in behavioural and stress-related traits, and by adaptive we mean phenotypic change favouring traits in the same direction as selection. Using 7 years of monitoring of urban and forest great tits, we first tested for differences in exploratory behaviour, aggressiveness and breath rate, between both habitats. Second, we performed habitat-specific analyses of selection on the three former traits using (a) reproductive success and (b) survival estimated via capture-mark-recapture models, as fitness estimates, to determine whether shifts in these behavioural and stress-related traits were aligned with patterns of ongoing selection. We found that urban birds displayed higher exploratory behaviour and aggressiveness, and higher breath rate, compared to forest birds. Selection analyses overall revealed that these shifts were not adaptive and could even be maladaptive. In particular, higher handling aggression and higher breath rate in urban birds was associated with lower fitness. Higher exploration scores were correlated with lower survival in both habitats, but higher reproductive success only in forest males. Overall, differences in patterns of selection between habitats were not consistent with the phenotypic divergence observed. Taken together, these results highlight that phenotypic shifts observed in cities do not necessarily result from new selection pressures and could be maladaptive. We hypothesise that divergences in behavioural traits for urban birds could result from the filtering of individuals settling in cities. We thus encourage urban evolutionary scientists to further explore the adaptive potential of behavioural traits measured in urban habitats (a) by replicating this type of study in multiple cities and species, (b) by implementing studies focusing on immigrant phenotypes and (c) by measuring selection at multiple life stages.

INTRODUCTION

Human‐induced environmental change is a world‐wide phenomenon taking place at an unprecedented rate. These changes include habitat degradation, climate change, increased presence of invasive species, pollution and over‐exploitation of resources, and affect biodiversity in a pervasive way (Pelletier & Coltman, 2018). All species living on our planet have faced environmental changes during their past evolutionary history, yet the present extent and rate of human‐induced alterations represent unprecedented novel challenges. Although a wide range of traits can be concerned, the first response of animals towards rapid environmental changes often implies behavioural shifts on a multitude of facets in response to different stimuli (reviewed in Tuomainen & Candolin, 2011). For instance, behavioural shifts can affect foraging strategies, mating success or predator avoidance behaviours (Lowry et al., 2013), hence they are closely linked to individual fitness. Consequently, behaviour is pervasively considered as a crucial factor determining how animals cope with environmental changes and new selective pressures (Baldwin, 1896; Mayr, 1963). In addition to behaviour, traits that are related to stress responses may play an equally important role in wild animal tolerance and adaptation to human‐induced environmental changes. Indeed, anthropogenic environmental change often results in highly altered new habitats with high levels of perturbation (Sih et al., 2011) and living in such habitats might thus require higher stress tolerance or else induce higher levels of chronic stress that can have deleterious effects on fitness. For instance, vertebrates can cope with environmental challenges by way of activation of the HPA axis which regulates the secretion of glucocorticoid hormones (Sapolsky et al., 2000). However, while some acute stress response occurring immediately after exposure to a stress stimuli (e.g. acute release of glucocorticoid) can promote survival (Landys et al., 2006; Wingfield, 2006), prolonged or repeated exposure to stressors can induce chronic stress that greatly impact fitness (French et al., 2007; Lendvai et al., 2007). Hence, a shift in how organisms handle acute and/or chronic stressors could be of major importance in facilitating adaptation to novel environmental stressors.

Urbanisation is one of the most extreme human‐induced environmental change for biodiversity (Newbold et al., 2015). Urbanised areas combine multiple perturbations on relatively limited areas of land, including deforestation, fragmentation, artificialisation, pollutions and high level of human disturbance (Grimm et al., 2008). While some species are unable to persist in these deeply modified habitats, others manage to tolerate the new conditions or even to thrive on the new resources. Behavioural plasticity has been identified as a key predictor of the successfulness of a species in urban habitats, in particular the capacity to innovate when facing new opportunities (Lowry et al., 2013). Furthermore, at an intraspecific level, some behavioural profiles (also referred as personalities, coping styles or temperament; Réale et al., 2007) can be better suited to cope with novel environmental challenges than others. In particular, is has been suggested that bolder, more exploratory and more aggressive individuals are particularly well suited to the urban habitat (Møller, 2008). In line with these hypotheses, numerous studies in mammals and birds have found that urban individuals indeed, displayed on average bolder (Evans et al., 2010; Prosser et al., 2006; Uchida et al., 2016), less careful (Chapman et al., 2012), more exploratory (Martin & Réale, 2008) and more aggressive (Evans et al., 2010; Scales et al., 2011) behaviours than their forest counterparts. Regarding stress‐related traits (i.e. traits reflecting a physiological response to stress), patterns of responses are less consistent across species and traits (Sepp et al., 2018). Indeed, some studies found higher (e.g. Torné‐Noguera et al., 2014), lower (e.g. Abolins‐Abols et al., 2016) or equal (e.g. Senar et al., 2017) levels of acute response (using breath rate as a proxy) to an induced stress in urban individuals, and this is also the case for other measures of stress response (heterocyte to lymphocyte ratio: Powell et al., 2013; basal corticosterone level: Fokidis et al., 2009; Sepp et al., 2018).

Despite numerous evidence in vertebrates for divergence in behavioural traits and stress responses between urban and rural populations, little is known about the evolutionary implications of such urban‐linked shifts (but see Lambert et al., 2021, Table 1 for a list of studies ‘suggesting’ behavioural adaptations to urbanisation). In particular, to our knowledge no study ever investigated if urban linked shifts in aggressiveness and exploration behaviour or in breath rate, emerged in response to habitat‐specific selective pressures, even though they are largely assumed adaptive in the urban habitat (e.g. Lowry et al., 2013). Behavioural particularities observed in urban habitats could result from adaptive responses to new environmental conditions, such as accepting new food resources or avoiding human disturbance, but they also could be maladaptive and not confer any advantage to individuals. The majority of urban environmental conditions are relatively recent compared to environmental conditions experienced during a species evolutionary history. As a result, phenotypic reaction norms of individuals might not be suited respond in an optimal way to these environmental conditions that were never encountered during species evolutionary histories, and shifts might result in maladaptation (Sih, 2013; Tuomainen & Candolin, 2011). In addition, in such artificial environments cues exploited could be decoupled from the true quality, reinforcing these maladaptive responses (potentially leading to ecological/evolutionary traps, Demeyrier et al., 2016). Maladaptive responses might be particularly common for stress‐related traits in urban dwellers. For instance, in waterbirds, fear of humans perceived as predators can result in temporary or permanent desertion of the nest that can heavily impact individual fitness (Carney & Sydeman, 1999). In the current context of ongoing massive erosion of biodiversity, it is now more than ever necessary to understand demographic and evolutionary dynamics associated with these phenotypic shifts.

TABLE 1

Repeatability of handling aggression, breath rate and exploration score, estimated with rpt() R function using the best model selected (see Table 3), with associated standard error (SE), 95% confidence intervals (CI) and p‐value

Trait	Number of measures	Number of individuals with replicated measures (at least 2)	Repeatability	SE	CI	p
Handling aggression	851	138	0.383	0.056	0.276–0.495	7.2 × 10−10
Breath rate index	874	146	0.440	0.055	0.34–0.547	2.06 × 10−14
Exploration score	850	139	0.523	0.050	0.428–0.627	2.98 × 10−13

TABLE 3

Model selection for handling aggression, breath rate index and exploration scores as a function of habitat, year, sex, age, capture rank, temperature, habitat × year, habitat × sex and habitat × age. All models included individual as random effect and HA and BRI also included manipulator identity as random effect. Models were ranked according to their AICc values, models in black are equally supported (ΔAICc < 2), detailed effects sizes are presented in Table S6.

Models	df	Deviance	AICc	ΔAICc
Handling aggression
habitat + year + sex + habitat × sex	12	3,599.54	3,623.80	0.00
habitat + year + sex + temperature + habitat × sex	13	3,599.46	3,625.70	1.96
habitat + year + sex + capture rank + habitat × sex	13	3,599.54	3,625.80	2.04
habitat + sex + habitat × sex	7	3,612.68	3,626.20	2.37
habitat + year + sex + temperature + capture rank + habitat × sex	14	3,599.46	3,627.80	4.00
Breath rate
habitat + year + temperature + capture rank	12	2,797.74	2,822.10	0.00
habitat + year + temperature + capture rank + habitat × an	17	2,789.33	2,824.10	1.97
habitat + year + sex + temperature + capture rank	13	2,797.73	2,824.10	2.05
habitat + year + age + temperature + capture rank	13	2,797.73	2,824.20	2.05
habitat + year + sex + temperature + capture rank + habitat × sex	14	2,797.06	2,825.60	3.45
Exploration score
habitat + temperature + capture rank	6	4,803.23	4,815.30	0.00
habitat + capture rank	5	4,805.64	4,815.70	0.38
habitat + sex + temperature + capture rank	7	4,801.96	4,816.10	0.76
habitat + temperature + capture rank	5	4,806.03	4,816.10	0.77
habitat + temperature + capture rank + habitat × sex	8	4,799.95	4,816.10	0.79
habitat + sex + capture rank	6	4,804.33	4,816.40	1.10
habitat + sex + capture rank + habitat × sex	7	4,802.34	4,816.50	1.14
habitat + sex + temperature + habitat × sex	7	4,802.78	4,816.90	1.58
habitat + sex + temperature	6	4,804.82	4,816.90	1.59
habitat + age + temperature + capture rank	7	4,802.85	4,817.00	1.65
habitat	4	4,802.97	4,817.00	1.68
habitat + age + capture rank	6	4,805.20	4,817.30	1.96

Our study takes place in a context where differences in behaviour and stress response between urban and forest populations have been previously reported for the focal species, but where eco‐evolutionary implications of these differences remain mostly speculative. Our main objective is to understand to which extent urban shifts in behavioural and stress‐related traits are adaptive. In this study, we focus on three personality traits: handling aggression, exploration behaviour and breath rate under constraint. Handling aggression reflects aggressive behaviour in response to manipulation by humans and distress behaviour (Senar et al., 2017). Exploratory behaviour measures how individuals explore a novel environment and is often used as a proxy for risk‐taking behaviour (e.g. Nicolaus et al., 2012). Breath rate is a non‐invasive measure of primary stress response whereby rapid breath rate is indicative of higher stress response (Carere & van Oers, 2004; Krams et al., 2014). These three traits are expected to be correlated with each other following the pace‐of‐life syndrome hypothesis postulating that personality and physiological traits might have coevolved with life‐history traits (Réale et al., 2010). Individuals are positioned along a slow–fast pace‐of‐life gradient with fast individuals (displaying higher aggressiveness, more exploratory behaviour and higher breath rate) at one end and slow individuals (displaying lower aggressiveness, lower exploratory behaviour and lower breath rate) at the other.

Using a pair of forest and urban populations of great tits Parus major, we first quantified divergence between urban and forest birds for the three studied traits and estimated survival in each habitat using capture–mark–recapture (CMR) models. Since the focal traits can display between‐individual covariance among behaviours—and thereby form so called ‘behavioural syndromes’ (Sih et al., 2004)—we also investigated correlations between each pair of traits. We subsequently performed selection analyses to determine whether phenotypic divergences were aligned with patterns of ongoing linear and quadratic selection in each habitat. Since personality traits might both affect reproduction and survival and not necessarily in the same direction (Dingemanse & Réale, 2005; Smith & Blumstein, 2008), we performed selection analyses using both survival (controlling for imperfect recapture rates using CMR modelling) and reproductive success as fitness proxies. Based on previously reported results (see above), we predicted that urban birds would be more aggressive, more exploratory and have a higher breath rate under constraint compared to forest birds. Following mainstream theories on urban behavioural adaptation (e.g. Lowry et al., 2013), we also predicted that high aggressiveness and more exploratory behaviour would be under positive—directional—selection in the urban habitat, but not higher breath rate.

MATERIALS AND METHODS

Data collection

Population monitoring

A pair of forest and urban populations of great tits Parus major were monitored during the breeding period in two study sites in France. The forest of La Rouvière (43°40′N, 3°40′E) is a Mediterranean forest dominated by downy oaks Quercus pubescens, and with 10% of holm oaks Quercus ilex, located c. 20 km North‐West of Montpellier city. Since 1991, great tit breeding events in nest‐boxes have been monitored in this forest. Nest‐box number has fluctuated between 51 and 92 across years. In the city of Montpellier (43°36′N, 3°52′E) 203–223 nest‐boxes, placed both in urban parks and on streets, have been monitored since 2012 with similar protocols as in the forest. For detailed maps of nest‐box positions see Figure S1.

During the breeding season, nest‐boxes were visited at least weekly, to record brood development from nest building to fledging, providing records of reproductive success for each pair. Adults were captured inside nest‐boxes when feeding 9–15 days old nestlings, measured and ringed with unique metal rings provided by the French CRBPO. Nestlings were measured and ringed at 15 days old (for more details on monitoring protocol see Caizergues et al., 2018).

All experimental protocols were approved by the local ethics committee for animal experimentation of Languedoc Roussillon (CEEA‐LR, 05/06/2018) and regional institutions (bylaw issued by the Prefecture no. 2012167‐003). Captures were performed under personal ringing permits delivered by the CRBPO (Centre de Recherche par le Baguage des Populations d'Oiseaux) for the Research Ringing Program number 369.

Behavioural assays

While breeding monitoring and chick ringing started in 2012 in the city, parental captures started in 2013, and behavioural trials were performed from 2014 onwards. In this study, the last year of data included in the analyses is 2019. Following a parent capture in a nest‐box, each breeding bird was submitted to this sequence of events: it was scored for handling aggression, isolated in a cloth bag for 5 min (resting), measured for breath rate, kept in an acclimatisation compartment adjacent to the main open‐field cage, scored for exploration rate in the open‐field during 4 min, extracted from the open‐field cage, ringed (if not already ringed) and measured for morphological traits. Time between capture and release took on average 20 minutes, and in any case always <30 minutes.

First, right after capture, each bird was tested for handling aggression (HA) by the capturer (see detailed protocol in Dubuc‐Messier et al., 2017 & Table S1). The bird was handled with one hand by the capturer; while facing outward from the manipulator it was nagged at with a finger of the other hand (see Figure S2A) following a standard procedure. After 15 s of observation the capturer attributed a handling aggression score ranging from 0 (unreactive bird) to 3 (aggressive bird striking the handler each time and spreading its wings and tail) with increments of 0.5 (see detailed scoring in Table S1, Charmantier et al., 2017; note that the link between HA and aggressive behaviour in nature requires further exploration). Right after the HA test, the bird was put to rest in a cloth bag.

Second, after 5 min of resting, the bird was removed from the cloth bag, held still on its back by the handler (see Figure S2B) who counted its breath rate. From 2013 to 2016 breath rate was estimated as the number of chest moves during 30 s, while since 2017 the protocol was updated to measure the time to complete 30 chest moves. Measures from 2013 to 2016 were thus converted to approximate the time needed to complete 30 chest moves to obtain the same scale measures. Each count was replicated twice and the handling breath rate was averaged over the two measures. The protocol had a small but significant effect on the breath rate measure (see Appendix S1) and was thus taken into account in the analyses by using residuals of the regression between breath rate and protocol, hereafter referred as breath rate index (BRI). Breath rate is often used as a proxy of acute stress response (Carere & van Oers, 2004; Krams et al., 2014) whereby a smaller BRI (hence higher breath rate) reflects a stronger stress response.

Third, the bird was submitted to a ‘novel environment’ test using an experimental open‐field cage built following Stuber et al. (2013; see Figure S3). This experiment is classically used to estimate exploration behaviour in great and blue tits (Dubuc‐Messier et al., 2017; Stuber et al., 2013). After measuring breath rate, the bird was placed in an acclimatisation compartment right next to the main open‐field cage (Figure S3B, 3) for 2 minutes before being released in the exploration room (Figure S3B, 1). The bird's behaviour was recorded for 4 min. Videos were analysed with the software BORIS (Friard & Gamba, 2016) to count the number of flights and hops during the 4 min. The number of flights and hops r‐squared transformed was used as a proxy of exploratory behaviour, hereafter called exploration score (ES), following previous studies (Dingemanse et al., 2002).

Statistical analyses

Data analyses were performed in four steps: (a) quantifying repeatabilities, between trait correlations and phenotypic divergence between forest and urban great tits for the three traits of interest; (b) estimating habitat‐specific survival probability; and finally estimating the strength and direction of selection on the three studied traits via (c) reproductive success and (d) survival. Analyses of phenotypic divergence and reproductive selection analyses were performed using R software (R Core Team, 2018). Selection analyses via survival were performed using E‐SURGE (Choquet, Rouan, & Pradel, 2009).

Phenotypic divergence, repeatability and between‐individual correlations

Differences between urban and forest great tits in the three behavioural traits described above (HA, BRI, ES) were explored using linear mixed models (r package lme4, Bates et al., 2015). Fixed effects for the full model included: habitat (forest vs. urban), year (categorical variable), sex, age (1 year old (=yearling) vs. 2+ years old (=adult)), rank of capture (i.e. how many times the focal bird was captured until the current capture, ranges from 1 (for first capture) to 6) and air temperature, as well as habitat by year and habitat by sex interactions. For BRI, body mass and hour of the day were also included in the full model since they are known to be highly linked with metabolism. Individual ring number was included as random effect to account for the non‐independence of repeated individual measures, and manipulator identity was also included as random effect for handling aggression and breath rate models.

Model selection was performed using the dredge() function from the MuMIn r package (v1.43.17, Barton, 2020): starting from the complete model all possible models combining fixed effects were automatically run and ranked based on their quasi Akaike criterion corrected for sample size (AICc). Model averaging of the best models (ΔAICc < 2) were performed to obtain average estimate for each effect using the model.avg() function.

Repeatability of each trait was estimated using the rpt() function from the rptR r package (Stoffel et al., 2017), with the best linear mixed model selected using the process described above.

To test for an association between traits forming a ‘syndrome’ we estimated between‐individual correlations across the three focal traits (Dingemanse & Wright, 2020) in each habitat. To measure the strength of a syndrome, it is important to correct for within‐individual correlation or other sources of association between traits when estimating between‐trait correlations (Cleasby et al., 2015). Hence, we estimated between‐individual correlations using bivariate linear mixed models (Dingemanse et al., 2010) computed with MCMCglmm (Hadfield, 2010) and following Ferrari et al. (2013). Models included two of the traits as variable responses, significant fixed effects that were present in the best models selected in the between‐habitat divergence test described above, as explanatory variables, and individual ID as random effect. These models allowed to estimate the variance of each trait and the covariance between two traits ‘1’ and ‘2’, and to decompose (co)variances into individual (VI1, VI2 and CovI1,2), and residual (VR1, VR2 and CovR1,2) (co)variance components for the random intercepts. Correlations between two traits were calculated as CovX1,2/√(VX1 × VX2) (Dingemanse & Dochtermann, 2013), where X represents the component analysed (between‐individual or residual) (for detailed protocol see Appendix S2).

Hence, we calculated two types of correlations. First, the between‐individual correlation, RI1,2, represents effects responsible for the consistency in the correlation between the two traits across individuals and over time; strong RI1,2 across the three traits would be indicative of a behavioural syndrome. Second, the residual correlation RR1,2 represents the within‐individual correlation influenced by the traits' phenotypic plasticity across an individual's lifetime, as well as by trade‐off between traits (Dingemanse & Dochtermann, 2013). Finally, as it is usually reported in studies of behavioural syndrome, we calculated the total phenotypic correlation between each pair of traits to allow comparison with previous studies.

Analyses of survival

Based on capture–recapture data collected during seven breeding seasons (2013–2019), we estimated individual survival and recapture probabilities across habitats (urban vs. forest), ages (yearling vs. adult) and sexes following Lebreton et al. (1992). Survival analyses were performed using a CMR framework with the E‐SURGE software. This framework is based on individual histories of capture (i.e. series of ‘1’ and ‘0’ indicating whether a focal bird was (re)captured or not each year) from which is estimated a survival probability (Φ) while taking into account imperfect detectability estimated via a recapture probability (p). Prior to fitting the CMR model, we tested whether the data met the Cormack–Jolly–Seber's (CJS) assumptions (homogeneity of capture and survival probabilities and independence between individuals), and tested for over‐dispersion with the goodness‐of‐fit tests using U‐CARE (v2.3.2, Choquet, Lebreton, et al., 2009). We initially tested for a time (i.e. year) effect on survival (Φ) and recapture (p) probabilities, but since it was never retained in the best models, it was removed early on during analyses in order to limit over‐parametrisation.

Then, we included effects of habitat, age and sex, using the capture–recapture dataset based on captures of parents between 2013 and 2019. All combinations of models including habitat × age, habitat × sex, habitat, age, sex or i (=constant) effects on the two parameters Φ and p were tested, and model selection was performed using the quasi Akaike criterion (QAICc).

Selection analyses

Via reproductive success

We estimated selection operating via reproductive success for the three focal traits. Since our study system offers little access to lifetime reproductive success, we opted to focus on annual reproductive success via the annual number of fledglings. The male and the female of each pair shared the same annual reproductive success, hence to avoid replicated measures within a single model, males and females were analysed separately.

Following the classic Lande and Arnold (1983) approach, first we quantified directional selection acting on each trait by calculating linear selection differentials. To identify the target of natural selection when studying multiple traits, we also estimated multivariate linear selection gradients β after controlling for indirect selection on the other traits. Linear selection differentials were estimated for each trait following Equation (1), and linear selection gradients were estimated on all traits together following Equation (2):

In both equations, the strength and direction of selection is given by the slope (β or β ) of the regression between relative fitness (ω) and the standardised (per sex) value of a trait (x or x ), α represents the intercept and ε the error. Interaction terms (habitat × x ) indicate whether selection differs between habitats.

Model selection was performed across all models, starting from the complete model described in equations 1 and 2 to a null model (=constant), using AICc. Once more, all models within two points of AICc were considered equivalent and estimates were obtained by a model averaging procedure. All linear mixed models included individual identity and year as random effects to control for repeated measurements on the same individuals and variation across years. To test for correlational selection (i.e. interaction term between traits), the null model of reference contained the additive effect of the three traits and habitat as fixed effect (see details in Appendix S3).

The aim of this study was to explore how directional, hence linear selection, was aligned with the documented phenotypic divergence between forest and city birds. However, nonlinear selection could also differ between the forest and the urban habitats, and thereby shape differences in trait variances. Although differences in variance are not a focal objective here, we estimated quadratic selection differentials and gradients to present a complete picture on comparing selection. Quadratic selection differentials (γ) and selection gradients (γ for each trait x ) as well as correlational selection gradients (γ for traits x and x ) were estimated following these equations:

Where, similarly as in Equations (1) and (2), α represents the intercept, β and β represent respectively the linear selection differential and the linear selection gradient for each trait x , and ε represents the error. While a positive estimate of quadratic selection differential/gradient indicates disruptive selection, a negative estimate is indicative of stabilising selection. More details about methods and results on quadratic selection differentials can be found in Appendix S3.

Via survival

To estimate selection acting on behavioural traits via survival we tested whether an individual's probability of survival from 1 year to the next was associated with the value of the trait, based on capture–recapture data combined with behavioural data collected between 2014 and 2019. As previously described, we tested whether the data met CSJ assumptions, tested for over‐dispersion with the goodness‐of‐fit tests using U‐CARE and performed survival analyses with E‐SURGE. We performed three sets of models, one on each of the traits studied in order to preserve statistical power. In our dataset, most individuals presented more than two measures for each trait. Hence, we used the best linear unbiased predictor (BLUP) for each individual as phenotypic value (e.g. Bergeron et al., 2013). BLUPs were extracted from linear mixed models on repeated measures of each trait, including age, year, hour of the day, air temperature and rank of capture as fixed effects. In the model for BRI, protocol of measure was also included as fixed effect to correct for the bias induced by the change in protocol since 2017 (see Appendix S1).

In this set of analyses, because the number of possible models highly increased with each added effect, we first selected the appropriate model structure for recapture probability (p), starting from the full model p (hab×sex+trait+time) and sequentially removing each effect from the full model. The model with the lowest AICc was selected (for a detailed list of the five best models see Table S2). We then ran models with the same procedure for survival probabilities, starting from the full model Φ (time+habitat×age+hab×sex+hab×sex×trait) with p following the best structure previously selected. A habitat × sex × trait interaction was interpreted as selection on the trait varying between habitats and sexes. Similarly, a habitat × trait interaction indicated that selection on the trait varied between habitats. In addition, to test for potential quadratic selection acting on the traits we also added the trait2 and habitat × trait2 effects in the best 10 models. In the end, models were ranked according to the AICc criterion and all models within a ΔAICc < 2 were considered equivalent.

RESULTS

Repeatable non‐correlated traits

We confirmed that all three traits were repeatable (Table 1). Interestingly, almost all between‐individual phenotypic correlations RI1,2 between pairs of traits, measured in each habitat, had credibility intervals overlapping 0 (Table 2), suggesting that the three focal traits do not covary in a behavioural syndrome (detailed estimates of (co)variance components for each trait or pair of traits are available in Tables S4 and S5). Hence selection acting on one of the focal traits is unlikely to induce indirect selection on the other traits.

TABLE 2

Between individual (RI1,2), residual (RR1,2) (above the diagonal) and total (below the diagonal) correlations between behavioural and stress‐related traits with their lower and upper credibility intervals (RX1,2 = CovX1,2/√(VX1 × VX2)) in urban and forest environments

Phenotypic divergence between urban and forest habitats

When testing for phenotypic differences between habitats, model selection revealed a habitat effect for the three studied traits (Table 3). Urban individuals displayed on average higher handling aggression scores, and this difference was more pronounced in males (meanurban‐males = 2.62 ± 0.23, meanforest‐males = 1.74 ± 0.11; Figure 1a) than in females (meanurban‐females = 2.09 ± 0.14, meanforest‐females = 1.67 ± 0.12; Figure 1a).

FIGURE 1

Divergence between urban and forest great tits on three behavioural and stress‐related traits (a) handling aggression (presented for each sex), (b) breath rate index and (c) exploration score (M ± SE estimated from models accounting for effects of year, sex, age, rank of capture, temperature, habitat × year, habitat × sex, (and body mass and hour for BRI) see Table S6). Note that the breath rate measured corresponds to the time for 30 breaths hence lower values of breath rate index imply a higher breath rate.

Two equivalent models were selected when modelling BRI (Table 3), both including a habitat effect whereby urban individuals had higher breath rate (i.e. a lower index) than forest birds (meanurban = −0.730 ± 0.100, meanforest = −0.32 ± 0.27; Figure 1b).

Model selection for ES revealed 12 equivalent models (Table 3) all including a habitat effect. Urban individuals were more active explorers than forest ones (ES meanurban = 5.88 ± 0.37, meanforest = 3.54 ± 0.72; Figure 1c). Sex, temperature, capture rank and habitat × sex effects were present in numerous models, suggesting they may affect ES. Since so many models were found equally supported, we do not draw any strong conclusion regarding these effects. Detailed effect sizes obtained using model averaging are presented in Table S6.

Habitat‐specific survival

On the overall capture–recapture dataset, a goodness‐of‐fit test conformed to the Cormack–Jolly–Seber's assumptions (χ 2 = 30.6059, p = 0.246, ĉ = 1.177; see Table S3 for detailed GOF test results). Deviances and AICc were corrected for the ĉ value in all models (ĉ = χ 2/df).

In the analysis aiming to estimate survival and recapture for each sex, age and habitat, model selection revealed two models with substantial support (delta AICc < 2, Table 4). For both models, survival probability depended on habitat and age in interaction, with forest yearlings showing a substantially higher survival probability than the three other groups (Table 5).

TABLE 4

Model selection for survival probability in relation to habitat, age, sex (N = 938, full dataset from 2013 to 2019) and Handling Aggression (N = 881) or Breath Rate Index (N = 668) or Exploration Score (N = 661) measured between 2013 and 2019 for the overall dataset and 2014 and 2019 for other datasets. Models are ranked according to QAICc values, and only the best models (i.e. lowest QAICc) are presented for each trait. Equally supported models (ΔQAICc < 2) are represented in bold. For all models (except in the global analysis), the probability of recapture is considered constant (see Table S2).

Models	# Par.	Deviance	QAICc	ΔQAICc
Global dataset
Φ hab.age p .	5	1,689.36	1,441.71	0
Φ hab.age p hab	6	1,689.16	1,443.55	1.85
Φhab.age+sex p.	6	1,689.35	1,443.72	2.01
Φhab.age psex	6	1,689.36	1,443.73	2.02
Φhab.age phabitat+sex	7	1,689.15	1,445.57	3.86
Handling aggression
Φ HA+hab.age+hab.HA 2 p .	8	1,410.61	1,426.61	0
Φ HA+HA 2 +hab.age p .	7	1,412.80	1,426.80	0.16
Φhab+HA+HA 2 p.	5	1,419.33	1,429.39	2.64
Φage+HA+HA 2 p.	5	1,421.68	1,431.73	5.00
Φage+HA+hab.HA 2 p.	4	1,421.66	1,433.73	6.99
Breath rate index
Φ hab.age+hab.BRI 2 p .	7	1,251.68	1,265.79	0
Φ hab.age+BRI+BRI 2 p .	7	1,252.52	1,266.63	0.84
Φhab.age+BRI+hab.BRI 2 p.	8	1,251.66	1,267.80	2.01
Φhab+BRI+BRI 2 p.	5	1,260.52	1,270.58	4.78
Φhab.age p.	5	1,265.65	1,275.71	9.92
Exploration score
Φ ES+ES 2 +hab.age p .	7	1,255.10	1,269.22	0
Φ ES+hab.age p .	7	1,257.55	1,269.63	0.41
Φhab+age p.	5	1,261.50	1,271.56	2.34
Φhab+ES+ES 2 p.	5	1,262.93	1,272.99	3.78
Φhab+ES p.	4	1,264.99	1,273.03	3.82

TABLE 5

Estimates of survival (Φ) and recapture (p) probabilities in forest and urban great tits (yearlings = 1 year olds, adults = 2 years and older). Estimates are obtained from model averaging of the equally supported best models presented Table 4.

Parameter	Age class	Forest		Urban
		Est.	95% CI	Est.	95% CI
Φ	Yearling	0.616	0.488−0.731	0.435	0.355–0.523
	Adult	0.413	0.333–0.497	0.454	0.403–0.505
p	All	0.727	0.632–0.786	0.716	0.637–0.786

Reproductive selection

When estimating reproductive selection on HA, BRI and ES in females using number of fledglings as a fitness proxy, we found no evidence for linear selection on HA nor on ES (Table 6). BRI was under positive significant linear selection in females, that is, females with a lower breath rate had a higher reproductive fitness. However, the difference in selection between urban and forest females for BRI was not significant, as suggested by the habitat × BRI interaction non included in any of the best models (Table S7).

TABLE 6

Reproductive linear selection gradients acting on handling aggression (HA), breath rate index (BRI) and exploration score (ES) in urban and forest great tits, estimated using model averaging for the best equally supported models presented in Table S8. Equally supported models (ΔQAICc < 2) are represented in bold.

Sex	Habitat	Sample size	Linear selection gradients
			HA	BRI	ES
Females	Forest	103	−0.022 ± 0.095	0.132 ± 0.046	0.077 ± 0.103
	City	276	0.078 ± 0.116		0.011 ± 0.119
Males	Forest	93	−0.060 ± 0.050	0.021 ± 0.064	0.230 ± 0.100
	City	247		0.011 ± 0.058	−0.034 ± 0.117

In males, we found no evidence for directional reproductive selection on HA and BRI (Table 6) but positive reproductive selection favouring males with higher ES in the forest but not in the city. This difference in selection between forest and city was significant (habitat × ES interaction included in the two best models see Table S7).

When looking at nonlinear selection gradients (see Appendix S3 and Table S8 for detailed results on selection differentials), quadratic terms were present in most of the best models (Table S9), however, estimates were significant only for BRI females and males and ES in males, suggesting slight stabilising selection acting on female BRI (Table 7), diverging selection acting on male BRI (γ/2 = 0.070 ± 0.031, Table 7), whereby fathers with low or high BRI produced more fledglings; as well as stabilising selection on male ES (γ/2 = −0.088 ± 0.060, Table 7), suggesting that intermediate values of exploration scores were favoured. Finally, there was no evidence for correlational selection acting on the focal traits.

TABLE 7

Reproductive quadratic selection gradients and correlational selection acting on handling aggression (HA), breath rate index (BRI) and exploration score (ES) in urban and forest great tits, estimated using model averaging for the best equally supported models presented in Table S9. Equally supported models (ΔQAICc < 2) are represented in bold.

Sex	Habitat	Sample size	Quadratic selection gradients (γ/2)			Correlational selection
			HA2	BRI2	ES2	HA × BRI	HA × ES	BRI × ES
Females	Forest	103	0.030 ± 0.055	−0.051 ± 0.038	0.015 ± 0.067	−0.098 ± 0.046	0.081 ± 0.045	−0.044 ± 0.046
	City	276		−0.056 ± 0.031	−0.015 ± 0.069
Males	Forest	93	0.018 ± 0.064	0.070 ± 0.031	−0.088 ± 0.060	−0.018 ± 0.048	0.033 ± 0.051	−0.006 ± 0.052
	City	247

Viability selection

Goodness‐of‐fit tests conformed to the Cormack–Jolly–Seber's assumptions for each trait‐specific dataset (ES dataset: χ 2 = 29.222, p = 0.084, ĉ = 1.461; BR: χ 2 = 28.681, p = 0.094, ĉ = 1.434; HA dataset: χ 2 = 24.7, p = 0.213, ĉ = 1.235). Deviances and AICc were corrected for the ĉ values in all sets of models (ĉ = χ 2/df).

The best supported models for the recapture probability are presented in Table S2: in each case, the model with the lowest AICc estimated p constant. For the three traits, we found two equally performing models (Table 4, ΔAICc < 2). The habitat × age interaction was present in each set of selected models, in accordance with results from the global model (Section 3.3).

For all three traits, the best model included a negative quadratic trait factor suggesting that all three traits were under stabilising selection (Figure 2a–c). When estimating viability selection on HA, the best model included HA as additive effect as well as habitat × HA2 interaction; however, the second best and equivalent model (Table 4, ΔAICc < 2) did not include this interaction while it retained a quadratic HA2 effect. Taken together these results suggest that high and low values of HA are counter‐selected in both habitats, yet possibly more so in the urban habitat (Figure 2a). For breath rate, BRI2 effect was in interaction with habitat in the lowest AICc model (Table 4) providing again some evidence for stronger stabilising selection in the city (Figure 2b). Regarding exploration score, the linear and quadratic dimensions of selection were present in the best model but not in interaction with habitat, revealing that in both habitats birds were experiencing similar stabilising selection, particularly counter‐selecting low exploration behaviour (Figure 2c). However, the quadratic term was not present in the second best selected model and Figure 2c illustrates the weak stabilising force of selection.

FIGURE 2

Survival rate of forest (green) and urban (black) great tits in relation to Handling Aggression (a), Breath Rate Index (b) and Exploration Score (c). Dashed lines represent yearling individuals and solid lines represent adult birds of 2 years or more. Individual values of survival are calculated using the following equation ϕ ind = (1/(1 + exp(−β 0 + β 1 × x ind + β 2 × x ind 2)) with x ind the individual value of a trait and β estimates provided by the CMR models containing the trait effect and the lower AICc (Table 4, β 1 = linear selection differential, and β 2 = quadratic selection differential).

DISCUSSION

In this study, we investigated differences in behavioural traits and a stress‐related trait, between great tits breeding in a forest versus an urban habitat, and explored whether the observed divergences were aligned with patterns of reproductive and viability selection ongoing in each habitat. The three traits studied were all repeatable yet not strongly correlated, suggesting they did not co‐evolve in a behavioural ‘syndrome’. Overall, we confirmed previous findings that urban individuals were on average more aggressive, had a higher breath rate when handled by humans, and were faster explorers than their forest counterparts (Figure 1; Charmantier et al., 2017). In addition, urban yearling great tits had lower annual survival probability compared to forest yearlings (43% vs. 62%), while adult birds had similar survival in both habitats (40% vs. 41%, Table 5). Selection analyses provided evidence for stabilising viability selection and both linear and quadratic reproductive selection acting on the three traits studied. Importantly in the context of divergent phenotypes, there was evidence for a divergence in selection between the forest and urban habitats in two cases only: (a) stabilising viability selection favouring intermediate levels of handling aggression and breath rate was stronger in the urban habitat (Figure 2a,b), and (b) higher exploratory behaviour was associated with higher reproductive success in forest but not urban males (Table 6). Overall, the patterns of selection found here are not aligned with a divergence in phenotype.

As predicted, we found phenotypic divergence between urban and forest habitats for the three studied traits. Higher aggressiveness and exploratory behaviour are in line with previous results and studies in urban vertebrate populations (Charmantier et al., 2017; Evans et al., 2010; Martin & Réale, 2008). Regarding breath rate, previously reported results were more contrasted. Our finding that urban individuals displayed higher breath rate is similar to some studies (Torné‐Noguera et al., 2014) but contrasts with others (no difference: Senar et al., 2017; lower: Abolins‐Abols et al., 2016), suggesting that the shift observed could be due to another factor than urbanisation. Shifts in behavioural and stress‐related traits in urban habitats are often regarded as evolutionary adaptations in response to novel urban environmental conditions (Lambert et al., 2021; Sepp et al., 2018), a hypothesis in line with findings of urban–rural genetic divergences associated with genes involved in aggressiveness and exploratory behaviours (in particular SERT and DRD4, Van Dongen et al., 2015; Riyahi et al., 2017; Mueller et al., 2020). However, Lambert et al. (2021) reported recently that only six urban studies convincingly demonstrated an evolutionary adaptation to urbanisation, and none of these six examples concerned animal behaviour. One of the key elements missing in most studies reviewed by Lambert and colleagues is the adaptive nature of the observed divergence in phenotype. Here we provide one of the first studies investigating the link between individual behaviour and fitness in an urban context to test for the adaptive nature of an urban‐specific shift in behaviour (see also Halfwerk et al., 2019).

Correlations of behavioural traits are often observed in animals and referred as ‘behavioural syndromes’ or even ‘pace‐of‐life syndromes’ (POLS) if life‐history traits are involved (Réale et al., 2010). POLS represent (co)variation of personality and life‐history strategies across individuals. At two ends of a POLS continuum, fast individuals display low survival, high reproductive outputs, precocious reproduction and bold behaviours, while slow individuals show high survival, low reproduction and shy personalities. In this study, urban birds displayed features of a faster pace‐of‐life: higher exploratory behaviour, higher aggressiveness, higher breath rate and lower survival at yearling stage. These signatures of a faster urban life contrast with previous findings on reproductive traits associated with a slower pace‐of‐life in the same population (i.e. smaller clutches, Caizergues et al., 2018), and with a meta‐analysis showing slower POLS in urban birds (Sepp et al., 2018, but note this meta‐analysis did not consider behavioural traits). Note also that no correlation was found across the three studied traits, which contrasts with previous findings (Carere & van Oers, 2004) and with the prediction that these traits co‐evolve.

Our selection analyses revealed that shifts in avian behavioural and stress‐related traits in the urban habitat were not aligned with patterns of ongoing directional or quadratic selection. Regarding behaviour, we detected stabilising selection acting on HA via survival and on ES via both survival and reproduction, counter‐selecting extreme aggressive and exploratory behaviours in both habitats. While viability stabilising selection acting on HA was stronger in the urban context (yet with a similar phenotypic optimum), reproductive linear selection acting on ES in urban males was relaxed compared to strong selection favouring faster forest explorers. Taken together, these results suggest that the phenotypic divergence between forest and city great tits in both HA and ES does not result from habitat‐specific patterns of ongoing selection, and that the higher values of these traits observed in urban habitats are not adaptive. However, it is possible that spatio‐temporal variation in natural selection within the city masks the complexity of selection on these characters. First, higher aggression or exploratory behaviour might be advantageous only in some parts of the city (for instance in highly urbanised areas where food is scarcer and competition potentially higher) but not in others. If there is such spatial variation in selection, estimating selection at a fine scale within the city might be necessary to apprehend the adaptive nature of high aggression and exploratory behaviour in urban birds. Second, recent work by Corsini et al. (2020) suggests that selection on great tit morphology might act mainly during the early stages of life. If strong early‐life selection also impacts behavioural traits, our analysis focused on adult behaviour could potentially miss an important stage of selection occurring before individuals start to reproduce.

Breath rate under constraint is linked to stress physiology and is used as a proxy of acute stress response, that is, the short‐term response to a stressful event, in contrast with chronic stress referring to a long‐term response (Dantzer et al., 2014). Higher breath rate shows higher short‐term stress response to a stressor. In opposition to what is expected in case of adaptation to urban habitats, we found that urban birds displayed higher breath rate in response to handling and thus a higher acute stress response than forest individuals, while high breath rate was associated with reduced reproductive success in females and reduced survival in both sexes (Figure 2b). Note, however, that we did not measure glucocorticoid secretion levels under stress, and therefore we cannot conclude with certitude that higher breath rate was associated with higher physiological stress response. In addition, we detected stabilising viability selection acting on BRI that might be stronger in the urban habitat as suggested by the habitat × BRI2 interaction present in one of the two best selected models (Table 4; Figure 2b). Surprisingly, however, when considering reproductive success, while we detected slight stabilising selection on females BRI, we on the contrary found disruptive selection in males favouring extreme values of BRI in both habitats. Taken together these results suggest that higher breath rates in urban great tits are maladaptive. Note that previous studies found that acute stress responses associated with urbanisation could differ depending on population, species and traits studied (Powell et al., 2013). In particular, it would be interesting to explore how birds handle chronic stressors present in urban habitats from a physiological point of view (Iglesias‐Carrasco et al., 2020), especially because inadequate long‐term physiological responses might have deleterious repercussions on the birds' quality of life and ultimately their fitness, and thus might be under stronger selection than acute‐stress response.

Interestingly, across the three traits we found one case of relaxed selection (ES in males) but two instances of stronger stabilising selection (HA and BRI) urban habitats. The later results contrast with recent meta‐analyses pointing towards relaxed selection following anthropogenic disturbances (Fugère & Hendry, 2018; Lahti et al., 2009) and increased morphological variation in urban tit populations (Thompson et al., 2021). Whether urbanisation results generally in stronger or weaker selection is still highly debated and insufficiently explored empirically. While novel disturbances are typically expected to lead to stronger selection in the city (Alberti et al., 2017), urban features such as reduced predation or access to supplementary food (Lahti et al., 2009) could result in relaxed selection (Branston et al., 2021). It is likely that urban‐related changes in selection will depend both on the trait studied as well as on the proxy of fitness used, especially since a trait can undergo opposite selection pressures across the life cycle (see e.g. BRI results with both disruptive and stabilising selection depending on the fitness parameter studied). Hence future research should aim at comprehensive work that will disentangle the different factors driving stronger or weaker selection pressures in the cities.

The result of the non‐adaptive nature (i.e. not aligned with ongoing selection) of behavioural and stress response shifts linked with urbanisation raises numerous questions regarding the origins of such shifts and the implications for populations. First, despite the fact that behavioural shifts in urban habitats do not confer better fitness outputs compared to forest habitats in this case study, they could be advantageous to colonise urban habitats or tolerate the new environmental conditions. Indeed, bolder great tits are known to disperse more and on longer distances than shyer ones (Dingemanse et al., 2003) and could thus be more prone to colonise urban habitats. Second, behavioural shifts can result from plasticity in response to the environmental conditions experienced by individuals. For example, some personality traits can emerge from endocrine stress physiology (Baugh et al., 2017) and experiments demonstrated that stress stimuli inflicted to individuals could induce a reduction of their neophobic behaviours in house sparrows (Gormally et al., 2018). Hence, some behavioural shifts can emerge from a habituation process, in contrast with other traits such as aggressive behaviour which are consistent across lifetime (Cavalli et al., 2018), and phenotypic flexibility could be responsible for phenotypic shifts between habitats without necessarily providing fitness benefits to individuals.

The plastic versus genetic origin of the behavioural and stress response differences between urban and rural environments remains largely debated (Minias et al., 2018; Riyahi et al., 2017). While we have shown here that the three focal traits are repeatable, repeatabilities of 38% to 52% leave an important margin for plastic responses in these traits. In particular, a plastic response in the focal traits could be initiated during early life stages. For instance, hormonal levels experienced in mothers and transmitted to eggs can affect offspring personality (Rokka et al., 2014). In addition, DNA methylation could also play a role in the behavioural shifts observed in urban habitats (Caizergues, Le Luyer, et al., 2022; Riyahi et al., 2015). In any case, determining the contribution of genetics and plasticity in behavioural and stress responses shifts observed in the urban habitat is an important future challenge, and would greatly benefit from experimental designs such as cross‐fostering or common garden experiments (Diamond et al., 2017).

This study is we believe, the first to formally test for a link between avian behavioural as well as stress‐related traits, and fitness components, in a context of urbanisation. Our results suggest that, in opposition with expectations formulated in the literature, urban great tits display a faster pace‐of‐life with more exploratory behaviour and enhanced aggressiveness that may not be adaptive. In addition, urban individuals display higher breath rates even if this phenotype is counter‐selected, revealing maladaptation. However, further investigations and complementary analyses will be necessary to conclude on general patterns of selection ongoing in urban compared to natural habitats. In particular, replication in other populations will be a crucial step to understand whether the results found here are generalisable. In addition, selection analyses as performed here are data hungry while molecular tools can offer great opportunities to detect molecular adaptation with smaller datasets. We thus encourage urban evolutionary biologists to combine both quantitative genetic and genomic approaches to fully understand evolutionary processes ongoing in urban habitats (Perrier et al., 2020).

AUTHORS' CONTRIBUTIONS

A.E.C. collected the data, carried out the statistical analyses and wrote the manuscript; A.G. and A.C. conceived and funded the study, collected the data and contributed to the writing of the manuscript; S.P. collected the data and R.C. participated in the survival data analyses. All authors gave approval for publication.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Data S1

Click here for additional data file.