The Loss of Landscape Ecological Functionality in the Barcelona Province (1956–2009): Could Land-Use History Involve a Legacy for Current Biodiversity?

Abstract

Could past land uses, and the land cover changes carried out, affect the current landscape capacity to maintain biodiversity? If so, knowledge of historical landscapes and their socio-ecological transitions would be useful for sustainable land use planning. We constructed a GIS dataset in 10 × 10 km UTM cells of the province of Barcelona (Catalonia, Spain) for 1956 and 2009 with the changing levels of farming disturbance exerted through the human appropriation of photosynthetic net primary production (HANPP), and a set of landscape ecology metrics to assess the impacts of the corresponding land-use changes. Then, we correlated them with the spatial distribution of total species richness (including vascular plants, amphibians, reptiles, birds and mammals). The results allow us to characterize the main trends in changing landscape patterns and processes, and explore whether a land-use legacy of many complex agroforest mosaics maintained by the intermediate farming disturbance managed in 1956 could still exist, despite the decrease or disappearance of those mosaics before 2009 due to the combined impacts of agroindustrial intensification (meaning higher HANPP levels), forest transition (meaning lower HANPP levels) and urban sprawl. Statistical analysis reveals a positive impact of the number of larger, less disturbed forest patches, where many protected natural sites have been created in 1956–2009. However, it also confirms that this result has not only been driven by conservation policies and that the distribution of species richness is currently correlated with the maintenance of intermediate levels of HANPP. This suggests that both land-sharing and land-sparing approaches to biodiversity conservation may have played a synergistic role owing to the legacy of complex land cover mosaics of former agricultural landscapes that are now under a serious threat.

1. Introduction

Europe has experienced a long-lasting and widespread forest transition on the steep terrain [1,2,3,4,5], whereas urban-industrial facilities and linear infrastructures are taking ever more land on the plains. In 2012 woodland was the largest land use, covering 34% of total continental area, followed by arable land and permanent crops (25%), and pastures and mosaics (17%), while artificially sealed built-up areas and infrastructures occupy around 4% [6]. From 1990, the total forest area has remained relatively stable in Europe, but forests continue to expand in some regions, such as the Mediterranean [7]. There are authors who praise this forest transition for some of its environmental benefits, such as carbon sequestration and lower human pressure over large mountain areas where species richness can refuge [8,9,10], while others have raised concerns about other negative impacts on biodiversity [11,12,13,14,15].

These contrasting views on forest transition are tightly linked to the ongoing debate between the “land sparing” and “land sharing” approaches to biological conservation [16]. The land-sparing view argues that intensifying agroindustrial production in some areas may allow setting others aside to nature conservation [17,18]. However, land-sharing proponents point out that having a set of isolated natural reserves is not enough for biodiversity conservation, whereas spreading a wildlife-friendly farming in a complex land-matrix can maintain a greater farm-associated species richness at the landscape level that compensates for the decline of species richness at the plot level exerted by farming disturbance [19,20,21]. Furthermore, the land cover diversity of complex agroecological land matrices can also provide a much-needed ecological connectivity to prevent the isolation of protected sites, thus helping to increase biodiversity at the regional extent. According to the intermediate disturbance hypothesis, landscape heterogeneity can contribute to a dynamic biodiversity peak at intermediate levels of ecological disturbance, through the interaction among land cover diversity, ecosystem complexity, and dispersal abilities of the colonizing species that escape form the most disturbed patches and are sheltered in less disturbed ones [16,21,22,23,24,25].

Indeed, studies conducted in the Barcelona area have shown that the two approaches can jointly explain current species richness distribution [26,27]. Accordingly, both land-sparing and land-sharing policies can be combined in a sound conservation strategy aimed at enhancing all ecosystem services [28]. Larger patch units kept more undisturbed can provide refuge for populations of wild species that activate their dispersal abilities in areas more disturbed by farming, whereas agroforest landscape mosaics can provide ecological connectivity among refuge areas. Furthermore, from a metapopulation perspective having larger landscape patches kept less disturbed would not only act as sanctuaries, but also as source populations that can recolonise other fragmented and disturbed units, establishing a dynamic network of habitats that help maintain the populations of many species [29,30,31].

