Evaluating temporal turnover in avian species richness in a Mediterranean semiarid region: Different responses to elevation and forest cover

When studying the effects of global change on biodiversity, it is far more common for the effects of climate change and land‐use changes to be assessed separately rather than jointly. However, the effects of land‐use changes in recent decades on species richness in areas affected by climate change have been less studied. We assess the temporal turnover in species richness of an avian community between a historical period and a modern one as a consequence of global change.


| INTRODUC TI ON
Understanding long-term changes in distribution patterns of species caused by global change is one of the most pressing tasks in ecology and biodiversity conservation (Jetz et al., 2007;Powers & Jetz, 2019).Climate change causes range shifts towards cooler temperatures, which may lead to shifts along the latitude or uphill shifts along elevation (i.e.mountaintop extinctions and lowland range expansions (Freeman et al., 2018;Tayleur et al., 2015)).In addition to climate change, global change is caused by land-use change directly generated by the intensification or abandonment of human activities (Newbold, 2018;Sage, 2020;Sirami et al., 2017).The effects of climate change in areas experiencing land-use change may generate different species distribution patterns, which may have a potential effect on biodiversity although such interactive effects have been less studied (Newbold, 2018), and potentially having even greater impacts on biodiversity in combination than alone (Cabral et al., 2023;Neff et al., 2022).
Therefore, when assessing the effects of climate change in areas undergoing land-use changes, altitude may be a proxy for changes in temperature and precipitation due to climate change (Freeman et al., 2018), although species and their associated vegetation types may change their elevation ranges synchronously (when both change in the same direction) or not (when they change in different directions; Santos et al., 2015).These patterns may generate different range shifts and should be evaluated for species conservation.
Forest ecosystems host a large part of the biodiversity worldwide, and they have been heavily modified in many countries mainly due to two processes with opposing effects: agricultural abandonment and deforestation (Foley et al., 2005).On the one hand, the agricultural abandonment has caused an expansion and regrowth of the natural areas due to succession, increasing the agroforestry and scrubland areas (Ameztegui et al., 2016;Whytock et al., 2018), whereas traditional farmlands and the species inhabiting them are vanishing.Therefore, these processes of passive habitat rewilding (i.e.natural revegetation) in agricultural areas may generate range shifts in species associated with agroecosystems, which are suffering from widespread declines at the European level (Eglington & Pearce-Higgins, 2012;Gaüzère et al., 2020;Reino et al., 2018).On the other hand, deforestation has caused habitat loss and fragmentation for forest species (Laurance et al., 2000), which may generate a landscape barrier for range shift of communities affected by climate change (Marjakangas et al., 2023).Therefore, the study of the range shifts of animal species in forest areas undergoing landscape transformations (agricultural abandonment vs. deforestation) and located along elevational gradients may fill a priority knowledge gap.Studied jointly, these effects can help understand the impact of global change (Chen et al., 2011;Newbold, 2018;Whytock et al., 2018).
The effects of global change have been evaluated at different spatial scales and extents, from global (Jetz et al., 2007), through continental (Gregory et al., 2009;Princé & Zuckerberg, 2015), country-level (Reino et al., 2018), regional (Clavero et al., 2011) and down to the local scale (Wilbanks & Kates, 1999).Evaluation at the regional scale allows to identify in greater detail changes in the landscape and movements along altitudinal gradients of species (Palacio et al., 2020).Moreover, biodiversity responses to habitat creation was found to depend on local-and landscape-scale factors that interact across time and space (Whytock et al., 2018), so there is an increasing interest in understanding how land-use change is changing the local diversity, since land-use change is an important determinant of the functioning of ecosystems (Newbold, 2018;Newbold et al., 2020).Regarding the temporal scale of global change effects, evidence is mostly based on projections and modelling assuming different global change scenarios in a future context (Reino et al., 2018;Titeux et al., 2017;Triviño et al., 2011).Only a handful of studies have been able to evaluate or compare historical data with modern data over a long time interval (e.g.Daskalova et al., 2020;Freeman et al., 2018;Marjakangas et al., 2023;Moritz et al., 2008;Palacio et al., 2020).
Such studies are particularly valuable since they provide 'hard' information about the temporal trajectories of ecosystems rather than probabilistic forecasts only, and allow actually testing for the presence of trends as well as investigating causal hypotheses about the mechanisms underlying the observed changes (Cheng et al., 2019).In order to address this, it is necessary to have and maintain historical data with which to compare modern data, as well as maintaining the same field design for both periods in order to make an accurate comparison (Tingley & Beissinger, 2013;Yackulic & Ginsberg, 2016).Generally, the studies to understand biodiversity changes have assessed selected species, usually threatened ones or species of conservation concern (Gómez-Catasús et al., 2018;Moreno-Zarate et al., 2020), as well as specific taxa (Cheng et al., 2019;Eglington & Pearce-Higgins, 2012).
However, analyses of entire communities are much rarer, due to the difficulty of conducting fieldwork that samples all the species of the studied taxonomic groups (Freeman et al., 2018;Palacio et al., 2020;Santos et al., 2015).These patterns are complex to analyse, so the estimate of the species richness can be a measure of ecosystem functioning and a key parameter to study global change effects on communities (Newbold, 2018).Nevertheless, the temporal dynamics of biodiversity can be decoupled from rates of change, so it must be complemented with other measures of temporal turnover such as the richness-based species-exchange ratio (hereafter species turnover rate; Hillebrand et al., 2018).
bird community, climate change, global change, historical data, modern data, multispecies occupancy models, richness-based species-exchange ratio, species richness Semiarid Mediterranean ecosystems represent an ideal scenario to study the effects of global change, where different species are found at the limit of their Eurasian distribution and others are endemic to these areas (Esteve et al., 2015).These ecosystems at the edge of their distribution limit allow for surprising changes, with species appearing and disappearing more frequently.Moreover, this area offers the possibility to assess community dynamics in relation to biogeographical limits (the aridity margin) of the biome's distribution, which are areas considered to be highly sensitive to climate change (Guiot & Cramer, 2016).In addition, the Mediterranean region has many mountains with considerable altitudinal gradients, which allows us to evaluate the movement of species along this gradient (Clavero et al., 2011) and the consequent changes in local species richness (Maphisa et al., 2019).Finally, the region is located in a climatic gradient generating a biogeographical ecotone between the Mediterranean conifer-sclerophyllous-broadleaf forest (hereafter, pine forest) biome and the Mediterranean scrub (hereafter, scrubland) biome (Esteve et al., 2013;Walter, 1981).These Mediterranean ecosystems are suffering both habitat revegetation due to agricultural abandonment and anthropization processes in the last decades with legacy effects on biodiversity (Jiménez-Franco et al., 2022), so this area represents an ideal scenario to evaluate the influence of land-use change and its effects on species richness (Clavero et al., 2011;Newbold et al., 2020).
In this study, we assess the temporal turnover in avian species richness in a semiarid Mediterranean region at a regional scale (11,317 km 2 ) between a historical period (1991)(1992) and a modern one (2012 and 2017) as a consequence of global change.We used bird community surveys for this purpose, since bird species play a vital role in key ecological processes, including seed dispersal, pollination, pest control, nutrient cycling and soil formation (Cannon et al., 2019) and have been used as bioindicators in multiple scenarios of global change (Gregory et al., 2009).Specifically, our aims are: (1) to analyse the relationships between changes in the bird species richness and key environmental variables (elevation and forest cover) along a comparison time span of 21-25 years; and (2) to evaluate the temporal turnover in the avian species richness by mapping species richness for each period, as well as the species turnover rate between the two periods.We used bird community survey data from the historical period considering multiple 1-km transect bird surveys as spatial replicates for larger cells interpreted as 'sites' (Jiménez-Franco et al., 2019), and resurveying them in modern time using the same protocol.This data set allows us to fit two hierarchical multispecies occupancy models for each period, obtaining species-specific estimates of occupancy probability in relation to environmental covariates that are interpreted as surrogates for global change effects (Clement et al., 2016), as well as estimates of average occupancy probability of the entire bird community (Zipkin et al., 2009).We obtain the turnover in species richness between two periods as the difference in species distribution (Kéry & Royle, 2016).Although this application of community occupancy models between different temporal periods has already been explored (Tingley & Beissinger, 2013), the use of spatial replicates with occupancy models is a suitable option to implement this modelling framework (Jiménez-Franco et al., 2019;Karanth et al., 2011;Martínez-Martí et al., 2016;Srivathsa et al., 2018).Numerous studies have demonstrated that detection probability commonly varies among species, over time and among ecological covariates, and there may be seriously biased inferences when this variability is ignored (Guillera-Arroita, 2017; Kellner & Swihart, 2014;Zuberogoitia et al., 2020).The design of this study should be able to identify patterns of change at the community level between different periods that are robust to patterns of imperfect detection.Specifically, we expect that the bird community can have different responses due to the joint evaluation of drivers of global change, the climate change and the land-use change.On the one hand, species distributions may be displaced to higher elevations because of climate change and pine forest decay in lower areas, leading to greater species richness there.On the other hand, we also expect to observe that the climate change effects may be weaker than other drivers of global change, such as the land-use changes.For example, coastal areas (with lower elevations) may increase their bird species richness due to the generation of different secondary habitats (e.g.green areas from urbanizations, passive habitat rewilding due to agricultural abandonment), which are potential for forest bird species.

