Land Use Changes Threaten Bird Taxonomic and Functional Diversity Across the Mediterranean Basin: A Spatial Analysis to Prioritize Monitoring for Conservation

Land use changes rank among the highest threats to biodiversity, but assessment of their ecological impact is impaired by data paucity in vast regions of the world. For birds, land use changes may mean habitat loss or fragmentation, changes in resource availability, and disruption of biotic interactions or dispersal pathways. As a result, avian population sizes and assemblage diversity decline in areas subjected to urbanization, agricultural intensification, and land abandonment worldwide. This threat is especially sensitive in hotspots such as the Mediterranean basin, where avifaunas of several biogeographic origins meet, encompassing numerous endemic taxa, and ecological specialists with low resilience to habitat modifications. Here, we correlated several facets of bird taxonomic and functional diversity to a fine-grained land-use change classification, in order to identify priority areas in need for enforced protocoled bird sampling in a conservation prospect. For this, we computed the species richness, functional richness, originality and specificity of 211 bird assemblages based on bird extent-of-occurrence data for 279 species and 10 ecological traits. We used a spatialized regression model to correlate bird diversity patterns with bioclimatic gradients and land use change between 1992 and 2018, assessed from an unsupervised clustering on 2 km resolution data. We showed that species-rich bird assemblages are subjected to agricultural intensification, while functionally diverse assemblages are mainly undergoing desertification and land abandonment. Unfortunately, most of these changes occur in areas where protocoled bird surveys with sufficient spatial and temporal resolution are lacking. In light of these results, we urge for the setting of bird monitoring programs targeted mainly on parts of North-Africa and the Levant, in order to allow a region-level evaluation of the threat posed by recent land use changes on the exceptional avifaunistic diversity of the basin. Fostering such regional-scale evaluations of congruences between human threats and centers of diversity is a necessary preliminary step for a pragmatic response to data deficiencies and ultimately setting appropriate responses to avoid the collapse of avian assemblages.

Land use changes rank among the highest threats to biodiversity, but assessment of their ecological impact is impaired by data paucity in vast regions of the world. For birds, land use changes may mean habitat loss or fragmentation, changes in resource availability, and disruption of biotic interactions or dispersal pathways. As a result, avian population sizes and assemblage diversity decline in areas subjected to urbanization, agricultural intensification, and land abandonment worldwide. This threat is especially sensitive in hotspots such as the Mediterranean basin, where avifaunas of several biogeographic origins meet, encompassing numerous endemic taxa, and ecological specialists with low resilience to habitat modifications. Here, we correlated several facets of bird taxonomic and functional diversity to a fine-grained land-use change classification, in order to identify priority areas in need for enforced protocoled bird sampling in a conservation prospect. For this, we computed the species richness, functional richness, originality and specificity of 211 bird assemblages based on bird extent-of-occurrence data for 279 species and 10 ecological traits. We used a spatialized regression model to correlate bird diversity patterns with bioclimatic gradients and land use change between 1992 and 2018, assessed from an unsupervised clustering on 2 km resolution data. We showed that species-rich bird assemblages are subjected to agricultural intensification, while functionally diverse assemblages are mainly undergoing desertification and land abandonment. Unfortunately, most of these changes occur in areas where protocoled bird surveys with sufficient spatial and temporal resolution are lacking. In light of these results, we urge for the setting of bird monitoring programs targeted mainly on parts of North-Africa and the Levant, in order to allow a region-level evaluation of the threat posed by recent land use

