Long-term dynamic of nestedness in bird assemblages inhabiting fragmented landscapes

Nestedness is a common pattern of species assemblages in fragmented landscapes. The spatial pattern and ecological drivers of nested communities have been widely explored, but few studies investigated their long-term variability. To investigate the variability of nestedness and species-specific fragment occupancy of forest birds in a fragmented landscape affected by environmental changes over 16 years. Data were obtained from the Monitoring Program of Breeding Birds in Lombardy (northern Italy). For two study periods (1997–2001, 2013–2017), we tested for overall nestedness and for sites and species nestedness independently using the NODF metric. We tested for nestedness drivers (variable selection on multiple linear regression models) and evaluated the effect of species ecological traits on fragment occupancy changes (variable selection on multiple linear regression model). The community showed a significant nestedness driven by both selective extinction and selective colonization in both study periods. Sites nestedness was significant in the second study period only. Over 16 years the effect of distance from source areas was completely lost and only local isolation conditions drove selective colonization in the second study period. Between the two study periods, we discovered a general occupancy decrease of interior species and a significant occupancy increase of generalist species characterized by large size and large minimum area requirement. Nestedness drivers of the investigated community significantly varied over time probably because of both environmental and demographic changes. Long-term studies are crucial to explore spatial pattern changes and to address management strategies for species conservation in fragmented landscapes.