| Study area and species
This study was carried out in the region of Murcia (SE Spain; Figure 1), an area of 11,317 km 2 that has a semiarid Mediterranean climate, and where during the last five decades  the annual temperature averaged 16.7°C, with a tendency to increase by 0.135°C per decade (Garrido et al., 2015).There is a wide climatic gradient in the study area which allows it to host several types of ecosystems that represent an ecotone between Mediterranean and arid subtropical: semi-desert areas, Mediterranean scrub and coniferous forest (Esteve et al., 2015).In compliance with the European Birds and Habitats Directives, 22 Special Protection Areas for birds have been designated in the region (Abellán et al., 2011).Nowadays, as a result of the climatic gradient and the different ecosystems, the breeding bird community of the study area consists of around 120 species (Calvo et al., 2017).
This study focuses on bird species inhabiting Mediterranean forest ecosystems in the region of Murcia, where forests are dominated by one tree species, the Aleppo pine (Pinus halepensis), a conifer that may reach up to 22 m in Mediterranean areas (Mitsopoulos & Dimitrakopoulos, 2013), with open areas made up of Mediterranean shrub of the 'garrigue' type.A total of 73 avian species were recorded during the historical surveys (Table S1), most of them passerines, and the families best represented being Sylviidae, Turdidae, Fringillidae and Paridae.These four families encompass more than a third of the forest bird community (31 species), which indicates different states of forest maturity and include representatives of relevant trophic and functional guilds (granivores, frugivores, seed dispersers, etc.), as well as of different zoogeographical origins, for example, boreal versus Mediterranean species (Blondel et al., 2010).For the modern period, we recorded a total of 77 species, which are the same as those in the historical period, but with four species typical for the forested areas being newly present: Bucanetes githagineus, Regulus ignicapilla, Streptopelia decaocto and Sylvia atricapilla.'Forest birds' is used here in a broad sense, including species with preferences ranging from strictly wooded areas to scrubland of low height and/ or cover.