INTRODUCTION
Land use changes range amongst the most pressing threats to biodiversity (Sala et al., 2000;Foley et al., 2005;Millennium Ecosystem and Assessment, 2005;Montanarella et al., 2018). Urban expansion, the replacement of agro-pastoral mosaics by intensive industrial agriculture and land abandonment decrease habitat availability and connectivity worldwide (Newbold et al., 2015;Tilman et al., 2017). These dynamics generally tend to decrease taxonomic and functional diversity at multiple spatial scales, with strong variations across taxa and regions (Jetz et al., 2007;Titeux et al., 2016;Newbold, 2018). Identifying where and how land use changes threaten high-diversity areas is essential for the priorization of conservation efforts to halt the loss of natural heritage. However, the lack of temporal ecological data prevents the assessment of biodiversity trends in vast portions of the world, including hotspots such as the Mediterranean basin.
Rather than a uniform drop in biodiversity, land use changes trigger a decline in ecological specialist species, paralleled with an increase in generalist species ("winner-loser dynamics, " McKinney and Lockwood, 1999;Devictor et al., 2008). The determinants of these dynamics are multiple and interact together, but can roughly be summarized in a neat loss of habitat extent and suitability, disruption of connectivity through habitat fragmentation and direct disturbance or pollution (e.g., pesticides, Meehan et al., 2011). Besides a direct impact on species' local survival and demographic rates, these processes may interact with sensitivity to other threats, notably climate change (Mantyka-Pringle et al., 2012;Jantz et al., 2015) and invasion by alien species (McKinney, 2006;Sol et al., 2017). At an ecosystem level, land use changes may in turn alter metapopulation dynamics and biotic interaction networks, sometimes resulting in the loss of essential ecosystem services (Hooper et al., 2012;Hautier et al., 2015;Oliver et al., 2015).
Uncontrolled land exploitation and the lack of strong environmental policies increases pressure on biodiversity in places undergoing rapid human and/or economic growth, which are often insufficiently covered by biodiversity datasets (Amano and Sutherland, 2013;Meyer et al., 2015;Siddig, 2019;Nori et al., 2020). The Mediterranean basin presents a particularly contrasted pattern, with numerous standardized monitoring programs allowing an in-depth understanding of biodiversity dynamics in its north-western part (mainly Spain, France and Italy, Gregory et al., 2005), but a near-absence of similar initiatives along its southern and eastern shores. Apart from sparse and local surveys usually performed within limited timespan and with low sample sizes, most biodiversity knowledge in north Africa and the Levant therefore relies on opportunistic data that can hardly be used to estimate trends or rates of change (Siddig, 2019). This deficiency is especially worrying since just as human presence, south Mediterranean ecosystems are either coastal or limited in surface due to the border imposed by the Sahara desert.
Urban sprawl and artificialization of soils on the coastline and in periurban areas occur at a rapid pace in Mediterranean countries, particularly along the southern coast (Fao and Bleu, 2018). This trend involves pollution, habitat destruction, and overexploitation of yet limited natural resources. Fast urbanization and artificialization of fertile coastal areas reduce even more arable lands, which are yet particularly scarce, leading to the intensification of irrigated and rainfed agriculture and overexploitation of pastures (Puigdefábregas and Mendizabal, 1998;Caraveli, 2000;Voltz et al., 2018). Along the northern coast of the basin, which is also the best studied, these processes decrease biodiversity and trigger compositional changes in vertebrate species assemblages (Concepción and Díaz, 2010;Sokos et al., 2013;De Solan et al., 2019). For instance, the growth of intensive olive monocultures, at the expense of traditional agroecosystems such as winter cereals, extensive pastures and low-input olive farming (Stoate et al., 2001), led to a strong reduction in biodiversity (Santos and Cabral, 2004;Siebert, 2004). Although the effects strongly depend on the context and on the specific environmental boundary conditions (Queiroz et al., 2014), the conversion of former mountainous pastoral fields into vast extents of forest can lead to major losses in bird diversity. This was observed in several case studies in southern France and Spain (Suárez-Seoane et al., 2002;Blondel et al., 2010;Queiroz et al., 2014). While land abandonment might be beneficial for some generalist species, it reduces habitat quality and availability for open-habitats species (Preiss et al., 1997;Sirami et al., 2006). It also reduces landscape heterogeneity, and leads to the establishment of a secondary vegetation cover which may reveal unsuitable to many forest species (Suárez-Seoane et al., 2002;Poyatos et al., 2003;Geri et al., 2010). Given these intense modifications of landscapes, the Mediterranean basin therefore urgently needs an overview of recent land use changes and of their consequences on biodiversity.
Priorities must be set in order to overcome the limited resources available for south-Mediterranean biodiversity surveys, the political and economic obstacles preventing the emergence of a strong research network around the basin, and the urgency of addressing ecological threats related to land use changes (Lavorel et al., 1998). In the absence of longitudinal datasets allowing the estimation of temporal ecological dynamics, a first necessary step is to identify the degree of overlap between land use change and biodiversity, as assessed from the best possible sources of spatial data (Tilman et al., 2017). Specific surveys can then be set where potential conflicts are identified (e.g., high levels of land conversion toward intensive agriculture or urban areas) in a framework that optimizes the use of available sampling effort. Because taxa are unequally known, this approach is mostly realistic for the best studied groups, such as birds, for which most species are described and their distributions relatively well delineated in expert-validated maps (BirdLife International Nature Serve, 2012). The results of such studies are necessarily uncertain because species' range margins are unequally accurate, spatial resolution is usually coarse and species assemblages are inferred from the overlap of species' range maps rather than observed from field surveys (Herkt et al., 2017). In spite of these limitations, these distribution data remain the only option for a first evaluation of priority areas to set up dedicated monitoring.
In the present study, we aimed at identifying where recent land use change overlap with areas of high bird diversity across the Mediterranean basin. Because biodiversity and its threats are organized into multiple, partly orthogonal facets, relying on taxonomic diversity alone (e.g., species richness) is usually insufficient (McGill et al., 2006;Bellard et al., 2012). Functional diversity, or the diversity of traits within an assemblage, is now a classical way to incorporate species' ecological characteristics into biodiversity assessment (Devictor et al., 2010;Cadotte et al., 2011;Carmona et al., 2016). This approach rests on the assumption that ecological traits (species' characteristics that describe their relationships to their biotic and abiotic environment, sensu largo) are directly related to ecosystem processes and services (Tilman et al., 1997;Newbold et al., 2016). Irrespective of changes in species richness, functional diversity varies over regional scales under the influence of climatic gradients, distribution of resources and anthropogenic imprint, revealing the imprint of niche-related processes and evolution on biodiversity patterns (Oliveira et al., 2016;Barnagaud et al., 2017;Le Provost et al., 2020).
We predicted that the most intense land use changes would be associated with increasing urbanization and agricultural intensification in the semi-arid agroecosystems of the basin (southern shore and Iberian Peninsula) and in parts of the north-west, and, inversely with land abandonment mainly in the northern hills and mountains. These regions are climatic transition zones between warm deserts and coastal areas in the south, and between coastal areas and continental climates in the north; they have formed biodiversity reservoirs since the last glacial maximum (Blondel et al., 2010). Consequently, we expected that the strongest patterns of anthropogenic land use change would be associated with high bird diversity.
Although a land use typology already exists in the Mediterranean area (Malek and Verburg, 2017), it is static and therefore omits temporal dynamics. We thus constructed our own typology of land use changes from which we derived continuous gradients, which we correlated with bird taxonomic and functional diversity. We used these results to infer types of land use dynamics and areas of the Mediterranean basin where systematic bird monitoring should be set in priority to quantify the effective ecological imprint of land use change and limit their effects wherever possible.

