New Evidence on the Linkage of Population Trends and Species Traits to Long-Term Niche Changes

Simple Summary: Temporal changes of ecological niches in relation to both population trends and species traits (life-history and ecological characteristics) are a poorly investigated topic in animal ecology even though they can provide crucial information for investigating how species respond to ongoing environmental changes. To ﬁll this gap, we assessed niche changes from 1992 to 2017 for 71 common breeding birds in Northern Italy, using several niche metrics and relating them with both population trends and species traits. Our results showed that 32% of the species tended to change their niche in relation to the accessible environmental conditions. We found a tendency to adjust the niche towards warmer thermal conditions. This process may depend on several drivers, which act differently across species. Species undergoing a niche expansion showed increasing populations, while species tending to retain their niche showed a tendency to be associated with decreasing populations. We also found evidence for a non-random association between niche temporal changes and species traits. Investigating these relationships is a pivotal point to implement effective conservation strategies for wildlife. Abstract: Despite the assessment of long-term niche dynamics could provide crucial information for investigating species responses to environmental changes, it is a poorly investigated topic in ecology. Here, we present a case study of multi-species niche analysis for 71 common breeding birds in Northern Italy, exploring long-term niche changes from 1992 to 2017 and their relationship with both population trends and species traits. We (i) quantiﬁed the realized Grinnellian niche in the environmental space, (ii) compared variations in niche breadth and centroid, (iii) tested niche divergence and conservatism through equivalency and similarity tests, (iv) calculated niche temporal overlap, expansion and unﬁlling indices, and (v) investigated their association with both population changes and species traits. Results supported niche divergence (equivalency test) for 32% of species, although two-thirds were not supported by the similarity test. We detected a general tendency to adjust the niche centroids towards warmer thermal conditions. Increasing populations were positively correlated with niche expansion, while negatively correlated with niche overlap, albeit at the limit of the signiﬁcance threshold. We found moderate evidence for a non-random association between niche changes and species traits, especially for body size, clutch size, number of broods per year, inhabited landscape type, and migration strategy. We encourage studies correlating long-term population trends and niche changes with species traits’ information and a speciﬁc focus on cause-effect relationship at both the single and multiple-species level.


Introduction
The quantification of the realized niche of a species [1] (sensu Grinnell [2][3][4][5], hereafter niche) and how it varies across space and time has received increasing attention in ecological Birds 2022, 3 150 studies, with important implications in evolution, biogeography, and conservation [6][7][8][9][10][11][12][13][14][15]. Many studies have focused on niche evolution in related taxa [16][17][18][19][20][21][22][23], or on niche quantification of alien species to assess their potential as invaders in non-native ranges [24][25][26][27][28][29][30][31][32]. However, niche quantification and change in avian species remain poorly investigated aspects, especially if compared to other research topics (e.g., population trend assessment). As well as the general tendency in ecology, in ornithology most studies on niche also have focused on evolutionary aspects or on invasive species. Moreover, several studies have focused on niche-tracking under climate change by a correlative modeling approach (e.g., species distribution models [33]) in order to predict future changes in species distributional range [34][35][36][37][38]. Both the niche-tracking concept and its application in predictive species distribution modeling assume niche conservatism [5], which implies that in the future a species will occupy the same environmental conditions that it occupies today. Several studies have supported niche conservatism on both ecological and evolutionary timescales [6,39,40]. However, these studies focused on conservation of the fundamental niche, while a separate issue is whether the realized niche is also conserved. Some studies have highlighted that species may not retain their realized niche between the native and invaded range [41,42] and that variations may be expected under climate changes [43].
Very little attention has been given to within-species niche changes in the same geographic area over time. Assessing the change of the niche of a species in a relatively long-term period (e.g., decades) could allow measuring the environmental plasticity of the species in adapting to suboptimal environmental conditions and in colonizing new habitats, as well as the tendency in moving the niche boundaries or position towards different combinations of environmental gradients. This information may be combined with population trends to shed light on mechanisms that regulate species capability in responding to environmental changes, which can be pivotal to plan effective conservation actions. We ascertained a lack of studies on this potential relationship (but see [44,45]), although long time-series data coming from standardized bird monitoring programs (e.g., North American Breeding Bird Surveys [46]) exist. Ralston et al. [44] found a positive association between changes in climatic niche breadth and population trends. Ralston et al. [45] also showed that increasing species tended to show greater levels of climatic niche expansion compared to declining species, while declining species had significantly greater climatic niche unfilling compared to increasing species. Notwithstanding the relevance and the novelty of their results, these two studies did not consider the habitat characteristics (e.g., land cover variables) in niche quantification, thus omitting important scenopoetic variables-i.e., variables that are invariant to species [47]-constituting the realized Grinnellian niche. Another important question is how temporal changes in niche occupancy (i.e., variations of the density of occurrence corrected by the environmental availability in the multidimensional niche space) are related to population trends. The inference on the cause-effect relationship between niche changes and population trends needs to be carefully evaluated [45]. Indeed, it is reasonable to expect that a change (both an expansion and an unfilling) in a species niche could be the reason for the observed population trend (e.g., [48] for the isotopic niche). However, it is also reasonable that the variation in a species niche may be the result of a change in the species population abundance (due to niche-independent potential factors, e.g., epidemic, interspecific competition, density-dependent factors), which could have led to a niche reshaping with expansion, unfilling, centroid shift, variations in breadth or density of occurrence.
Beyond the relevance of a species-specific approach, in the past decade the increasing attention on a trait-based approach allowed highlighting emergent ecological patterns in groups of birds sharing similar life-history and ecological traits [49][50][51][52][53]. The relationship between the ecological niche and species traits has been mainly explored under an evolutionary perspective ( [54][55][56]), but there is a lack of studies focusing on ecological aspects (but see [57,58]). Revealing such types of ecological signal would point out a non-random evolution of niche changes, suggesting the existence of shared ecological pressures within groups of species sharing similar characteristics.
Birds 2022, 3 152 by a restricted pool of expert surveyors, thus limiting the presence of inter-observer detection bias. For the purpose of this study, we retained the most common species [53] with a frequency of occurrence higher than or equal to 1% in both T1 and T2 (number of species = 71).

