Changes in thermal niche position and breadth of bird assemblages in Spain in relation to increasing temperatures

Animal communities around the world are responding to climate change by altering their taxonomic composition, mainly through an increase in the colonisation rate of warm‐dwelling species and the local extinction of cold‐dwelling ones. We assessed whether the taxonomic composition of bird assemblages in peninsular Spain has changed in accordance with the recent increase in temperature. We also evaluated the role of species thermal affinities and population dynamics on these changes.


| INTRODUC TI ON
Climate change is among the main drivers of global change.
Increasing temperatures, changing precipitation patterns, or the increased frequency of extreme meteorological events are having effects on the populations of many animal (e.g.Herrando et al., 2019;Pearce-Higgins et al., 2015) and plant species worldwide (Gottfried et al., 2012).A host of different species-specific responses have already been recorded, such as changes in population size and reproductive success (Both et al., 2006;Stephens et al., 2016), shifts in species distribution-poleward and toward higher elevations- (Chen et al., 2011) and phenological changes, such as earlier reproduction and migration (Cohen et al., 2018;Radchuk et al., 2019).Recently, increasing interest has led to assessing whether those species-specific responses translate into coherent community patterns, and signs of climate change on entire animal (Lindström et al., 2013;Oliver et al., 2017) and plant assemblages (Feeley et al., 2020;Gottfried et al., 2012) have been sought.
The search for signs of climate change within species assemblages has most often revolved around responses to the worldwide rising temperatures.It has mostly been found a pattern of thermophilisation in line with an increased colonisation rate of warmdwelling species and local extinction of cold-dwelling ones across taxa and ecosystems (e.g.temperate butterflies and birds : Devictor et al., 2012a, Lehikoinen et al., 2021;mountain vegetation: Gottfried et al., 2012;marine biota: Burrows et al., 2019.However, response intensity varies locally (Gaüzère et al., 2015) and regionally (Foden et al., 2019;Tayleur et al., 2016) according to community taxonomic composition, and the environmental context.
The development of indices that reflect species' thermal affinities has proven useful to detect the response of species and assemblages to rising temperatures.The Species Thermal Index (STI), averages the temperature and the species experience in their areas of distribution (Devictor et al., 2008;Lindström et al., 2013).Thus, the index can be thought of as the position of a species in this single dimension of its niche.Most often STIs are calculated averaging mean temperatures throughout the species' ranges, though it could be argued that indexes based on maximum and minimum temperatures may actually be more informative of a species' tolerance limits.This is because species may compensate for temperature changes around their optima through physiological and behavioural adaptations, with limited energetic costs and consequences in fitness.In contrast, energy demand and fitness decrease sharply beyond critical temperatures (du Plessis et al., 2012).Species' thermal tolerance limits are asymmetric and species' responses to temperature change at warm versus cold edges (Araújo et al., 2013;Massimino et al., 2015), driving different responses of assemblages to minimum or maximum temperatures.
There is also evidence that the geographical variation of ambient climate is related to thermal tolerance (Khaliq et al., 2014).The indexes are built under the implicit assumption that species cannot change their thermal tolerances over short study periods, which may not be appropriate for taxa with short generation times (e.g.bacteria).These nuanced views of thermal niches suggest using a set of indices based on minimum, maximum and range-wide temperature variation to provide a more complete description of species' thermal niches than that based on simple averages.
The thermal niche of a species assemblage can be summarised with the Community Thermal Index (CTI; note we use the word 'community' here following common use, although we adhere to the term 'assemblage' following Fauth et al. (1996)).This index is estimated by averaging the STIs of the species occurring in an assemblage.
Assemblages with high CTIs host, on average, more thermophilic species than those with low CTIs.A change in CTI, when calculated from occurrence data, reveals that the assemblage has experienced compositional changes beyond those due to changes in species relative abundances.Therefore, the assemblages' CTIs increases with time that has been recorded by previous studies match the distribution shifts that would be expected if the species distribution followed rising temperatures (Chen et al., 2011;Gaüzère et al., 2015;Lindström et al., 2013).
It is commonly recognised that drivers of global change other than temperature increases need to be considered when assessing shifts in species distribution and community composition (Gaüzère et al., 2020;Oliver & Morecroft, 2014).Among them, rapid humaninduced habitat transformations have a major influence on bird assemblages (Oliver et al., 2017).However, habitat has only been included in a few CTI studies, and when the conditioning role of habitat and land-use dynamics has been studied, a considerable geographical variation has been found (Clavero et al., 2011;Gaüzère et al., 2020).Indeed, climatic and habitat niches might be associated, which further complicates telling apart their influence on aggregate indicators.Several studies described that cold-climate species in the Western Palearctic favour closed habitat (Barnagaud et al., 2012;Seoane et al., 2023) and, consequently, bird assemblages in forested areas generally have lower CTIs (Clavero et al., 2011;Gaüzère et al., 2020;but see Oliver et al., 2017).This suggests that an increase in forest cover might lead to a decrease in CTI, in contrast to the effect of rising temperatures.
The Mediterranean Basin is a recognised hotspot of threatened biodiversity in the Palearctic (Myers et al., 2000).Within this area, the Iberian Peninsula stands out due of its climatic and habitat heterogeneity, which provides interacting ecological gradients that might condition the ability of the species to follow their thermal niches.This large region (~0.6 M km 2 ) features both Mediterranean and Atlantic climates, and thus environments range from predominantly warm and arid conditions in the southeast to generally colder biodiversity monitoring, bird atlas, climate change, ecological niche, habitat change, population dynamics, thermal index and more humid conditions in the north.The oceanic influence (or lack thereof) further intensifies variations in climate, with large thermal amplitudes toward the inner part of this region.Climate forecasts warn that the Iberian Peninsula is particularly exposed and vulnerable to climate change within the already responsive Mediterranean region (Diffenbaugh et al., 2007;Giorgi, 2006).Regarding change in habitats, total forest area and forest density have increased in the region as has been the case in most of Western Europe (see for instance Ameztegui et al., 2021;Barbier et al., 2010).
Our main aim here is to evaluate whether the composition of bird assemblages has recently changed in agreement with the temperature rise experienced in continental Spain over the last two decades, considering both the modulating effects of habitat change and different aspects of species' thermal niches.For that purpose, we used community thermal indices (CTIs) for bird assemblages recorded in 10 × 10 km UTM grid squares during two periods.We took advantage of two major citizen science-based projects joined by thousands of volunteers throughout the whole country whose goal was to describe species distribution.These were the second (Martí & Del Moral, 2003) and (the most recent) third (Molina et al., 2022) breeding bird atlases of Spain.We expected a generalised CTI increase in the more recent atlas as a result of the northward shifts of the species' distribution driven by rising temperatures.However, we also expected the antagonistic effects of increased forest cover and temperature, that might attenuate the change in CTIs (Barnagaud et al., 2013).
To deepen previous studies on thermophilisation of communities, we further explored different formulations of the thermal indices and identified the species and population dynamics most responsible for assemblage changes.We expected that indices based on minimum or maximum temperatures across ranges, rather than those based on averages, would be more sensitive to local temperature changes because they are closer to the species' tolerance limits, where adaptive thermoregulation is compromised (Boyles et al., 2011).We also anticipated an interplay between populations dynamics (extinction, colonisation) and species' thermal preferences, so that both the local colonisation of warm-dwelling species and the local extinction of cold-dwelling species are responsible for changes in CTIs (Gaget et al., 2018).Finally, we expected the indices would show an increase of thermal generalist species if the pattern of community composition homogenisation with climate warming and habitat change that has been found elsewhere is at play in our study system (Clavero & Brotons, 2010;Davey et al., 2012;Gaget et al., 2020;Gaüzère et al., 2015).

| Bird data
We obtained species distributions from the last two Spanish breeding bird atlases, developed during two periods approximately 15 years apart (1998-2002in Martí & Del Moral, 2003;and 2014and -2018and in Molina et al., 2022)).These were nationwide citizen science-based projects enhanced with additional information and professional sampling where needed to fill information gaps.Volunteers committed to different sampling schemes that differ in their degree of difficulty.Most volunteers joined non-quantitative surveys, and the data here will be analysed as such.The atlases include the occurrence of all species recorded during the main breeding period (spring) in a grid of approximately 5300 10 × 10 km UTM squares covering the whole country.We followed a two-step process for square selection.First, we excluded records from islands and North African territories (we thus only considered records from the Iberian Peninsula).Secondly, we selected 2082 grid squares previously reviewed by local experts and national coordinators and considered to be well-surveyed (Figure 1).Grid squares where less than 10 species were registered were excluded.