Study Area
Defining and mapping the Mediterranean basin from a biological point of view has been a subject of debate among biogeographers for more than a century. There are no sharp borders with neighboring regions and many factors have to be considered, including vegetation, climate, latitude, and elevation. A consensus however exists among ecologists, historians, and geographers, which all agree that the unity and the specificity of this region is provided by its climatic pattern of hot, dry summers and humid, cool, or cold winters. Within the Mediterranean Basin itself, climate also changes with rising altitude in mountain ranges and when traveling from west to east. Overall, a sharp gradient exists between the colder, wetter northwestern and northeastern quadrants of the basin and the hotter, more arid, southeastern and southwestern parts in North Africa and the Levant. At its outer limits, the area borders a wide diversity of climatic and biogeographic zones, including boreal forests in central and northern Eurasia, the vast steppe regions of central Asia and northwestern Africa, and the hot subtropical deserts of northeastern Africa.
The extent of the studied area was thus defined according to the limits of the Mediterranean biome, classified as "Mediterranean forests, woodland and scrub" (Dinerstein et al., 2017). The region covers 27 countries (Figure 1). We defined a grid with a 1 • × 1 • resolution (approximately 110 km × 110 km at the equator) on a cylindrical equal area projection, which is considered as the optimal resolution for managing the bird data we used (Kissling et al., 2012). We then selected only the cells covered at least by 30% the Mediterranean biome, leaving us with 211 pixels.

Bioclimatic Data
Bird distributions are the outcomes of both short-and long-term biogeographical dynamics that arise under the joint effects of climate, habitats and biotic processes operating simultaneously over multiple timescales. We therefore characterized the bioclimatic gradients influencing bird assemblage diversity by a combination of current and historical climatic data. We retrieved Worldclim global-scale climatic data at 30 s spatial resolution 1 (Fick and Hijmans, 2017). We calculated the annual median temperature based on the monthly mean temperature and recorded temperature seasonality, annual precipitation, precipitation seasonality, and winter, spring, summer, and autumn precipitation, for the present  and for the Last Glacial Maximum (about 22,000 years ago). We used as well Worldclim elevation data at 30 s spatial resolution.
We upscaled the past and present bioclimatic variables at 1 • × 1 • spatial resolution by calculating the median of each variable based on the values contained in each pixel. We subsequently synthetized these variables in a Primary Component Analysis (PCA) through the ade4 R package (Dray and Dufour, 2007;R Core Team, 2020). Our rule of thumb for all multivariate analyses performed in this study was to keep the smaller number of axes so that the sum of eigenvalues are at least 90%, leaving us in this case with four axes. The final bioclimatic variables are thus the coordinates of each 1 • × 1 • pixel on each of the four axes.

Land Use Change
We retrieved 1992 and 2018 land cover maps ESA CCI Land Cover time-series at 300 m spatial resolution 2 (Esa, 2017). We only retained pixels with centroids falling into the limits of the studied area for the analysis. We then aggregated the initial 36 land use classes into 11 classes in order to limit the imprint of high-resolution land use changes incompatible with the grain of the study (Supplementary Table 1). Since restraining our analyses to the grain of available bird data could have hidden meaningful co-structures between bird diversity and land use changes at finer grains, we resampled the land use maps from their 300 m original spatial resolution at 2, 30, and 90 km and calculated the percentage of each land use category in these larger pixels. We then retrieved the amount of change in each land use class by subtracting the land use values of 1992 to those of 2018. As the objective is to focus on change, we discarded the pixels that do not display any change between the two dates.
In order to synthetize land use change patterns into interpretable trajectories, we performed an unsupervised clustering on land use change data. We chose partitioning around medoids (PAM; Kaufman and Rousseeuw, 2005), known to be more robust than k-means due to its lower sensitivity to noise and outliers (Han et al., 2001;Park and Jun, 2009;Arora et al., 2016). We used the sampling-based approach of PAM proposed in the algorithm Clustering Large Applications (CLARA; Kaufman and Rousseeuw, 2005) to reduce computing time and RAM storage problems (R package cluster; Maechler et al., 2019). We specified the number of clusters through graphical displays of silhouettes (Rousseeuw, 1987), which assess graphically whether objects lie well within the cluster, and give the user an appreciation of the relative quality of the clustering. Then, we estimated the optimal number of clusters so that it maximizes the average silhouette over a range of possible values (Factoextra R package; Kassambara and Mundt, 2017).
We performed separate cluster analyses on the land use change data at 2, 100, and 350 km in order to test the sensitivity of the method to the chosen scale (Supplementary Figure 1). To match bird data resolution, we calculated the percentage of each cluster in 1 • × 1 • pixels, which we synthesized through a fourdimensional correspondence analysis (CA) through the ade4 package (Dray and Dufour, 2007;R Core Team, 2020).