Figure 1.
Main land covers in the study area and locations of the sampling sites. Land covers are represented according to 3rd-level of DUSAF land-cover map [59] with some mergers (see Table  S1). 10-km square grid (WGS84, UTM 32N) is superimposed on sampling sites. X and Y coordinates are in km.

Habitat and Climatic Variables
Species niche was quantified by using 17 land-cover variables, four topographic variables and four climatic variables. Land-cover variables were recorded on the field as fractional cover within a 250-m circular buffer around each point count and were classified according to the 3rd-level of DUSAF digital map [59]. Some variables were merged because they were strongly underrepresented in the dataset and indicated similar types of land covers (Table S1). Land-cover variables included continuous urban matrix (e.g., dense urban areas, industrial areas, infrastructures; C110), discontinuous urban matrix (e.g., farmsteads; C112), arable lands (e.g., maize, wheat, horticulture; C211) paddy fields (C213), permanent orchard plantations (olive groves, vineyards, orchards; C220), wood plantations (C224), meadows and pastures (C231), broadleaved forests (C311), coniferous forests (C312), mixed forests (C313), grasslands (C321), shrublands (C322), no tree vege- Figure 1. Main land covers in the study area and locations of the sampling sites. Land covers are represented according to 3rd-level of DUSAF land-cover map [59] with some mergers (see Table S1). 10-km square grid (WGS84, UTM 32N) is superimposed on sampling sites. X and Y coordinates are in km.