| Field work
An extensive monitoring programme was conducted in the entire region of Murcia with the aim of characterizing breeding bird communities and their temporal turnover in species richness.This field work was developed as a basis for assessing the state of the region's forest heritage and the effects of global change on its biodiversity (Esteve, 1991;Esteve et al., 2015), concentrated in two periods: a 'historical period' (year 1991 and 1992) and a 'modern period' (year 2012 and 2017), that is, 21-25 years of time span (calculated as the minimum-maximum number of years between historical and modern data; Moritz et al., 2008).The survey protocol was characterized by transects of 1 km length, which were distributed across all the forested areas in the study area in a random fashion and were surveyed once during a breeding period (May-July; Hernández & Barberá, 1997).Each transect survey was conducted by walking and recording the number of each species detected (by sight or song/ call) along it (Tellería, 1986).A total of 377 1-km transects were conducted during the historical period as part of a forestry plan in the region, mainly developed by V. Hernández, G. G. Barberá, A.
Giménez and M. A. Esteve (Hernández & Barberá, 1997).Then, we conducted a new avian survey of 198 transects with a similar protocol during the modern period, developing each transect individually by M. León-Ortega or J. Martínez-Ródenas.The locations of the transects between both periods have different spatial coordinates, due to the different sampling effort and their random distribution within the forest systems studied, also responding to long-term forest habitat change between periods (Figure 1).