Bird Distributions
We retrieved extent-of-occurrence vector layers for all bird species occurring in the study area (n = 279) from Birdlife International 3 (BirdLife International Nature Serve, 2012) and intersected them with the 1 • × 1 • grid defined previously to obtain 211 contiguous species assemblages. These data are the most exhaustive possible and are reasonably smoothed for survey effort and other biases through expert advice, but they tend to be over conservative in poorly studied areas and thus to underestimate species' distributions, especially close to range margins (Herkt et al., 2017). This issue may affect the composition of species assemblages at the Saharan edge and in the Levant, although we have currently no reason to believe that it might induce a directional shift in assemblage-level diversity indices at the spatial resolution considered. Furthermore, because of the absence or near-absence of atlas data or opportunistic records usable for our purpose in these parts of the survey area, extent-of-occurrence maps are the least-worst solution for an overview of bird assemblages all around the Mediterranean basin.

Ecological Traits
A major source of heterogeneity in functional diversity indices lies in the choice of ecological traits (Calba et al., 2014). Since we did not aim at testing hypotheses on specific aspects of birds' ecology, we considered traits that represent the widest possible range of independent ecological strategies while avoiding trivial correlations (Villéger et al., 2008), as well as any circularity with responses to land use (hence, we discarded traits expressing habitat selection). We considered four biometric traits (mean body length, mean wing length, mean bill length, mean body mass) from Storchová and Hořák (2018). These measurements correlate with key aspects of birds' life history that affect their responses to environmental changes, including reproductive strategy, home range, and demography. We also included three traits directly associated with reproduction (mean clutch size, mean brood number per year, and mean age of first breeding) and two categorical traits associated with strategies in response to resource use and seasonality (territoriality and propensity to seasonal migration). Eventually, we characterized diet, which is key to multiple aspects of species' life history and ecosystemic functions (Kissling et al., 2012;Sekercioglu et al., 2016;Barnagaud et al., 2019), with the proportions of use for 10 types of preys (Wilman et al., 2014). 4

Functional Diversity Indices
Since functional diversity indices are inherently dependent on the dimensionality of the functional space (i.e., the smallest possible convex hull containing all the species of the regional assemblage), we first summarized the trait matrix in a principal coordinate analysis on a Gower distance matrix weighted for the non-independence of diet proportion data (Pavoine et al., 2009). We selected six principal coordinate axes, the minimum dimensionality providing a suitable representation of the initial trait dissimilarity among species (mean square error below 2%, method described in Maire et al., 2015). Not only may the range of traits be affected, but also species' functional originality and redundancy, which affect the resilience of ecosystem processes to disturbance (Mouillot et al., 2013). For instance, the conversion of a semi-arid area into agricultural fields may remove desert specialists but maintain generalist species and promote invasion by synanthropic species, resulting in no neat change in taxonomic diversity but profound modifications of functional diversity. Conversely, changes in forest practices or type of crops within an ancient agricultural area may lead to species replacements without deep changes in ecological functions. Therefore, we computed three complementary indices in addition to species richness to quantify the functional diversity of each bird assemblage based on this functional space (for a full description of these indices, refer to Mouillot et al., 2013;computations implemented in R using scripts available at http://villeger. sebastien.free.fr/Rscripts.html): -Functional richness (Fric), expressed as the proportion of the total functional space filled by a given assemblage, the closest trait-based equivalent to species richness differentiating assemblages that consist of few or numerous different traits; -Functional specialization (Fspe), the scaled average distance of species to the centroid of the assemblage in the functional space, a measure of whether an assemblage mainly consists in generalist (central) or specialist (marginal) species; -Functional originality (Fori), the scaled average pairwise distance of species within an assemblage in the functional space, as a measure of functional redundancy.
We discarded several other possible indices based on their similarity with these three or because their biological interpretability was limited in our context. All indices (and especially functional richness) were correlated to some extent with species richness. We therefore computed null distributions for each index by drawing randomly 99 species compositions for each assemblage while preserving observed species richness. We then calculated standardized effect sizes as the difference between observed indices and their null mean, divided by their null standard deviation.