Habitat and Climatic Variables
Species niche was quantified by using 17 land-cover variables, four topographic variables and four climatic variables. Land-cover variables were recorded on the field as fractional cover within a 250-m circular buffer around each point count and were classified according to the 3rd-level of DUSAF digital map [59]. Some variables were merged because they were strongly underrepresented in the dataset and indicated similar types of land covers (Table S1). Land-cover variables included continuous urban matrix (e.g., dense urban areas, industrial areas, infrastructures; C110), discontinuous urban matrix (e.g., farmsteads; C112), arable lands (e.g., maize, wheat, horticulture; C211) paddy fields (C213), permanent orchard plantations (olive groves, vineyards, orchards; C220), wood plantations (C224), meadows and pastures (C231), broadleaved forests (C311), coniferous forests (C312), mixed forests (C313), grasslands (C321), shrublands (C322), no tree vegetation under evolution (C324), areas with sparse or absent vegetation (e.g., outcrops, debris deposits, glaciers; C330), wetland vegetation (C410), rivers and streams (C511), and natural and artificial lakes (C512). Topographic variables were derived from the Digital Elevation Model (DEM, 20-m ground-resolution, downloadable from http://www.geoportale.regione.lombardia. it/, (accessed on 21 September 2021)), and included sine (sin) and cosine (cos) of the aspect (adimensional), slope (grades) and elevation (m), all measured by averaging pixel values within a 250-m circular buffer around each survey site. Climatic variables were obtained from the European Climate Assessment & Dataset (ECAD) with a resolution of 0.1 degree [68]. They included the annual cumulative rainfall (rr), and the annual average (tavg), minimum (tmin) and maximum (tmax) temperature. For each period, after trimming the original grid by our study area, we spatially associated the value of each climatic variable to the corresponding point count, by averaging values within a period between the two years immediately preceding the first year of the time interval and the last year-i.e., 1990-1996 for T 1 and 2013-2017 for T 2 . This way, we obtained values representative of climatic conditions in the two periods, limiting the effects of climatic oscillation that might be detected during a single year and taking into account that birds could have exhibited a delayed response to climatic parameters [69].

Niche Quantification and Temporal Changes
In recent years, the quantification and the analysis of spatio-temporal niche changes benefited from novel methodological insights [72][73][74][75][76][77]. The most important progress concerns the use of the environmental space (E-space [72][73][74]76]) instead of the geographical space (G-space [33,78,79]). This way, no assumptions are required for the model-based approach [72,74]. The use of the E-space permits better assessment of the niche overlap because it takes into account the environmental availability and the analogy between ranges [72,80], and it overcomes some limitations on the inference of niche similarity [76].
Niche analysis was performed separately for the 71 bird species using the R package ecospat [81]. The E-space was assessed through an ordination technique using a PCAenvironment framework, retaining the two main axes resulting from the input of the 25 habitat and climatic variables presented in Section 2.2. The axes of the PCA maximize the ecological variance present in the study area for both periods. The PCA scores of species occurrence in each period were projected onto a grid of cells bounded by the minimum and maximum PCA scores of both periods-i.e., gridded E-space-and a kernel density function was used to estimate the density of occurrence of the species in each cell of the grid [74].
We evaluated the species niche breadth and centroid along each of the PCA axes for each period, by computing the variance and the median of the scores along the principal components, respectively. To detect a significant variation in the niche breadth and centroid between the two periods across all species, we performed the Wilcoxon rank sum test. To compare the overlap of the environmental niche of each species between T 1 and T 2 , we used the metric Schoener's D [72,82]. This index quantifies the amount of niche overlap by computing, for every grid cell in the E-space, the absolute differences in densities of occurrence-corrected by the prevalence of the environments in their range, i.e., niche occupancy-between the two periods. It varies from zero (no overlap) to one (complete overlap). Additionally, we calculated three niche dynamic indices between the two periods, namely the expansion, unfilling, and stability index (all varying from a minimum of zero to a maximum of one). The expansion index represents the E-space occupied by the species in T 2 only, the unfilling index indicates the E-space occupied by the species in T 1 only, and the stability index represents the E-space where the species occurred in both periods (T 1 ∩ T2) [8]. This decomposition can provide more information about the drivers of niche change between the two periods. We computed the three indices by limiting their calculation to the shared available E-space (intersection = 0), in order to exclude that the observed values were due to differences in the accessible E-space [76]. Since the stability index is complementary to the expansion index, it was not reported in the results in order to avoid redundant information. We evaluated the correlation among Schoener's D, and expansion and unfilling indices by the non-parametric Kendall test [83]. To test the statistical significance of niche overlap (Schoener's D), we ran equivalency and similarity tests, comparing the overlap of the two observed niches to the overlap of simulated niches. The equivalency test assesses whether two niches are more or less equivalent than expected by chance when randomly reallocating all occurrences between the two ranges. The similarity test assesses whether the overlap between observed niches in two ranges is different from the overlap between niches when their occurrence density grids are randomly shifted in the background environment [72,74,78]. To ensure a robust statistical inference, each test was performed based on 1000 simulations [84,85]. For both tests, we assessed the hypotheses of both niche divergence and niche conservatism-i.e., one-side tests-and for the similarity test both niches were randomly shifted. For both the equivalency and similarity tests, a significant p-value (≤0.05) means that the two niches are less (niche divergence) or more (niche conservatism) equivalent/similar than a random expectation.

Relationship between Niche Temporal Changes, Population Trends and Species Traits
To assess the relationship between niche temporal changes and population trends we performed a Kendall non-parametric correlation analysis [83] between each of the niche metrics (Schoener's D, expansion and unfilling index) and the species population trends (% change). Compared to the Spearman rho statistic, Kendall tau provides more accurate p-values when the sample dimension is small and there are many ties in the data [86], as in our case. We derived species trends from [53], where long-term trends (1992-2019) for the species considered in this study were assessed in the same study area, using data derived from the same bird monitoring program. We assigned a zero value to species with non-significant trends.
To investigate the relationship between niche metrics (Schoener's D, expansion and unfilling index) and species traits, we performed distinct analyses for continuous traits (12 life-history traits, namely body length, wing length, tail length, bill length, tarsus length, body weight, clutch size, incubation period, fledging period [87]; number of broods per breeding season [53,87]; annual fecundity [88]; dispersal ratio [89,90]; and six species specialization indices, namely foraging habitat, acquisition behavior, nesting habitat, foraging substrate, diet [53,91,92]; overall [53]) and categorical traits (four uncorrelated traits of interest [53], namely migration strategy [53], nest type [53,87], landscape type-i.e., habitat preferences at landscape scales [53]-, and diet [53,93]). For continuous traits, we performed a principal component analysis (PCA) using all 18 traits as input variables. This way, we reduced the dimensionality of traits into new uncorrelated components representing species life-history and ecological characteristics. Then, we analyzed the correlation between the main axes of the PCA-picked out through eigenvalues, scree plots and percentage of variance explained-and the niche metrics, which were input as supplementary quantitative variables. PCA was performed using the R packages FactoMineR [94] and factoextra [95]. For categorical traits, we performed a Kruskal-Wallis rank sum test separately for each combination of niche metric and trait. In the case of overall significance (p-value ≤ 0.05), we ran a Dunn's test (significant threshold = 0.05) for pairwise comparisons between trait levels with the Bonferroni adjustment [96]. See Table S2 for details and referencing literature about traits. All statistical analyses were performed in R software v. 4.0.2 [71].

Niche Quantification and Temporal Changes
In defining the E-space, the first two PCA axes explained 31.3% of the total variance. PC1 explained 24.2% of the variance, while PC2 7.1%. The first axis represents a thermal gradient, showing a strong negative correlation with elevation and a strong positive correlation with temperatures ( Figure 2, Table S3 and Figure S1a). The major contributions to PC2 were given by nine variables ( Figure S1b). Annual rainfall (rr), broadleaved forests (C311), mixed forests (C313), meadows and pastures (C231), and slope were positively correlated with this component, while arable lands (C211), grasslands (C321), shrublands (C322), and areas with sparse or absent vegetation (C330), showed a negative correlation ( Figure 2 and Table S3). The second niche component reflects habitat characteristics along an ecological gradient of open (variables negatively correlated with PC2) vs closed (variables positively correlated with PC2) habitats, despite the inclusion of meadows and pastures in the latter group. The strong contribution of annual rainfall to PC2 could reflect the major water requirements of natural habitats such as broadleaved and mixed forests. gradient, showing a strong negative correlation with elevation and a strong positive correlation with temperatures ( Figure 2, Table S3 and Figure S1a). The major contributions to PC2 were given by nine variables ( Figure S1b). Annual rainfall (rr), broadleaved forests (C311), mixed forests (C313), meadows and pastures (C231), and slope were positively correlated with this component, while arable lands (C211), grasslands (C321), shrublands (C322), and areas with sparse or absent vegetation (C330), showed a negative correlation ( Figure 2 and Table S3). The second niche component reflects habitat characteristics along an ecological gradient of open (variables negatively correlated with PC2) vs closed (variables positively correlated with PC2) habitats, despite the inclusion of meadows and pastures in the latter group. The strong contribution of annual rainfall to PC2 could reflect the major water requirements of natural habitats such as broadleaved and mixed forests. Niche breadths and centroids showed a large variability across species along each of the two PCA axes (Table S4). Despite differences within species, we detected a general change only for the niche centroid along PC1 between T1 and T2, representing a shift towards warmer climatic conditions (W = 1991, p-value = 0.031; Figure S2).
Niche quantification in T1 and T2, and niche temporal changes of the 71 species were reported in Figure S3. The whole extent and the centroids' position of the available Espace (solid lines and dashed arrows in Figures 3a-d and S3) were largely similar between T1 and T2, meaning that the background E-space in the two periods was comparable. Only a small portion in the bottom area of the available E-space in T2 did not match with the Espace in T1 (shrublands, grasslands and areas with sparse or absent vegetation at higher altitudes). Some species (e.g., Common Quail Coturnix coturnix, Western Yellow Wagtail Motacilla flava, and Cetti's Warbler Cettia cetti) tended to fill a small portion of the total available E-space, while others (e.g., Common Cuckoo Cuculus canorus, Grey Wagtail Motacilla cinerea, Eurasian Wren Troglodytes troglodytes, Common Chaffinch Fringilla coelebs) occurred in most of the available E-space. Regarding changes between T1 and T2 of the density of occurrence and the occupied E-space, different types of patterns can be identified ( Figure S3). Some species, such as the Song Thrush (Turdus philomelos, Figure 3a), ex-  Table S3 for variables' coordinates.
Niche breadths and centroids showed a large variability across species along each of the two PCA axes (Table S4). Despite differences within species, we detected a general change only for the niche centroid along PC1 between T 1 and T 2 , representing a shift towards warmer climatic conditions (W = 1991, p-value = 0.031; Figure S2).
Niche quantification in T 1 and T 2 , and niche temporal changes of the 71 species were reported in Figure S3. The whole extent and the centroids' position of the available E-space (solid lines and dashed arrows in Figures 3a-d and S3) were largely similar between T 1 and T 2 , meaning that the background E-space in the two periods was comparable. Only a small portion in the bottom area of the available E-space in T 2 did not match with the E-space in T 1 (shrublands, grasslands and areas with sparse or absent vegetation at higher altitudes). Some species (e.g., Common Quail Coturnix coturnix, Western Yellow Wagtail Motacilla flava, and Cetti's Warbler Cettia cetti) tended to fill a small portion of the total available E-space, while others (e.g., Common Cuckoo Cuculus canorus, Grey Wagtail Motacilla cinerea, Eurasian Wren Troglodytes troglodytes, Common Chaffinch Fringilla coelebs) occurred in most of the available E-space. Regarding changes between T 1 and T 2 of the density of occurrence and the occupied E-space, different types of patterns can be identified ( Figure S3). Some species, such as the Song Thrush (Turdus philomelos, Figure 3a), extended their niche in most of the whole E-space. Other species exhibited a directional expansion along one of the niche components. For example, the Great Spotted Woodpecker (Dendrocopos major, Figure 3b) showed an expansion (without remarkable niche unfilling) towards coniferous forests and scattered vegetation at higher altitudes. Conversely, other species mainly showed a shift with unfilling and expansion along the same dimension, such as the Eurasian Bullfinch Pyrrhula pyrrhula along PC1 towards warmer climatic conditions, lower altitudes and mixed/broadleaved forests (Figure 3c). Other species showed a relatively stable niche, with moderate areas of unfilling/expansion and non-directional patterns (e.g., the Eurasian Blackcap Sylvia atricapilla, Figure S3). Finally, 17 species (e.g., the Eurasian Wryneck Jynx torquilla, Figure 3d) exhibited an overall niche contraction.
pansion along one of the niche components. For example, the Great Spotted Woodpecker (Dendrocopos major, Figure 3b) showed an expansion (without remarkable niche unfilling) towards coniferous forests and scattered vegetation at higher altitudes. Conversely, other species mainly showed a shift with unfilling and expansion along the same dimension, such as the Eurasian Bullfinch Pyrrhula pyrrhula along PC1 towards warmer climatic conditions, lower altitudes and mixed/broadleaved forests (Figure 3c). Other species showed a relatively stable niche, with moderate areas of unfilling/expansion and non-directional patterns (e.g., the Eurasian Blackcap Sylvia atricapilla, Figure S3). Finally, 17 species (e.g., the Eurasian Wryneck Jynx torquilla, Figure 3d) exhibited an overall niche contraction. (d) Eurasian Wryneck. In each graph, the red area represents the E-space exclusively occupied in T2, the green area the E-space exclusively occupied in T1, and the blue area the E-space occupied in both T1 and T2. The color intensity of the filled area represents the density of occurrence of the species in T2. Green and red solid lines represent the whole extent of the available E-space in T1 and T2, respectively. Arrows indicate the centroids' shift from T1 to T2 for the available E-space (dashed arrow, where we did not detect any shift) and for the occupied E-space (solid arrow).
Across the 71 species, the median value of niche overlap between T1 and T2, assessed through the Schoener's D index, was 0.60 (range = 0.14-0.81). Most species showed low or moderate values of the expansion index (median = 0.04, range = 0-0.66), but the Dunnock (Prunella modularis), the Song Thrush, the Water Pipit (Anthus spinoletta), the Willow Tit (Poecile montanus), and the Eurasian Skylark (Alauda arvensis) experienced a large expansion (>0.40). Similarly, the median of the unfilling index was 0.06 (range = 0-0.32), but the (d) Eurasian Wryneck. In each graph, the red area represents the E-space exclusively occupied in T 2 , the green area the E-space exclusively occupied in T 1 , and the blue area the E-space occupied in both T 1 and T 2 . The color intensity of the filled area represents the density of occurrence of the species in T 2 . Green and red solid lines represent the whole extent of the available E-space in T 1 and T 2 , respectively. Arrows indicate the centroids' shift from T 1 to T 2 for the available E-space (dashed arrow, where we did not detect any shift) and for the occupied E-space (solid arrow).
Mallard (Anas platyrhyncos), the Common Kestrel (Falco tinnunculus), the Common Quail, the Eurasian Wryneck, the Melodius Warbler (Hippolais polyglotta), the Marsh Tit (Poecile palustris) and the Carrion Crow (Corvus corone) showed higher unfilling values (>0.25) (Figure 4 and Table S5). The Schoener's D resulted negatively correlated with both the expansion (Kendall tau = −0.34, p-value < 0.001) and the unfilling (Kendall tau = −0.27, pvalue = 0.001), and the expansion was negatively correlated with the unfilling (Kendall tau = −0.26, p-value = 0.001). Niche equivalency and similarity tests were summarized in Table 1. Twenty-three of the 71 species showed a niche equivalency significantly smaller than expected by chance (equivalency test: p-value ≤ 0.05 for the niche divergence hypothesis), suggesting the existence of differences in niches between T1 and T2. On the other hand, seven species (Common Cuckoo, Eurasian Wren, European Robin Erithacus rubecula, Common Blackbird Turdus merula, Lesser Whitethroat Curruca curruca, Eurasian Blackcap, Goldcrest Regulus regulus) significantly conserved their niche between the two periods (equivalency test: pvalue ≤ 0.05 for the niche conservatism hypothesis). The similarity test did not accept the niche divergence hypothesis for any species (p-value > 0.05). However, the niche conservatism hypothesis was accepted for 56 species (p-value ≤ 0.05 whose Schoener's D resulted to be among the 59 highest values, showing that most species retained their own niche over time when both occurrence density grids were randomly shifted into the background environment. Niche equivalency and similarity tests were summarized in Table 1. Twenty-three of the 71 species showed a niche equivalency significantly smaller than expected by chance (equivalency test: p-value ≤ 0.05 for the niche divergence hypothesis), suggesting the existence of differences in niches between T 1 and T 2 . On the other hand, seven species (Common Cuckoo, Eurasian Wren, European Robin Erithacus rubecula, Common Blackbird Turdus merula, Lesser Whitethroat Curruca curruca, Eurasian Blackcap, Goldcrest Regulus regulus) significantly conserved their niche between the two periods (equivalency test: p-value ≤ 0.05 for the niche conservatism hypothesis). The similarity test did not accept the niche divergence hypothesis for any species (p-value > 0.05). However, the niche conservatism hypothesis was accepted for 56 species (p-value ≤ 0.05 whose Schoener's D resulted to be among the 59 highest values, showing that most species retained their own niche over time when both occurrence density grids were randomly shifted into the background environment.

Species Schoener's D E (D) E (C) S (D) S (C)
Great    Regarding the relationship between niche metrics and continuous traits, PCA reordered them into four principal axes. Dim1 (40.5% of the total variance explained) represented a gradient in body size (positive correlation), Dim2 (19.4%) was related to most of the species specialization indices (positive correlation), Dim3 (10.3%) to clutch size (positive correlation), and Dim4 (7.9%) was positively correlated with the number of broods per year. Dispersal ratio was negatively correlated with both Dim3 and Dim4, and positively correlated with Dim1 (Table S6). Analyses revealed a weak relationship among the niche metrics and the four PCA dimensions (Table S6). Specifically, the unfilling index showed a slight positive correlation with Dim1 (0.23)-i.e., body size-and Dim3 (0.19)-i.e., clutch size- (Figure 6a), and a weak negative correlation with Dim4 (-0.17), thus demonstrating an association with species producing fewer broods ( Figure S4a,b). The expansion index resulted in a weak and negative correlation (−0.13) with Dim1-i.e., the higher the expansion, the smaller the body size (Figure 6a,b)-and a weak and positive correlation (0.10) with Dim4-i.e., the higher the expansion, the higher is the number of broods- (Figure S4a,b). The dispersal ratio showed a weak negative association with the expansion index (Figure 6a and Figure S4a) and a non-unidirectional relationship with the unfilling index (Figure 6a). The Schoener's D highlighted a slight negative correlation (−0.12) with Dim3 ( Figure 6a) and a slight positive correlation (0.10) with Dim2 (Figure 6b), the latter suggesting that specialist species tended to retain their niche more than generalist species. niche metrics and the four PCA dimensions (Table S6). Specifically, the unfilling index showed a slight positive correlation with Dim1 (0.23)-i.e., body size-and Dim3 (0.19)i.e., clutch size- (Figure 6a), and a weak negative correlation with Dim4 (-0.17), thus demonstrating an association with species producing fewer broods ( Figure S4a,b). The expansion index resulted in a weak and negative correlation (−0.13) with Dim1-i.e., the higher the expansion, the smaller the body size (Figure 6a,b)-and a weak and positive correlation (0.10) with Dim4-i.e., the higher the expansion, the higher is the number of broods- (Figure S4a,b). The dispersal ratio showed a weak negative association with the expansion index (Figures 6a and S4a) and a non-unidirectional relationship with the unfilling index (Figure 6a). The Schoener's D highlighted a slight negative correlation (−0.12) with Dim3 ( Figure 6a) and a slight positive correlation (0.10) with Dim2 (Figure 6b), the latter suggesting that specialist species tended to retain their niche more than generalist species.  Table  S2 for further details on species traits.

Relationship between Niche Temporal Changes, Population Trends and Species Traits
Regarding categorical traits, significant differences were detected for the unfilling index respect to migration strategy (Kruskal-Wallis χ 2 = 7.91, df = 2, p-value = 0.019), wherein Dunn's test highlighted a significant greater unfilling for long-distance compared to short-distance migrants (p-value = 0.020; Figure 7a). Moreover, all niche metrics showed Abbreviations' meaning for the continuous traits: len = body length; wing = wing length; tail = tail length; bil = bill length; tar = tarsus length; wei = body weight; clu = clutch size; bro = number of broods per breeding season; fec = annual fecundity; inc = incubation period; fle = fledging period; disr = dispersal ratio; SSI. = prefix for the specialization indices; fh = foraging habitat; acq = acquisition behavior; nes = nesting habitat; fs = foraging substrate; tr = diet; ov = overall. See Table S2 for further details on species traits.
Regarding categorical traits, significant differences were detected for the unfilling index respect to migration strategy (Kruskal-Wallis χ 2 = 7.91, df = 2, p-value = 0.019), wherein Dunn's test highlighted a significant greater unfilling for long-distance compared to short-distance migrants (p-value = 0.020; Figure 7a). Moreover, all niche metrics showed differences with respect to the landscape type (unfilling index: Kruskal-Wallis χ 2 = 10.24, df = 3, p-value = 0.017; expansion index: Kruskal-Wallis χ 2 = 19.74, df = 3, p-value < 0.001; Schoener's D index: Kruskal-Wallis χ 2 = 13.48, df = 3, p-value = 0.004). Specifically, pairwise comparisons showed a greater unfilling for farmland birds compared to species inhabiting several types of landscapes (group "several", p-value = 0.014; Figure 7b), and a greater expansion for woodland birds and species inhabiting natural-open habitats compared to farmland species (p-value = 0.010 and p-value < 0.001, respectively; Figure 7c). Finally, species of natural-open habitats showed a lower niche overlap (Schoener's D) compared to both the group of "farmland" and "several" (p-value = 0.044 and p-value = 0.002, respectively; Figure 7d). wise comparisons showed a greater unfilling for farmland birds compared to species inhabiting several types of landscapes (group "several", p-value = 0.014; Figure 7b), and a greater expansion for woodland birds and species inhabiting natural-open habitats compared to farmland species (p-value = 0.010 and p-value < 0.001, respectively; Figure 7c). Finally, species of natural-open habitats showed a lower niche overlap (Schoener's D) compared to both the group of "farmland" and "several" (p-value = 0.044 and p-value = 0.002, respectively; Figure 7d).  Table S2 for further details on species traits.  Table S2 for further details on species traits.

Niche Quantification and Temporal Changes
Niche temporal changes of 71 common breeding birds in a wide region of northern Italy over 26 years revealed a high complexity of patterns across species. The equivalency tests supported the niche divergence hypothesis for 23 species (32%). This means that their niche occupancy in the E-space differed between T 1 and T 2 , suggesting that these species changed their niche in relation to the accessible environmental conditions over time. Evidence for niche divergence requires two conditions [78,98]: (1) niche characteristics differ between ranges, and (2) these differences are greater than background environmental divergence. However, Rodriguez-Rodriguez et al. [99] and Cuervo et al. [100] argued that niches are already divergent when found less equivalent than expected by chance-i.e., significant equivalency test-and the niche conservatism hypothesis evaluated by the similarity test is not accepted. In our study, only seven species (Mallard, Common Kestrel, Eurasian Skylark, Dunnock, Song Thrush, Melodius Warbler, and Eurasian Bullfinch) satisfied this criterion, thus showing statistical evidence of niche divergence. For the other 16 of the 23 species, for which the equivalency test supported niche divergence (Table 1), we cannot confirm this hypothesis since the similarity test retained a larger niche similarity. We found strong evidence of niche conservatism (both the equivalency and the similarity tests supported the hypothesis) for the Common Cuckoo, the Eurasian Wren, the European Robin, the Common Blackbird, the Eurasian Blackcap, the Lesser Whitethroat, and the Goldcrest. The first five species occupied a large portion of the available E-space in both periods ( Figure S3), showing they were able to use a wide spectrum of climatic and habitat conditions. On one hand, it is possible they did not undergo great ecological pressures in the occupied E-space, on the other hand the environmental changes that occurred in the study area may have not significantly affected their niche because of a large range of environmental tolerance. Conversely, the Goldcrest is a forest specialist species and it is likely less prone to modify its habitat and climatic requirements even under ongoing environmental changes (it is actually the only woodland bird along with the Common Chiffchaff Phylloscopus collybita suffering a population decline in the study area [53]). Regarding the Lesser Whitethroat, we underline that the species showed a Schoener's D (0.56) just smaller than the overall median, and we did not exclude that the species rarity, especially in the first period (frequency of occurrence = 1%), could have affected the significance of the tests. Caution is in fact recommended when inference is conducted on rare species [44].
The niche dynamic indices, expansion and unfilling provided useful additional information on niche temporal changes. Differently from the Schoener's D, they assess the niche overlap regardless of the density of occurrence. We observed that 27% and 31% of the species experienced a niche expansion and unfilling greater than 10% (values > 0.1), respectively. Moreover, as highlighted by the correlation analysis between the two indices, species showing an expansion tended to not experience a niche unfilling, and vice versa. Furthermore, our findings suggested that decreasing niche overlap (Schoener's D) is associated alternatively with niche expansion or unfilling, underlining the role of occupying new habitats, or losing original habitats, in determining differences in the observed niches.
We found an overall shift of the niche centroids towards greater values along PC1. This means a general tendency to experience higher thermal conditions, but the relationship between habitat and climate, often acting in a synergistic way [101][102][103], should be considered. Our study area is characterized by an important altitudinal range (from 2 m to 4000 m a.s.l.), and habitat types, as well as anthropization intensity, are strictly associated with the elevational gradient and thus with temperatures. The characteristics of open and closed habitats differ according to the altitudinal range (e.g., open habitats are essentially wild in uplands while mainly anthropized in lowlands), and such differences must be taken into account to understand niche temporal changes. The general centroids' shift towards warmer climates can actually depend on several drivers that can be species-specific. Firstly, a species may undergo higher temperatures because of an increase of temperatures in the environment (we found an increase of the average temperature from 9.53 • C to 10.54 • C during the study period in the whole study area; data elaborated from the ECAD dataset [68]). Several studies underlined the ability of birds to track their climatic niche (e.g., [104,105]), while lags may be observed especially under rapid environmental change on a short timescale [34,43,106]. In the presence of such lag-effects, or the absence of nichetracking for several reasons (e.g., species thermal physiological tolerances, species dispersal ability, habitat availability and landscape connectivity), raising temperatures might result in a passive shift of niche centroids towards warmer climates. Species responses to this phenomenon may be idiosyncratic and difficult to predict [106], although some groups are expected to suffer more than others (e.g., mountain birds [107][108][109][110]). Secondly, a species may move towards lower altitudes, and therefore to warmer climates, in response to fine-scale changes in the habitat structure and characteristics (including anthropogenic ones as forestry management, agricultural practices, intensity of grazing and mowing). For example, we found that the Song Thrush, the Common Firecrest (Regulus ignicapilla), and the Eurasian Bullfinch showed a shift in density of occurrence from high-altitude coniferous to middle and low-altitude mixed and broadleaved forests ( Figure S3). Finescale changes, not detectable in our study, can directly affect the species occurrence and habitat suitability within a specific type of habitat [84,111]. In addition, modifications in habitat availability due to land uses changes might push a species to track the habitat component of its niche to find similar conditions. In our study, for example, the Mistle Thrush (Turdus viscivorus) moved its centroid towards high-altitude shrublands, maybe as response to an expansion of forests upwards and neo-colonization of wood at the expense of abandoned meadows in uplands [112]. We also remark that some species may have extended their niche towards lower altitudes independently from environmental pressures, but as response of density-dependent factors (e.g., increasing populations, [45] and see Section 4.2) leading to colonize new environments. For example, the Black Redstart (Phoenicurus ochruros), a species originally associated with open-natural habitats, showed a meaningful increment in density of occurrence in closed habitats (mixed and broadleaved forests) and in lowland open habitats (where it colonized human buildings).
Despite we are confident of consistency of our results, we acknowledge that some limitations could affect species niche quantification. Specifically, we could not account for species detection probability because data did not include multiple surveys in the same site. Detection probability may represent an issue when comparing niche across species [113]. However, intra-species niche comparisons are not supposed to be considerably affected by detection probability when comparing niche metrics of a species across different geographic areas [15] or over time. In fact, in the case of a standardized sampling design, the species habitat-specific detection function is supposed to be constant over time [114], and its effect on niche comparison may be considered negligible. Moreover, our data collection and its aggregation into three-year time intervals further reduced the potential noise arising from imperfect detection (see Section 2.1). Measuring the effect of detection probability in niche modeling is certainly an issue that deserves attention in future research.