| Species selection
We considered all breeding species recorded in the surveys except for those subjected to reintroduction programmes, non-indigenous, vagrants, accidental or pelagic species (n = 256, see Appendix S1, Table S1.1 in Supporting Information).The spatial distribution of the species within the first group is strongly influenced by human activities.We also considered vagrants and pelagic species not very informative of the relationship between terrestrial temperatures and occurrence.

| Species thermal indices
The species' thermal indices were calculated using the species' global distributions (BirdLife International and Handbook of the Birds of the World, 2019) and monthly average temperatures for 1970-2000 (Worldclim 2.0: Fick & Hijmans, 2017), to describe long-term thermal conditions in the distribution ranges.We selected the breeding period distributions, excluding occurrence areas derived from introductions.We compared the data of one quarter of these global distributions with that of other more geographically restricted sources, to informally assess the accuracy of the database (i.e.Brazil, 2009;eBird, 2021;Grimmett et al., 2016;Ollé & Trabalon, 2019;Shirihai & Svensson, 2018a, 2018b;Svensson et al., 2009).
We considered temperature data from April to July-the main breeding season for most species in our study area-at a spatial resolution of 30 arcseconds (~1 km 2 ), and extracted monthly mean temperatures, maximum temperatures during the warmest month (July), minimum temperatures during the coldest month (April), and the difference between the latter two (temperature range).
We obtained four thermal indices for each species (Species Thermal Index-STI) by combining the global species' distribution and the climate information.The STI1 shows the mean temperature of the breeding season (April-July) throughout the species' breeding distribution range.The STI2 is the average of the maximum temperatures above the percentile 95 in July and the STI3 is the average minimum temperature below the percentile 05 in April.These three indices provide a description of thermal affinities that consider niche position (STI1), and proxies for critical thermal maximum (STI2) and minimum (STI3).The fourth index (Species Thermal Range-STR) represents the average thermal range (April-July) throughout the distribution area and can be understood as species thermal breadth.
To our view, this multiple perspective on thermal niches is a useful but neglected avenue for studies on the thermophilisation of assemblages.