Statistical Analyses
We used a linear regression model to quantify the relationship between each functional diversity index or species richness (over the 211 grid cells i), with the j climatic axes (BC 1 -BC 4 ) and the k land use change axes (LUC 1 -LUC 4 ). All slope coefficients were considered as assemblage-level random effects to account for spatial non-stationarity in the effects of climate and land system on bird assemblages. The general form of the regression model was therefore: In Eq. (1), µ was either the mean of a Gaussian distribution (functional diversity indices) or the log-transformed expectancy of a Poisson distribution (species richness), α was an intercept and β and γ were random regression coefficients. We estimated parameters in a Bayesian framework with flat normal priors for α, ν, and ω (null mean and variance of 1,000) and flat uniform priors for √ 1/τ and √ 1/ψ (min = 0, max = 10). We ran three Monte-Carlo Markov chains of 20,000 burn-in iterations + another 20,000 for inference, thinned by 20. Convergence was satisfactory for all parameters based on Gelman and Rubin's statistics and visual inspection of chains for hyperparameters. We assessed fit by visual checks of the correlation between functional diversity indices and their values replicated by the models (median over 3,000 iterations); no major signal of lack of fit appeared (values aligned on a 1-1 line) apart from a slight underestimation of upper values in most indices.
The statistical analyses were conducted with R software version 3.5.3 (R Core Team, 2020), and JAGS version 4.1.0 (Plummer, 2003).

Bioclimatic Gradients
The first axis of the PCA describing Mediterranean bioclimatic gradients (BC 1 , Table 1) was characterized by median and spring precipitation (positive side) and median temperatures (negative side). We therefore interpreted this axis as a latitudinal gradient ranging from Northern rainy and temperate zones to dry and warm semi-arid areas in the southern bank of the basin (Figure 2A). The second axis BC 2 was characterized by winter precipitation and precipitation seasonality (high values in the negative side, low values in the positive side, Table 1). Positive values corresponded to areas under the influence of the Atlantic Ocean, the Mediterranean islands and the Albanian, Greek, and Levant western coasts ( Figure 2B). They were described by a high variation of precipitation rates throughout seasons and high precipitation regimes in winter ( Table 1). Negative values indicated low precipitation rates throughout the year along Spanish, French and Italian coasts and hinterland and semiarid Maghreb (Table 1 and Figure 2B). The third axis BC 3 was dominated by a variation of temperature seasonality (high values in the negative side, low values in the positive side, Table 1) along a coastal/inland gradient ( Figure 2C). The Aegean part of Greece, Turkey and the Levant did not display the same pattern, with an overall high temperature seasonality. The Levant coasts displayed a high variation of precipitations (BC 2 , Figure 2C and Table 1) and a high variation of temperatures through seasons (BC 3 , Figure 2C and Table 1). Conversely, western Mediterranean and Atlantic coasts were characterized by temperature stability. The semi-arid Maghreb, characterized by a low seasonality of precipitations (BC 2 , Figure 2B and Table 1), displayed however high yearly temperature variations (BC 3 , Figure 2C and Table 1). A joint reading of BC 2 and BC 3 therefore shows the spatial discrepancy between temperature and precipitation seasonality. The fourth axis was dominated by elevation variability ( Table 1).
The positive values highlighted areas where mountains contrast with the flatness of the shores, while the negative values characterized large plains or plateaus ( Figure 2D).