Relationship between Nyche Temporal Changes, Population Trends and Species Traits
The assessment of the relationship between niche metrics and population trends highlighted interesting patterns, partially confirming findings from Ralston et al. [45]. We found a positive correlation between population trends and niche expansion, while unfilling did not correlate with trends. Metapopulation theory predicts that in growing metapopulations with many individuals dispersing from local "source" populations, it is likely that local "sink" populations settled in suboptimal environments exist. On the contrary, in declining metapopulations it is probable that only the largest local populations persist in the most suitable environments, but not enough migrants are generated to colonize suboptimal environments and to support local "sink" populations [115][116][117]. However, we cannot exclude that a niche expansion due to stochastic phenomena, interspecific interactions with competitive exclusions, and/or fine-scale increase of habitat quality, could have pushed the species to extend the niche with a positive secondary effect on population size. In any case, some meaningful exceptions can be found for the positive relationship between niche expansion and population trends. We remark that the Eurasian Skylark, the bird undergoing the most negative population trend in Lombardy from 1992 to 2019 [53], showed one of the highest niche expansions during this time interval. In the first period (1992-1995-1996) [118][119][120][121], and the fact that the Eurasian Skylark retained the original core area may be due to a partial lag-effect for which the species remains present in unsuitable environments that were previously suitable [122]. Additionally, we do not exclude that the two local populations that seem to emerge may undergo different ecological pressures leading to divergent population dynamics between the two areas. This example highlights the importance of integrating demographic responses with niche temporal changes, especially to implement conservation efforts and actions for threatened species.
Ralston et al. [45] reported a significant relationship between niche unfilling and bird trends for the realized climatic niche in 46 birds breeding in North America. By also including habitat characteristics in our niche modeling, we did not find any association.
Overlooking the fact that our study was carried out in a different biogeographic area and at very smaller spatial scale, increasing the complexity of the niche multidimensional space can reveal different patterns. Embracing this complexity can be the key to a better interpretation of the niche-trend relationship. Niche shifts in the occupied E-space-i.e., unfilling followed by expansion-could also mask the actual relationship between unfilling and population trends, but the negative correlation emerged between the unfilling and the expansion indices suggests excluding this hypothesis. We also found a negative correlation between the niche overlap and population trends very close to the significance threshold (p-value = 0.053). This finding might suggest that species tending to retain their niche more than others (in terms of niche occupancy within the E-space), may be more prone to demographic declines. This static response may represent a drawback when changes occur in the environment or in interspecific relationships, leading to population decline in a long-term perspective. Finally, dispersal limitation, depending on habitat configuration but also on species-specific characteristics, may play a crucial role in metapopulation and niche changes [117], and integrating this parameter in the modeling framework may improve our ability to understand this complex relation.
Regarding our last question, results highlighted a moderate relationship between niche temporal changes and species traits. Overall, the unfilling index was the most related with continuous traits. Findings suggested that birds with larger body size are more related with niche unfilling compared to smaller species. Large-bodied species have generally smaller population sizes, lower reproductive rates, and larger home or geographic range requirements than small-bodied species; moreover, they typically occupy higher trophic levels [123]. All of these characteristics make these species more prone to extinction [124,125], and understanding their niche changes could provide important insights to enhance conservation strategies.
Interesting findings also emerged for the relationship between unfilling and the reproduction-related traits (number of broods per year, average clutch size, fecundity). Niche unfilling grew with increasing clutch size-average number of eggs laid within a single brood-and with decreasing number of broods per year. It may suggest that species investing much effort in increasing the number of broods, rather than producing more eggs per brood, might be more adaptable to environmental changes and may be able to expand their niche. However, annual fecundity-a product of clutch size and number of broods [88]-did not show a straightforward relationship with niche unfilling, and further studies should investigate the issue.
Dispersal ratio, a surrogate of species mobility and dispersal ability [89,90], showed that species with higher dispersal ability would not tend to expand the niche, but findings were inconclusive in respect to unfilling. We also have to consider that the dispersal ratio could be an inappropriate indicator of the dispersal ability [64,126], as well as that habitat configuration and landscape connectivity significantly affect the dispersal process [127], leading to idiosyncratic responses between species. Regarding the degree of species specialization (Dim2 of PCA), niche metrics were not linked with it, although a weak signal showed that specialist birds tend to retain their niche more than generalist ones. It supports the hypothesis that specialists require a narrow combination of environmental characteristics [128], and they are less prone to modify their ecological needs in response to environmental change [65].
Concerning the categorical ecological traits, neither nest type nor diet showed a relationship with niche metrics. Conversely, long-distance migrants-wintering in Sub-Saharan Africa-showed a higher niche unfilling compared to short-distance migrants-wintering in Europe or North Africa-and sedentary species, although statistically significant differences emerged only in the first comparison. Phenological mismatch [129,130] may represent a putative driver for long-distance migrants. Despite the fact that species can attempt to adjust to advancement in spring phenophases through behavioral plasticity, for example, by a reduction in time between arrival and breeding [131] or anticipating the breeding ground arrival [132,133], these adjustments may remain insufficient to track phenological shift [129,134,135]. Mismatches might lead, for example, to an asynchrony between food availability and species demand, generating spatially consistent directional selection on timing, which could promote rapid evolutionary change [135]. Moreover, prey availability may affect habitat suitability, as recently documented in raptors [136], leading to potential niche mismatch over time. Under climate and environmental changes, long-distance migrants may be outcompeted by sedentary species having similar niches, due to the residents' ability to better judge the onset of the breeding season [137]. Such effects could affect long-distance migrants more than short-distance ones, which may be more capable of anticipating breeding ground conditions due to the proximity of their wintering areas to the breeding areas [138]. In addition, niches could differ between wintering and breeding grounds [58], with separate dynamics over time also due to different drivers acting in wintering and breeding areas [139]. Exploring long-term niche dynamics throughout several phases of the annual cycle (breeding, migration, wintering) represents a key point to look at for future research.
We also found that species inhabiting natural-open landscapes (which are all birds inhabiting mountain areas) and woodland birds underwent a higher niche expansion in respect to farmland species. Overall, birds of natural-open habitats extended the niche along the open-closed habitat gradient, within the limit of their ecological requirements, as well as towards greater thermal conditions (see Section 4.1 for discussion on this point) and, as in the case of the Black Redstart, towards lowlands. Within woodland species, the direction of expansion showed a large interspecific variability. Some species extend the niche towards different type of forests (e.g., the Eurasian Bullfinch), other towards multiple directions of the available E-space (e.g., the Eurasian Nuthatch Sitta europaea) or along boundaries of the niche (e.g., the Common Chiffchaff). Forest species populations are increasing in Lombardy [53]. This can further support the hypothesis that population demography may act as a driver to affect niche expansion over time. On the contrary, farmland birds experienced a greater niche unfilling compared to species inhabiting several types of landscapes. Farmland birds are negatively influenced by intensification of agricultural practices [140]. They could suffer fine-scale environmental changes (e.g., grazing intensity, mowing, use of pesticides, land consolidation) affecting habitat suitability, as well as a reduction of available habitats because of the reduction of agricultural lands that occurred in our study area over the considered period. Finally, the greater niche overlap of species inhabiting several types of landscapes and farmland birds compared to species of naturalopen landscapes may reflect a larger environmental tolerance for the first group and more niche fidelity for the second group, leading, in the latter case, to an overall niche unfilling when environmental changes occurred.