| Community thermal indices
We calculated a set of CTI for the assemblage of bird species in each of the 10 × 10 km UTM grid squares of each of the breeding bird atlases.We obtained four different CTIs: CTI1, CTI2, CTI3 and CTR.The first three were calculated as the average of the STI1, STI2 and STI3 of the species present in the assemblage, respectively.CTI1 is the most widely considered index in thermophilisation studies.However, it reflects the thermal optimum of the species, so it might be less sensitive to temperature changes due to the species' potential for phenotypic plasticity and microevolutionary adaptations.CTI2 and CTI3 are based on proxies for critical thermal limits and might be more responsive to filtering effects of increasing temperatures.A high CTI2 indicates a predominance of species tolerant of high temperature, while high CTI3 indicates species sensitive to low temperatures are common within the assemblage.Comparing the temporal change of the indexes and their responses to the temperature gradient allows to assess whether temperature have a larger effect on cold-tolerant species (CTI3) or on warm-tolerant species (CTI2).The CTR (community thermal range) is based on the average temperature range of the species (STR) that make up the assemblage, and thus informs on the average niche breadth (Gaget et al., 2020).An increase in CTR would suggest a trend favouring thermal generalist species.We calculated CTIs for each of the four-year periods covered by the atlases, which could potentially absorb the lag in community response to temperature rise that has been reported in other studies based on yearly CTIs (Devictor et al., 2012a).
Finally, we considered the potential influence of local temperatures and habitat on CTI changes.Local temperatures were obtained from Chelsa (v.2.1.,Karger et al., 2017), averaging data for each fiveyear sampling period in each square.We used the CORINE Land Cover Accounting Layers built for years 2000 and 2018, to link forest cover with the community indices for the first and second sampling periods, respectively (EEA, 2019; Table 1).