Land Use Change Patterns
The unsupervised clustering of land use change variables, at 2 km resolution, resulted in an optimal number of nine clusters. Using Highlighted values: three major values (in red) and three minor values (in blue). The percentages are the amount of variance explained by each of the four axes (BC 1 -BC 4 ). LGM = Last Glacial Maximum. The "axis denomination" row synthetizes the aspects described by the axis, in its positive side (+) and in its negative side (-).
the 30 and 90 km resolutions instead did not substantially alter the results of bird-land use change relationships (Supplementary Figure 1). Therefore, we only described the results based on the nine 2 km resolution clusters. Cluster C1 gathered the pixels with very low rates of change, whatever the land use category ( Table 2). We thus chose to discard a posteriori the concerned pixels from subsequent analyses, as our objective was to focus on broad patterns. Cluster C2 brought together the areas undergoing a strong urban sprawl, eroding periurban rainfed agricultural lands on coasts, large valleys and plains (French Rhône valley, Italian Padania plain, Figure 3) and around cities (e.g., Madrid, Barcelona). Clusters C3, C4 and C6 concerned pixels that experienced forest regain after the abandonment of mountainous traditional agricultural lands and pastures; the three clusters described the same process, but with an increasing magnitude (Table 2 and Figure 3). Cluster C7 grouped pixels under the influence of agricultural intensification, corresponding to the sprawl of rainfed and irrigated croplands over natural areas with sparse vegetation, mostly in semi-arid agroecosystems of the High Plateaus of the Atlas, in the Spanish regions of Murcia and Andalucía, and at the Turkish-Syrian border and the Alep region (Figure 3). Cluster C5 reported an increase in pastures, shrublands, and on the mixed class with cropland and natural vegetation, at the expense of homogeneous forest cover ( Table 2). In contrast to cluster C7, the increase of irrigated and rainfed croplands conveyed by Cluster C5 mostly characterized northern countries, and was associated to an increase in pastures and shrublands and to a decrease of forests (Figure 3). The pixels concerned by clusters C8 and C9 underwent two antagonistic processes, mostly occurring in pre-Saharan lowlands (Table 2 and Figure 3). Cluster C9 highlighted a desertification process, with bare areas spreading over sparse natural vegetation particularly in Algeria and Tunisia, while cluster C8 pointed out a regain of natural vegetation on bare soils (mostly in Morocco, and in Algeria in lower proportions ( Table 2 and Figure 3).
The CA used to upscale the 2km land use change clusters at the resolution of bird data was structured by four axes ( Table 3). The first axis LUC 1 highlighted a longitudinal dichotomy in land use change processes. The Northern bank and the Maghreb coastline (negative values) were mostly characterized by the abandonment of agricultural activities in mountains and the regain of forests (clusters C3, C4, and C6). The positive values characterized agricultural intensification (cluster C7), sparse natural vegetation regain over bare soils (cluster C8), and desertification (cluster C9) in the semi-arid agroecosystems of Maghreb and Spain (Table 3 and Figure 4). LUC 2 and LUC 3 traduced the interaction between these three latter processes. Agricultural intensification (cluster C7) was located in the positive side of the two axes, and was opposed to desertification (cluster C9) in the negative side of LUC 2 and to sparse vegetation regain on bare areas (cluster C8) in the negative side of LUC 3 (Table 3 and Figure 4). LUC 4 was related to the replacement of former agricultural lands by urban areas on the positive side, and by forests in the negative side (Table 3). This axis expressed the dichotomy between soil artificialization on coasts and periurban areas, and the abandonment of agricultural activities in inland rural and mountainous areas (Table 3 and Figure 4).

Bird Diversity Patterns
Species and functional richness displayed opposite patterns in the Western part of the Mediterranean area, with the highest numbers of species but the lowest functional richness (Figures 5A,B). The East (Balkans, Turkey, and Levant) presented less obvious patterns, being a transition zone without any barrier effect (Figures 5A,B). Functional specialization and originality (Figures 5C,D) demonstrated similar patterns, with homogenous and ordinary assemblages in the north, and specific and original species in the southern part of the Iberic peninsula, Turkey, the Levant and semi-arid zones of Maghreb. Functional diversity indices were moderately inter-correlated (Pearson's R 2 FRic,FSpe = 0.62; R 2 Fric,FOri = 0.57; R 2 FOri, Fspe = 0.59).

Effects of Climate and Land Use Change on Birds
Regression parameters on bioclimatic variables revealed generally high effects on all dimensions of bird diversity (Figure 6). The third axis, related to temperature seasonality patterns, had the lesser impact. As observed on the maps of diversity indices (Figure 5), species richness and functional diversity were opposed along The values in red correspond to the three major positive values by cluster; the values in blue correspond to the three minor negative values by cluster. The percentages in the cluster column correspond to the percentage of pixels contained in the cluster.  Highlighted values: three major values (in red) and three minor values (in blue). The percentages are the amount of variance explained by each of the four axes (LUC 1 to LUC 4 ). The "axis denomination" row synthetizes the aspects described by the axis, in its positive side (+) and in its negative side (−). bioclimatic gradients (Figure 6, BC 1 ). Species richness was higher in cooler and wetter areas, in the mountains, and in areas with low seasonality, while warmer and more arid area, lowlands and areas with high precipitation seasonality had higher functional diversity (Figure 6, BC 2 , BC 4 ). Temperature seasonality (BC 3 ) seemed to impact positively functional richness and negatively functional originality. Due to its high correlation with the first bioclimatic Axis (R 2 BC1, LUC 1 = 0.51) we did not incorporate the first land use change Axis (LUC 1 ) in the model to avoid collinearity among predictors. The confidence intervals of land use change variables were larger than those of bioclimatic variables, which suggests more spatially variable coefficients. However, as for bioclimatic gradients, the maps of local effects did not display any clear pattern that could have suggested regional variations in bird-land use change relationships (Supplementary Figures 2-5). High species richness and low functional richness were located in areas experiencing agricultural intensification (mostly in the Northern Atlas and in Andalucía) (Figures 5, 7, LUC 2 , LUC 3 ). Conversely, desertification occurred in a few welldelineated zones with high functional richness and functional specialization (Figure 7, LUC 2 ). Revegetation of former semiarid bare areas was associated with high functional diversity and low species richness in a few locations of the High Atlas (Figures 5, 7, LUC 3 ). Conversely, agricultural abandonment and its subsequent reforestation corresponded with high species richness and high functional originality and specialization (Figure 7, LUC 4 ), for instance in former refuges or transition areas such as the Balkans, Turkey, or the Hispanic southern steppes (Figure 5). FIGURE 6 | Average effect of bioclimatic variables on bird diversity (ν). FOri, Functional Originality; FRic, Functional Richness; FSpe, Functional Specialization; RS, Species Richness. Whiskers correspond to 95% credibility intervals, boxes correspond to the interquartile distance and dots represent median estimates. P-values are computed as the number of MCMC iterations above zero relative to chain length. BC 1 (North (+)/South (-) climatic gradient), BC 2 (Low (+)/high (-) precipitation seasonality), BC 3 (Low (+)/high (-) temperature seasonality) and BC 4 (High (+)/low (-) topographic variability) refer to the Primary Components Analysis axes for bioclimatic variables.

DISCUSSION
Increasing pressure on land threatens biodiversity throughout the Mediterranean basin. This justifies large-scale conservation responses that are impaired by the lack of adequate monitoring data. Here, we identified regions where intense land use change dynamics threaten centers of bird diversity. We showed that while desertification, shrub encroachment, and agricultural land abandonment overlap with high bird functional diversity, agricultural intensification is associated with species-rich areas. Our results therefore reveal associations between land use change and multiple facets of bird diversity that are well delineated in space. In light of these conclusions, we propose directions to orient bird monitoring in order to respond efficiently to land use change throughout the basin. In addition, we show the value of combining high-definition land use change assessments with multiple facets of bird diversity to answer deficiencies in biodiversity sampling.
Our land use change modeling approach aimed at spatializing in a spatially continuous, homogenous, reproducible, and data-driven framework the broad ongoing processes in the Mediterranean area, which have been already extensively identified in the literature through small-scale and scattered case studies. Nevertheless, the land use data currently available over large spatial scales are not a sufficient material to represent ongoing landscape dynamics in the Mediterranean area, particularly with respect to the agricultural dimension. An appropriate estimation of agricultural intensification should be supported by data on the inputs, the type of cultivated crops, the size of the parcel, and many more dimensions that would enable going a step further land use change, toward a land system change assessment (Verburg et al., 2013). Although some of these information have been at best modeled and gridded worldwide at reasonable spatial resolutions (the MAPSPAM Spatial Production Allocation Model 5 for example), they are static and thus do not enable capturing dynamics. Furthermore, fine-grained land use data with an adequate temporal dimension are scarce and local, based on agricultural censuses that cannot be extrapolated over large areas. These data are, in particular, inexistent or inaccessible for most of the North-African, Levant and Balkan countries, where unfortunately stakes to birds are the highest. Therefore, bird monitoring protocols should systematically be accompanied with standardized, fine grained land use assessments. 5 https://www.mapspam.info/ We used static, large-scale distributional bird maps, which can only be stacked to assess the potential diversity of bird assemblages under strong assumptions on data quality and representativeness (Herkt et al., 2017). The absence of any temporal dimension, the coarse spatial grain, and the lack of data on the actual composition of local assemblages impair the possibility to assess the actual imprint of land use change on bird diversity, a consequence of the Wallacean shortfall outside Europe (Hortal et al., 2015). In particular, most North-African regions and vast parts of the Levant completely lack suitable surveys. Our study should thus be seen as a call for enforced effort for bird monitoring in these areas, and a guideline to target specific areas where strong land use change dynamics could alter exceptionally high levels of diversity.
Our results highlight several areas where intense land use change, toward either intensification or abandonment, may particularly threaten taxonomically or functionally diverse bird assemblages. Detecting deficiencies in bird surveys and anticipating future monitoring programs require a precise delineation of these areas, for which we suggest a few directions below. In particular, a vast species-rich area undergoes agricultural intensification in the Iberian Peninsula (southeastern Spain and the northeastern Pyrenean foothills). This area is already well surveyed by opportunistic and standardized monitoring programs, which have identified a combination of land abandonment, agricultural intensification and urbanization that triggers an overall reduction in bird diversity and biotic homogenization (Clavero and Brotons, 2010;Lasanta and Vicente-Serrano, 2012;Nainggolan et al., 2012;Sokos et al., 2013). These results should drive attention on the possibility of a similar biodiversity decrease in less-surveyed areas where we detected large clusters of agricultural intensification, which encompass the Moroccan Rif and Medium Atlas, as well as their northern lowlands in Meknes hinterland. These areas host unique bird assemblages matching continental and thermophile birds, including North African endemics (such as Levaillant's Green Woodpecker Picus vaillantii). Unfortunately, large scale bird monitoring programs are completely lacking there, and any assessment of land use effects on birds would need to rely on sparse surveys and opportunistic data.
Conversely, the desertification process, known to be subsequent to forest clearing and overgrazing in semi-arid steppe regions of North Africa and the Middle East (Nasr, 2003;Blondel et al., 2010), was correlated more with functional diversity rather than with species richness. This overlap concerned a large slick on the Libyan and Tunisian coasts, the high plateaus of the Saharan Atlas in the center-North of Algeria and isolated places such as parts of the Moroccan Drâa River, the Northeastern Cyrenaic coast around Derna, the West Bank and Northwestern Syria. In some of these places, the semi-arid habitats undergoing desertification were located close to patches where sparse vegetation colonizes bare areas with high bird functional diversity, typically in the Algerian Saharan Atlas and in southern Morocco. The closing of semi-arid open habitats could be either the result of an abandonment of extensive pastoral practices, or of active reforestation programs aiming toward soil erosion and runoff prevention that are in progress in North Africa and the Middle East (Blondel et al., 2010;Fao and Bleu, 2018). Associations of desertification and vegetation regain were mostly located along the border between the Saharan desert and the Mediterranean biomes, which host assemblages matching desert specialists and scrub-related bird species. These patterns should therefore incite to set up long-term bird monitoring schemes in these specific locations, where unique transitional bird assemblages may be jeopardized by even small changes in vegetation structure.
The border of species' distributions are naturally fluctuating due to demographic processes and high resource stochasticity (Sexton et al., 2009), which implies that the patterns we highlight here may not be of major concern for bird conservation unless monotonous trends of desertification and vegetation regains persist over the long term. Since the feasibility of a bird monitoring program is currently limited in much of the Saharan desert, efforts could be concentrated in southern Morocco, which is regularly visited by ornithologists. Similarly, our results show that vast extents of agricultural land abandonment in the Balkans correlate with high functional diversity areas. The abandonment of traditional land-use practices has been shown to be one of the most prevalent factors of biodiversity loss in the Mediterranean region, triggering habitat homogenization that leads to the loss of specialist open-habitat species (Sirami et al., 2006;Clavero and Brotons, 2010;Sokos et al., 2013). Land abandonment in the Balkans thus deserves specific attention since this area, which has acted as a climatic refuge in several past glaciations, is one of the main centers of European biodiversity (Griffiths et al., 2004). The biogeographic transition area located in western Turkey undergoes similar levels of land abandonment and should therefore also be subjected to closer bird monitoring.
Interestingly, species richness was opposed to all facets of functional diversity along most of the investigated land use change axes. Studies showing incongruences between taxonomic and functional diversity have warned against relying on species richness as a surrogate of all facets of biodiversity (Devictor et al., 2010;Lyashevska and Farnsworth, 2012). Our results are fully in line with the recommendation of considering simultaneously multiple taxonomic and trait-based indices (McGill et al., 2006;Cadotte et al., 2011), and add several interesting insights. Most importantly, the recurring opposition of species and functional diversity reveals that divergent land use dynamics could favor different facets of bird communities to the detriment of each other. For instance, shrub encroachment on bare lands in functionally rich bird assemblages may lead to the loss of typical desert species in northern Africa (such as several lark or wheatear species), while creating suitable conditions for generalist species (e.g., finches). This biotic homogenization dynamics is already well known in Europe and North America (La Sorte and McKinney, 2007;Le Viol et al., 2012) but has been little documented in North Africa. This lack of knowledge is specifically worrying since areas undergoing shrub encroachment are typically located at mid/high-elevations in the Atlas mountains, to which some species are restricted due to their affinities with cold climate (e.g., the North African subspecies of horned lark Eremophila alpestris atlas, del Hoyo et al., 2004). For such species, vegetation regain may inflate the effect of climate change as their altitudinal range narrows while growing vegetation decreases habitat suitability. Similar issues may arise in relation with desertification and agricultural intensification in other areas, enforcing an urgent need for monitoring programs to settle suitable actions for the preservation of these rangerestricted species or taxa.
Another point to note is that all the facets of functional diversity considered in our study were associated with similar land use types, although they were only moderately correlated and were not fully congruent in space. This is likely the consequence of strong biogeographic and climatic imprints on the distribution of the species and traits we considered, which therefore implies that land use change impacts on bird assemblages could be modulated by these large-scale constraints. Consistent with this interpretation, we did not detect any residual pattern in the spatial distribution of bird-land-use correlations in models accounting for bioclimatic variables. Hence, the correlation among the diversity indices used in our study should incite to incorporate the complementary effects of climate and land-use in models predicting the composition of bird assemblages (Dale, 1997;Williams and Newbold, 2020).

CONCLUSION
Our study revealed that centers of bird diversity are located in areas subjected to strong land use changes throughout the Mediterranean basin. Two serious issues arose behind this result, which should be addressed urgently through dedicated monitoring. First, the centers of bird taxonomic and functional diversity in the basin diverge under the interacting effects of long-term anthropogenic dynamics, biogeographic heritage and climatic gradients. Furthermore, land use changes exhibit complex patterns that are not fully correlated with climate. As a consequence, we showed that the land use dynamics that may threaten taxonomic diversity differ from those that may affect functional diversity; this calls for a better integration of spatial non-stationarity in the investigation of land use change impacts on bird diversity. Second, faunistic data are deficient on wide extents of the Mediterranean basin, where ongoing anthropogenic processes could already have impacted bird assemblages. In the near-absence of surveys setting a before-impact reference, the diversity of these assemblages prior to the initiation of current land use dynamics will remain unknown. Because a systematic sampling of birds is unrealistic in most of the basin, we identified data-deficient areas where protocoled bird monitoring is both urgent and feasible. This prioritization approach, albeit suboptimum, will probably be the most efficient way to evaluate the impacts of land use change in data-deficient regions (Hortal et al., 2015).
Land use changes have been identified as a major driver of biodiversity collapse globally (Newbold, 2018). Facing this challenge, our study warns that vast extents of land are being modified at a high pace without any possibility to assess the ecological consequences. In response, the level of threat that land use change imposes on species assemblages has to be evaluated at a much finer resolution than currently available global data. Sustainable conservation responses will, in turn, directly depend on the elaboration of efficient monitoring based on this initial evaluation. We therefore hope that our study will help setting monitoring priorities for birds in the Mediterranean basin and stimulate similar approaches in other biodiversity hotspots.

AUTHOR CONTRIBUTIONS
JF conceived the idea, led the study, and performed land use change and bioclimatic typologies. J-YB provided the bird data, computed bird diversity indices, and supervised the ecological part of the study. JP and EW constructed the statistical models. JF and J-YB led the writing. JF, J-YB, EW, and JP contributed substantially to the manuscript at all stages. MD supervised the land use dynamics assessment and reviewed the writing. AB supervised the work package and reviewed the writing. All authors contributed to the article and approved the submitted version.