Diversity of European habitat types is correlated with geography more than climate and human pressure

Abstract Habitat richness, that is, the diversity of ecosystem types, is a complex, spatially explicit aspect of biodiversity, which is affected by bioclimatic, geographic, and anthropogenic variables. The distribution of habitat types is a key component for understanding broad‐scale biodiversity and for developing conservation strategies. We used data on the distribution of European Union (EU) habitats to answer the following questions: (i) how do bioclimatic, geographic, and anthropogenic variables affect habitat richness? (ii) Which of those factors is the most important? (iii) How do interactions among these variables influence habitat richness and which combinations produce the strongest interactions? The distribution maps of 222 terrestrial habitat types as defined by the Natura 2000 network were used to calculate habitat richness for the 10 km × 10 km EU grid map. We then investigated how environmental variables affect habitat richness, using generalized linear models, generalized additive models, and boosted regression trees. The main factors associated with habitat richness were geographic variables, with negative relationships observed for both latitude and longitude, and a positive relationship for terrain ruggedness. Bioclimatic variables played a secondary role, with habitat richness increasing slightly with annual mean temperature and overall annual precipitation. We also found an interaction between anthropogenic variables, with the combination of increased landscape fragmentation and increased population density strongly decreasing habitat richness. This is the first attempt to disentangle spatial patterns of habitat richness at the continental scale, as a key tool for protecting biodiversity. The number of European habitats is related to geography more than climate and human pressure, reflecting a major component of biogeographical patterns similar to the drivers observed at the species level. The interaction between anthropogenic variables highlights the need for coordinated, continental‐scale management plans for biodiversity conservation.