| Statistical analyses
We built mixed-effects models with spatial correlations to address the change in the CTI between the study periods.The response variable in each model is thus the index calculated per square in each of the sampling periods.The model's fixed component included (in this order): geographic location (latitude and longitude, included as orthogonal cubic polynomials to allow for large non-linear geographical trends), forest landcover and temperature of the square as continuous variables, and the sampling period as a two-level factor (Table 1).The continuous explanatory variables were standardised.We included the interaction between forest cover and period, and the interaction between temperature and period, to respectively assess whether the relationship of forest cover and temperature with the community index varied between sampling periods.The model included We followed this procedure using each of the thermal indices as response variables.Our aim was to evaluate whether the way we characterised the average thermal affinity of the species in a community (CTI1, CTI2, and CTI3) determines the community response to climate change that we can detect, and whether, as we expected, bird communities have generally gained species with wider thermal ranges and lost species with narrower thermal ranges (with CTR).
We used a variance partitioning approach based on linear models as an approximation to the (marginal) relative contribution of the predictors on the thermal indices to show their joint and independent effects.We considered three groups of explanatory variables, with (i) temperature, (ii) forest cover to account for habitat and (iii) longitude and latitude to explain non modelled geographical variation.For each period, we built three models including just one of the three groups of predictors, three models including two of them, and a saturated model which included all three of them.Using set theory, we assigned the amount of variation attributable to each predictor separately; the joint contribution of two of them (i + ii; i + iii; ii + iii), and the effects of the three predictors which are inseparable from each other (i + ii + iii).The final amount of variation of each category results from averaging the scores obtained for both study periods.
We repeated this procedure with each of the thermal indices.
We assessed the contribution of each species to either increases or decreases in community thermal optimum, and the extent to which these specific contributions are attributable to distinct population dynamics (colonisation, extinction, or persistence).To do so we used a simple jackknife procedure (as in Princé & Zuckerberg, 2015;Tayleur et al., 2016; but see Gaüzère et al., 2019 for a solution generalisable to more complex analytical scenarios).Firstly, we calculated the CTI for each grid square and period after removing one species at a time.Secondly, these jackknifed thermal indices for each grid square and period were subtracted from the whole community index in that sampling unit and period.Thus, for example, a positive result identified a warmdwelling species that contributed more than the assemblage's average species to reach a high community thermal index.By subtracting the value obtained for each grid in the second period from that of the first period, we obtained an index of the contribution of a species to the change in the CTI of a grid.Lastly, to associate these contributions to population dynamics, we multiplied the matrix of subtractions by three dichotomic matrices showing which square was colonised by the species (if recorded in the second sampling period but not in the first), which one experienced a local extinction (if recorded in the first sampling period but not in the second), and which one suffered no changes in occupation (if either presence or absence was recorded in both periods).The sum of the values for each species in these matrices allowed us to sort species in terms of their relative contribution to the community thermal optimum due to recent colonisations, extinctions or persistence.Species in the first and fourth quartiles of species thermal index distribution were classified as warm-or cold-dwelling species, respectively, while the rest were classified as intermediate (see a detailed description of this analysis in Appendix S2 in Supporting Information).Then, we modelled the relative specific effect on the population process (colonisation vs. extinction), the classified species thermal index (warm-dwelling, cold-dwelling and intermediate species) and their interaction, to provide a descriptive summary of the species' effects.
The spatial analyses (obtaining the species' thermal indices and model predictors) were carried out using qgis 3.10.Subsequent statistical analyses were performed in R v4.1.2.(R Core Team, 2021) TA B L E 1 Explanatory variables for models of community thermal indices in bird assemblages in peninsular Spain.These were obtained for a grid of 10 × 10 km UTM squares overlayed on the study area.

Temperature
Standardised mean annual temperature of each square in each sampling period obtained from Chelsa v.2.1 at a horizontal resolution of 30 arc sec (Karger et al., 2017).This dataset is based on downscaled air temperature two meters above the ground modelised from the data collected from many sources (mainly weather stations, weather balloons, aircraft, ships, and satellites).Mean annual temperatures were averaged for 1998-2002 and 2014-2018, and assigned to the first and the second sampling periods, respectively (in degrees Celsius; note raw Chelsa data is in Kelvin × 10 degrees)

Period
Each of the sampling periods considered (1998-2002; 2014-2018) Forest cover*period Interaction between forest cover and sampling period Temperature*period Interaction between temperature and sampling period a The CORINE Accounting Layers constitute an effort to control the methodological differences between the CORINE versions.
and specialised packages ('spaMM' v. 4.1.20: Rousset & Ferdy, 2014, to build the spatial mixed models).Thermal indexes are robust to particularities of input datasets (Devictor et al., 2012b;Lindström et al., 2013), which in combination with a qualified bird dataset (see above) leads us to believe there is no room for systematic biases.Accordingly, the community thermal optimum (CTI1) increased slightly (on average 0.02 ± 0.0073°C).This index increased in parallel with local temperatures (β temperature = 0.65 ± 0.020) and showed a higher gradient in the first sampling period (β temperature*perio d = −0 .07± 0.025), which suggests that cooler grid squares experienced a greater rise in the index (Figure 2).The index was lower in the grid squares with larger forested areas (β forest = −0.12± 0.010), and this relationship did not change significantly between the study periods.CTI values increases between sampling periods were conspicuous for grid squares with below-average temperature at any given forest cover (see Appendix S3, Figure S3.1 in Supporting Information).
Lastly, forest cover was negatively associated to all thermal indexes, which implies that forest habitats tend to be occupied by colddwelling species.

| Community thermal ranges
The range of species thermal tolerance within an assemblage decreased with temperature and to a lesser extent with forest cover (Table 3, and see Figure S3.4 in Supporting Information).The separate and joint effects of temperature and geographic location explained most of the CTR variability, whereas habitat was less important (Figure 3).We note that thermophilic species also tend to have narrower thermal breadths in our dataset (Appendix S1, Figure S1.1).However, our model was only able to explain a small fraction of the variability of the index.

| Species contribution
The species that have contributed the most to increase the commu-  4, Figure 4).
On the other hand, cold-dwelling species contributed to a lesser extent to the generalised increase of CTIs mainly through local extinctions-thus disappearing from a large part of their former distribution (5.3°C [3.7-6.9],Table 4, Figure 4).