However, either opposing or combining them, these broad approaches to nature conservation require more researches on how biodiversity is maintained, improved, or lost through land use changes [32,33,34,35,36,37,38]. Industrialization of agriculture through the “green revolution” has been a major driver of biodiversity loss worldwide [39]. This has been particularly true for farm-associated biodiversity [40], as shown by the collapse of common farmland birds populations experienced all over Europe since the 1960s, which contrasts with the recovery of formerly endangered species of rare specialist birds as a result of the successful nature protection policies developed [41,42,43]. In Catalonia, a decline of bird populations that live in agroforest mosaics has also been reported (http://www.ornitologia.org/ca/), as well as butterflies and Mediterranean orchids that inhabit the open spaces of those landscape mosaics [44,45].

All these examples highlight that agriculture, and the cultural landscapes it produces, can either decrease or increase species richness depending on the way it is practiced [46]. Overcoming the global food-biodiversity dilemma requires advanced knowledge on how species richness is maintained or lost in different land-use patterns, according to the level and the spatiotemporal pattern of ecological disturbances caused by farming [34,47,48]. If we want to ensure ecosystem services in the future, better operative criteria and indicators are needed on when, where, and why the energy throughput driven by farming increases or decreases the mosaic pattern of cultural landscapes in a way that affects its capacity to maintain biodiversity [49]. To address this big societal challenge requires more landscape history research in different bioregions and territories [26,45,50,51].

This calls for multidisciplinary studies of human-nature interaction and coproduction in agroecosystems [46] through socioecological integrated models, like the energy-landscape integrated analysis (ELIA) and the intermediate disturbance-complexity model (IDC) [52,53,54]. Similar to the IDC model, this work relies on the human appropriation of net primary production (HANPP) and a set of landscape ecology metrics as indicators to study to what extent the loss of habitat differentiation in increasingly homogenous landscapes, because of the abandonment of traditional agroforest mosaics, has had negative impacts for biodiversity conservation even considering the positive impacts of creating larger forest units not affected by farming.

To that aim, we assess the land-use changes in the Barcelona province in 1956–2009 from a landscape ecology perspective and test their impact on the current species richness distribution. One underlying hypothesis is the existence of a historical land-use legacy [55,56,57]: Could the mosaic pattern of the former agricultural landscapes that existed in 1956 (i.e., before industrial farming was spread in cultivated flatter areas and forest encroachment took place following cropland abandonment in steeper areas) still explain current locations of species richness in the Barcelona province? What role has the traditional land-sparing conservation policies played in that outcome? This kind of land-use historical legacies of past landscape configurations which are still contributing to biodiversity conservation at present has been demonstrated in other researches that open an interesting field of study [55,56,57]. Section 2 describes the methodology, Section 3 presents the main results, Section 4 discusses them and Section 5 presents the main conclusions of this preliminary study, its limitations and pending tasks, pointing out some hypothesis to follow up.

2. Case Study and Methods

The dependent variable to be explained is the spatial distribution of the total biodiversity recorded in the Biodiversity Data Bank of Catalonia (BDBC; http://biodiver.bio.ub.es/biocat) in the 48 cells of 10 × 10 km of the Catalan province of Barcelona in Spain (Figure 1 and Figure 2). Following the intermediate disturbance-complexity (IDC) modelling [26,50,58] to study the interaction between certain levels of farming appropriation of the biomass produced by the photosynthesis in the study area, and the complexity of land-cover patterns of cultural landscapes to assess their capacity to maintain species richness [58], we consider the following explanatory variables: 1) The different amount of farming disturbance exerted on the territory through the human appropriation of net primary production (HANPP); and 2) landscape ecology metrics that assess the capacity of the land cover patterns and ecological processes to provide differentiated habitats for biodiversity maintenance: Shannon-Wiener index, polygon density, edge density, largest patch index, effective mesh size, and ecological connectivity index.

All these GIS data were accounted for each of the 48 UTM 10-km cells within the borders of the Barcelona province (Figure 1 and Figure 2). The landscape structure of these cells was taken from the 2009 Land Cover Map of Catalonia (Spain), reclassified into a set of principal land cover categories to compare with the 1956 land cover map digitised and photointerpreted from the aerial photograph made by the US Army (www.creaf.uab.es/mcsc/). We also used the land cover map of 1993 digitised from the satellite image, although the lower quality of land cover photointerpretation prevents us from using it in most statistical analysis.

Then we examined whether the farming-driven land use and cover change (LUCC) experienced in the province of Barcelona from 1956 to 2009 can explain the current distribution of species richness of vascular plants, amphibians, reptiles, birds, and mammals registered by the BDBC. We mapped through GIS the total species richness accounted in the 48 cells of 10 × 10 km (Figure 1), to perform statistical analyses on the relationships between the variations of HANPP and landscape ecology metrics in 1956–2009 with the locations of total species richness currently observed in each cell [59]. Finally, in order to assess the role played by nature conservation policies and their interaction with the ongoing land use changes, a specific correlation analysis was made for the same study area.

2.1. HANPP as a Measure of Farming Ecology Disturbance

The ecological disturbance exerted by farming is accounted by means of the human appropriation of net primary production (HANPP), calculated according to the standard method set forth by Haberl et al. [60]:

$$HANPP=\frac{HANP{P}_{harv}+HANP{P}_{LUC}}{NP{P}_0}$$

where HANPP_harv is the photosynthetic net primary production (NPP) appropriated through harvest, HANPP_luc is the change in the NPP due to farming-induced land use changes, and NPP₀ is the potential NPP. HANPP_luc is defined as the difference between the NPP of the potential (NPP₀) and actual (NPP_act) vegetation cover (this is, the whole net biomass produced within a year):

$$HANP{P}_{luc}=NP{P}_0-NP{P}_{act}$$

Thus, HANPP was calculated for each land-cover category and UTM cell of 10 × 10 km as the weighted sum of fixed coefficients multiplied by the proportion of the land area occupied by each land cover type. To estimate these HANPP values it is necessary to assess the different levels of NPP annually produced in each land type in the study area, and the amounts of biomass harvested by farming. NPP₀ values have been taken from Haberl et al. [60]. Harvest ratios appropriated from each land type were assessed according to the agricultural statistical sources available for Barcelona province in each time point, transformed into carbon using the conversion factors, the residue/product losses, and the unharvested biomass ratios given by Guzmán et al. [61].

2.2. Landscape Ecology Metrics to Assess Land Cover Diversity and Fragmentation

Several metrics were used to differentiate the (positive) effect of land cover diversity from the (negative) effect of land cover fragmentation on biodiversity [62,63,64]. Landscape heterogeneity was assessed through the Shannon–Wiener index (H’) used as a measure of land cover equi-diversity in each UTM cell:

$$H^{\prime }=-{{\displaystyle \sum }}^{}\left(p_i·ln p_i\right)$$

where p_i is the proportion of the land unit occupied by each type of vegetal cover.

The largest patch index (LPI) measures the area of the largest polygon in each cell. The polygon density (PD) indicates the number of polygons in each cell. Edge density (ED) is the sum of the polygon perimeters in each cell. Effective mesh size (MESH) is the sum of the areas of the polygons squared, divided by the size of the study area, as a measure of the inverse of fragmentation.

$$MESH=\frac{{{\displaystyle \sum }}^{}{A}_i^2}{{{\displaystyle \sum }}^{}{A}_i}$$

where A_i is the area of each polygon (measured in square km).

2.3. Landscape Ecology Metrics to Assess Ecological Connectivity

After having assessed the landscape structure (patterns) with the above set of metrics, we calculated another index to account for the landscape ecological functionality (processes). The ecological connectivity index (ECI) assesses the functionality of the land matrix according to its capacity to maintain the horizontal flow of energy, matter, and information that sustain biodiversity through the ecosystems’ landscape patterns and trophic chains [65]. The ECI used is based on a simplification of the original methodology proposed by Marull and Mallarach [63] following Lindenmayer and Fischer [66]. It relies on defining a set of ecological functional areas (EFA) considered as focal habitat patches to be connected, and a GIS computational model of cost-distance of displacement that includes the effect of modelled anthropogenic barriers (urban areas and infrastructures) considering the type of barrier, the range of distances, and the kind of land-use involved. Then the model applies the CostDistance function using two databases: a “source” surface for each type of EFA and an “impedance” surface resulting from applying the effects of barriers to the potential affinity matrix. The result is a cost-distance adapted to each type of EFA. By calculating the value of the sums of cost-distances adapted, this computational model of ecological connectivity brings about a normalized range that varies from zero to ten:

ECI = 10 – 9 [ln (1 + x_i)/ln (1 + x_t)³]

where x_i is the value of the sum of the cost-distance by pixel and x_t the maximum theoretical cost distance. This ECI helps to emphasize the role played by all sorts of landscape units (forest, pasture, and cropland) to keep up ecological connectivity [67].

All the above indicators and indices were accounted to test the variation in the impacts of farming disturbance and landscape ecology patterns and processes in 1956 and 2009 on the spatial distribution of the total species richness observed in 2009, in order to check whether their current locations maintain the land-use legacy of former farming disturbance through the agroforest landscape mosaics.

2.4. Statistical Analysis of the Relationship between Farming Disturbance and Landscape Ecology Metrics in 1956 and 2009 with Total Species Richness in 2009

We first identified the highest differences among the values of each variable in the two time points by performing two-tailed tests using the Bonferroni p-values adjustment for all pairwise comparisons, with equal variances with a significance level of 0.05. Then we performed a negative binomial regression analysis [59], as our data are over-dispersed count variables with a Poisson distribution. The negative binomial regression analysis was carried out for the total species richness measured in 2009 as dependent variable, considering as explanatory variables the variation in HANPP and the values of land-cover metrics (H’, LPI, ED, EMS, PD, and ECI) in 1956 and 2009. By regressing the differences on the values for HANPP and land-cover metrics between 1956 and 2009 we can grasp how certain dynamics on landscape patterns and processes carried out in each cell of the grid affected the current spatial distribution of total species richness.

Despite data in 1993 are useful to show the main trends on the overall land uses, and to analyse some landscape processes and changes in HANPP values, the procedure used to build the map did not allow for comparing landscape ecology metrics with the detail reached in the maps for 1956 and 2009. Therefore, we limited to 1956 and 2009 the landscape ecology assessment of most of these historical land-use changes, and the statistical analysis of these indicators with biodiversity data.

Finally, based on the data of protected areas through conservation, and their spatial overlapping with the other variables assessed in the 10 × 10 km² cells, we analysed the Pearson correlation of these natural sites with the significant variables of the previous statistical analysis.

3. Results

The spatial distribution in the 48 cells of 10 × 10 km of the Barcelona province of the total species richness registered by the BDBC is shown in Figure 1.

Table 1. Land use and cover change (LUCC) in the Barcelona province (1956, 1993, and 2009).

Land Covers	1956		1993			2009
	ha	%	ha	%	1956 = 100	Ha	%	1956 = 100
Forest	195,526.4	40.7	223,060.6	46.5	114.1	233,357.3	48.6	119.3
Scrubland and pastures	88,710.9	18.5	69,751.0	14.5	78.6	68,691.5	14.3	77.4
River corridor and wetlands	1,759.7	0.4	1,755.3	0.4	99.7	2,097.0	0.4	119.2
Cropland	173,140.0	36.1	135,288.2	28.2	78.1	104,359.8	21.7	60.3
Unproductive	8,246.7	1.7	10,688.9	2.2	129.6	10,106.3	2.1	122.5
Road and rail networks	2,246.6	0.5	3,509.4	0.7	156.2	7,487.2	1.6	333.3
Urban area	10,369.7	2.2	35,946.7	7.5	346.7	53,900.8	11.2	519.8
Total	480,000.0	100.00	480,000.2	100.0		480,000.0	100.0

Source: Our own accounted with GIS from the digital land cover maps explained in the text.

presents the data obtained on the land-use and cover changes experienced in the Barcelona province by accounting with GIS the digital maps of 1956, 1993, and 2009. These maps are shown on the left side of Figure 2.

The average values of the HANPP indicators and landscape ecology metrics in the whole study area were accounted with GIS in the digital land cover maps seen on the left side of Figure 1, in most cases only for 1956 and 2009 (

Table 2. Comparative analysis of the landscape ecology metrics applied in the Barcelona province in all the 10 × 10 km² sample cells in 1956, 1993, and 2009.

Landscape Ecology Metric and HANPP Values	1956 (A)	1993 (B)	2009 (C)
Polygon Density—PD (number of polygons)	3,081.19 -	-	3,786.98 A
Edge Density—ED (km)	31.51 -	-	31.58 -
Largest Patch Index—LPI (ha)	1,964.82 -	-	2,096.68 -
Effective Mesh Size—EMS (km2)	256.49 C	-	121.51 -
Shannon-Wiener Index (H’)	0.589 C	-	0.504 -
Ecological Connectivity Index—ECI	6.81 BC	5.65 -	5.23 -
Net Primary Production actual—NPPact (TM C year−1 ha−1)	69.08 -	92.69 A	86.04 A
Net Primary Production harvested—NPPharv (TM C year−1 ha−1)	24.97 -	35.82 AC	25.73 -
Human Appropriation of NPP—HANPP (%)	60.87 BC	48.78 -	50.42 -

Note: Results of two-tailed tests with equal variances with a significance level of 0.05. For each significant pair, the key under the category (A, B, C) indicates when a value of the variables considered is, for the specific period shown in a column, statistically different from the other two dates shown in the other columns. These results appear below the category, and they have been adjusted for all pairwise comparisons using the Bonferroni test. For 1993 only HANPP values and indicators of landscape processes at the cell level were considered. Source: Our own, calculated from the data and methods explained in the text.

). This data includes the pairwise statistical results obtained by comparing each variable in 1956, 1993, and 2009 through two-tailed tests, in order to assess which values had the highest differences in these years.

Negative binomial regression analysis [60] of the capacity to explain the distribution in the 48 UTM cells of 10 × 10 km of the total species richness observed in the province of Barcelona in 2009 was done, considering as independent variables the historical variation from 1956 to 2009 in the levels of farming disturbance exerted through HANPP, and the land cover patterns and processes assessed through all the landscape ecology metrics accounted. The incidence rate ratio (IRR) value informs on how a change in the independent variables (landscape patterns and processes) affects the dependent variable (total biodiversity). Values over 1 show positive relation while values under 1 shows negative relation (

Table 3. Negative binomial regression analysis for the total biodiversity (2009), considering the variation (Δ) in human appropriation of net primary production (HANPP) and land-cover metrics listed in Table 2 (H’, LPI, ED, EMS, PD, and ECI) in all the10 × 10 km sample cells of the Barcelona province in 1956–2009.

Total Biodiversity	IRR	Std. Err.	Z	P > |z|
Δ HANPP	1.132	0.064	2.20	0.028 (**)
Δ H’	0.802	0.047	−3.73	0.000 (*)
Δ LPI	1.127	0.055	2.45	0.014 (**)
Δ PD	1.197	0.130	1.65	0.099
Δ ED	0.996	0.151	−0.03	0.977
Δ EMS	1.058	0.038	1.60	0.111
Δ ECI	0.990	0.045	−0.22	0.827
cons.	813.297	27.236	200.10	0.000

Note: IRR (incidence rate ratio) is the estimated rate ratio for a one unit increase in a standardized test score, given that the other variables are held constant. * = level of significance at 10%, ** = level of significance at 5%, *** = level of significance at 1%; likelihood-ratio test of alpha = 0; chibar2 (01) = 3902.05; Prob ≥ chibar2 = 0.000. Source: Our own, calculated from the data and methods explained in the text.

).

In order to better understand how the debate between the two main conservation policies fit our previous results, we performed a correlation analysis between the total surface areas under nature protected figures and the statistically significant variables obtained in the negative binomial regression analysis (∆HANPP, ∆H’, and ∆LPI) in their variation from 1956 to 2009 (

Table 4. Pearson correlation between total surface area under nature conservation (m²), and the variation (Δ) in human appropriation of net primary production (HANPP) and land-cover metrics (H’, LPI), in all the 10 × 10 km sample cells of the Barcelona province in 1956–2009.

Surface Area under Conservation Figures	Pearson Correlation	P > |z|
Δ HANPP	−0.029	0.845
Δ H’	−0.195	0.184
Δ LPI	0.370	0.010

Source: Our own, calculated from the data and methods explained in the text.

).

4. Discussion

4.1. Land Use and Cover Change (LUCC) from 1956 to 2009

From 1956 to 2009 the province of Barcelona underwent the three main land use and cover changes (LUCC) experienced all over Europe [1,2,3,4,5,6]: simplification of agricultural landscapes, forest encroachment, and urban sprawl (Table 1 and Figure 2). Cropland abandonment in 40% of the cultivated area in 1956 (68,780 ha lost up to 2009), together with 23% reduction of pastureland and scrub (20,019 ha lost), has resulted in 19% of forestland growth (37,831 ha more) and 387% increase in the land taken by urban development and infrastructures (48,772 ha more). In addition to the 40% reduction, cropland area also changed its composition by losing many complex landscape mosaics because of the disintegration among cropping, livestock rearing, and forestry [12,13,14,54]. This feature is better assessed through the landscape ecology metrics (Table 2) discussed in the next Section 4.2 and Section 4.3.

Urban sprawl has been mainly concentrated in the flatter areas of the Barcelona Metropolitan Region (BMR) near the coast, at the expense of cultivated areas in some of the most fertile soils, while cropland and pastureland abandonment led to forest encroachment mainly in steeper lands of mountain areas further from the coast, as shown on the left side of Figure 2. On the right side, the ecological connectivity maps resulting from these land cover changes show a general reduction of ECI values particularly (but not only) in the BMR. Rural areas in the Northern part of the Barcelona province have also experienced a decrease in ECI values, showing that even moderate expansions of built-up areas and infrastructures may have relevant negative impacts on the landscape ecology connectivity through their barrier effect [11,32,67].

4.2. Trends in the Net Primary Production and its Human Appropriation from 1956 to 2009

The average results in the 48 UTM cells of all the variables accounted in the Barcelona province are the result of a series of LUCC changes that in many cases have experienced opposite trends (Table 1 and Figure 2). However, despite these contrasting changes the overall outcome is meaningful. HANPP values decreased from 61% of photosynthetic biomass production (NPP) in 1956 to 49–50% in 1993–2009, meaning that cropland and grazing abandonment of steep areas has predominated over the land-use intensification carried out in flat lands by industrial farming and livestock fattening in feedlots. There has also been a rise in the NPP from 69 to 93 TM C/ha/year from 1856 to 1993, followed by a small decrease up to 86 TM C/ha/year in 2009, which largely explains the overall HANPP reduction despite the greater amount of harvest intensity in the land that remained cultivated.

The average NPP increase was, in turn, a combined result of spontaneous afforestation following rural abandonment and cropland intensification through synthetic fertilization, irrigation, and use of high-yielding industrial seeds. While the biomass harvested in newly forested areas declined, it largely increased in industrial cultivated areas, leading to a slight increase in the average NPP harvested from 25 TM C/ha/year in 1956 to 36 TM C/ha/year in 1993, followed by a subsequent decrease to 26 TM C/ha/year in 2009. The land occupied by urban-industrial development at the expense of intensively cultivated flat lands also contributed to the last decreases experienced both in the photosynthetic NPP production and harvested biomass from 1993 to 2009.

4.3. The Loss of Landscape Complex Mosaics from 1956 to 2009

The Shannon-Wiener Indices (H’) of land cover equi-diversity in the 48 UTM cells show a net decrease from 0.59 in 1956 to 0.50 in 2009, which means a reduction of landscape heterogeneity. How this LUCC homogenization has been combined with landscape fragmentation is a more complicated issue. Largest patch indices (LPI) have slightly grown from 1965 ha in 1956 to 2097 in 2009. This mainly expresses the greater extent of larger forestland units created by forest regrowth, although the greater extent of crop monocultures also contributed to this. In fact, when selecting only those cells already dominated by forest in 1956, the average increase in LPI has been 245 ha, while the rest decreased by 72 hectares on average. These LPI indicators, combined with lower Shannon-Wiener Indices, point to a decrease or disappearance of former agroforest landscape mosaics. The steady decrease in the mean values of Effective Mesh Size (EMS) and the increase in Polygon Density (PD) throughout the period corroborate the growing landscape fragmentation, which took place at the same time as a general trend toward land cover homogenization.

The loss of agroforest land-use mosaics replaced by larger and more homogeneous land covers has gone hand in hand with a greater landscape fragmentation mainly driven by urban sprawl and the spread of linear infrastructures that exert a strong barrier effect. This LUCC combination has been detrimental to ecological connectivity. The normalized ECI values, which range from zero to ten, decreased 23% from 6.8 in 1956 to 5.7 in 1993 and 5.2 in 2009 on average.

All in all, HANPP values and landscape metrics taken together confirm that the combination of land-use intensification and fragmentation of some flat areas with land-use abandonment in the rest of the territory has led to a loss of complex agroecological land cover mosaics and the creation of larger forestland patches in stepper lands in a similar way as experienced in the rest of the Mediterranean, and all over Europe as well [32,33,36].

4.4. Statistical Impact on Total Species Richness in 2009 of Landscape Homogenization and the Increase of Larger Forest Units from 1956 to 2009

The negative binomial regression analysis provides the more relevant answer to the research question raised in this article—i.e., whether there has been a land-use legacy either of former agroforest mosaics currently disappearing (to test for the land-sharing approach), and of the increasing presence of larger forestland units (to test for the land-sparing approach) on the current species richness distribution in the Barcelona province. The explanatory variables for total biodiversity are the change in HANPP values (farming disturbance) and in the landscape ecology metrics (Table 3). The results give clues in two directions.

HANPP results show that the increase of farming disturbance has had a positive effect on the distribution of total species richness in 2009. The average change of HANPP in the whole study area from 1956 to 2009 is −10.6%, and the lowest biodiversity observed in 2009 is statistically associated with those areas of lower values of HANPP, where afforestation has been the usual LUCC. This implies a land-sharing approach as an explanation.

At the same time, however, the increase of LPI values in 1956–2009 is also significantly correlated with the spatial distribution of total biodiversity observed in 2009 (Table 3), while there is a negative significant correlation with the changes in spatial land cover equi-diversity H’. These results mean that there is a problem of landscape fragmentation and confirm the positive contribution provided by larger forest areas to help maintain high species richness. All this implies the land-sparing approach to nature conservation as an explanation.

Therefore, we can conclude that both land-sharing and land-sparing LUCC changes from 1956 to 2009 have contributed to the spatial distribution of total species richness observed in 2009. This suggests that they should be combined into future policies of nature conservation, a result similar to those obtained by applying ELIA and IDC models to Catalonia and the province of Barcelona [26,27,44,45].

However, these are preliminary results that are limited by the lack of biodiversity data in 1956, and by the coarse grid of 10 × 10 km cells we had to use because of the available data of species richness for 2009. The existence of metapopulation dynamics of many species could also affect the results obtained, and would require performing in the future the same type of analysis at different scales using finer data on specific species.

4.5. The Impact of Past Conservation Policies on Landscape Functioning

The nature conservation policies implemented through nature protected sites have only affected significantly the increases on LPI values, while they have had no impact in changing HANPP and land cover H’ values (Table 4). This is particularly important, because it shows that these traditional land-sparing policies do not guarantee adequate conditions for biodiversity conservation in the whole land matrix in the long term [58]. Land-sharing improvements are also needed, mainly through more integrated and sustainable agriculture, livestock rearing, and forestry.

5. Conclusions

After constructing a GIS dataset to assess from a landscape ecology point of view the main LUCC in 48 UTM cells of 10 × 10 km in the Barcelona province from 1956 to 2009, we obtained two relevant set of statistical results by testing to what extent the spatial distribution of total species richness currently observed correlates with the changes experienced from 1956 to 2009 in the levels of farming disturbance (HANPP), land cover equi-diversity (H’), grain size (LPI), fragmentation (EMS), ecotony (PD, ED), and ecological connectivity (ECI). On the one hand, these results have checked the positive role for biodiversity maintenance played by traditional farming that combined intermediate levels of HANPP with land cover diversity (H’) giving rise to heterogeneous landscapes where agroforest mosaics predominated. However, this result contrasts with the significant positive contribution provided by larger afforested areas following cropland abandonment in mountain areas, where many natural protected sites were created during 1956–2009. These contrasting results call for the need to combine land-sharing and land-sparing approaches to nature conservation, and to reinforce the former once the latter seems to have reached a limit either in its functional role or its spatial extent.

On the other hand, the results suggest the existence of a historical land-use legacy of former agricultural landscape mosaics that might still affect the spatial species richness distribution in 2009, despite the worrying decreases these mosaics have experienced with the land-use changes from 1956 to 2009. Although preliminary, they seem to suggest that the location of species richness across the landscape preserves to some extent a memory of the former agroforestry mosaics that existed in 1956 and have been decreasing or even disappearing.

We acknowledge that this hypothesis requires more research based on better biodiversity data analysed with a higher spatial resolution. But our preliminary results allow us to conclude that this type of research deserves to be developed in future research. We have offered some evidence on the main variables that can explain this legacy and need to be carefully analysed, as well as on the existence of complementarity between land-sparing and land-sharing improvements in biodiversity conservation. Metapopulation dynamics could also help explaining the synergistic result of combining larger landscape patches that remain less disturbed, and capable of offering refuge areas, together with more heterogeneous land use mosaics with higher levels of anthropogenic disturbances but able to maintain high ecological connectivity.

The possibility that there exist historical land-use legacies in the coevolution between cultural landscapes and species richness has become a very important issue for biodiversity conservation policies. If so, promoting restoration of agroforest mosaics in more complex agroecological landscapes would enhance farm-associated biodiversity by recovering those habitats previously fragmented or lost with the negative impacts of industrial farming, withdrawal of extensive livestock rearing, and abandonment of multifunctional use of forests. Therefore, wildlife-friendly farming could contribute to biodiversity conservation by supplementing nature-protected areas mainly created in isolated forest remnants. All this requires advances in the landscape agroecology research here explored, to study at different scales, and with higher spatial resolution, the existing relationships between farming management and farm-associated biodiversity considering different taxa and specific species.