3
Vol:. (1234567890) 2015; Matthews et al. 2015;Betts et al. 2017;Chase et al. 2020). Investigating how the destruction and fragmentation of natural habitats affect animal communities is of fundamental importance for biodiversity conservation (Baselga et al. 2015). Traditionally, most studies assumed that species richness was a proper metric to measure community changes in fragmented landscapes (Banks-Leite et al. 2012), but during the last decades the study of nested subset patterns in communities has become more and more popular (Fischer and Lindenmayer 2005;Heino et al. 2009;Ulrich 2009).
A community occupying multiple discrete sites is considered nested when species assemblages occurring at relatively species-poor sites are proper subsets of species assemblages occurring at progressively more species-rich sites (Patterson and Atmar 1986). This nonrandom pattern suggests that ecological mechanisms structure the community in a predictable way (Lindenmayer and Fischer 2013). Nestedness of communities inhabiting islands or habitat fragments can be caused by four main mechanisms: (1) selective extinction, (2) selective colonization, (3) habitat nestedness, and (4) passive sampling (Patterson and Atmar 1986;Worthen 1996). If a community ordered by the area of occupied sites is significantly nested, selective extinction is the dominant mechanism and the community corresponds to a set of predictable extinction sequences (Lomolino 1996;Wright et al. 1997). Under this hypothesis, fragment area is the main driver of nestedness because species loss follows gradients of species sensitivity to habitat size, i.e., species with large area requirements and small range size will disappear first when fragments' area is reduced (Wright et al. 1997;Watling and Donnelly 2006;Li et al. 2019). Selective colonization may produce nestedness because species have different dispersal ability and sites are characterized by different degrees of isolation (Patterson and Atmar 1986): strongest dispersers (e.g., species with high body mass and dispersal ratio) occur in both less and highly isolated sites, whereas poorer dispersers (e.g., species with low body mass and dispersal ratio) occur in less isolated sites only (Heino et al. 2009). Habitat nestedness is the dominant ecological mechanism when the nestedness of species assemblages is due to their reliance on habitats that in turn have a nested distribution (Calmé and Desrochers 1999;Honnay et al. 1999;Tan et al. 2021). Finally, passive sampling can be a major driver of nestedness if sites are simply more likely to be occupied by the most abundant species (Cutler 1994;Higgins et al. 2006). Despite the identification of nestedness drivers should involves analyses that account for both site characteristics and species ecological traits ), few studies analyzed their combined roles in generating nestedness (Dardanelli and Bellis 2021).
Distribution patterns in nature are known to change over time (Soininen et al. 2018) making essential studies that analyze the temporal dynamics of animal community characteristics, such as nestedness (Wilby and Shachak 2000;Klaassen and Nolet 2008;Tunberg and Krång 2008). Luther et al. (2020) asserted that repeated assessments are crucial for investigating how habitat fragmentation affects nestedness dynamics both in the short and in the long-term. In fact, the response of animal communities to environmental changes occurring both at local (e.g., changes in habitat quality) and landscape (e.g., habitat loss and fragmentation) scale can be often noticeable over long periods due to system inertia, which imperil the observation of short-time responses (Lindenmayer and Fischer 2013). While several studies have focused on the short-term variability of nestedness (e.g., Patterson 1990;Fernández-Juricic 2002;Bloch et al. 2007;Murgui 2010;Seoane et al. 2013;Zhou and Chu 2014; De la Hera 2019), very few studies have investigated the long-term variation in nestedness causality and results were divergent (Sebastián-González et al. 2010).
In this study, we investigated the long-term spatial dynamics of bird assemblages inhabiting forest fragments in a wide area of northern Italy. The region was originally covered by large broadleaved forests and has been historically fragmented by human activities; the more recent and impacting was the deforestation during World Wars (Gallinaro 2004). From the second half of the twentieth century forest cover partially increased and then stabilized in the last decades, but severe land cover changes occurred in the surrounding landscapes. To evaluate if these changes could have affected the bird community composition in the study area, we compared the spatial patterns of bird assemblages in two study periods (1997-2001 and 2013-2017) 16 years apart (considering the central year of the study periods). In this time interval, forest cover did not change in the study area, but the landscape matrix underwent different land cover and land use changes. Specifically, urban areas increased (+ 2.8%) at the expense of arable lands (comparison of digital land cover maps ERSAF 2002ERSAF , 2017, whereby an intensification of agricultural practices occurred (Lassini et al. 2007). This process could have reduced the connectivity of the study area for forest birds by reducing matrix permeability due to the increase in monocultures, use of chemical fertilizers and pesticides, mechanization, and irrigation and to the loss of semi-natural elements (e.g., hedgerows) (Pellegrini and Fernández 2018). Specifically, in the study area, about 16% of continuous hedgerows were substituted by discontinuous ones during the considered time interval (comparison of digital hedgerow layers ERSAF 2002ERSAF , 2017. Discontinuous hedgerows are typically characterized by poorer vegetation structure and lower vegetation density and floristic diversity and are often unsuitable as connective elements (Dondina et al. 2016). We hypothesized that the decrease of matrix permeability in the study area might have affected the distribution of dispersal-limited species (such as interior species with low dispersal abilities) altering the spatial pattern of the whole community.
To detect possible community changes during the considered time interval, for both study periods, we investigated whether a significant nested subset pattern exists and examined three possible driving ecological mechanisms: selective extinction, selective colonization, and the statistical artefact of passive sampling. We did not test the habitat nestedness hypothesis because the habitat composition of fragments in our study area was relatively homogeneous (see "Materials and methods" section). Moreover, we explored possible changes of fragments occupancy by individual species during the study period and evaluated if they were driven by species ecological traits. Finally, the results emerged from the analyses were discussed to provide insights for a conservation-oriented management of the landscape in the study area.

Study area
The study was carried out in an area of c. 15,300 km 2 in Lombardy (northern Italy). The area is composed of three main sub-regions: the Prealps in the North, the lowland (below 300 m a.s.l.) in the central part, and the Apennines in the South (Fig. 1). The area is crossed from West to East by the Po River and from North to South by four main rivers, namely the Ticino, Adda, Oglio and Mincio. Broadleaved forests cover c. 12.5% of the study area in both study periods (ERSAF 2002(ERSAF , 2017. Continuous forests occupy the hilly and mountain areas of the Prealps and the Apennines and the banks of the Ticino River in the western part of the lowland. The remaining part of the lowland counts several small residual forest fragments (94.2% of which are smaller than 10 ha; 1999: 94. 4% ERSAF 2002: 94.1% ERSAF 2017 scattered in the intensively cultivated lowland area.

Bird data
Bird data were acquired from the Monitoring Program of Breeding Birds in Lombardy (Bani et al. 2009). Under this project, surveys were carried out almost continuously from 1992 to 2021 at sampling units randomly extracted each year at least 500 m apart (Tirozzi et al. 2021). Data were collected through a standardized method based on the 10-min unlimited distance point count technique (Blondel 1981;Fornasari et al. 1998): surveys were carried out each year during the breeding season (10 May to 20 June), from sunrise to 11.00 a.m. only in good weather conditions (Bani et al. 2009).
For this study, we retained all point counts performed in forest fragments or within 250 m from their boundaries. We selected this threshold distance because it corresponds to the maximum detection distance for most of the species surveyed . For our analyses, we only considered forestdependent species, i.e., species breeding in forest habitats. Forest-dependent species were then classified into three groups characterized by a different forest specialization: generalist (species linked to forest remnants but tolerant towards anthropogenic matrix), interior (species intimately linked to forest remnants), and edge species (species that require transitional forest habitats) (we followed the classification proposed by Bani et al. 2006 revised based on Cramp and Simmons 2006). We discarded edge species from the dataset because their response to forest fragmentation is the opposite of that of interior and generalist species  and their inclusion in the dataset could mask possible nested subset patterns.
To investigate the long-term spatial dynamics of bird assemblages we selected the point counts performed between 1997 and 2001 (first study period) and between 2013 and 2017 (second study period). The bird community considered for each study period was composed of the whole set of generalist and interior forest-dependent species detected during the 5-year time interval to avoid the effect of interannual variability and to obtain a more balanced and robust sample (see Rocchia et al. 2018 for a similar approach).
Since the project did not include multiple surveys within the same season, species detection probability cannot be explicitly considered. Nevertheless, the large amount of data used for this study and the aggregation of data into 5-year time intervals overcome the potential limit of imperfect detection, lowering the noise generated by stochasticity in species detection Bani et al. 2019). Moreover, 10-min point counts have been indicated as an effective survey technique to detect a sufficiently high percentage of the species present in a site (Fuller and Langslow 1984;Bani et al. 2019).

Landscape data
We used the land-use digital maps available for the central year of each study period (1997-2001: DUSAF 1.1 update to 1999ERSAF 2002-2017: DUSAF 5 update to 2015ERSAF 2017. For each study period, we selected forest fragments Fig. 1 Location of the study area in Lombardy (northern Italy) and the study area with broadleaved forest cover in green. Forest layer refers to DUSAF 5 (ERSAF 2017) (second study period). Neglectable forest cover changes (+ 2%) have occurred between the first and second study period (DUSAF 1.1 update to 1999ERSAF 2002DUSAF 5 update to 2015ERSAF 2017 belonging to one of the following categories: broadleaved forests (DUSAF class: 311), mixed forests (DUSAF class: 313) and broadleaved reforestations (DUSAF class: 2242). In ArcGIS 10 (ESRI 2011), we merged all fragments that were < 25 m apart. This measure corresponds to the typical width of a secondary road (easily crossed even by the less mobile bird species), whose digitalization often leads to artificial subdivisions of forest fragments ). To control for a possible sample size effect, we only retained forest fragments with a point count density higher than 1 point/500 ha. We thus obtained 844 and 543 forest fragments for the first and second study period, respectively (Fig. S1 in Online Resource 1). In both study periods, the composition of forest patches was highly homogeneous (94% and 84% of forest plots merged to obtain forest fragments was composed of broadleaved forests in the first and second study period, respectively) making the variability related to the forest type very low. The pool of forest fragments selected for the two periods was not the same because sampling units were randomly extracted each year but, given the large number of fragments selected for each study period, we are confident that they are both representative of the bird communities in the study area (see Canedoli et al. 2017 for a similar approach).

Nestedness analyses
For each study period, we built an occurrence and an abundance matrix by associating to each fragment (rows) the presence (1) or absence (0) of each species (columns) and the maximum number of detected individuals, respectively.
We used the occurrence matrix to test for nestedness pattern in bird assemblages in each study period. The degree of nestedness was quantified through the NODF metric (Nestedness metric based on Overlap and Decreasing Fill) (Almeida-Neto et al. 2008). NODF is based on two properties: decreasing fill (i.e., standardized differences in occurrence matrix row and column fills) and paired overlap (i.e., the overlap of presences in two adjacent columns), and is largely considered the most robust nestedness index (Matthews et al. 2015). To determine whether the observed NODF score for the two time periods was significantly higher (i.e., nested) than those expected by chance (null hypothesis), we compared the NODF score of the maximally packed matrix with the NODF scores of 500 highly constrained random community matrices using the tswap permutation method (preserving both row and column totals) (Miklós and Podani 2004). Nestedness was estimated for both the whole occurrence matrix and independently for species (NODF among rows, NODFr) and sites (NODF among columns, NODFc). The analyses were carried out by using the oecosimu function with nestednodf method in the 'vegan' package (Oksanen et al. 2019) in R (R Core Team 2020).

Nestedness drivers
To test for the effect of selective extinction and selective colonization in structuring the nested subset pattern in the two time periods, we investigated the statistical relationship between the fragments and species rank in the maximally packed matrix with ranked fragments metrics and species ecological traits, respectively.
We selected six metrics to characterize the spatial configuration of forest fragments: Area, defined by fragment surface (ha); Shape index, accounting for the edge effect and calculated as the ratio between patch area (ha) and perimeter (m); DNSA, accounting for source area effect and defined by the distance (m) of the fragment from the nearest source area (i.e., the nearest forest fragment with an area higher than 1000 ha, Dondina et al. 2017); Isolation FA, accounting for forest isolation and defined by the total amount of forest and poplar plantations cover in a 1-km buffer around the fragment (ha); Isolation NF, accounting for the archipelago effect (Simberloff and Abele 1982;Gotelli 2008) and defined by the number of forest and poplar plantation fragments in a 1-km buffer around the fragment; and Isolation NH, accounting for landscape permeability (Mortelliti et al. 2010) and defined by the number of continuous hedgerows in a 1-km buffer around the fragment. We jointly considered forests and poplar plantations to calculate Isolation FA and Isolation NF metrics because poplar plantation can play the role of "soft" matrix elements (sensu Lindenmayer and Fischer 2013) for forest birds (Chiatante et al. 2019). We considered the number of continuous hedgerows in the surrounding of fragments because evidence regarding the role of these landscape elements in increasing connectivity for forest birds exist (Bennett et al. 2004;Ernoult et al. 2006).
The gradient of fragment metrics was comparable for the two study periods (Fig. S2 in Online Resource 1). We selected five life-history and ecological traits commonly associated with bird extinction proneness or dispersal ability: Ecological group (interior/ generalist); Range size (km 2 ) (obtained from Birdlife International species factsheets BirdLife International 2021); MAR, minimum area requirement (ha), calculated as the minimum fragment area occupied by each species (see Tan et al. 2021); Body mass (g); Dispersal ratio [for each species: mean wing length (mm)/ cube root of mean body mass (g), Fischer and Lindenmayer 2005;Li et al. 2019]. Body mass and mean wing length were obtained from Storchová and Hořák (2018). Species-specific life-history and ecological traits are reported in Table S1 in Online Resource 1.
As regards to the fragment metrics, we considered the selective extinction as a driver of nestedness if the relationship between fragment rank and fragment log (Area) and/or Shape index was significant, while we considered the selective colonization as a driver of nestedness if the relationship between fragment rank and log (DNSA) and/or Isolation FA and/or Isolation NF and/or Isolation NH was significant. Similarly, as regards to the species life-history and ecological traits, the selective extinction was considered as a major driver of nestedness if the species rank was significantly related to the Ecological group (interior species are expected to be more extinction prone than generalist, Dondina et al. 2017) and/or the Range size (Dardanelli and Bellis 2021) and/or the MAR (Tan et al. 2021), while the selective colonization was the major driver if the species rank was significantly related to the Body mass and/or the Dispersal ratio (Schoener and Schoener 1984;Cook and Quinn 1995;Henle et al. 2004;Jenkins et al. 2007) and/or the Ecological group. The dispersal ability of a species, in fact, depends on the combined effect of the species' intrinsic movement ability (linked to body mass and dispersal ratio) and the species' forest specialization degree, with generalist species able to cross disturbed landscapes and interior species generally limited to well preserved woodlands and reluctant to cross heavily human-modified areas .
To investigate the possible drivers of nestedness, we ran a variable selection procedure on two full multiple linear regression models for each study period.
In the first full model, the fragments rank in the maximally packed matrix was the response variable and log (Area), Shape index, log (DNSA), Isolation FA, Isolation NF, and Isolation NH were the independent covariates. In the second full model, the species rank in the maximally packed matrix was the response variable and the Ecological group, Range size, MAR, Body mass, Dispersal ratio, and species Order (insert as a control variable) were the independent covariates. Before running the full models, we standardized all the continuous independent variables and checked for pairwise Pearson's correlation coefficient between covariates (correlation coefficient < 0.7; see Dormann et al. 2013). Variables' selection was carried out by iterative addition and deletion of variables following an Information-Theoretic Approach (Anderson et al. 2000(Anderson et al. , 2001Burnham and Anderson 2002) based on the Akaike Information Criterion (AIC) . We considered only the best model, i.e., the model with the lowest AIC. Finally, we checked the diagnostic plots of final models to assess residuals' distribution normality. Models' development and variables' selection were carried out through the 'MASS' package (Venables and Ripley 2002) in R.
To evaluate whether the nestedness of bird assemblages in forest fragments was caused by passive sampling, we used the Coleman's (1981) random placement model (Calmé and Desrochers 1999;Wang et al. 2012;Li et al. 2019). According to Coleman et al. (1982), the expected species richness in a fragment depends on the relative area of the fragment and the relative species abundances in the study area. The model was performed on the transpose of the abundance matrix with species in rows and fragments in columns. The model was run through the 'sars' package (Matthews et al. 2019) in R. Model fit was assessed through a diagnostic plot of fragment area (log transformed) against species richness, alongside the values predicted by the model (see Wang et al. 2010). The model should be rejected if more than a third of the observed data fall beyond one standard deviation from the expected curve (Wang et al. 2010;Matthews et al. 2019).

Changes of fragments occupancy
To assess whether changes of fragments occupancy by species occurred between the two study periods were driven by species ecological traits, we ran a variable selection procedure on a multiple linear regression model following an Information-Theoretic Approach and the AIC. The response variable in the full model was the number of position changes in the species rank in the maximally packed matrices obtained for the two study periods (the rank of the species in the first study period minus the rank of the species in the second study period), while the independent covariates were Ecological group, Range size, MAR, Body mass, Dispersal ratio, and species Order (continuous independent variables were standardized). We verified that the residuals of the final model had a normal distribution by checking diagnostic plots.

Results
From 1997 to 2001 (first study period), 1139 point counts were performed and 26 forest-dependent bird species (excluding edge species) were found (23% interior species, 77% generalist species, Table S1). The mean number of species detected at each point was 4.735 (SE = 0.071), while the mean number of species detected within each forest fragments was 5.091 (SE = 0.992). From 2013 to 2017 (second study period), 941 point counts were performed and 25 forest-dependent bird species (excluding edge species) were found (20% interior species, 80% generalist species, Table S1). The mean number of species detected at each point was 5.302 (SE = 0.080), while the mean number of species detected within each forest fragments was 5.917 (SE = 0.142). The species detected in the two study periods were the same, except for the lesser spotted woodpecker (Dryobates minor) not detected in the second period.

Nestedness drivers
As regards to the fragment metrics, the best model selected for the first study period showed an adjusted R 2 of 0.310 and revealed significant relationships between fragments rank in the maximally packed matrix and Shape index, log (DNSA), Isolation FA, and Isolation NF, and a marginally significant relationship between fragments rank and Isolation NH (Table 1; Fig. 2). Conversely, the best model selected for the second study period showed an adjusted R 2 of 0.165 and revealed significant relationships between fragments rank in the maximally packed matrix and Isolation FA and Isolation NH and a marginally Table 1 Results of the best models selected for the first (1997)(1998)(1999)(2000)(2001) and second (2013-2017) study period (analyses performed on bird assemblages in forest fragments, Lombardy, Italy) Response variable: fragments rank in the maximally packed matrix (study period I: N fragments = 844; study period II: N fragments = 543). Independent variables: Area, Shape index, DNSA (distance from the nearest source area), Isolation FA (area covered by forests and poplar plantations in a 1-km buffer around the fragment), Isolation NF (number of forest and poplar plantation fragments in a 1-km buffer around the fragment), Isolation NH (number of continuous hedgerows in a 1-km buffer around the fragment)  Fig. 2). As regards to the species life-history and ecological traits, the best model selected for both the first (adjusted R 2 = 0.499) and the second (adjusted R 2 = 0.513) study period revealed a significant relationship between species rank in the maximally packed matrix and the Ecological group [net of the effect of the control variable (species Order), Table 2; Fig. 3]. Moreover, the best model selected for the second study period showed a marginally significant relationship between species rank in the maximally packed matrix and the MAR [net of the effect of the control variable (species Order), Table 2; Fig. 3]. The random placement model showed that the nestedness of bird assemblages was not due to passive sampling in both study periods. Specifically, 84.95% and 88.40% of the observed data fall beyond one Fig. 2 Relationship between fragment rank in the maximally packed matrix resulting from nestedness analyses performed on bird assemblages in forest fragments (Lombardy, Italy) and fragment metrics: a log-transformed area, b Shape index, c log-transformed DNSA (distance from the nearest source area), d Isolation FA (area covered by forests and poplar plantations in a 1-km buffer around the fragment), e Isolation NF (number of forest and polar plantations fragments in a 1-km buffer around the fragment), and f Isolation NH (number of continuous hedgerows in a 1-km buffer around the fragment). In black: first study period (1997)(1998)(1999)(2000)(2001). In red: second study period (2013)(2014)(2015)(2016)(2017). Solid lines are fitted values plotted on the response scale. Dashed lines are 95% confidence intervals of the fitted values. Only relationships included in the best multiple regression model for each period (step AIC selection) are plotted standard deviation from the expected curve in the first and second study period, respectively (Fig. 4), rejecting the passive sampling hypothesis.

Changes of fragments occupancy
The species rank position changes observed between the two study periods ranged from a maximum relegation of 8 positions (− 8) for the European robin (Erithacus rubecula) and the common chiffchaff (Phylloscopus collybita), a small interior and generalist passerine species respectively, to a maximum advancement of 8 positions (+ 8) for the Eurasian sparrowhawk (Accipiter nisus), a large generalist raptor species (Table S2 in Online Resource 1). One of the six interior species detected in the first study period disappeared in the second period and other three interior species showed valuable relegations (European robin: position change in species rank = − 8, Eurasian wren Troglodytes troglodytes: position change in species rank = − 7, Eurasian nuthatch Sitta europaea: position change in species rank = − 4). Other two small generalist passerine species characterized by large MAR showed relevant position relegations (common chiffchaff: position change in species rank = − 8, Eurasian blue tit Cyanistes caeruleus: position change in species rank = − 4). Conversely, five generalist species characterized by large sizes (and therefore by a good dispersal ability) and a general large MAR showed valuable position advancements (Eurasian sparrowhawk: position change in species rank = + 8, and European honey buzzard Pernis apivorus, Eurasian jay Garrulus glandarius, common wood pigeon Columba palumbus, Eurasian green woodpecker Picus viridisc all with a position change in species rank = + 4).
The analysis developed to investigate if species ecological traits affected the number of position changes in the species rank selected a best model with an adjusted R 2 of 0.310 and revealed significant relationships between the number of position changes and the Ecological group and MAR (Table 3; Fig. 5).

Discussion
The present study was aimed at investigating the variability of nestedness and nestedness drivers in Table 2 Results of the best models selected for the first (1997)(1998)(1999)(2000)(2001) and second (2013-2017) study period (analyses performed on bird assemblages in forest fragments, Lombardy, Italy) Response variable: species rank in the maximally packed matrix (study period I: N species = 26; study period II: N species = 25). Independent variables: Ecological group (two level factor), Range size, MAR (minimum area requirement), Body mass, Dispersal ratio, Order (five-level factor). Ecological group (Generalist) and Order (Accipitriformes) were the reference category of the Ecological group and Order factor, respectively Few studies have investigated the nestedness variability, and the variability of its ecological drivers, over the long-term and results were contradictory (Sebastián-González et al. 2010). For instance, Heino et al. (2009) did not find any temporal differences of nestedness in a community of stream insects. Azeria et al. (2006) found consistency in nestedness pattern in bird assemblages over a longtime interval, but some noise in the system (changes in occupancy pattern) were documented. Conversely, Azeria and Kolasa (2008) and Sebastián-González et al. (2010) found significant annual changes in nestedness and its causes in invertebrate and bird assemblages, respectively.

Fig. 3
Relationship between species rank in the maximally packed matrix resulting from nestedness analyses performed on bird assemblages in forest fragments (Lombardy, Italy) in the two study periods (1997-2001 and 2013-2017) and species life-history and ecological traits. a Species ecological group (light grey = generalists, dark grey = interior species), and b species minimum area requirement (MAR). Solid lines are fitted values plotted on the response scale. Dashed lines are 95% confidence intervals of the fitted values. Only significant (p < 0.05) and marginally significant (p < 0.1) relationships included in the best multiple regression model for each period (step AIC selection) are plotted. **p < 0.01 Fig. 4 Comparison of observed data and expected values under the random placement model in two study periods (1997-2001 and 2013-2017) (analyses performed on bird assemblages in forest fragments, Lombardy, Italy). Expected values (black line) and associated standard deviations (± 1 SD; gray lines) are shown. Black dots represent observed species richness Table 3 Results of the model selection performed to investigate if species ecological traits drive the number of position changes in the species rank in the maximally packed matrices obtained for the two study periods (1997-2001 and 2013-2017) (analyses performed on bird assemblages in forest fragments, Lombardy, Italy) Only results of the best model are shown. Response variable: number of position changes in the species rank in the maximally packed matrices obtained for the two study periods. Our results showed that, overall, the nestedness of bird assemblages did not significantly change during the studied time interval, but its drivers have considerably varied over time.
In the first study period (1997)(1998)(1999)(2000)(2001), sites did not result significantly nested. However, fragments rank in the maximally packed matrix resulted significantly associated to fragment metrics revealing a not negligible action of both selective extinction and colonization. The key role of extinction in shaping community assembly has been documented in many studies (Dardanelli and Bellis 2021;Ruzzier et al. 2021). A reasonable trigger for extinctiondriven spatial patterns is the faunal relaxation process, typical of extremely fragmented forest landscapes (Brooks et al. 1999;Ferraz et al. 2007). Our results suggested that, in the investigated landscape, interior species that prefer regularly shaped fragments had the highest extinction vulnerability and disappeared first from forest fragments. The edge effect was thus the main process driving the selective extinction in the first study period, supporting the hypothesis that the effect of fragment size on bird communities may generally be mediated through edge effect (Banks-Leite et al. 2010).
In the first study period, selective colonization seems to be another not negligible driver of the spatial pattern of the investigated bird community. Models revealed that interior species disappeared first from forest fragments far from source areas and highly isolated from surrounding fragments, because of the lack of nearby forest or poplar plantation fragments and, secondly, of landscape connectivity elements such as hedgerows. It is not surprising that colonization shapes the spatial patterns of a bird community. In fact, assemblages of organisms with strong dispersal abilities, like birds, often exhibit a colonization-driven nested pattern as they can easily move between neighboring fragments (Sebastián-González et al. 2010). However, ecological traits generally associated with dispersal ability such as body mass and dispersal ratio were not significantly linked to species rank. As suggested by other authors (Dardanelli and Bellis 2021), dispersal ratio and body mass were probably not appropriate indicators of dispersal ability for birds in the study of nestedness drivers.
In the second study period, sites resulted significantly nested. The increase of urbanization and agriculture intensification and the decrease of the number of connectivity elements, such as continuous hedgerows, led to a general increase of the matrix hostility in the investigated landscape during the considered time interval. The decrease of matrix permeability probably led to a tightening of the isolation conditions between forest fragments, making the landscape more similar to an islands-ocean system and strengthening the effect of fragment metrics on bird communities . Specifically, the effect of the extinction process is evident in the second study period. Models revealed that interior species with large MAR were more extinction prone in small fragments. Together with the results obtained for the first study periods, this evidence confirmed the importance to conserve large and regularly shaped Fig. 5 Relationship between the number of position changes in species rank in the maximally packed matrices obtained for the two study periods (the rank of the species in the first study period minus the rank of the species in the second study period) and species life-history and ecological traits (analyses performed on bird assemblages in forest fragments, Lombardy, Italy). a Species ecological group (light grey = generalists, dark grey = interior species), and b species minimum area requirement (MAR). Solid lines are the fitted values plotted on the response scale. Dashed lines are the 95% confidence intervals of the fitted values. Only significant (p < 0.05) relationships included in the best multiple regression model (step AIC selection) are plotted. *p < 0.05 fragments to safeguard viable local populations of interior species sensitive to fragment area.
Conversely, only local isolation conditions drove the selective colonization in the second study period, while the effect of distance from source areas was completely lost. This evidence can be partially explained considering the results of the analysis of fragments occupancy changes occurred between the two study periods. Specifically, 67% of the interior species detected in the first study period disappeared or showed considerable position relegations during the investigated time interval. This result is aligned with other studies (e.g., Massimino et al. 2015;Gregory et al. 2019) that discovered a general population decline of specialist forest species in Europe. In our study area, according to the changes of nestedness drivers occurred during the considered time interval, these species probably disappeared from small and locally isolated fragments regardless of whether they were near or far from one or more source areas. In fact, the increase of matrix hostility occurred in the investigated landscape during the study period could had made even short distances from source areas impossible to face by internal species typically reluctant to cross highly disturbed areas. The interruption of the dispersal flow from the source areas probably caused the local extinction of internal species characterized by high MAR within small fragments relatively near to the source areas, and it can also be the cause of why the effect of fragment area is evident only in the second study period. Contrariwise, the importance of local isolation conditions in driving selective colonization in the second study period suggested that interior species persist when the arrangement of forest fragments, poplar cultivations and hedgerows form archipelago conditions where local populations support each other through the process of the "internal colonization" (Simberloff and Abele 1982;Gotelli 2008).
The disappearance of interior species from forest fragments relatively close to source areas may have been due, as well as to the increase of matrix hostility, even to alteration of the internal characteristics of woodlands due to management actions (not directly considered in this study). Structural and floristic woodland alterations are known to have a strong impact on the occurrence of small interior species characterized by strict ecological requirements (Dondina et al. 2015) and their effect can be evident with consistent delay due to the phenomenon known as "extinction debt" (Lindenmayer and Fischer 2013).
Another possible cause of the disappearance of internal species from small fragments where resources are typically limited, could have been the competitive exclusion due to the arrival of large generalist species. In fact, multiple generalist species characterized by large sizes (and therefore by a good dispersal ability) and a general large MAR showed valuable position advancements between the two study periods. We expect that demographic changes played a leading role in driving the observed change of species distribution. In fact, generalist forest species showed stable or increasing population trends in Europe during the very last decades (Gregory et al. 2019). Specifically, a significant positive population trend for the period 1992-2019 has been recently detected for both the Eurasian jay, the common wood pigeon, and the Eurasian green woodpecker in Lombardy (Tirozzi et al. 2021). When species show positive demographic trends, they often tend to enlarge their range by colonizing sub-optimal habitats, such as small forest fragments. The arrival of large generalist species in small fragments could have also caused the competitive exclusion of smaller generalist species, explaining the position relegation observed for some of these species, namely the common chiffchaff and the Eurasian blue tit. We can expect that, if competitive exclusion among generalist species were to continue in the future, it could generate a species checkerboard pattern (Heino et al. 2009) and a general decrease of the degree of nestedness within the bird community in the study area (Gotelli and McCabe 2002).

Conclusion
In conclusion, our study demonstrated that significant changes of nestedness drivers and species occupancy within fragments can occur over the long-time. Caution must thus be posed when community analyses are performed in a single snapshot of time or over short-time periods. Only long-term study, in fact, can detect community changes due to environmental and/or demographic modifications and possible extinction delay due to the phenomenon of "extinction debt". From a management point of view, the detected importance of selective extinction as a nestedness driver confirmed the importance to conserve large and regularly shaped fragments to safeguard viable local populations of interior species in the long-term. Additionally, our study pointed out the importance of maintaining and/ or enhancing the permeability of the matrix to safeguard the landscape connectivity that sustain dispersal flows from source areas even for species characterized by a low dispersal ability. In general, the protection of residual forest fragments should be integrated in functional ecological networks (composed of source areas, archipelagos of both large and regularly shaped and small and irregularly shaped fragments, and a "soft" matrix) aimed to conserve viable and resilient forest bird metapopulations in highly exploited and dynamic landscapes.