| INTRODUC TI ON
The need to preserve dynamic ecosystems under changing climates and increasing anthropogenic pressure challenge traditional conservation approaches that are based on the current distribution of species. An application-oriented way forward may lie in protecting those landscape elements that support the coexistence of many species. Indeed, the habitat (or ecosystem) approach for conservation has been recently highlighted by the IUCN as a necessary step for conservation. Under this view, habitat diversity is a complex, spatially explicit measure of biodiversity (Bunce et al., 2013), which has proven to be a prominent driver for species diversity of a variety of taxa at the landscape scale (Alsterberg et al., 2017;Dianzinga et al., 2020;Gibb et al., 2020;Keppel et al., 2016;Kerr & Packer, 1997).
According to the EU Habitats Directive (Council Directive 92/43/EEC), the term "habitat" refers to an environmental unit defined by specific abiotic and biotic factors. Although alternative definitions exist (Davies et al., 2004;Drakou et al., 2011;Hall et al., 1997;Kearney, 2006;Mitchell, 2005;Yapp, 1922), this formulation provides a pragmatic operational tool for characterizing landscape elements of conservation priority.
Term "habitat" in the context of the EU Habitats Directive has a particular meaning, which deviates from the autecological species-related concept "habitat" in ecology. As an "environmental unit" that includes species assemblages and site conditions, the term "habitat" used in this context is closer to the concept of ecosystems (sensu Keith et al., 2013), even if some units are rather defined by mere plant communities (reflected in phytosociological terminology) (e.g., 6190 "Rupicolous pannonic grasslands (Stipo-Festucetalia pallentis)") and others are classified based on abiotic site conditions (e.g., 8240 "Limestone pavements," 8320 "Fields of lava and natural excavations"), places (e.g., 8310 "Caves not open to the public"), or geographical units (e.g., 1150 "Coastal lagoons," 1620 "Boreal Baltic islets and small islands"). Some units are characterized by vegetation structures (e.g., 5400 "Phrygana") while others by typical species (e.g., 6160 "Oro-Iberian Festuca indigesta grasslands").
For the operational meaning of habitats, fluxes of energy and matter, and the processes that are forming an ecosystem, are not considered, even if these might be important (e.g., carbon sequestration, evapotranspiration). The concept of habitat is thus more focused on ecological compartments with a main emphasis on vegetation. Being aware of this inconsistency, using habitat types as indicated by EU Habitats Directive, provides a standardized and legally established tool for monitoring and assessment complex units, which is crucial for nature conservation.
Habitat diversity can be measured as the number of different habitats in a given area (Hortal et al., 2009;Triantis et al., 2006)-herein referred to as "habitat richness." Habitat richness can be monitored in situ or by remote sensing techniques Radeloff et al., 2019;Tuanmu & Jetz, 2015). Moreover, the accessibility of habitat distribution data is steadily growing, often provided in the form of maps, which may sometimes be proxies forand the only available information on-the distribution of specific groups (e.g., plants and invertebrates).
Habitat richness reflects environmental conditions and can be used as an explanatory variable for modeling the distributions and abundances of species or communities (Heidrich et al., 2020;Leclère et al., 2020). Many studies have focused on the factors regulating the spatial variation in species richness (Brown & Lomolino, 2005;Field et al., 2009;Gaston, 2000;Howard et al., 2020;Quintero & Jetz, 2018;Rosenzweig, 1995). Species richness is typically correlated with variables such as climate (Gao & Liu, 2018;Thuiller et al., 2005), latitude (Gaston, 2000(Gaston, , 2007Hillebrand, 2004), topographic models, and boosted regression trees. The main factors associated with habitat richness were geographic variables, with negative relationships observed for both latitude and longitude, and a positive relationship for terrain ruggedness. Bioclimatic variables played a secondary role, with habitat richness increasing slightly with annual mean temperature and overall annual precipitation. We also found an interaction between anthropogenic variables, with the combination of increased landscape fragmentation and increased population density strongly decreasing habitat richness. This is the first attempt to disentangle spatial patterns of habitat richness at the continental scale, as a key tool for protecting biodiversity. The number of European habitats is related to geography more than climate and human pressure, reflecting a major component of biogeographical patterns similar to the drivers observed at the species level. The interaction between anthropogenic variables highlights the need for coordinated, continental-scale management plans for biodiversity conservation.

K E Y W O R D S
anthropogenic impact, biodiversity conservation, environmental predictors, European habitat directive, habitat richness, terrain ruggedness index heterogeneity (Hortal et al., 2009), and anthropogenic pressure Malavasi et al., 2016). Contextually, the habitat amount hypothesis (Fahrig, 2013) predicts that species richness in equalsized sample sites should increase with the total amount of habitat.
Similarly, the number of species inhabiting a region can be explained by the number of vegetation types as a surrogate of habitat diversity (Jiménez-Alfaro et al., 2016). However, the relative importance of habitat amount and its spatial configuration (e.g., fragmentation, connectivity, or perimeter/area ratio) on biodiversity patterns is still subject of debate (Fahrig, 2021;Saura, 2021aSaura, , 2021b. Despite the importance of habitat richness and the large amount of spatial data available for many terrestrial habitats, there is a knowledge gap on the mechanisms that determine the spatial patterns of habitat richness, especially at continental scales. To our knowledge, no research has been done on the underlying factors associated with habitat richness at continental level. In Europe, biodiversity conservation policy focuses on habitat protection (e.g., Council Directive 92/43/EEC), and achieving measurable improvements of the conservation status of Natura 2000 habitats is one of the main targets of the 2030 Biodiversity Strategy. As mandatory biodiversity monitoring, every EU member country is obliged to provide on a regular basis (every six years-art. 17 Habitat Directive) habitat distribution maps on a 10 km × 10 km grid map. Thus, disentangling the role of environmental factors determining habitat richness is a key challenge both for basic understanding of biodiversity distribution and for guiding conservation strategies (Mücher et al., 2009;Zhang et al., 2020).
To investigate how habitat richness is correlated with environmental factors across the EU, we used habitat distribution data from the third report of the EU Habitats Directive (EEA, 2020a).
Specifically, we seek to answer the following questions: (i) how do bioclimatic, geographic, and anthropogenic variables affect habitat richness? (ii) Which factor is the most important? (iii) How do interactions among these variables influence habitat richness and which combinations produce the strongest interactions?

| Habitat richness
We used habitat distribution maps of 222 terrestrial habitats of community interest from the EU 2007-2012 reporting period, obtained from the European Environment Agency (EEA, 2020a). Distribution of the protected habitat type in Europe is based on the standard grid (10 km × 10 km) provided by EEA for habitat monitoring (EEA, 2013), and we assume it is representative of the whole ecosystem diversity. From a total of 231 habitat types, 9 marine habitats were excluded according to the appendix 2 of the "Lists of existing marine Habitat types and Species for different Member States" (European Commission, 2021). We calculated EU habitat richness (thereafter "habitat richness") by summing up the number of habitat type reported in each cell of the above-mentioned map provided by the EEA. All the EU countries that contributed to the third report were included, except Greece due to lack of habitat reporting there.
Despite being composed of equal-area cells (i.e., 10 km × 10 km), the EEA grid has many cells located in coastal areas, often reducing the terrestrial surface within the cell. It has been widely reported that both species richness and habitat richness increase as a function of area (Lomolino, 2000;Rosenzweig, 1995;Triantis et al., 2012).
Because of the effect of area on habitat richness, we applied a modified log-log power function (Arrhenius, 1921;Triantis et al., 2012) to normalize habitat richness values based on within-cell land area. This involved adding 1 to habitat richness in order to have a normalized value of 0 when no habitat was present (instead of having minus infinite) and the absolute value of the denominator in order not to have inconsistent values when area was between 0 and 1 km 2 . can be found in Hoffmann et al. (2018). To test the assumption that habitat richness is a proxy of biodiversity conservation status, we calculated Pearson's r correlation coefficient (Pearson, 1931) to assess the correlation between habitat richness and species richness.

| Environmental variables
We investigated the relationship between habitat richness and three groups of environmental variables: bioclimatic, geographic, and anthropogenic. Multicollinearity among variables was assessed through Pearson's r. Within variable pairs holding Pearson's r > .7 (see Figure S1), the one judged to have less ecological importance was discarded from model building (Dormann et al., 2013;Elith et al., 2006).
We extracted the 19 bioclimatic variables from WorldClim (Fick & Hijmans, 2017) at 10 km × 10 km resolution. Through variable selection (see Brandt et al., 2017), we ensured balanced representation of temperature and precipitation (both mean and seasonal variation-see Table 1 for a summary of the explanatory variables).
Geographic variables, such as latitude and longitude, strongly affect species and habitat richness at different spatial scales (Drakou et al., 2011). Moreover, northing and easting (i.e., latitude and longitude, respectively, as geographic Cartesian coordinates) may be used to account for spatial autocorrelation. For these reasons, we included these two factors as explanatory variables in all the models. In order to capture geographic heterogeneity (Dufour et al., 2006), we included the terrain ruggedness index (TRI-Riley , 1999). TRI is defined as the mean of the absolute differences in elevation changes (Riley et al., 1999), and it is highly positively re-  (Table 1). Landscape fragmentation was extracted from EEA (EEA, 2011), and it was calculated based on the Effective Mesh Density, which is a measure of the degree to which movement between different parts of the landscape is interrupted by fragmentation geometry (FG) (Jaeger, 2000). The more FGs fragment the landscape, the higher the effective mesh density and hence the frag-

| Data analysis
To test the significance of habitat richness as proxy of biodiversity, we first analyzed the relation between habitat richness and the richness of the species listed in the Annex species of the Birds and Habitats Directives, as reported by the single countries, by using Pearson correlation coefficient.
We then used generalized linear models (GLMs), generalized additive models (GAMs), and generalized boosted models (GBMs) to investigate how habitat richness responds to the selected climatic, geographic, and anthropogenic variables at the EU scale.
To identify the interactions to be included in the models, we ran all the possible combinations of variable pairs using the glmulti function from the glmulti package (Calcagno & de Mazancourt, 2010).
To assess the best settings for the GBM models, all the possible combinations of three different numbers of trees (10,000, 15,000, and 20,000), four interaction depths (3, 5, 7, and 9), three shrinkages (0.01, 0.1, and 0.5) and three bag fractions (0.65, 0.8, and 1) were employed. We selected the combination with the lowest root mean square error (RMSE).
We run autocovariate models to filter out spatial dependence (Harisena, 2021). These models are usually based on autocovariate estimation directly on the response variable. Crase et al. (2012) developed a related procedure, in which autocovariates are quantified from the model residuals, rather than the raw data. This leads to an autocovariate that captures only the variance not explained by explanatory variables (Fletcher & Fortin, 2018). Moreover, models including autocovariates typically provide unbiased estimates of fixed effects, as demonstrated by Bardos et al. (2015). To fit autocovariate models, we calculated autocovariates on the model residuals and then we used these covariates in a new model. These autocovariates TA B L E 1 Summary of explanatory and response variables. In bold variables selected for model building were calculated using the autocov_dist function in the spdep package (Bivand & Wong, 2018).
Variable importance for GLMs and GAMs was calculated excluding explanatory variables one by one from the models. Then, the contribution of each variable was assessed using the difference in deviance (D 2 ) between the model with and without that variable (Márcia Barbosa et al., 2013).
To produce response curves, we used the modified inflated response curve (Zurell et al., 2012) for abundance models. Interaction curves were produced by setting the explanatory variables to their mean values.
In addition, we calculate Pearson's r correlation for each candidate explanatory variable (1st and 2nd order polynomials) in order to show the bivariate relationship between normalized habitat richness and the entire set of environmental variables ( Figure S3).

| RE SULTS
Geographic variables showed similar effects across the different models, having the greatest cumulative contribution (Figure 2).
Northing and Easting were the most important variables, with about 25% of relative influence each. Terrain ruggedness index (TRI), with about 20% of relative influence, was the third most important variable affecting habitat richness.
Climatic variables had widely different relative influence, with annual mean temperature having far more influence (ca. 15%) than the other three climatic variables taken together (annual temperature range, total annual precipitation, and precipitation seasonality; Figure 2). Finally, anthropogenic variables had only a minor contribution (ca. 5%) (Figure 2).

| Effects of environmental variables on habitat richness
Both mean annual temperature and annual precipitation had a positive effect on normalized habitat richness, which also increased linearly with TRI ( Figure 3). Habitat richness diminished toward eastern regions of the EU, while latitude showed an initial negative effect on habitat richness, it then slightly increased. The interactions between geographic and bioclimatic variables have shown a strong impact on habitat richness (Figure 4). The highest values of habitat richness were observed at low latitudes, where annual precipitation was moderately high. The positive trends observed for mean annual temperature and mean annual precipitation were considerably stronger in more rugged cells (higher TRI in Figure 4).
Other interactions among bioclimatic variables did not show any remarkable trend. The landscape fragmentation index on its own showed a positive effect on habitat richness. In contrast, population density index did not reveal any clear pattern (Figure 3), though it has interesting interaction with landscape fragmentation values; in particular, high landscape fragmentation affects positively habitat richness at low population but negatively at high population density ( Figure 4).

| Model performance
The explained deviance of the GLMs without autocovariate was 0.22, whereas after accounting for spatial autocorrelation, the explained deviance reached 0.59 (Table 2)

| DISCUSS ION
Our analysis shows that the relationship between habitat richness and the richness of species with conservation concern is positive and monotonic. Indeed, it is widely reported that a high number of habitats support a high number of species (Hortal et al., 2009) and vice versa (Jiménez-Alfaro et al., 2014). This confirms that habitat diversity at continental scales can be used as a pragmatic proxy for species richness and is thus a useful tool to assess the status of biodiversity conservation (Kallimanis et al., 2008;Pyšek et al., 2002;Saura, 2021a). We also found that geographic variables are the most relevant variables to shape habitat richness, while the effects of bioclimatic and anthropogenic variables were less evident but still significant.

F I G U R E 3
Relationships between explanatory variables and normalized habitat richness. The density of (10 km × 10 km) grid cells is indicated by hexagonal binning using the viridis color scale (varying from high density in yellow to low density in violet). Gray lines represent the 100 inflated response curves averaged across the three models used: generalized linear models (GLMs), generalized additive models (GAMs), and boosted regression trees (BRTs). Red lines are the median value, violet lines are the mean value of the inflated response curves. Geographic variables: NORTH = northing, EAST = easting, TRI = terrain ruggedness index. Bioclimatic variables: BIO_1 = mean annual temperature, BIO_7 = temperature annual range, BIO_12 = total annual precipitation, BIO_15 = precipitation seasonality. Anthropogenic variables: FRAG_IND = landscape fragmentation index, POP_DENS = human population density

| Geographic variables
Geographic variables showed the strongest association with habitat richness. We found that a major predictor of large-scale habitat richness was latitude, which is considered a proxy for other environmental variables (e.g., solar radiation and productivity; Archibald et al., 2010, Qian & Ricklefs, 2011. This variable is widely used as a predictor of species richness and diversity (Gaston, 2007;Hillebrand, 2004). In particular, our findings support the general tendency of biodiversity to decrease from lower to higher latitudes (Fine, 2015; F I G U R E 4 Surface plots show the interactions among the explanatory variables, x-and y-axis represent pairs of explanatory variables and z-axis is the magnitude of the interaction on the response variable. Only interactions above the given threshold (|z| = 0.3) are displayed. Geographic variables: TRI = terrain ruggedness index, NORTH = northing and EAST = easting. Bioclimatic variables: BIO_1 = mean annual temperature, BIO_7 = temperature annual range, BIO_12 = total annual precipitation, BIO_15= precipitation seasonality. Anthropogenic explanatory variables: FRAG_IND = landscape fragmentation index, POP_DENS = human population density TA B L E 2 Explained deviance (D 2 ) and root mean square error (RMSE) of the models with and without autocovariate ("RAC")  , 1984;Stevens, 1989). At smaller scales (i.e., national), Drakou et al. (2011) report a positive relationship with habitat richness for both latitude and longitude. In our study, a weak decrease of habitat richness was observed for longitude, probably due to lower habitat richness of the eastern countries, due to continentality and likely as a result of varying completeness in reporting the full lists of habitat types for some SE countries.

MacArthur
Habitat richness was also positively correlated with topographic complexity. The importance of environmental heterogeneity in controlling biodiversity is widely recognized in ecological theory (Hjort et al., 2015;Huston, 1994;Marini et al., 2011;Stein et al., 2014). TRI can be considered one of the main factors contributing to explain habitat and species richness (López-González et al., 2015;Stein et al., 2014;Tews et al., 2004), and its effect is expected to increase with spatial grain (Stein et al., 2014). Indeed, topographically complex areas offer a larger number and more different local conditions than topographically simple areas, leading to a higher number of habitat types packed into the same area (Stein et al., 2014). TRI contributes to species richness not only by providing an abundance of niches in space but also by offering relatively stable niches in time (Davies et al., 2007;Irl et al., 2015;Thuiller et al., 2006). Several studies showed an increase of species richness in relation to the increase of surface complexity (Cramer & Verboom, 2017;Farwell et al., 2021) and highly heterogeneous areas supporting more species than areas of lower heterogeneity (Hortal et al., 2009;Kallimanis et al., 2008;Rosenzweig, 1995).

| Climate variable
Temperature and precipitation variables were the second most important group in explaining habitat richness. It is widely reported that climatic variables are considered as main drivers of broad-scale patterns in species richness (Grytnes & McCain, 2007;Thuiller et al., 2004;Vetaas et al., 2019;Xu et al., 2014). We observed a positive association of habitat richness with annual precipitation as reported also for aquatic, coastal, and forest EU habitats (Drakou et al., 2011).
Precipitation is very unevenly distributed across time (among and within years) and space in Europe (Rajah et al., 2014;Zolina, 2012).
Due to climate changes, an increase in precipitation variability is expected in the near future, and this phenomenon could negatively affect both species and habitat diversity (Adler & Levine, 2007;Pearson & Carroll, 1998). Moreover, not only precipitation but other factors such as potential evapotranspiration (Adhikari et al., 2019) or soil water availability (Daws et al., 2002;Vetaas & Ferrer-Castán, 2008) should be considered to better predict species and habitat distributions.
Among climate variables, mean annual temperature was the most strongly associated with habitat richness. Positive correlation among plant species richness and temperature is widely reported (Diogo et al., 2020;Gottfried et al., 2012), but this relationship may be reversed when water availability is limited (Pausas & Austin, 2001).
Annual temperature range has a slight positive effect on habitat richness with a final slight decrease, highlighting the secondary role of temperature range in shaping species and habitat distribution (Austin & Niel, 2011). Interactions among climatic variables did not produce any strong effects as compared to the single bioclimatic variables. Instead, the interactions of mean annual temperature and mean annual precipitation with TRI suggest a positive effect on habitat richness of environmental heterogeneity with both the increase in temperature and precipitation.

| Anthropogenic variables
Population density is considered one of the main proxies defining human pressure on nature (Sanderson et al., 2002;Venter et al., 2016). Among the anthropogenic variables considered, the LFI showed the strongest association with habitat richness, as already found by other authors (e.g., Jaeger, 2000). This anthropogenic factor breaks ecological interrelations between the habitat patches and decreases their ability to provide various ecosystem services (Jaeger, 2000).
The decrease of habitat richness along with increasing human pressure was not as strong as expected. Indeed, environmental factors are typically more important for explaining species richness than human impacts (Howard et al., 2020). Moreover, the weaker influence of anthropogenic variables was probably due to the coarse resolution used in this study (Curtis et al., 2018;Niemiec et al., 2018;Woodbridge et al., 2020).
On the other hand, the synergistic interaction between fragmentation and population density had a strongly negative influence on habitat richness, also affecting habitat conservation (Ewers & Didham, 2006). This finding suggests that the interaction of anthropogenic variables (Newbold et al., 2015) could be of greater importance than climate interactions (Holman et al., 2017;Lehsten et al., 2015). Thus, land-use modification in the near future should be planned in order to decrease landscape fragmentation and increase habitat connectivity.  -Leite et al., 2020) and biodiversity loss (Campagnaro et al., 2018;Cardillo et al., 2004;Chase et al., 2020;Davies et al., 2006;Leclère et al., 2020;Pacifici et al., 2017). Landscape transformation and habitat degradation outside protected areas may also contribute to landscape fragmentation leading to isolated "islands" with low connectivity (Chase et al., 2020;Keeley et al., 2018). The present study highlights how biodiversity policies of the EU, such as Habitats Directive, have a central role not only in biodiversity conservation (Gameiro et al., 2020) but also providing continental scale data which are fundamental to investigate biodiversity patterns, as already demonstrated by Hoffmann et al. (2018) in relation to the European protected area network. This finding has important implications for conservation planning through the use of European habitat inventories originally created for reporting regional status of the Natura 2000 network.

| CON CLUS IONS AND IMPLIC ATIONS FOR CONS ERVATION PL ANNING
Knowing the inherent vulnerability of some habitats could aid decisions regarding European conservation priorities and could form the basis for a biodiversity conservation roadmap (Arlidge et al., 2018).

ACK N OWLED G M ENTS
We would like to thank the Biodiversity and Macroecology Group (BIOME) and Richi from Zapap for the continuous support during the project. Research contributing to this study was funded by the project "Development of a National Plan for Biodiversity Monitoring" (Italian National Institute for Environmental Protection and Research -ISPRA). BIOME Group was partially supported by the H2020 SHOWCASE (Grant agreement No 862480) and by the H2020 COST Action CA17134 'Optical synergies for spatiotemporal sensing of scalable ecophysiological traits (SENECO)'.

CO N FLI C T O F I NTE R E S T
No conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available Alessandro Chiarucci https://orcid.org/0000-0003-1160-235X