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

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.

INTRODUCTION

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 (Jung et al., 2020; 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 for—and 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, 2007; Hillebrand, 2004), topographic heterogeneity (Hortal et al., 2009), and anthropogenic pressure (Liu et al., 2018; Malavasi et al., 2016). Contextually, the habitat amount hypothesis (Fahrig, 2013) predicts that species richness in equal‐sized 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, 2021a, 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?

MATERIALS AND METHODS

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 km2.

Relations with species richness

Species richness data were obtained from the distributions of species reported in the Annexes of the Birds Directive (Council Directive 2009/147/EC; Annex I to V) and the Habitats Directive (Council Directive 92/43/EEC; Annex II, IV, and V) reported by EU member states for 2007–2012 (EEA, 2020a, 2020b). Further details 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).

TABLE 1

Summary of explanatory and response variables. In bold variables selected for model building

Group	Acronym	Description	Type	Data source
Bioclimatic	BIO_1	Mean Annual Temperature	Numeric	Fick and Hijmans (2017)
Bioclimatic	BIO_4	Temperature Seasonality	Numeric
Bioclimatic	BIO_7	Temperature Annual Range	Numeric
Bioclimatic	BIO_12	Annual Precipitation	Numeric
Bioclimatic	BIO_15	Precipitation Seasonality	Numeric
Bioclimatic	BIO_17	Precipitation of Driest Quarter	Numeric
Anthropogenic	FRAG_IND	Landscape Fragmentation Indicator	Raster	EEA (2011) and Jaeger (2000)
Anthropogenic	STREET_LENGTH	Total street length	Shapefile	Meijer et al. (2018)
Anthropogenic	STREET_DENSITY	Total street length divided by the area of a 10 km x 10 km cell	Shapefile	Meijer et al. (2018)
Anthropogenic	POP_DENS	Population density	Raster	EEA (2009) and Gallego (2010)
Geographic	TRI	Terrain Ruggedness Index	Derived from DEM	EU‐DEM (2021)
Geographic	NORTH	Northing		/
Geographic	EAST	Easting		/
Response	HAB_RICH	Habitat richness	Count	/
Response	NORM_HAB_RICH	Habitat richness normalized	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 et al., 1999). TRI is defined as the mean of the absolute differences in elevation changes (Riley et al., 1999), and it is highly positively related to elevation (Amatulli et al., 2018). To obtain the Europe‐wide TRI, we used a freely available, 20‐m resolution digital elevation model (DEM) for the EU (EU‐DEM, 2021). We resampled the DEM to 200 m × 200 m pixel size using the arithmetic mean for aggregation. Then, for each 10 km × 10 km grid cell, we summed all the 2500 values obtained for TRI to obtain the topographic heterogeneity within each cell.

Anthropogenic factors and their interaction with bioclimatic and geographic variables strongly affect animal and plant biodiversity at global terrestrial scale (Mantyka‐Pringle et al., 2012, 2013). Three ‘anthropogenic’ variables were considered: landscape fragmentation index (LFI), population density, and total street length and street density (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 fragmentation (EEA, 2011; Jaeger, 2000; Moser et al., 2005; Schmiedel & Culmsee, 2016). Population density was also extracted from EEA (EEA, 2009). Both landscape fragmentation and population density were upscaled from 1 km2 to 100 km2 resolution, aggregating cells by median values due to their skewed distributions. Data on street length were extracted from the “Global biodiversity model for policy support” (GLOBIO—Meijer et al., 2018). Total street length was calculated by combining the length of all street types within each cell (highways, primary, secondary, tertiary, and local roads). Street density was calculated dividing the total street length by the area of each cell.

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 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).

Data processing was performed with R 3.6.3 (R Core Team, 2020), using the following packages: glmmulit (Calcagno & de Mazancourt, 2010), gbm3 (Hickey, 2016), ggplot2 (Wickham, 2016), ggeffect (Lüdecke, 2018), plot3D (Soetaert, 2019), tidyverse (Wickham et al., 2019), dplyr (Wickham et al., 2020), lhs (Carnell, 2020), patchwork (Pedersen, 2020), mgcv (Wood, 2020), and GGally (Schloerke et al., 2021).

RESULTS

Among the member countries of the European Union (EU, excluding Greece and including the former member United Kingdom), habitat richness (or habitat type richness; see Methods) per 10 km × 10 km cell (N = 43,240 cells in total) ranged from 0 to 43 (Figure 1a). While normalized habitat richness (i.e., habitat richness corrected for actual cell area) ranged from 0 to 1.5 (Figure 1b).