| DISCUSS ION
Our results show that bird assemblages in peninsular Spain have undergone compositional changes during the breeding season, associated with the increasing temperatures during the last two decades.
This reorganisation was due to an increase in the prevalence of thermophilic species, suggesting a larger colonisation of warm-dwelling species and, to a lesser extent, a more frequent local extinction of species with a greater affinity for cold climates.However, the assemblages' response intensity differs depending on how the thermal affinity of the species is summarised.
The thermal optima of the bird assemblages increased between sampling periods, with an annual increase around 0.0012°C year −1 (considering the 16.5 years between the centre of both sampling periods).This result is smaller than, but in line with studies based on species occurrence in other regions of Europe (0.0044°C year −1 in France (Devictor et al., 2008); 0.005°C year −1 in Sweden (Tayleur et al., 2016)) and similar to an increase of 0.0028°C found in a region -Catalonia-within our study area in an abundance based study (Devictor et al., 2012a).In summary, bird assemblages are acquiring a large proportion of thermophilic species in Western Europe and most likely elsewhere (Devictor et al., 2012a;Lehikoinen et al., 2021).However, the way communities respond to increasing temperatures is not consistent throughout regions (Devictor et al., 2012a;Lehikoinen et al., 2021;Tayleur et al., 2016) Note: Models included the geographic locations (as third-degree orthogonal polynomials of longitude and latitude), forest cover (according to CORINE accounting layers), temperature (according to CHELSA) and sampling period (1998-2002 vs. 2014-2018).Indices are CTI1 for thermal optimum, CTI2 for thermal maximum and CTI3 for thermal minimum.p-values are provided after log-likelihood ratio tests between nested models, respecting marginality (note we provide values for whole terms for polynomials and factors).λ: variance attributable to spatial random effects for each sampling period.φ: residual variance.An estimate of r 2 is provided as the r 2 of a linear regression of actual observations on marginal model predictions.
niche position based on thermal optima are more common, here we explored those based on thermal maxima and minima, all of which seem to be positively related to temperature.Thus, they reflect that communities found in warmer environments harbour, on average, more thermophilic species than those in colder environments.
However, the index based on species thermal minima across species ranges of distribution (CTI3), which describes species tolerances to low temperatures during the breeding season, is the most responsive to temperature variations.In this regard, it might be considered the best indicator of climate change perturbations in communities among those studied here.
The greater link between the thermal minima community index and local temperatures might be due to species' thermal minima (STI3) regulating animal distributions more than thermal maxima (STI2) (Araújo et al., 2013).This is consistent with previous evidence that the cold-edge distribution limit of birds is more determined by physiological temperature constraints (Massimino et al., 2015) and with the fact that realised niches-on which STIs are based-are more informative of species' critical minimum temperature than critical maximum across animal taxa (although the evidence is weaker for birds: Araújo et al., 2013).Ecological studies have long speculated that a suite of biotic processes (such as competition or facilitation) and abiotic factors (e.g.precipitation, water availability) are of greater relative importance in defining the distribution at the warm-edge and thus the temperature maxima at which species are found (Hawkins et al., 2003;MacArthur, 1972).In line with this, precipitation patterns seem to influence bird population responses to climate change more the closer they are to the equator (Pearce-Higgins et al., 2015).However, temperature changes have also been shown to result in distribution shifts-directly or indirectly (by altering biotic interactions)-at the species warm edge (Cahill et al., 2014;Princé & Zuckerberg, 2015), which further complicates a mechanistic interpretation of changes in assemblage species composition.
In contrast to the indexes of thermal position, average thermal breadth decreases with local temperature.This suggests that thermal niches are more constrained in the warmer localities of our study area; in other words, the colder environments of our study area might not limit the occurrence of thermophilic species as much as the warmer environments limit the occurrence of cold-dwelling species.However, we were not able to detect the pattern of community homogenisation between sampling period we hypothesised, with assemblages gaining more generalist species due to climate warming or habitat change (Davey et al., 2012;Gaüzère et al., 2015).
Indeed, we found that warm-dwelling species, with increasing populations and expanding ranges, tend to have narrower thermal niches (Figure S1.1; Gaget et al., 2020).It has been suggested that thermal specialists have higher upper critical temperature limits than thermal generalists (Boyles et al., 2011), and to our view, this may explain their success in a scenario of increasing ambient temperature.
Our results show a role of habitat in the average thermal niche of assemblages, which points to antagonistic interacting effects of climate warming and land-use change (Gaget et al., 2020).The areas with higher forest cover tend to harbour fewer thermophilic species than those with lower forest cover.Although the effect is limited, it is consistent with previous findings arguing that community composition integrates both habitat and climatic species' niches, in a non-hierarchically way (Barnagaud et al., 2013).
Forest assemblages may tend to have low CTIs because forests have microclimate conditions that moderate temperature variations and favour cold-dwelling species, or because they signal the biogeographical history of species (Barnagaud et al., 2013;Zellweger et al., 2020).However, it should be noted that forests and woodlands are more widespread in the cooler and wetter North and North-western corners of our study area, in contrast to the warmer and dryer South and East, which complicates separating the pure effects of habitat and temperature.Our analysis of species' contribution suggests that habitat changes might have had some effect on recent changes in community thermal optima (see Table S2.3, Figure S2.3).Among the species that have con- Other land-use changes may have also contributed to the increase of community thermal optima in our study area.This is the case of farmland and steppe environments, which have undergone an intense transformation in recent decades in the Iberian Peninsula, as a result of rural abandonment and agricultural intensification (Guerrero et al., 2010;Traba & Morales, 2019).Although an increase in bird species diversity has been detected in some farmland environments, this is due to the colonisation of non-farmland-dependent species, and masks a sharp decline in specialised farmlanddependent species (Díaz et al., 2022;Filippi-Codaccioni et al., 2010).
Species from this latter group are generally warm-dwelling species (Clavero et al., 2011), so their decline might also decrease thermal optima or lessen their increase (Oliver et al., 2017).Lastly, we underline that our analysis is based on species occurrence data gathered at a coarse spatial resolution.We speculate that we might have found a stronger response of communities to increasing temperatures if we had had species abundance data at a finer resolution because abundance is more sensitive to impacts than occurrence (Waldock et al., 2022; but see Princé & Zuckerberg, 2015).
However, obtaining high-quality data at high resolution using large-scale citizen-based projects may not be feasible, in particular in sparsely populated, remote or difficult-to-access areas (Galván et al., 2021;Loss et al., 2015).Indeed, choosing sampling requirements wisely facilitates public involvement in large-scale research (Brown & Williams, 2019;Carrascal & Del Moral, 2021).Here, we show that a simple, yet structured citizen-science project has provided useful information and a warning regarding changes in species distributions.
The current study warns about the rapid and directional compositional changes that bird assemblages have experienced in a narrow lapse of time (little more than a decade) despite the counteracting effects of increasing local temperatures and habitat transformations: warm-dwelling species have generally gained a higher prevalence to the detriment of cold-dwelling ones.Such impacts of climate change highlight the importance of monitoring programmes with the widest possible spatial coverage, where citizen-science projects are most useful to cope with the huge sampling effort needed.It should also be noted that, despite often overlooked, our