| Resurvey field protocol and environmental covariates
We considered the 1-km transects of each bird survey as nested within a 'site', that is, as spatial observations replicated within the site, thus enabling estimation of detection probability in an occupancy model (Hines et al., 2010;Srivathsa et al., 2018).We F I G U R E 1 (a) Location of the study area (Region of Murcia, SE Spain), the distribution of 1-Km spatial transects for the bird surveys in the historical (1990,1991) and modern period (2012,2017), as well as the environmental covariates: (b) elevation, (c, d) percentage of forest cover for the two time periods.defined sites as cells in a grid that covers our study area (Jiménez-Franco et al., 2019).We selected cell size to be 2 × 2 km, based on the species in our bird community (most of them passerines with small home ranges; Jiménez-Franco et al., 2019;Rechetelo et al., 2016), which may be suitable for evaluating the effects of environmental covariates on bird species occupancy at a regional scale (Kéry et al., 2013;Lipsey et al., 2017).Therefore, a total of 377 and 198 1-km spatial transects were conducted and grouped in 226 and 139 sites for the historical and modern period respectively (Table S2).
Elevation of the study area was obtained from a digital elevation model map (http:// centr odede scarg as.cnig.es) (mean (minmax) = 503.2(1-1847) m).The percentage of forest cover (FOREST) was estimated from the CORINE Land Cover map (http:// centr odede scarg as.cnig.es) for the year 1990 and 2012 for the historical and modern period (we considered no relevant changes during the same sampling period), respectively, (mean = 11 and 15% for the historical and actual period respectively).All maps had an original resolution of 100 × 100 m and were rescaled to a grid size of 2 × 2 km to match the spatial resolution used for our bird surveys.GIS analyses were carried out with the raster package (Hijmans, 2016) in R 3.6.2(R Core Team, 2019).

| Statistical analyses
We fitted two Bayesian multispecies occupancy models (Dorazio et al., 2006;Dorazio & Royle, 2005), one for each period (historical and modern) to estimate species richness and occurrence correcting for false-negative errors and to evaluate the influence of elevation and forest cover on bird richness and occurrence.These models are an extension of the single species site occupancy model (MacKenzie et al., 2002;Tyre et al., 2003), whereby a hierarchical structure combines community-and species-level attributes within a single analytical framework based on spatial replicates (Jiménez-Franco et al., 2019).The model for occurrence is specified as z(i, k) ∼ Bern i,k where i,k is the probability that species k occurs at site i.The true occurrence is imperfectly observed, and we define the detection model for the detection frequency Y i,k for species k at site i as Y(i, k) ∼ Binomial p ik ⋅ z ik , J(i) , where p i,k is the detection probability of species k for the jth spatial replicate at site i, given that species k is in fact present at site i, and J(i) is the number of spatial replicates (i.e.transects) in cell (site) i.We incorporated site-specific environmental characteristics into probability of occupancy by using covariates (Kéry & Royle, 2009, 2016;Zipkin et al., 2009).Linear and quadratic effects of elevation (ELEV) and percentage forest cover (FOREST) were included.Both covariables were standardized by centring and division of the result by the SD (bird dataset can be found in Jiménez-Franco et al., 2023).
In our model, we defined the probability of occupancy as follow: The inverse-logit of 0,k is the expected occurrence probability for species k at a site with 'average' habitat characteristics (i.e. with zero values of all normalized covariates).The coefficients from 1,k to 4,k are the effects of the elevation (linear and squared) and of the percentage of forest cover (linear and squared), for species k respectively.All species-specific parameters (i.e.0,k , 1,k , 2,k , 3,k , 4,k ) are modelled hierarchically as draws from separate normal distributions with hypermeans and hypervariances that are estimated from the data.These hyperparameters can be interpreted as describing the average species and the heterogeneity among species in the response of occupancy to these covariates (Kéry & Royle, 2016).
We assumed that detection probabilities varied depending on the species but were not influenced by survey characteristics: logit p k ∼ Normal lp ,2 lp .The two models for each period were fitted using JAGS (Plummer, 2003), run from R 3.6.2(R Core Team, 2019) through package jagsUI (Kellner, 2015), using vague priors, three chains, 15,000 iterations and a burn-in of 5000 iterations, with a thin rate of 2 (see R and JAGS code in Appendix S1).
Convergence was assessed by examining the Rhat values for each parameter estimate (Brooks & Gelman, 1998).We present posterior means and standard deviations for point estimates and the Bayesian analogue to a standard error.We also computed species richness maps for each period considering the environmental covariates of models (elevation and forest cover) for all the study area, as well as the richness difference between both periods.Moreover, we complemented these richness estimates with a measure of temporal turnover, the richness-based species-exchange ratio (turnover rate), SERr (Hillebrand et al., 2018).
This index is a measure of the proportional exchange of species between both periods, that is, the sum of immigrations and extinctions as a fraction of the total number of species, which is the complement of Jaccard's similarity index, a commonly used metric in biodiversity change studies (Hillebrand et al., 2018).The SERr is quantified as: SERr = (S imm + S ext )/S tot , where S imm is the number of species immigrating (newly recorded from the previous period), S ext is the number of species extinct (lost from the previous period) and S tot is the total number of species across both samples.Such a presence-absence-based SER quantifies the gross change in species composition (S imm + S ext ) rather than the net change (Δrichness = S imm − S ext ) on a closed scale between 0 and 1, where 0 means all species persist and 1 means all species change.Maps were produced with the raster package in R.Moreover, we compare the models for both periods by evaluating the estimates of detectability and occupancy for each species, as well as the estimates of regression coefficients for predictor variables.A linear correlation was performed among species estimates of detectability and occupancy, comparing both time periods.Since sites of field surveys were in general different, we do not analyse the differences for sites specifically.
The mean observed avian richness was 15.3 species (14.6-16.1,95% CI) for the historical period and 10.9 species (10.3-11.4,95% CI) for the modern period.However, the estimated mean avian richness was much higher than observed, being 28.2 (27.4-28.9;mean, 95% CRI) and 26.7 (26.1-27.3)for the historical and modern period respectively (Figure 2).The mean of observed and estimated species richness for each cell site is shown along the gradient of the environmental covariates for all sites (Figures S1 and 3 respectively).
Regarding the bird richness in relation to elevation, the estimates were higher for the intermediate levels of elevation in the historical period (Figure 3).However, avian species richness had a homogeneous trend for the modern period, being higher than the historical period in both lower and higher elevations, and lower in medium elevations (Figure 3).Moreover, species richness estimates were slightly higher for intermediate values of forest cover (20-60%), following the same trend for both periods (Figure 3).
The modelled maps for avian species richness showed different patterns for each period (Figure 4), with species richness being high mainly in the mountain ranges in the historical period, while in the modern period, the species richness increased both in low and high elevational gradients.The highest temporal turnover with an increase in avian richness was found in the northeastern area and southeastern area (coastal areas; Figure 4), coinciding both with the highest and lowest elevations of the study area (Figure 1b).
Moreover, the richness-based species-exchange ratio (SERr), used as a measure of temporal turnover rate, had values higher than 0.5 for the study area, being higher in the southeastern area (coastal sections) and lower in the northeastern areas (Figure 5).
The estimates of community occupancy in relation to environmental covariates showed similar patterns to avian species richness for each period (Figure S2).However, the community mean occupancy probability showed an optimum at a medium elevation in the historical period, but this parameter had a monotonous relationship with elevation for the modern period (Figure S2).Moreover, mean occupancy probability was slightly higher with a 20-40% of forest cover for both periods (Figure S2).The results of community occupancy models for each species showed that estimates of occupancy probability varied greatly among them (Figure S3).These estimates of occupancy were weakly correlated between the two time periods (r = .68,p < .001),as opposed to the naïve estimate, which was highly correlated between both periods (r = .85,p < .001).Similar results were shown for detectability, where the estimates of detectability varied among species and the values were weakly correlated between both periods (Figure S3; r = .67,p < .001).Fringilla coelebs was the species with the highest increase in occupancy (Figure S4; change of mean probability occupancy of 0.66, Table S1), and Luscinia megarhynchos was the species with the highest decrease in occupancy (Figure S4; change of mean probability occupancy of −0.66, Table S1).Moreover, Petronia petronia and Coturnix coturnix were the species with a highest increase and decrease in detectability, respectively (change of mean probability of detection of 0.46 and −0.48, respectively; Table S1), presumably an effect of changes in their mean abundance.Species-specific predictions of occupancy probability in relation to elevation and forest cover showed different trends.They were fairly variable between the historical and the modern period and among species (Figure S5).

| Species richness estimates by using hierarchical models
Our results showed that both the observed and estimated bird species richness decreased during the last 21-25 years of time span for the surveyed study area.These results are consistent with the fact that more bird species decreased in occupancy rather than increased.
Furthermore, comparing the decrease in observed bird species richness between the two periods, we obtained a 28% decrease, while there is an 8% decrease in estimated richness when accounting for detectability.This difference in the species richness reduction between the two periods was smaller when detectability is taken into account.That is, the difference between observed and estimated species richness was greater in the modern period.This result suggests that species are less frequent with respect to their historical distribution and are expanding into their new areas.Therefore, sampling units with low densities (due to rarity or recent expansion) have increased, and for approximately the same sampling effort per transect, non-detections will be more frequent in the modern period, and consequently mean observed richness will be more underestimated.After accounting for detectability, our results showed that the estimated species richness increased proportionally more in the modern period, although this species richness was still lower than in the historical period.We highlight that in changing conditions such as those in our study, hierarchical models provide a better estimate F I G U R E 2 Bird species richness in the historical and modern period, considering both the observed in survey data and estimated with the multispecies occupancy models for each time period.

| Temporal turnover in species richness
The decrease in bird species richness is spatially heterogeneous in the study area, affecting intermediate altitudes, which are the majority among the forest systems studied.However, at both ends of the elevational gradient, species richness increased slightly between the two time periods (Figure 3a).This spatial pattern of species richness increase corresponds to two distinct zones of the study area.
First, the species richness increased with elevation within the most elevated sections (NW in the map of Figure 1), which may be related to climate change (Freeman et al., 2018;Santos et al., 2015;Figure 1c,d).Second, we also observed an increase in bird species richness in the southeastern part of the study area, although since it is associated with lowland areas, there must be more complex anthropic habitats (Sirami et al., 2017).In other words, the increase in intensive agriculture (and specifically irrigated arboreal crops), with higher penetration in forest areas (Martínez & Esteve, 2002), as well as of urban gardens with trees, can generate habitats structurally similar to pinewoods, while providing additional water resources and food in peripheral areas to forests that attract forest bird species (Alderson & Sander, 2022;Pagaldai et al., 2021;Zamora-Marín et al., 2021).The forest interface becomes more diffuse, creating  (Gaüzère et al., 2020;Newbold, 2018).
Our results also showed a very high degree of species replacement (a high SERr index) in general for the whole study area, but with two opposite patterns of increase in species richness, associated with the two previously indicated syndromes and their drivers.
On one hand, there is an increase in species richness at high elevations (northwest of the study area) with a lower species turnover rate, which suggests an addition of new species with less replacement in the bird community, due to the effects of climate change in forested areas (Gil-Tena et al., 2007).On the other hand, there is an increase in species richness at low elevations (southeast of the study area) with a higher species turnover rate.This second rates of species turnover (greater in the lower areas).Therefore, the changes in species richness did not correspond to the rates of species turnover in the study area, especially in the southeastern part of the study area.This result supports the argument that biodiversity changes have a complex coupling with species richness trends (Hillebrand et al., 2018), which may instead be influenced by different global change drivers.Moreover, range shift of avian community in a regional spatial scale may also suffer the effects of landscape fragmentation and landscape barriers that decrease the species movements (Marjakangas et al., 2023).

| Forest cover and species richness
Our results showed a homogeneous and similar pattern of bird species richness in relation to forest cover for both periods, although with recent average values lower than those present in the previous historical period.Although our results showed that forest cover can be assimilated in practice rather to a 'control variable', we consider that attention should be paid to it in more detailed studies, since the study area has undergone land-use changes that can directly affect forest areas.For example, agricultural abandonment and natural succession leading to passive habitat recovery that can increase pine forest areas, replacing scrublands and agricultural areas (Queiroz et al., 2014).In parts of our study area this increase in forest cover responds to the replacement of scrubland (covering gaps within forested areas) and to the spontaneous revegetation processes occurring in abandoned agricultural areas, all benefiting forest bird species (Bowen et al., 2007;Gil-Tena et al., 2007).In addition, 13% of the forests are in decline, with severe decay, and as a result they have opened in clearings, increasing their internal heterogeneity (Esteve-Selma et al., 2003).This indicates that the closed forest has been patched with areas of scrub, as an indication of climatic imbalance.historical (1990-1991) and modern (2012-2017).This index quantifies the gross change in species composition (S imm + S ext ) rather than the net change (Δrichness = S imm − S ext ) on a closed scale between 0 and 1, where 0 means all species persist and 1 means all species are exchanged.
All this means that this variable (forest cover) may lose weight in the distribution of species.Previous studies have shown that changes in land-use affect species richness (Palacio et al., 2020), with forest bird populations, and hence species richness generally increasing, while they decrease in agricultural and scrubland environments (Cabodevilla et al., 2022;López Bermúdez et al., 2016).In this line, forest expansion due to agricultural abandonment (Lasanta- Martínez et al., 2005) may be beneficial for forest specialist species (Reino et al., 2018) by providing food resources and structural cover, although species richness may also be influenced by species functional traits (Lisón et al., 2022).

| Species-specific changes
Species-specific model estimates showed that forest birds species responded differently over time, that is, some species increased the probability of occupancy between the historical and modern periods, while other species decreased in that probability between these same periods (Figure S4).However, in the net balance there is a dominance of species that show a decrease in probability of occupancy.Furthermore, the effects of elevation and forest cover were different for each species and also for each period (Figure S5).
It should be highlighted that the effects on the bird community are complex and more detailed studies are needed to interpret the decline of particular species.Note that due to these complex patterns we performed an additive model rather than an interaction model, since our study covariates are quantitative and the interaction variables might be very hard to interpret at both the species and community level.Abundance and distribution studies for each species could clarify trends over time.We suggest that future studies should thoroughly evaluate the species-specific changes in areas affected by the combined effects of the main drivers of global change in order to elucidate their specific distributions and conservation measures.

| CON CLUS IONS
First, we must highlight the relevance and opportunity of the designs that consider detectability, in the study of situations of great dynamism in the assemblages of bird communities.In these periods of change it is more important than ever to consider detectability for occupancy estimates and consequently species richness, and failure to do so can result in biased estimates and erroneous conclusions.
Our study design has allowed us to implement a multispecies framework in which to compare past bird surveys with modern surveys conducted under equivalent designs, and also to estimate parameters that govern changes in species presence/absence in relation to environmental variables.Our results explored the temporal turnover of species richness of the regional forest bird ensemble over a period of species richness in relation to the observed data.Our results suggest that in these periods of change, it is more important to consider detectability for estimates of occupancy and consequently species richness(Devarajan et al., 2020).Failure to do so can result in biased estimation and erroneous conclusions(Kellner & Swihart, 2014;Tingley & Beissinger, 2013).Numerous studies have demonstrated that detection varies among species, over time and among habitats, and there may be serious consequences when this variability is ignored (Guillera-Arroita, 2017).For example, failure to correct for imperfect detection may result in bias in estimated relationships with ecological covariates, estimates of species distribution or abundance that are inaccurate or mask trends, and improper selection of indicator species(Banks-Leite et al., 2014).It is also important to highlight that in order to compare the temporal turnover from historical data, and to compare occupancy and detectability among years, we should consider the field design established, so future surveys must be conducted following a protocol similar to the previous one(Jiménez-Franco et al., 2019;Yackulic & Ginsberg, 2016).
drivers involved: it could be due to less exposure to climate change at lower elevations (temperatures on the Mediterranean coast are tempered by the presence of the sea, thus moderating temperature changes) and a higher rate of land-use change with the creation of F I G U R E 3 Relationships between the number of bird species (estimated richness, Nsite) and the covariates of each sampled site: (a) elevation; (b) scrub-forest land cover (%), that represents a gradient with the two land uses, where 0% is the scrub land use up to 100% representing the pine forest land use.Each point represents the richness of each cell site surveyed (n = 226 and n = 139 for the historical and modern bird survey respectively).Lines represent smoothing splines and grey shaded area denotes a 95% credible interval.Comparison of aggregation schemes shown in different colours: historical (blue) and modern (red).F I G U R E 4 Estimated bird species richness for the study area (Region of Murcia, SE Spain), considering both periods: (a) historical (1990-1991) and (b) modern (2012-2017).(c) Temporal turnover of avian species richness between two periods.Positive values indicate an increase in bird species, whereas negative values indicate a decrease in bird species between the historical and modern period.
an area less sensitive to climate change due to the water and food resources provided by these border environments, resulting in two patterns: species richness increase in high-altitude due to the effects of climate change; species richness increase in low-altitude areas with increased species turnover rates due to the effects of land-use change.According to our hypotheses, these results can indicate a dual pattern of global change due to the jointly effects of climate change and land-use change pattern may be related to the land-use changes in the peripheral areas forest systems in recent decades, such as the expansion of irrigated agricultural areas, agroforestry and urban areas.These landscape changes have conformed heterogeneous landscapes with moderate human influence, which may favour biodiversity(Estrada-Carmona et al., 2022).Considering our bird community, this would have resulted in the attraction of urban (such as house sparrows Passer domesticus) or ecotone species (many fringillids), as well as some naturally expanding species (like Streptopelia decaocto) favouring urban-forest interfaces, and others that would benefit from the supply of water resources (disappeared or reduced in forest areas).All these species would end up coexisting in these interface zones, as well as in the immediate forest systems (where they are recorded).The coastal and pre-coastal lowland areas of the south and southwest of the region were originally considered impoverished in some species that were in practice introduced in favour of irrigated agroforestry systems (Streptopelia turtur, Turdus merula).At the same time, some species from more arid habitats would have been lost, but the net result has been the same in the two previously cited syndromes, an increase in richness in the two extreme sectors of the altitudinal gradient, but with different

F
Richness-based species-exchange ratio, SERr for the study area (Region of Murcia, SE Spain), considering both periods: 21-25 years of time span in relation to global change, specifically to the joint effects of climate change and land-use change.We observed that elevation had an influence on the species richness, while forest cover behaves more as a control variable, since species richness, even declining, did not change its relationship with forest cover between time periods.We also found two different patterns of increase in bird species richness associated with different drivers of global change: in elevated areas with low turnover rate, mainly in response to the effects of climate change, by which new species can colonize elevated areas of forest systems; and in low-lying areas with a high turnover rate, associated with changes in land use related to human anthropization.Our results are therefore interpreted as a consequence of the combined effects of climate change and land-use change, which may be of great importance for assessing species changes and predicting population trends in the context of global change.The biodiversity patterns observed in our study area are a reflection of the pressures on an area affected by global change, and forest systems can serve as an open-air laboratory in which birds are a valuable indicator.