FIGURE 1

Maps of habitat richness obtained by overlapping habitat distribution maps with the standard grid of 10 km × 10 km cells provided by the European Environmental Agency (EEA) for habitat monitoring. Greece was excluded due to the lack of habitat reporting (lack of habitat reporting delivery of Article 17 data in 2013). Panel a shows the habitat richness. Panel b shows normalized habitat richness

Both maps show a heterogeneous distribution of habitat richness, peaking along the main mountain chains of Southern Europe and Central Europe as well as the Baltic area (Figure 1). Moreover, the habitat richness was positively correlated (r = .34, p < .001) with the richness of the species listed in the Annex species of the Birds and Habitats Directives as reported by every country (see Figure S2).

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.

FIGURE 2

Boxplots showing percentages of relative influence and cumulative contribution of single and grouped explanatory variables in accounting for habitat richness. The contributions of variables were estimated across the GLM, GAM, and BRT models. Panel a: the relative influence of each explanatory variable. Panel b: the cumulative contribution for each group of the explanatory variables selected. 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

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).

FIGURE 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

FIGURE 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

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). The root mean square error (RMSE) indicated large differences among the explained deviance by the GLM and GAM on one side and the BRT on the other side, with the former models showing 0.27 and 0.26, respectively, and the latter 0.85. By adding the autocovariates to the models, GLM and GAM showed substantial improvements (reducing the RMSE to 0.2 in each case), whereas the BRT did not.

TABLE 2

Explained deviance (D 2) and root mean square error (RMSE) of the models with and without autocovariate (“RAC”)

Model	D 2	RMSE
GLM	0.22	0.27
GLM with RAC	0.59	0.20
GAM	0.27	0.26
GAM with RAC	0.59	0.20
GBM	NA	0.85
GBM with RAC	NA	0.85

DISCUSSION

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.

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; MacArthur, 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.

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.

CONCLUSIONS AND IMPLICATIONS FOR CONSERVATION PLANNING

The evaluation of the variables associated with habitat richness here performed for the first time at the European continental scale revealed that geographic and climatic variables are more influential than anthropogenic variables for explaining habitat richness distribution. Despite not showing the strongest correlations, we found that human activities were indeed relevant in controlling the distribution of habitat richness, as the interaction among anthropogenic variables had a strong negative effect on habitat richness. Thus, for environmental management, it is important to consider the cumulative effect of interactions between natural and anthropogenic variables. Indeed, increasing human populations, long‐term land cover changes, and pressures from invasive alien species have all been linked to habitat transformation (Banks‐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. 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).

CONFLICT OF INTEREST

No conflict of interest.

AUTHOR CONTRIBUTIONS

Marco Cervellini: Data curation (lead); Investigation (lead); Methodology (equal); Project administration (lead); Writing – original draft (lead); Writing – review & editing (equal). Michele Di Musciano: Formal analysis (lead); Investigation (lead); Methodology (equal); Writing – original draft (lead); Writing – review & editing (equal). Piero Zannini: Formal analysis (lead); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Simone Fattorini: Investigation (lead); Writing – review & editing (equal). Borja Jiménez‐Alfaro: Investigation (lead); Writing – review & editing (equal). Emiliano Agrillo: Investigation (equal); Writing – review & editing (equal). Fabio Attorre: Investigation (equal); Writing – review & editing (equal). Pierangela Angelini: Investigation (equal); Writing – review & editing (equal). Carl Beierkuhnlein: Investigation (equal); Writing – review & editing (equal). Laura Casella: Investigation (equal); Writing – review & editing (equal). Richard Field: Investigation (equal); Writing – review & editing (equal). Jan‐Christopher Fischer: Investigation (equal); Writing – review & editing (equal). Piero Genovesi: Investigation (equal); Writing – review & editing (equal). Samuel Hoffmann: Investigation (equal); Writing – review & editing (equal). Severin D. H. Irl: Investigation (equal); Writing – review & editing (equal). Juri Nascimbene: Investigation (equal); Writing – review & editing (equal). Duccio Rocchini: Investigation (equal); Writing – review & editing (equal). Manuel Steinbauer: Investigation (equal); Writing – review & editing (equal). Ole R. Vetaas: Investigation (equal); Writing – review & editing (equal). Alessandro Chiarucci: Conceptualization (lead); Investigation (equal); Writing – review & editing (equal).

Figure S1

Click here for additional data file.

Figure S2

Click here for additional data file.

Figure S3

Click here for additional data file.

Figure S4

Click here for additional data file.

Supplementary figure captions

Click here for additional data file.