Conclusions
Even though the study of long-term species niche dynamic could provide crucial information for investigating species responses to environmental changes, it is a poorly investigated topic in ecology. Our results provided evidence of the complexity of patterns in niche temporal changes across 71 birds breeding in a temperate region. Results supported niche divergence (equivalency test) for 32% of species, although two-thirds were not Birds 2022, 3 166 supported by the similarity test. We highlighted a general tendency to adjust the niche centroids towards warmer thermal conditions, which may depend on several drivers acting in an idiosyncratic way between species. Embracing different niche metrics, as well as tests for niche inference, is essential to understand niche changes.
To date, a totally overlooked aspect is the potential relationship between niche changes and population trends (but see [44,45]) and species traits (but see [57,58]). We found a positive association between niche expansion and increasing populations, and a signal that species retaining their niche over time could be more prone to undergo population declines. Moreover, we found moderate evidence for a non-random association between niche changes and species traits, especially for body size, clutch size, number of broods per year, inhabited landscape type, and migration strategy.
We encourage the implementation of studies incorporating long-term population trends and niche dynamics with a focus on cause-effect relationship. Moreover, we highlight the importance of enhancing knowledge about the role of species traits in determining niche changes over time, but also working at the single-species level due to the specificity of responses. Implementing the modeling framework by including the potential variations of niche throughout several phases of annual cycle (breeding, migration, wintering), bionomic variables (biotic interactions and resource-consumer dynamics [47,141]), and dispersal dynamics represents a crucial point. We believe that this research roadmap could provide meaningful insights to improve conservation plans aimed at preserving threatened species and biodiversity.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/birds3010011/s1, Table S1: Land-cover variables used in PCA (E-space), Table S2: List of species traits and referencing literature, Table S3: Coordinates of each variable respect to the PCA-axes (E-space), Figure S1: Contribution of each variable (%) to PC1 (a) and PC2 (b) (E-space), Table S4: Niche breadth and centroid for each species along the PCA-axes (E-space) in T 1 and T 2 , Figure S2: Notched box plots of niche breadth and centroid for each PCA-axis (E-space) in T 1 and T 2 , Figure S3: Niche quantification in the E-space for each species in T 1 (1992-95-96) and T 2 (2015- [16][17] and niche changes between the two periods, Table S5: Schoener's D, expansion and unfilling index for each species, Table S6: Coordinates of the 18 continuous species traits and of niche metrics respect to the four retained PCA axes (trait space), Figure S4  Data Availability Statement: Restrictions apply to the availability of data. Part of the data was obtained from the General Directorate for Agriculture of the Lombardy Region and part was collected by authors. Full dataset is available at the Department of Earth and Environmental Sciences, University of Milan-Bicocca, Piazza della Scienza 1, 20126, Milan, Italy, upon request to the corresponding author, after the permission from the General Directorate for Agriculture of the Lombardy Region.