F
I G U R E 1 Location of the study area and sampling units (open circles) considered for the analyses of community thermal indices of breeding bird assemblages in peninsular Spain (n = 2082 10 × 10 km UTM grid squares, Geographic Coordinate System: WGS 84).Circles are overlaying (a) long-term average summer temperature (April to July) for the study area during 1970-2000, and (b) mean temperature increase from 1998-2002 to 2014-2018 sampling periods.(c) shows the thermal niche optimum (CTI1) increase between both sampling periods at the sampled grid squares.correlated random effects by fitting a Matérn correlation function to the geographic centres of each square.Matérn functions are useful to describe latent spatial terms, and we fitted a distinct function per period to account for possible differences in the spatial autocorrelation structure between periods (preliminary models built with a common Matérn function for both sampling periods showed weak but significant spatial autocorrelation).
Figure3).The forest cover contributed significantly to explain the variation, but its overall influence was relatively minor.
nity thermal optimum index (CTI1) are warm-dwelling species that extended their breeding ranges and colonised new squares after the first sampling period.Examples are the Eurasian Collared Dove (Streptopelia decaocto), Zitting Cisticola (Cisticola juncidis), Red-rumped Swallow (Cecropis daurica) and Pallid Swift (Apus pallidus) (see a comprehensive description of the contribution of each species in Appendix S2).Overall, the colonisation of warm-dwelling species had the largest impacts on the index (9.4°C[95% CI: 7.8-11.1]on average, Table Examples of this group of species are the Common Sandpiper (Actitis hypoleucos), Goshawk (Accipiter gentilis), Tree Sparrow (Passer montanus) and Skylark (Alauda arvensis).Nevertheless, we have also identified some species that experienced changes contrary to our expectations based on temperature F I G U R E 3 Venn diagrams showing variation partitioning of joint and shared effects of temperature, geographic trend, and forest cover on community thermal indices of bird assemblages in peninsular Spain.Indices are CTI1 for thermal optimum, CTI2 for thermal maximum, CTI3 for thermal minimum and CTR for community thermal range.increase.These are warm-dwelling species that suffered frequent local extinctions (−6.9°C [−8.5 to −5.3] on average) and cold-dwelling species that colonised new grid squares (−6.4°C [−8.0 to −4.8] on average).Examples of the former are the Barn Owl (Tyto alba), Eurasian Thick-knee (Burhinus oedicnemus), Western Black-eared Wheatear (Oenanthe hispanica) and Rufous-tailed Scrub-robin (Cercotrichas galactotes).Examples of the latter are the Lesser Spotted Woodpecker (Dryobates minor), Eurasian Nuthatch (Sitta europaea), Golden Eagle(Aquila chrysaetos) and Songthrush (Turdus philomelos).These are fairly common species with different migratory status (from resident to long-distance migrants).Correlates of species' traits to speciesspecific contributions to thermal indexes were not obvious and we do not explore them here.
tributed the most to decreasing the thermal position indices in our study area are some enlightening examples of cold-dwelling forest species such as the Lesser Spotted Woodpecker (Dryobates minor), Nuthatch (Sitta europaea), Song Thrush (Turdus philomelos) and Grosbeak (Coccothraustes coccothraustes), whose effect was mainly through colonisation of new grid squares.Forest area and density have grown in Spain and other European countries in recent decades (Korhonen & Stahl, 2020; Ministry of Ecological Transition, 2018), and populations of forest bird species have increased accordingly (Escandell & Escudero, 2021; Seoane &Carrascal, 2008).These forest species had a notable effect on the change in community thermal optima mainly through colonisation of new grid squares.They had an effect opposite to the generalised thermal index increase.However, even though we only consider forest area change as a proxy for forest cover variation, other changes such as those related to forest structure, tree species composition and patch connectivity are more difficult to quantify but can lead to the colonisation of forest species and to the local extinction of species linked to open habitats(Martínez- Jauregui et al., 2016).
Coefficient table for the spatial (Gaussian) linear mixed models explaining community thermal index values in bird assemblages in peninsular Spain, in 2082 10 × 10 km UTM grid squares.
, which provides considerable spatial variation on any account of climate-driven changes in assemblage descriptors.Community thermal indices are collective metrics that describe the species average thermal niche position in assemblages, and these might conceivably be based on different perspectives on temperature (Arnan et al., 2015; Gaget et al., 2020).While descriptors of TA B L E 2 Coefficient table for the spatial (Gaussian) linear mixed model explaining community thermal ranges (CTR, average thermal breadth) in bird assemblages in peninsular Spain, in 2082 10 × 10 UTM grid squares.Coefficient table for the linear model summarising the relative species-specific effects on community thermal index CTI1 in bird assemblages in peninsular Spain, in 2082 10 × 10 km UTM grid squares.Effects were modelled on population processes (colonisation vs. extinction of grid squares between 1998-2002 and 2014-2018 sampling periods), species thermal indices (warm-dwelling, cold-dwelling and intermediate species), and their interaction.