This is an Accepted Manuscript of the article Guadilla-Sáez, Sara, et al. "Biodiversity conservation effectiveness provided by a protection status in temperate forest commons of north Spain." Forest Ecology and Management 433 (2019): 656-666, available online at: https://doi.org/10.1016/j.foreco.2018.11.040. ©2018. This manuscript version is made available under the CC-BY-NC-ND 4.0 license. Biodiversity conservation effectiveness provided by a protection status in temperate forest commons of north Spain Sara Guadilla-Sáeza,*, Manuel Pardo-de-Santayanab, Victoria Reyes-Garcíaa,c, Jens-Christian Svenningd,e a Institut de Ciència i Tecnologia Ambientals, Universitat Autònoma de Barcelona, Barcelona, Spain b Departament of Biology (Botany), Universidad Autónoma de Madrid, Madrid, Spain c Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain d Center for Biodiversity Dynamics in a Changing World (BIOCHANGE), Department of Bioscience, Aarhus University, Aarhus, Denmark e Section for Ecoinformatics and Biodiversity, Department of Bioscience, Aarhus University, Aarhus, Denmark *Corresponding author. E-mail address: sara.guadilla@uab.cat Abstract: The establishment and maintenance of protected areas is the backbone of global conservation strategies to halt biodiversity loss. However, despite the more than 200,000 legally designated protected sites worldwide, the rate of species extinction has not decreased, for which some debate the real effectiveness of protected areas to preserve biodiversity. Using data from tropical areas, many studies have attempted to test the effectiveness of protected areas by comparing species richness in protected and neighbouring unprotected sites, without reaching a consensus. Here, we extend this line of research with data from temperate deciduous forests inside and outside Picos de Europa National Park and Biosphere Reserve (N Spain). Specifically we compare data from mixed broadleaved woodlands, beech forests (Fagus sylvatica L.) and Pyrenean oak (Quercus pyrenaica Willd.) forests. We conducted botanical inventories and recorded ecological data from 25 0.2-ha concentric plots distributed in forest commons inside the reserve and from other 25 similar plots established in neighbouring not protected forest commons. Data were used to construct a set of ecological indicators and evaluated using modelling methods. We found no significant differences in species composition between plots in protected and non- protected forest commons, likely due to the similar management criteria applied in both land uses. We found less active management outside the protected area, which helps to maintain stands in a semi-natural state. In contrast, we observed the presence of silvicultural treatments inside the protected area, although these treatments were non-intensive, promoting vegetation composition associated to late-successional ecosystems. We only detected significant differences between plots inside and outside the protected area when relation between species richness was analysed with reference to forest habitat type. Precisely, plots of beech forests inside Picos de Europa were more homogenous than plots outside the protected area, which may indicate that management practices inside the protected area do not favour tree species diversity. Non-intensive silviculture management in beech forests inside Picos de Europa seems to promote the presence of the dominant tree species Fagus sylvatica L., which in the absence of perturbations is characterized by conforming monospecific vegetation communities. Overall, our results do not support the idea that protected areas hold more biodiversity than surrounding forest commons. Conservation treatments applied in protected areas should promote the presence of species associated to disturbances, particularly in stands tending to homogeneous species composition at late- successional stages, as this may enhance their resilience under the current rapid global changes. Keywords: Anthropogenic disturbances; biodiversity indices; protected area; species richness; temperate deciduous forests. 1. Introduction With 234,468 Protected areas (PA) already established worldwide (IUCN, UNEP-WCMC, 2017), site designation rate is considered one of the most remarkable conservation successes of the twentieth century (Gaston et al., 2008). If this trend continues, the goal of 17% coverage for terrestrial and inland waters by 2020 under Aichi Biodiversity Target 11 would be achieved (Gannon et al., 2017). However, in numerous cases, site designation has followed an exclusionary approach –referred to as ‘neo-protectionist’ or ‘fortress conservation’ approach (Wilshusen et al., 2002; Brockington, 2002)– resulting in the displacement and dispossession of communities residing in the newly protected site, and often led to contested actions (Laudati, 2010; Mahapatra et al., 2015). Moreover, even when less controversial ‘conservation-centric’ initiatives have been applied to the establishment of a PA, such as protected sites created in partnership with local people, conflicts with local residents have still arisen, especially when the establishment of a PA have prevented local users from the management of the surrounding natural resources (West et al., 2006). As the economic basis of many indigenous peoples and local communities is closely dependant on the goods obtained from neighbouring natural areas, particularly forests (Angelsen et al., 2014), when restrictions to local use are applied without providing suitable alternative livelihood options, struggles are likely to appear (Mahapatra et al., 2015), constituting a significant shortcoming to PA conservation efforts (Andrade and Rhodes, 2012). There is evidence of better biodiversity conservation outcomes from PA management strategies integrating local economic activities than from strictly conservation PA management regimes (Oldekop et al., 2016), for which many conservationists nowadays embrace a more integrative perspective for the establishment and management of PA (Shultis and Heffer, 2016). Notions such as social-ecological systems, resilience, or cultural landscape are holistic approaches which 2 consider people as part of their surrounding environment rather than mere passive users of landscape biophysical components. People modify their living and adjacent territories, sometimes causing the depletion of natural resources, but sometimes coevolving with nature in the benefit of a sustainable use and promotion of a dynamic mosaic of ecosystems at the landscape level. Along these lines, scholars increasingly advocate for considering PAs as 'social-ecological systems' (Hirschnitz-Garbers and Stoll-Kleemann, 2011; Cumming et al., 2015; Mathevet et al., 2016), a conceptual framework that contemplates both the social and ecological aspects of the system as equally important (Berkes, 2017). Under the social-ecological systems approach, the social and the natural systems are indeed coupled subsystems that co-evolve, which implies that societies are able to adapt to perturbations in the environment and vice versa. Given the rapid global changes occurring nowadays, the co-evolving capacity offered by the social-ecological systems approach brings to light the idea that societies have the opportunity to face environmental challenges without compromising long-term sustainability of ecosystems (Berkes et al., 2003). The growing scholarly emphasis on conservation efforts outside the physical boundaries of PAs focuses on reconcile management practices from land uses, such as farming or forestry, with biodiversity conservation (Kremen, 2015). Research on practices that may be both favourable to biodiversity and economically profitable has significantly increased attention to the potential that community-based resource management may bring to foster a sustainable use of ecological resources (Xu and Melick, 2007; Larson et al., 2016). Through practices such as clearing, livestock grazing, or swidden agriculture, humans have modified landscapes for millennia, partially replacing the ecological functions that megaherbivores used to play in shaping vegetation structures of terrestrial biomes (Sandom et al. 2014; Bocherens, 2018). While not all of these small and intermediate-scale disturbances have a positive effect in preserving biodiversity, overall, management practices applied by communities seem to enhance biodiversity through the creation of a mosaic of ecosystems (Agrawal and Gibson, 1999; Guèze et al., 2015). Aiming to further understand the conservation outcomes that may result from community resource management, here we study a community-based regime that can be considered a social- ecological systems: forest commons. Forest commons are characterized by having clearly defined boundaries and legal enforceable property rights and by providing resources to a variety of social groups that are usually involved in their management (Aryeetey et al., 2012). The importance of forest resources to support rural household economies strongly engaged local communities in the monitoring of their surrounding woodlands to prevent mismanagement and overexploitation, resulting in the implementation of management techniques that allowed the co- existence and long-term maintenance of diverse forest uses and habitats (Parrotta and Trosper, 2012; Kirby and Watkins, 2015). Nowadays, forests commons constitute 18% of global forest area and appear to contribute significantly to biodiversity conservation (Chhatre and Agrawal, 2008). In this work, we analyse the role of forest commons in long-term biodiversity persistence and the effectiveness of protected areas in this type of community-based regime. We do so by comparing a set of ecological and anthropogenic features observed in forest commons plots inside and 3 outside a PA classified as an IUCN category II (National Park), a very restrictive protection category with regard to human activities (Gray et al., 2016; Hewitt et al., 2016). Our study has three main goals. First, we test whether plots inside and outside the PA differ in their ecological characteristics (i.e. topography, edaphic factors), and how this relates to species richness. We based our null hypothesis on the general assumption that we will find the same tree species abundance and evenness in plots inside and outside the PA. The second goal of this research is to analyse the effect of the human intervention on species distribution in plots under protected and unprotected sites, an analysis performed by linking the variables measuring human disturbances (i.e. plot isolation, silvicultural systems) with tree species occurrence. Based on previous case studies highlighting the association of anthropogenic disturbance with species richness in human-dominated landscape (Guèze et al., 2015; Mod et al., 2016), our hypothesis is that anthropogenic disturbances can induce changes in species composition that would result in more heterogeneous species assemblages of the studied forest communities. Finally, the third goal of this paper is to study the relationship between tree species composition, and particularly species diversity, and a) spatial distribution and b) forest management approach. For this goal we compare three different forest habitat types occurring inside and outside the study PA, in an attempt to quantify the conservation outcomes resulting from the protection status between habitats of the same land use type, i.e. temperate deciduous forests. 2. Methods 2.1 Study site We conducted the study in the Liébana valley (57,500 ha), a wide depression located in the southwest of the Cantabria region, in northern Spain. High elevation differences characterize the region, with altitude ranges from 330 m to 2600 m above sea level. Liébana is surrounded by hills, a geographical feature that results in relative geographical isolation with the neighbouring areas and a high number of habitat types of unique ecological value. The bottom of the valley presents a Mediterranean microclimate with less rainfall than the rest of the Cantabria region; at higher altitudes we find an Atlantic climate (Rescia et al., 2008). Liébana’s mean temperature varies from 7.9 to 20.8 C and the annual average rainfall varies from 700 to 1500 mm in the mountainous parts (ETSIM, 1978). Topographic and climatic differences result in a very heterogeneous landscape with a wide array of vegetation types. From a land tenure perspective, common property is a key feature in the Liébana region, with almost 80% of its territory under this regime (Arbeo, 2012). More than three quarters of the woodlands in Liébana are forest commons (Ministerio de Medio Ambiente, 1997–2007). The area has a long tradition of human intervention, as reflected in the large number of local ordinances regulating forest uses since the 15th century (Pérez-Bustamante and Baró, 1988). Management activities have shaped the structure of the local ecosystems and their species composition, at the same time that they have allowed the persistence of high valuable habitats, to the point that 60% 4 of Liébana’s forest commons are currently under some category of protection, including the Picos de Europa National Park (hereafter Picos de Europa) in the northwest of the region. Picos de Europa was the first Park designated in Spain, in 1918, following the North American conservationist model that advocated for the preservation of wilderness areas by preventing human interventions. It also applied a State-led management (Hirschnitz-Garbers and Stoll- Kleemann, 2011; González, 2015). Since its designation, conflicts with local communities in Liébana’s neighbouring region arose due to the limitations that the protection status enforced on local uses (e.g., hunting, wood extraction) (Voth, 2007). Conflicts have shaped the negative perception of Liébana’s rural communities regarding the extension of Picos de Europa National Park in 1995, which included three municipalities of Liébana, conforming the second largest site of the Spanish National Parks and notably including human settlements within its boundaries. Although the management criteria established in 1995 aimed at making compatible local traditional uses –like livestock herding or fuelwood collection– with biodiversity conservation (Royal Decree No. 640/1994 of Spanish Government), confrontation with local population, particularly livestock farmers, resulted in the revoke of the National Park management plan in 2005 due to a legal action taken by local communities within Picos de Europa buffer zone (Spanish Government, 2005). The lack of a valid management plan, specific to the area, difficults the regulation of traditional uses. However, and despite the lack of a specific plan, several conservation initiatives are being undertaken inside the national park in response to national and regional environmental legislation, including the recovery plans for the endangered species brown bear (Ursus arctos) and Cantabrian capercaillie (Tetrao urogallus cantabricus) conservation. Altogether, this situation makes of Picos de Europa a very suitable site to analyse the effects of human activities and conservation management strategies on species diversity. 2.2 Description of the forest communities studied To test the ecological differences between stands inside and outside the PA, we studied the ecological features of the three more abundant habitats in the forest commons inside Picos de Europa, which are also present outside the PA: 1) Mixed broadleaf woodlands, defined as forest with a variable mixture of at least two native broadleaf species accounting for ≥70% of the plot forest cover; 2) Beech forests or forest with Fagus sylvatica L. as the dominant tree species, accounting for ≥70% of the plot forest cover; and 3) Pyrenean oak forests, or forest with Quercus pyrenaica Willd. as the dominant tree species accounting for ≥70% of the plot forest cover. 2.2.1 Mixed broadleaf woodlands Mixtures of broadleaf species cover most of the territory. Despite some human influence, mixed broadleaf woodlands are considered naturally distributed particularly in areas of contact between 5 different tree formations. The predominant species in the mixed broadleaf woodlands located in upper altitudes of Liébana are Quercus pyrenaica, Fagus sylvatica, Ilex aquifolium L., and Crataegus monogyna Jacq. At lower altitudes, the predominant tree species are conformed by Castanea sativa Miller, Fraxinus excelsior L., Tilia cordata L., Tilia platyphyllos Scop., Quercus ilex L., and Corylus avellana L. (ETSIM, 1978). 2.2.2 Beech forests (Fagus sylvatica L.) Forests of Fagus sylvatica are the predominant tree species formation in the Liébana valley and, generally, also the ones with the best conservation status and regeneration rates. These woodlands are shade-tolerant and occupy lands from 600 m to 1700 m in shadow slopes, and from 700 to 1300 m in sunny slopes. Beech woods form close, dense stands, where only shade- tolerant species can grow and where the competing tree species need to take advantage of clearings resulting from felling, browsing animals, or fires. As a result, F. sylvatica only appears in combination with other species in boundary areas with other forest types, or in areas where it displaces other species, as in the case of Quercus pyrenaica. A sparse understorey could accompany beech forests, composed by Ilex aquifolium, Crataegus monogyna, Sorbus aria L. Crantz, and Corylus avellana. However, most frequently, the accompanying species are Pteridium aquilinum (L.) Kuhn in Kerst and Rubus sp., partly due to the increased light following human interventions such as clear cuts (ETSIM, 1978; Godefroid et al., 2005; Kelemen et al., 2012). 2.2.3 Pyrenean oak (Quercus pyrenaica Willd.) Forests of Quercus pyrenaica occur widely in the Liébana valley, with the species also occurring as a shrub in combination with Corylus avellana, Crataegus monogyna, Erica and Ulex sp. in mixed broadleaf woodlands. Q. pyrenaica is tolerant to a wide range of site conditions, occurring from 400 to 1300 m in altitude on sunny slopes and from 300 to 1200 on shady ones. Despite its frequency as shrub strata, it appears most often as a tree when the forest is cleared, i.e. when human intervention is high (Tárrega et al., 2006), accompanied by Corylus avellana, Crataegus monogyna, Erica arborea L., Daboecia cantabrica (Hunds) C. Koch and Rubus sp. Occasionally, it also appears in combination with Prunus spinosa L., Erica vagans L., Calluna vulgaris L. Hull, leguminous species like Ulex europaeus L., Genista florida L., Cytisus sp., or Quercus ilex (ETSIM, 1978). 2.3 Local history of forest commons management As mentioned, the area has a long tradition of forest commons management. The invasion of Germanic tribes into the Iberian Peninsula in the fifth century brought the concept of woodlands 6 collective property to the northwest areas of the peninsula, where forest natural resources were used by local communities (Aranda, 1996). Some centuries later, during the Christian Reconquest (eighth to fifteenth century), the communal regime spread to other parts of the peninsula as a strategy followed by medieval kings to promote the settlement of Christian populations that would displace Muslims from the newly gained territories. Within this process, the crown, ultimate owner of the common lands, granted land concessions –including woodlands– to the new settlers, who organized in village councils or concejos to collectively manage and use natural resources (Behar, 1983; Pardo and Gil, 2005). Forest commons, or montes comunales in Spanish, were the most common type of tenure regime in Spanish woodlands until the nineteenth century, when Europe’s transition to capitalism gave rise to the establishment of a new liberal framework that initiated a process of privatisation of communal resources. Despite the long tradition of forest commons in Spain, the liberal framework did not recognize communal ownership and just distinguished between the public and private proprietorship of lands and goods. Under that political context, forest commons were classified as public properties and their management transferred to Spanish State Forest Administration (Beltrán, 2015; Guadilla-Sáez et al., 2017). In the Liébana region, forest commons belong to municipalities, while their use corresponds exclusively to local communities, with the regional forestry administration exerting its influence by monitoring access to grazing, firewood, and other forest goods (Balboa, 1999; Pérez-Soba and Solá, 2004). Remarkably, some forest-related stewardship customs still persist today in the study area, such as the neighbourhood councils or juntas vecinales, minor local bodies which have replaced former village councils and have legal rights to regulate the use of forest commons (Law No. 6/1994 of Cantabrian Regional Government). As a result, the regional forestry administration management intervention needs the approval from the neighbourhood councils before being applied in forest commons. However, and despite the regional legal enforcement of these minor local bodies role in the management of their common lands, the Picos de Europa Board of Trustees or Patronato –i.e. the participatory body aimed to integrate society to management activities and to promote further implications of local residents (Law No. 30/2014 of Spanish Government)–, does not include neighbourhood councils representatives. 2.4 Data collection We used a GIS procedure to randomly select plots in forest patches with, at least, 70% of tree canopy cover according to the Third Spanish National Forest Inventory (Ministerio de Medio Ambiente, 1997–2007). In each forest type, we located the centre of the plots according to a systematic sampling design in the intersection of a 125 × 125 m fishnet grid created with ArcGIS version 10.3.1 (ESRI, 2015). From the total possible labels, we selected 50 values using the ‘randbetween’ option of MS Excel and stored their spatial coordinates in a global positioning system (GPS) device. These values were taken as the centre of each plot. We inventoried 50 circular plots of 25 m radius (0.2 ha). Half of the plots were in forest commons inside Picos de 7 Europa (Fig.1, white-shaded area) and the other half were in forest commons located outside Picos de Europa and not affected by any other formal category of protection (Fig.1, grey-shaded area). Fig.1. Map of the study area, illustrating the location of the 50 plots across the Liébana region (Cantabria, Spain). We collected data during 2015 and 2016. We recorded a set of ecological parameters in each of the plots that included their general characteristics, the dominant tree species of the stand, stand structure properties, and distribution of the ground vegetation cover (see appendix B). We inventoried the dendrometric characteristics of the 10 adult trees closest to the centre of the plot, where an adult tree were defined as a tree with more than 3 meters height or with a minimum diameter at breast height (1.30 m) of 7.5 cm. To determine the abundance of each tree species, we identified all tree stems rooted within a sub-plot of 15 meters radius (aprox. 700 m2 in area). We also quantified the topographic variable slope for every plot using a SUUNTO clinometer and collected a surface soil sample to later analyse pH and texture parameters in the laboratory by using a glass electrode in a suspension of 1:10 soil:distilled water. As additional monitored field measurements, we inventoried other edaphic variables such as soil texture, organic matter thickness, and stoniness. We also recorded the presence of silvicultural treatments such as clear cuttings, brush removal, thinning, or ground improvements within 25 m radius (0.2 ha). We measured several anthropogenic variables associated with plot disturbance, including plot accessibility and linear distance from the plot centre to the nearest village and nearest path (i.e. unpaved roads or trails), calculated in desktop using ArcGIS version 10.3.1 (ESRI, 2015) in high 8 resolution aerial photographs derived from the Spanish Aerial Ortophoto National Plan (PNOA, 2015). We also measured population density of the nearest towns using data from the Spanish National Statistics Institute (INE, 2018). 2.5 Ecological indices Species presence indicators are frequently used to monitor effectiveness of a particular forest management treatment in biodiversity conservation (Canadian Council of Forest Ministers, 1997). For this research, we calculated a set of tree species composition indices to examine the heterogeneity in the composition of the studied communities, including the Shannon index (𝐻′ ), a separate measure of evenness (𝐽′ ) from the standardization of the Shannon index, Species Richness (𝐷𝑀𝑛 ), the complement of the Simpson index (1 − 𝐷), and the reciprocal form of the Berger-Parker index (1/𝑑) (Table 1). Although most studies have focused on the numerical richness to compare tree species diversity between different ecosystems (Hui et al., 2011), surrogate measures, such as stand structure indicators, are increasingly being used to provide a measure of biodiversity in forest communities (Pommerening, 2002). Following this trend, in this study, we have considered the Clark and Evans Index of aggregation (𝑅) to define the distance between neighbouring trees in a forest spatial structure unit (Neumann and Starlinger, 2001). We have also included the Uniform angle (contagion) index (𝑊𝑖 ), which tests the regularity of the distribution pattern of the trees, and the species complement for the Mingling index (1 − 𝑀𝑖 ) to test the heterogeneity of species among nearest neighbouring trees (Aguirre et al., 2003; Pommerening, 2002). The spatial unit considered for the estimation of the stand structure indices was the group size of the four nearest neighbouring trees to the reference tree, which is considered the optimum group size for evaluating spatial attributes (Hui et al., 2011). Table 1. List of the diversity and stand structure indices used to analyse the ecological data. Index Formula Definition Species composition indices where pi is the proportion of individuals found in Shannon (Magurran, 𝐻′ = − ∑ 𝑝𝑖 ln 𝑝𝑖 the i species referred to the total number of 2004) individuals. where 𝐻𝑚𝑎𝑥 is the maximum level of diversity Evenness (Elliott et al., 𝐻′ 𝐻′ 𝐽′ = = possible within a given population and 𝑆 is the 1997) 𝐻𝑚𝑎𝑥 ln 𝑆 total number of species. where 𝑆 equals the number of different species Richness (Magurran, 𝑆 𝐷𝑀𝑛 = represented in the sample, and 𝑁 is the total 2004) √𝑁 number of individuals in the sample. 9 Simpson (Magurran, where 𝑝𝑖 is the proportion of individuals in the ith 𝐷 = ∑ 𝑝𝑖2 2004) species. where 𝑁𝑚𝑎𝑥 is the number of individuals of the Berger-Parker 𝑁𝑚𝑎𝑥 𝑑= most abundant species and 𝑁 refers to the total (Magurran, 2004) 𝑁 number of individuals. Stand structure indices where 𝑟𝑖 is the distance from the reference tree to 1 𝑛 Clark-Evans (Vorčák et ∑𝑖 𝑟𝑖 its nearest neighbour, 𝑛 is the number of trees on 𝑅= 𝑛 𝑃𝑙 the sample plot and 𝑃𝑙 the area of the sample plot al., 2006) 0.5 × √ 𝑛 in square meters. where 𝑤𝑖𝑗 is 1 if the angle with the jth Uniform angle 𝑛 1 neighbouring tree is lower than the defined (contagion) index 𝑊𝑖 = ∑ 𝑤𝑖𝑗 𝑛 standard angle, and 𝑤𝑖𝑗 equals 0 otherwise. 𝑛 is 𝑗=1 (Pommerening, 2002) the number of trees on the sample plot. where 𝑣𝑖𝑗 is 1 if the jth neighbouring tree is not 𝑛 Species Mingling (Hui et 1 the same species as the ith reference tree, and 𝑀𝑖 = ∑ 𝑣𝑖𝑗 al., 2011) 𝑛 𝑣𝑖𝑗 equals 0 otherwise. 𝑛 is the number of trees 𝑗=1 on the sample plot. 2.6 Statistical analyses 2.6.1 Indices computed We estimated tree species presence and stand structure indicators for each forest habitat type and subjected the resulting data to the Shapiro-Wilk normality test (Table A.1). We used a two- sample t-test for testing for differences between forest habitat types inside and outside Picos de Europa for the normally distributed data and a two sample Wilconxon rank-sum (Mann-Whitney) test for the non-normally distributed data (Table A.2 to Table A.4). Significance levels (p-values) were adjusted to the number of tests carried out using a standard Bonferroni correction (Zuur et al., 2007). To facilitate the graphical comparison between stands of the three different forest habitats inside and outside the PA, we standardized correlation coefficients, representing in a common scale values obtained. We did all analysis using the STATA software version 13.1 (StataCorp, 2013). 2.6.2 Model selection We considered two explanatory environmental variables in our model: (1) plot topography and (2) soil characteristics. Due to their relevance for the research, we also included as explanatory variables the level of human intervention observed in the plots, including variables that measure 10 (3) plot isolation –i.e. distance to human infrastructures–, and (4) presence of anthropogenic disturbances –i.e. forest management practices–. As dependent variable, we used the abundance of individuals per tree species estimated by the Shannon index (𝐻′ ), a simple formula widely used to measure species richness in which the higher the value of 𝐻′ , the greater the species diversity of the studied system (Elliot et al., 1997; Magurran, 2004; Zuur et al., 2007). We evaluated the relative importance of the variables that measure environmental characteristics, plot isolation, and anthropogenic disturbances in explaining tree species distribution by using a Wald test analysis of 𝐻 ′ against different models combining the explanatory variables. To select the explanatory variables to be included in the final model, we checked the normality of the variables through a correlation matrix and applying the Shapiro-Wilk normality test (see Table A.5). Based on those results, we selected non- parametric statistics for several of the explanatory variables. To determine relations between the explanatory variables, we applied Pearson’s correlation test for the normally distributed variables and Spearman rank correlation test for the variables not-normally distributed (Table A.6). For categorical variables, we used a Pearson’s chi-squared correlation. After running the regressions, we removed multiple variables due to collinearity, identified by tolerances values approaching 0.1 in these predictors. The first step for the model selection consisted in fitting the global models of each set of explanatory variables to the data, examining the goodness-of-it of each model with a χ2 statistical test (Johnson and Omland, 2004). Due to our reduced sample (50 plots), models only included up to four variables at a time, as recommended by Harrell et al. (1996), representing all possible combinations of the variables (excluding interactions) with no model including more than one variable from each general category of variables (i.e. plot topography, soil characteristics, plot isolation, and anthropogenic disturbances). In total, we analysed a total of 479 models. We fit each model to the data by using an Ordinary least squares method and performed model comparison of all possible models by using the Akaike information criterion (AIC) for selecting the best set of explanatory variables in describing the variation of species diversity based on the minimum AIC score (Burnham and Anderson, 2002; Quinn and Keough, 2002; Johnson and Omland, 2004). We calculated the AIC difference (∆AICi) and Akaike weights (𝑤𝑖 ) for the ten best ranked models fitting the data to assess the statistical level of support for a given model (Table A.7). Due to the low AIC differences between models, we considered a subset of models with ∆AICi < 4 to estimate the relative importance of individual explanatory variables (Burnham and Anderson, 2002). We then repeated the analysis using Richness index as dependent variable. 3. Results 3.1 Protection status and ecological and human disturbance variables A one-way ANOVA analysis of variance test showed significant differences in the explanatory and dependent variables measured in forest commons inside and outside Picos de Europa (Table 2). 11 Mean values of longitude and latitude showed significant difference in plots inside (mean ± SD: 368965.8 ±5080.0) and outside (357567.4 ±4456.5) Picos de Europa, arguably because of the geographical location of the PA in the northwest part of the Liébana region. Our results suggest that there is a significantly thicker layer of organic matter in the soils of the plots outside (7.1 ±2.5) than in the soils of plots inside the PA (4.8 ±2.8). Additionally, there is a lower presence of surface stones in the plots outside the PA. Regarding the anthropogenic variables, plots inside Picos de Europa were more distant to towns (2225.6 ±907.6) than plots outside the PA (817.2 ±434.1), likely due to the reduced number of villages that have their municipality boundaries inside the protected area. Nearest town population density presented an average of 66.6 and 37.6 inhabitants, respectively, a situation explained by the presence of two very touristic villages inside Picos de Europa, Espinama and Pido. We also found significant differences in the implementation of silvicultural practices, with presence of some type of silvicultural systems –regeneration felling, forest cover improvement and/or ground improvements– in 80% of studied plots inside the PA and only in 52% for plots outside it. Specifically, plots inside Picos de Europa had a higher presence of cover improvement treatments (80%) than plots outside (48%). We did not detect significant differences in the presence of regeneration felling or ground improvement practices. The census of the 50 sample plots led to the identification of a total of 14 families and 17 tree species, with an average number of 2.52 (SD = 1.7) and 3.8 (SD = 1.5) species per plot in plots inside and outside Picos de Europa, respectively. We observed significant differences between species abundance distributions in Shannon index means across plots. Specifically, the mean Shannon index was lower in plots inside Picos de Europa (0.45 ±0.5) than in plots outside it (0.86 ±0.4). We did not find significant differences for the Richness index, which also showed lower values in plots inside (0.32 ±0.2) than outside Picos de Europa (0.47 ±0. 3). Table 2. Mean and standard deviation (SD) (or percentage for categorical variables) of variables considered, by location in relation to Picos de Europa. The One-way ANOVA analysis of variance compares the means in plots outside and inside the protected area, in bold significant differences based on their p-values (** p-value ≤ 0.01; * p-value ≤0.05) . Pool Outside Inside Mean Mean Mean Variable Code Specification p-value Na N (±SD) (±SD) (±SD) Explanatory variable Topography 363266.6 368965.8** 357567.4** Longitude UTMY Degrees .000 25 25 (±7450.6) (±5080.0) (±4456.5) 4779509 4778188* 4780830* Latitude UTMX Degrees .042 25 25 (±4620.2) (±4004.1) (±4890.2) 46.9 Slope SLO Percentage .276 24 48.9 (±12.5) 25 44.9 (±12.9) (± 12.7) Soil characteristics 12 pH PH Numeric scale 6.2 (±0.8) .177 24 6.0 (±0.7) 21 6.4 (±0.95) 1- Sandy 45.8% 43.5% 48.0% Texture TEX 2- Loam 22.9% .388 23 39.1% 25 8.0% 3- Clay 31.2% 17.4% 44.0% Organic matter OM Centimetres 6.0 (±2.9) .004 25 7.1** (±2.5) 25 4.8** (±2.8) 1- Without stones 26.0% 36.0%* 16.0%* 2- Low stony 32.0% 36.0%* 28.0%* Stoniness STO .043 25 25 3- Stony 18.0% 12.0%* 24.0%* 4- Very stony 24.0% 16.0%* 32.0%* Isolation 161.78 149.6 174 Distance to path DIST1 Meters .651 25 25 (±188.4) (±154.1) (±220.1) 1521.4 817.2** 2225.6** Distance to town DIST2 Meters .000 25 25 (±1000.9) (±434.1) (±907.6) Number of 52.7 37.6** 66.6** Town population POP .001 23 25 inhabitants (±31.7) (±30.6) (±26.3) Anthropogenic disturbances 0- No presence 32.0% 44.0% 20.0% Grazing GRA .071 25 25 1- Presence 68.0% 56.0% 80.0% Silvicultural 0- No presence 34.0% 48.0%* 20.0%* SILV .037 25 25 treatments 1- Presence 66.0% 52.0%* 80.0%* Regeneration 0- No presence 48.0% 60.0% 36.0% FEL .093 25 25 felling 1- Presence 52.0% 40.0% 64.0% 0- No presence 36.0% 52.0%* 20.0%* Cover improvement COV .018 25 25 1- Presence 64.0% 48.0%* 80.0%* Ground 0- No presence 96.0% 96.0% 96.0% GRO 1.00 25 25 improvement 1- Presence 4.0% 4.0% 4.0% Dependent variable Shannon 𝐻′ 0.64 (±0.5) .004 21 0.86** (±0.4) 25 0.45** (±0.5) Richness 𝐷𝑀𝑛 0.39 (±0.3) .064 21 0.47 (±0.3) 25 0.32 (±0.2) a Not available or unclear observations excluded from the analysis. 3.2 Ecological indices Although we observed higher species diversity in the mixed broadleaf woodlands plots inside Picos de Europa than in the plots outside, we did not find any significant differences in the ecological indices calculated (Fig.2 and Table A.2). In contrast, although protected and unprotected Fagus sylvatica plots were relatively similar in terms of spatial diversity, the richness and evenness was higher in plots outside (Fig.2 and Table A.3). Finally, we do not find any 13 Standardized value 0 -1 significant differences between Quercus pyrenaica plots inside and outside the protected area (Fig.2 and Table A.4). -2 ** * * A * A A A Shannon 4 Outside Inside 2 0 -2 -4 A B C A B C A B C A B C A B C A B C A B C A B C Shannon Evenness Richness 1- Simpson 1/Berger-Parker Clark-Evans 1- Mingling Contagion Fig. 2. Boxplots illustrating Z scores for the five species composition indicators (n=46) and the three stand structure ones (n=49) for mixed broadleaf woodlands (A), Fagus sylvatica forests (B), and Quercus pyrenaica stands (C). Asterisks indicate significant differences in means between protected and unprotected Fagus sylvatica forests based on their p-values after Bonferroni correction (** p-value ≤ 0.005; * p-value ≤ 0.025). 3.3 Model selection According to the AIC, the best-ranked model for explaining the evenness of tree species of forest commons in the Liébana valley is given by Eq. 1: 𝐻′ = (−0.002) × SLO + (−0.003) × PH + (−0.0002) × DIST2 + (−0.196) × 𝐹𝐸𝐿 + 1.287 (1) where 𝐻′ is the Shannon index value of the studied plots. This model had an AIC c on 47.75 which gave 𝑤𝑖 = 5.71% (Table A.7). The best-ranked model for explaining the species richness of forest commons in the Liébana valley is given by Eq. 2: 𝐷𝑀𝑛 = 0.028 × 𝑂𝑀 + (−0.0004) × 𝐷𝐼𝑆𝑇1 + 0.289 (2) where 𝐷𝑀𝑛 is the Richness index value for the studied plots. This model had an AIC c on 5.22 which gave 𝑤𝑖 = 2.33% (Table A.7). In both models Akaike weights have a very low value, indicating uncertainty that these models are the best approximating models to our data (Symonds 14 and Moussalli, 2011). Given this uncertainty, and following McCracken et al. (2015), we tested the dependent variables against the explanatory variables included in the subset of models with ∆AICi < 4, estimating the Akaike weights for these candidate explanatory variables and their relative importance across the subset models (Table 3). Table 3. Importance of environmental and anthropogenic variables considered across the 50 study plots. Importance was derived using Akaike weights (𝑤𝑖 ), reporting in bold the most important variables (with importance ≥0.4). See table 2 for variables definition. Topography Soil characteristics Isolation Disturbances Variable UTMY UTMX SLO PH TEX OM STO DIST1 DIST2 POP GRA SILV FEL COV GRO Response = Shannon index. Models where ∆𝐴𝐼𝐶 < 4 = 44 of 473. Importance .05 .06 .63 .44 .29 .04 .08 1 .08 .10 .32 .09 .08 Response = Richness index. Models where ∆𝐴𝐼𝐶 < 4 = 93 of 473. Importance .18 .17 .20 .78 .13 .40 .49 .02 .07 .11 .16 .12 .17 The analysis shows that the two abiotic variables slope of the plot (SLO) and soil acidity (PH), as well as the anthropogenic variable distance to nearest town (DIST2) are the explanatory variables most associated with tree species evenness. All these variables bear a negative association with the Shannon index thus suggesting that, in general, species evenness is higher in plots with lower slopes, less acid soils, and closest to human settlements. Presence of anthropogenic disturbances as regeneration felling practices (FEL) also show a negative association to species evenness, although the strength of this association is weak. In addition, the soil characteristic thickness of organic matter (OM) and the anthropogenic variables distance to nearest path (DIST1) and town (DIST2) seem to play an important role in the distribution of vegetation communities. Tree species richness was positively associated to organic matter thickness, indicating that many species prefer soils with a deep organic humus layer. In contrast, the Richness index was negatively associated to plot isolation variables, showing more species diversity in more accessible plots, or plots located closer to human settlements. 4. Discussion Our results generally suggest that tree species composition is less heterogeneous inside than outside the Picos de Europa National Park. Two main factors may explain this result: the dominance of monospecific Fagus sylvatica forest inside the protected area and the application of silvicultural systems oriented to promote the presence of beech forests inside the protected area. Our results also suggest that human intervention variables, particularly distance to nearest town, are more important drivers of tree species distribution and diversity in forest commons of the Liébana valley than the abiotic factors considered for the analysis. 15 4.1 Ecological variables The main goals of this work were to quantify differences in ecological characteristics between forest commons inside and outside Picos de Europa and to analyse how location in relation to the protected area relates to tree species distribution and diversity. We found that, compared to soils in the surrounding landscape, soils inside Picos de Europa are characterized by a higher abundance of stones and a thinner humus layer, two important soil parameters associated to plant diversity in forests (Cantero et al., 2003; French et al., 2008; Ren et al., 2012). Stoniness is negatively associated to species richness in the study area, a result in line with previous studies in Mediterranean environments (Ceacero et al., 2012). However, according to our best-ranked models, stoniness is not a relevant variable for tree species diversity. In contrast, topsoil organic matter content shows an important effect, particularly on the Richness index, showing a positive association with species composition. This result, which matches previous studies analysing the effect of organic matter mass in temperate deciduous forests in central Europe (Härdtle et al., 2003), may explain the lower value of the Richness index of plots inside Picos de Europa. However, topsoil organic matter is not associated in a significant way to the Shannon index, the species diversity indicator that actually differed between protected and unprotected plots. Our best-ranked models also show an association between species richness and pH soil acidity and plot slope. These results match findings observed in earlier studies analysing the effects of edaphic and topographic factors on species richness (Härdtle et al., 2003; French et al., 2008; Mod et al., 2016). Still, these two ecological factors are similar in all study plots, for which they do not help to explain differences on species composition between plots in the protected area and outside it. Finally, we observed lower Shannon index values in plots inside the PA, indicating a more heterogeneous landscape outside Picos de Europa. A possible explanation for this result is the dominance of Fagus sylvatica stands inside Picos de Europa, which contrast with the dominance of mixed stands in plots outside Picos de Europa. As beech forests are characterized for conforming monospecific stands (Krämer and Hölscher, 2009), this homogeneity of tree species might result in low values of the diversity Shannon index. Similar differences have been observed in other studies comparing ecological features of near-natural Fagus sylvatica stands and communities with less proportion of one dominant tree species (Hui et al., 2011). Overall, soil’s organic matter seems to be the only ecological factor studied contributing to explain the lower species inside the PA. Yet, soil organic matter is not associated to the Shannon index, the only indicator that significantly differed between plots inside and outside the PA. Thus, overall, our results do not show key association between ecological variables and species richness in Liébana forest commons. 16 4.2 Anthropogenic disturbances The second objective of this study was to quantify differences in human-dominated disturbances between forest commons under the protected status and forest commons outside it, and to analyse how these factors may influence tree species richness distribution in the study area. Best-ranked models for Shannon and Richness indices show a significant association between tree species composition and anthropogenic disturbances, particularly plot isolation and presence of silvicultural systems. The most important variables in our models are those measuring distance to the nearest path and village, and, with less direct importance, presence of regeneration felling treatments. These findings dovetail with recent research evaluating variables associated to species diversity, which also emphasize the importance of including anthropogenic disturbances in species diversity analyses, particularly in human-dominated landscapes (Guèze et al., 2015; Mod et al., 2016). Five aspects deserve further discussion. First, the results of the correlation analysis underscore the significant association between species richness and distance to the nearest human settlement, a finding previously reported in tropical forest of Bolivian Amazon (Guèze et al., 2015) and in pine and oak forests across Mexico (Silva-Flores et al., 2016). The finding, however, contrasts with one study in temperate oaks in Spanish Central Pyrenees showing no association between distance to nearest town and Quercus spp spatial distribution (Kouba et al., 2011). We argue that these differences in findings probably relate to variations in the socioeconomic characteristics in the study sites, which might have resulted in different impacts on the ecological system (Meyer and Crumley, 2011). Thus, during the second half of the twentieth century, the Pyrenean region suffered an important depopulation and a consequent abandonment of human activities, a phenomena that frequently results in the homogenization of the rural landscape mosaic as woody species colonized abandoned lands (Rotherham, 2013; Viedma et al., 2015; Lavorel et al., 2017). Although the Liébana region also suffered from depopulation during that time –i.e. population decreased by 54% from 1950 to 1981 (Reques, 1997)–, the increasing demand for food supply by the neighbouring industrialized areas favoured the specialization in livestock production of Liébana from the 1970s onwards (González, 2001). Traditional livestock farming operations such as hay making, pruning, or grazing on forest commons allowed the maintenance of a grassland- woodland mosaic, and the biodiversity dependent upon these practices, until the early 2000s. Timber harvest of native deciduous species like Fagus sylvatica continued until 1980s, although the twenty century saw a shifting to harvest non-native conifers from the commercial timber- producing plantations established by the State Forestry Administration intervention of public woodlands in Liébana. This shift favoured the co-occurrence of wild and synantrophic species (Ezquerra and Gil, 2004). Second, our modelling methods show a negative association between distance to nearest town and species diversity, i.e. forest species diversity is higher in areas close to human settlements. This association might be due to the human pressure potentially associated to higher accessibility, as Guirado et al. (2007) also noticed for periurban oak forests of north-eastern Spain, where 17 presence of recent human disturbances is associated to higher species richness, particularly synanthropic species. Other studies have also linked intensity of human disturbances to the presence of pioneer, non-native species in natural habitats (Battles et al., 2001; Gourlet-Fleury et al., 2013; Bauman et al., 2015). In our study, we observed Castanea sativa and Ilex aquifolium, species traditionally pruned for the provision of food (chestnuts) and winter fodder, in plots located less than 1000 meters away from towns. We also recorded plantations of the exotic Pinus radiata D. Don inside or around accessible plots located less than 650 meters away from human settlements. Inside these plots, we also recorded species untypical for the studied forest habitats, like Arbutus unedo L. and Pyrus sp. These results may provide further support to the idea that plot isolation does not necessarily result in an enriched forest habitat, as human intervention may sometimes increase the total number of species by introducing atypical species in a community (Helm et al., 2015). Third, this research also shows a relevant association between species diversity and distance to nearest path, including unpaved roads or walking trails trampling by foot, animals, or wheeled vehicles, but in which vehicle traffic is limited to forest rangers and local inhabitants use. Species diversity bears a negative association distance to nearest path; in other words, forest species diversity is higher in better accessible areas, a finding in line with a recent study analysing the effects of human path, trails, and roads on plant species richness (Root-Bernstein and Svenning, 2018). Nevertheless, distance to path does not significantly differ between plots inside and outside Picos de Europa, for which this anthropogenic disturbance does not assist in explaining differences on species richness between protected and unprotected sites. Fourth, our findings show that plots inside Picos de Europa tend to present more silvicultural treatments, particularly forest cover improvement treatments, than plots outside the park. Greater forestry operations inside the PA may be counterintuitive for practitioners considering the protection category of National Park as forest reserve where minimal intervention is applied to allow a continued succession and natural disturbances in the forests. In Spain, however, National Parks are actively managed by authorities who implement management strategies oriented to safeguard the natural systems that justified the PA designation (Law No. 30/2014 of Spanish Government). In Picos de Europa, these strategies center around silvicultural activities oriented to preserve the habitats that justified its designation, i.e. Atlantic Forest natural and semi-natural habitats such as Fagus sylvatica forest patches (Regional Forest Administration, pers. comm., September 2016). Although Picos de Europa National Park management plan was derogated in 2005 due to a legal action taken by local communities living in the buffer zone (Spanish Government, 2005), silvicultural operations are undertaken under the umbrella of national and regional legislation and European conservation initiatives like the LIFE+ Cantabrian Capercaillie Project (LIFE09 NAT/ES/000513). Thus, weeding and brushing out of trails have been carried out to facilitate the access of visitors to the National Park (Park ranger, pers. comm., September 2016). In addition, the conservation strategy followed by the National Park Administration to promote the presence of the endangered bird subspecies Cantabrian capercaillie (Tetrao urogallus cantabricus), considered an umbrella species in montane forest ecosystems (Blanco- 18 Fontao et al., 2011), includes silvicultural interventions on forest cover practices in an attempt to favour capercaillie’s habitat. These silvicultural practices included brush cutting and thinning applied to reduce canopy closure and to facilitate Vaccinium myrtillius L. growth, an important food source to capercaillie’s populations (Lakka and Kouki, 2009; Mikoláš et al., 2015). Accordingly, we found that plots inside Picos de Europa tend to present more silvicultural treatments, particularly forest cover improvement treatments, than plots outside the park. Overall, silviculture inside the reserve seems to affect tree composition by reducing woody species diversity. Finally, when studying the influence of silvicultural systems on species’ diversity, we find that the only variable with relative importance in our models is the presence of regeneration felling operations, which shows a negative association to the Shannon index. Our results mainly relate to the effects of selective felling, a low-intensity clear-cutting activity practice consistent in an individual-tree selection cutting 'that maintains or develops an uneven-aged forest structure over time' (Lexerød and Eid, 2006, p.503). Several reports have discussed the association of low intensity, close-to-nature silviculture to the presence of vegetation typical of late-successional stages (Battles et al., 2001; Saeki, 2007). On the one side, by benefitting the presence of late- successional ecosystems, selective felling contributes to a higher evenness of the first dominant tree species in the study area, Fagus sylvatica (Ministerio de Medio Ambiente, 1997–2007). On the other, when focusing on mixed broadleaf forests, our study brings deeper insights into the assumption that non-intensive silvicultural systems promotes uneven-aged stands that benefit species diversity by increasing vegetation composition associated to late-successional ecosystems of forest sites. Particularly, we observed Corylus avellana, as dominant tree species, and Crataegus monogyna, as accompanying tree species, in mixed broadleaf forests plots presenting selective felling practices. The latter finding is in line with published studies analysing the effects of selective felling in coniferous and deciduous species, which also consider selective felling as a management practice that favours biological diversity for those forest habitats (Atlegrim and Sjöberg, 2004; Martín-Alcón et al., 2015). In sum, silviculture inside the reserve seems to affect tree composition by reducing woody species diversity. Among the silvicultural operations considered, presence of regeneration felling practices is the only treatment with relative importance for species’ diversity, favouring the evenness of species associated to late-successional ecosystems. Remarkably, the anthropogenic variable distance to town is a key driver of species diversity in the study area, with a higher number of species in plots closer to human settlements. 4.3 Protection status and species diversity interaction in three forest habitats When comparing species diversity of three different temperate deciduous forest habitats occurring in forests commons inside and outside Picos de Europa,We observe similar tree species pattern for mixed broadleaf forests and Pyrenean oak forests regardless the protection 19 status. Our findings are in line with Gray et al. (2016) when comparing biodiversity between protected and unprotected areas with the same land use. We argue that the finding is explained by the fact that, despite differences in the protection status, both sites are managed by the Forest Administration of the regional government, who, in the absence of a valid management plan in the PA, applies similar forest management technical criteria inside and outside the protected area (Forest ranger, pers. comm., August 2016). Interestingly, we found differences in the ecological features of one forest habitat type, Fagus sylvatica stands, which have a significant more heterogeneous species composition outside Picos de Europa. Specifically, we observed significant differences for the Shannon, Evenness and Richness indices and for the reciprocal form of the Berger-Parker index (Fig.3), similar to the ones observed by Bilek et al. (2011) when comparing managed and unmanaged beech forests in Central Europe. The low value of Shannon index in the beech forests inside Picos de Europa indicate a single-layered Fagus sylvatica composition of these stands, in line with observations in pure beech forests of the Basque Country, a neighbouring region (Peña et al., 2011). Nevertheless, although a high portion of accompanying species may not be considered a natural pattern of beech stands, studies analysing the temporal and spatial dynamics of near-natural deciduous forests dominated by Fagus sylvatica show a forest cycle based on ‘gap-dynamics’, in which species composition oscillate due to small scale disturbances such as canopy openings (Wissel, 1992; Emborg et al., 2000). Briefly, by promoting the presence of monospecific beech stands without including species associated to the natural small-scale disturbances of these habitats, interventions carried out in beech forests inside Picos de Europa do not seem to contribute to higher stand diversification, plots inside the PA being poorer in terms of species than plots outside the PA. 5. Conclusion Results presented herein show more heterogeneity of forest communities outside than inside Picos de Europa National Park, arguably as silvicultural systems applied inside the designated area benefit the presence of single-layered Fagus sylvatica forests. Aiming to provide some management guidelines to favour biodiversity maintenance in this protected area, we recommend the inclusion of practices emulating natural disturbances of beech forests, promoting the appearance of pioneer light-demanding specialist species associated to these forest habitats. This recommendation implies a strategy to enhance the resilience capacity associated to Fagus sylvatica ecosystems, particularly relevant in the context of climate change in the study area, which might increase vulnerability to pests and phenological changes of forest habitats (OECC, 2012). Further research under scenarios of climate change is required to attain more integrative management recommendations. Within the range of management options for promoting tree species assemblages in this forest habitat type, our work brings into consideration the positive influence that human-induced 20 disturbances hold for increasing tree species richness and evenness. This option might be especially relevant in natural landscapes with long histories of forest use, such as the study site. Our study also reinforces the social-ecological systems approach to biodiversity conservation by providing an example of how the persistence of some stewardship customs and traditional uses of local communities in forest commons of the Liébana valley seem to result in the long-term maintenance of diverse natural habitats and their associated species. We recommend further research to identify which particular forest-related practices performed in forest commons outside Picos de Europa enrich community composition through habitat-specific species. Acknowledgments This research was supported by the Agència de Gestió d’Ajuts Universitaris i de Recerca AGAUR of the Government of Catalonia (2015FI_B00333) and by the Spanish Ministry of Economy and Competitiveness through the research project ‘Citizen Science and traditional agroecological knowledge: How to increase citizen's participation in the Spanish inventory of traditional knowledge related to biodiversity?’ (CSO2014-59704-P). This work contributes to ICTA 'Unit of Excellence' (MinECo, MDM2015-0552). We thank Peter Wilshusen and Maximiliem Guèze for their valuable comments on an earlier version of the text. We would also like to acknowledge and thank Miguel de Terán and Javier del Hoyo for their support during field work, and Picos de Europa National Park and Biosphere Reserve Administration and the Cantabrian Regional Government for collaboration and permission to conduct this study. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at Appendix B. Field plot form References Agrawal, A., Gibson, C.C., 1999. Enchantment and Disenchantment: The role of Community in Natural Resource Conservation. World Dev 27(4), 629–649. Aguirre, O., Hui, G., Gadow, K., Jiménez, J., 2003. An analysis of spatial forest structure using neighbourhood-based variables. For. Ecol. Manage. 183, 137–145. https://doi.org/10.1016/S0378-1127(03)00102-6 Andrade, G.S.M., Rhodes, J.R., 2012. Protected Areas and Local Communities: an Inevitable Partnership toward Successful Conservation Strategies? Ecol Soc 17(4), 14. http://dx.doi.org/10.5751/ES-05216-170414 21 Angelsen, A., Jagger, P., Babigumira, R., Belcher, B., Hogarth, N.J., Bauch, S., Börner, J., Smith- Hall, C., Wunder, S., 2014. Environmental Income and Rural Livelihoods: A Global-Comparative Analysis. World Dev 64, S12–S28. https://doi.org/10.1016/j.worlddev.2014.03.006 Aranda, G., 1996. La selvicultura en España hasta el siglo XIX. Ecología 10, 173-184. Arbeo, P. 2012. La sociedad económica de amigos del país de Liébana en el siglo XIX. Madrid: Bubok Publishing, S.L. Aryeetey, E., Devarajan, S., Kanbur, R., Kasekende, L. editors. 2012. The Oxford Companion to the Economics of Africa. New York: Oxford University Press Inc. Atlegrim, O., Sjöberg, K., 2004. Selective felling as a potential tool for maintaining biodiversity in managed forests. Biodivers. Conserv. 13, 1123–1133. https://doi.org/10.1023/B:BIOC.0000018148.84640.fd Balboa, X., 1999. La historia de los montes públicos españoles (1812-1936): Un balance y algunas propuestas. Historia Agraria 18, 95-128. Battles, J.J., Shlisky, A.J., Barrett, R.H., Heald, R.C., Allen-Diaz, B.H., 2001. The effects of forest management on plant species diversity in a Sierran conifer forest. For. Ecol. Manage. 146, 211– 222. https://doi.org/10.1016/S0378-1127(00)00463-1 Bauman, J.M., Cochran, C., Chapman, J., Gilland, K., 2015. Plant community development following restoration treatments on a legacy reclaimed mine site. Ecol Eng 83, 521–528. https://doi.org/10.1016/j.ecoleng.2015.06.023 Behar, R. 1983. The presence of the past: A historical ethnography of a Leonese village. PhD dissertation, Princeton University. Beltrán, F.J., 2015. Social and Environmental Filters to Market Incentives: The Persistence of Common Land in Nineteenth-Century Spain. J AGRAR CHANGE 15(2), 239–260. https://doi.org/10.1111/joac.12056 Berkes, F., Colding, J., Folke, C. editors. 2003. Navigating Social-Ecological Systems. Building resilience for complexity and change. United Kingdom: Cambridge University Press. Berkes, F., 2017. Environmental Governance for the Anthropocene? Social-Ecological Systems, Resilience, and Collaborative Learning. Sustainability 9(7), 2071-2105. https://doi.org/10.3390/su9071232 Bilek, L., Remes, J., Zahradnik, D., 2011. Managed vs. unmanaged. Structure of beech forest stands (Fagus sylvatica L.) after 50 years of development, Central Bohemia. For Syst 20(1), 122– 138. http://dx.doi.org/10.5424/fs/2011201-10243 22 Blanco-Fontao, B., Quevedo, M., Obeso, J.R., 2011. Abandonment of traditional uses in mountain areas: typological thinking versus hard data in the Cantabrian Mountains (NW Spain). Biodivers. Conserv. 20, 1133–1140. https://doi.org/10.1007/s10531-011-0016-1 Bocherens, H., 2018. The Rise of the Anthropocene since 50,000 Years: An Ecological Replacement of Megaherbivores by Humans in Terrestrial Ecosystems? Frontiers in Ecology and Evolution 6(3), 1–8. https://doi.org/10.3389/fevo.2018.00003 Brockington, D., 2002. Fortress Conservation. The preservation of the Mkomazi Game Reserve, Tanzania. Oxford: James Currey. Burnham, K.P., Anderson, D.R., 2002. Model Selection and Multimodel Inference: a practical information-theoretic approach. 2nd Edition. Springer-Verlag New York. Canadian Council of Forest Ministers, 1997. Criteria and indicators of sustainable forest management in Canada. Technical report. Ottawa, Canada. Cantero, J.J., Liira, J., Cisneros, J.M., Gonzalez, J., Nuñez, C., Petryna, L., Cholaky, C., Zobel, M., 2003. Species richness, alien species and plant traits in Central Argentine mountain grasslands. J. Veg. Sci. 14, 129–136. https://doi.org/10.1111/j.1654-1103.2003.tb02135.x Ceacero, C.J., Díaz-Hernández, J.L., del Campo, A.D., Navarro-Cerrillo, R.M., 2012. Interactions between soil gravel content and neighboring vegetation control management in oak seedling establishment success in Mediterranean environments. For. Ecol. Manage. 271, 10–18. https://doi.org/10.1016/j.foreco.2012.01.044 Chhatre, A., Agrawal, A., 2008. Forest commons and local enforcement. PNAS 105(36), 13286- 13291. https://doi.org/10.1073/pnas.0803399105 Cumming, G.S., Allen, C.R., Ban, N.C., Biggs, D., Biggs, H.C., Cumming, D.H.M., de Vos, A., Epstein, G., Etienne, M., Maciejewski, K., Mathevet, R., Moore, C., Nenadovic, M., Schoon, M., 2015. Understanding protected area resilience: a multi-scale, social-ecological approach. Ecol Appl 25(2), 299–319. https://doi.org/10.1890/13-2113.1 Elliot, K.J., Boring, L.R., Swank, W.T., Haines, B.R., 1997. Successional changes in plant species diversity and composition after clearcutting a Southern Appalachian watershed. For. Ecol. Manage. 92, 67–85. Emborg, J., Christensen, M., Heilmann-Clausen, J., 2000. The structural dynamics of Suserup Skov, a near-natural temperate deciduous forest in Denmark. For. Ecol. Manage. 126, 173–189. https://doi.org/10.1016/S0378-1127(99)00094-8 ESRI, 2015. ArcGIS Desktop: Release 10.3.1. Redlands, CA: Environmental Systems Research Institute. 23 ETSIM, 1978. Ordenación del paisaje. Estudios de planificación física. El Valle de Liébana. Madrid: Escuela Técnica Superior de Ingenieros de Montes. Ezquerra, F.J., Gil, L., 2004. La transformación histórica del paisaje forestal en la comunidad de Cantabria. Tercer Inventario Forestal Nacional. Madrid: Ministerio de Medio Ambiente. French, L.J., Smith, G.F., Kelly, D.L., Mitchell, F.J.G., O’Donoghue, S., Iremonger, S.F., McKee, A.M., 2008. Ground flora communities in temperate oceanic plantation forests and the influence of silvicultural, geographic and edaphic factors. For. Ecol. Manage. 255, 476–494. https://doi.org/10.1016/j.foreco.2007.09.014 Gannon, P., Seyoum-Edjigu, E., Cooper, D., Sandwith, T., Ferreira de Souza Dias, B., Pașca Palmer, C., Lang, B., Ervin, J., Gidda, S., 2017. Status and prospects for achieving Aichi biodiversity Target 11: Implications of national commitments and priority actions. Parks 23(2), 13– 26.Gaston, K.J., Jackson, S.F., Cantú-Salazar, L., Cruz-Piñón, G., 2008. The Ecological Performance of Protected Areas. Annu Rev Ecol Evol Syst 39, 93–113. https://doi.org/10.1146/annurev.ecolsys.39.110707.173529 Godefroid, S., Rucquoij, S., Koedam, N., 2005. To what extent do forest herbs recover after clearcutting in beech forest? For. Ecol. Manage. 210, 39–53. https://doi.org/10.1016/j.foreco.2005.02.020 González, M., 2001. Sociología y ruralidades. La construcción social del desarrollo rural en el Valle de Liébana. Madrid: Ministerio de Agricultura, Pesca y Alimentación. González, M., 2015. Denominaciones de Origen y Parques Nacionales: Actividad quesera en Picos de Europa. Revista de Humanidades 25, 57-84. Gourlet-Fleury, S., Beina, D. Fayolle, A., Ouédraogo, D.Y., Mortier, F., Bénédet, F., Closset-Kopp, D., Decocq, G., 2013. Silvicultural disturbance has little impact on tree species diversity in a Central African moist forest. For. Ecol. Manage. 304, 322–332. http://dx.doi.org/10.1016/j.foreco.2013.05.021 Gray, C.L., Hill, S.L.L., Newbold, T., Hudson, L.N., Börger, L., Contu, S., Hoskins, A.J., Ferrier, S., Purvis, A., Scharlemann, J.P.W., 2016. Local biodiversity is higher inside than outside terrestrial protected areas worldwide. Nat Commun 7, 12306. https://doi.org/10.1038/ncomms12306 Guadilla-Sáez, S., Pardo-de-Santayana, M., Reyes-García, V., 2017. The dismantling of forest commons in Spain. In: Proceedings of the XVIth IASC Conference Practicing the commons. Self- governance, cooperation and institutional change. Utrecht, the Netherlands, July 10-14. Available at: http://hdl.handle.net/10535/10373 Guirado, M., Pino, J., Rodà, F., 2007. Comparing the role of site disturbance and landscape properties on understory species richness in fragmented periurban Mediterranean forests. Landsc Ecol 22, 117–129. https://doi.org/10.1007/s10980-006-9009-y 24 Guèze, M., Luz, A.C., Paneque-Gálvez, J., Macía, M.J., Orta-Martínez, M., Pino, J., Reyes- García, V., 2015. Shifts in indigenous culture relate to forest tree diversity: A case study from the Tsimane’, Bolivian Amazon. Biol. Conserv 186, 251–259. http://dx.doi.org/10.1016/j.biocon.2015.03.026 Härdtle, W., von Oheimb, G., Westphal, C., 2003. The effects of light and soil conditions on the species richness of the ground vegetation of deciduous forests in northern Germany (Schleswig- Holstein). For. Ecol. Manage. 18, 327–338. https://doi.org/10.1016/S0378-1127(03)00091-4 Harrell, F.E., Lee, K.L., Mark, D.B., 1996. Multivariable prognostic models: Issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Stat Med 15, 361–387. Helm, A., Zobel, M., Moles, A.T., Szava-Kovats, R., Pärtel, M., 2015. Characteristic and derived diversity: implementing the species pool concept to quantify conservation condition of habitats. Divers. Distrib. 21, 711–721. https://doi.org/10.1111/ddi.12285 Hewitt, R., Pera, F., Escobar, F., 2016. Recent land cover changes in Spanish national parks and their surroundings. Cuadernos geográficos 55(2), 46–84. Hirschnitz-Garbers, M., Stoll-Kleemann, S., 2011. Opportunities and barriers in the implementation of protected area management: a qualitative meta-analysis of case studies from European protected areas. Geogr J 177(4), 321–334. https://doi.org/10.1111/j.1475- 4959.2010.00391.x Hui, G., Zhao, X., Zhao, Z., von Gadow, K., 2011. Evaluating Tree Species Spatial Diversity Based on Neighborhood Relationships. For. Sci. 57(4), 292–300. INE, 2018. Instituto Nacional de Estadística. http://www.ine.es/nomen2/index.do (accessed 6 February 2018). IUCN, UNEP-WCMC, 2017. The World Database on Protected Areas (WDPA). Cambridge (UK): UNEP World Conservation Monitoring Centre. http://www.protectedplanet.net/ (accessed 15 December 2017) Johnson, J.B., Omland, K.S., 2004. Model selection in ecology and evolution. Trends Ecol Evol 19(2), 101–108. https://doi.org/10.1016/j.tree.2003.10.013 Kelemen, K., Mihók, B., Gálhidy, L., Standovár, T., 2012. Dynamic Response of Herbaceous Vegetation to Gap Opening in a Central European Beech Stand. SILVA FENN 46(1), 53–65. https://doi.org/10.14214/sf.65 Kirby, K.J., Watkins, C. (Eds.), 2015. Europe’s Changing Woods and Forests: From Wildwood to Managed Landscapes. CABI, UK. 25 Kouba, Y., Alados, C.L., Bueno, C.G., 2011. Effects of abiotic and anthropogenic factors on the spatial distribution of Quercus faginea in the Spanish Central Pyrenees. Plant Ecol 212, 999– 1007. https://doi.org/10.1007/s11258-010-9880-0 Krämer, I., Hölscher, D., 2009. Rainfall partitioning along a tree diversity gradient in a deciduous old-growth forest in Central Germany. Ecohydrology 2, 102–114. https://doi.org/10.1002/eco.44 Kremen, C., 2015. Reframing the land-sparing/land-sharing debate for biodiversity conservation. Ann. N. Y. Acad. Sci. 1355, 52–76. https://doi.org/10.1111/nyas.12845 Lakka, J., Kouki, J., 2009. Patterns of field layer invertebrates in successional stages of managed boreal forest: Implications for the declining Capercaillie Tetrao urogallus L. population. For. Ecol. Manage. 257, 600–607. https://doi.org/10.1016/j.foreco.2008.09.042 Larson, L.R., Conway, A.L., Krafte, K.E., Hernandez, S.M., Carroll, J.P., 2016. Community-based conservation as a potential source of conflict around a protected area in Sierra Leone. Environ Conserv 43(3), 242–252. https://doi.org/10.1017/S0376892916000096 Laudati, A.A., 2010. The Encroaching Forest: Struggles Over Land and Resources on the Boundary of Bwindi Impenetrable National Park, Uganda. Soc Nat Resour 23, 776–789. https://doi.org/10.1080/08941920903278111 Lavorel, S., Grigulis, K., Leitinger, G., Kohler, M., Schirpke, U., Tappeiner, U., 2017. Historical trajectories in land use pattern and grassland ecosystem services in two European alpine landscapes. Reg Environ Change 17, 2251–2264. https://doi.org/10.1007/s10113-017-1207-4 Lexerød, N.L., Eid, T., 2006. Assessing suitability for selective cutting using a stand level index. For. Ecol. Manage. 237, 503–512. https://doi.org/10.1016/j.foreco.2006.09.071 LIFE09 NAT/ES/000513. Conservation of the Cantabrian Capercaille (Tetrao urogallus cantabricus) in its habitat in the Cantabrian Mountain range. European Commission, 2016. McCracken, M.E., Woodcock, B.A., Lobley, M., Pywell, R.F., Saratsi, E., Swetnam, R.D., Mortimer, S.R., Harris, S.J., Winter, M., Hinsley, S., Bullock, J.M., 2015. Social and ecological drivers of success in agri-environment schemes: the roles of farmers and environmental context. J Appl Ecol 52, 696–705. https://doi.org/10.1111/1365-2664.12412 Magurran, A.E., 2004. Measuring biological diversity. Oxford, UK: Blackwell Science. Mahapatra, A.K., Tewari, D.D., Baboo, B., 2015. Displacement, Deprivation and Development: The Impact of Relocation on Income and Livelihood of Tribes in Similipal Tiger and Biosphere Reserve, India. Environ Manage 56, 420–432. https://doi.org/10.1007/s00267-015-0507-z Martín-Alcón, S., Coll, L., Salekin, S., 2015. Stand-level drivers of tree-species diversification in Mediterranean pine forests after abandonment of traditional practices. For. Ecol. Manage. 353, 107–117. http://dx.doi.org/10.1016/j.foreco.2015.05.022 26 Mathevet, R., Thompson, J.D., Folke, C., Chapin, F.S., 2016. Protected areas and their surrounding territory: socioecological systems in the context of ecological solidarity. Ecol Appl 26 (1), 5–16. https://doi.org/10.1890/14-0421 Meyer, W.J., Crumley, C.L., 2011. Historical Ecology: Using what Works to Cross the Divide. In: Moore, T., Armada, L. (Ed.), Atlantic Europe in the First Millennium BC: Crossing the Divide. Oxford: Oxford University Press. Mikoláš, M., Svitok, M., Tejkal, M., Leitão, P.J., Morrissey, R.C., Svoboda, M., Seedre, M., Fontaine, J.B., 2015. Evaluating forest management intensity on an umbrella species: Capercaillie persistence in central Europe. For. Ecol. Manage. 354, 26–34. http://dx.doi.org/10.1016/j.foreco.2015.07.001 Ministerio de Medio Ambiente, 1997–2007. Tercer Inventario Forestal Nacional. Madrid: Dirección General para la Biodiversidad. Mod, H.K., Scherrer, D., Luoto, M., Guisan, A., 2016. What we use is not what we know: environmental predictors in plant distribution models. J. Veg. Sci. 27, 1308–1322. https://doi.org/10.1111/jvs.12444 Neumann, M., Starlinger, F., 2001. The significance of different indices for stand structure and diversity forests. For. Ecol. Manage. 145, 91–106. https://doi.org/10.1016/S0378-1127(00)00577- 6 OECC, 2012. Oficina Española de Cambio Climático. Boletín de la Red de Seguimiento del Cambio Global nº2. https://www.miteco.gob.es/en/red-parques-nacionales/red- seguimiento/rcg_boletin_02_tcm38-59615.pdf (accessed 30 October 2018). Oldekop, J.A., Holmes, G., Harris, W.E., Evans, K.L., 2016. A global assessment of the social and conservation outcomes of protected areas. Conserv. Biol. 30 (1), 133–141. https://doi.org/10.1111/cobi.12568 Pardo, F., Gil, L., 2005. The impact of traditional land use on woodlands: a case study in the Spanish Central System. J. Hist. Geogr 31, 390–408. https://doi.org/10.1016/j.jhg.2004.11.002 Parrotta. J.A., Trosper, R.L. (Eds.), 2012. Traditional Forest-Related Knowledge: Sustaining Communities, Ecosystems and Biocultural Diversity. Springer, Dordrech, New York. https://doi.org/10.1007/978-94-007-2144-9 Peña, L., Amezaga, I., Onaindia, M., 2011. At which spatial scale are plant species composition and diversity affected in beech forests? Ann. For. Sci. 68, 1351–1362. https://doi.org/10.1007/s13595-011-0147-2 Pérez-Bustamante, R., Baró, J. 1988. El gobierno y la administración de los pueblos de Cantabria. I Liébana. Santander: Diputación Regional de Cantabria. 27 Pérez-Soba, I., Solá, M.A., 2004. La tragedia de los comunales: Legalidad y realidad de los montes comunales en España. Estudios Agrosociales y Pesqueros 203, 187–232. PNOA, 2015. Plan Nacional de Ortofotografía Aérea. http://centrodedescargas.cnig.es/CentroDescargas/buscadorCatalogo.do?codFamilia=LIDAR (accessed 26 April 2016). Pommerening, P., 2002. Approaches to quantifying forest structures. Forestry 75(3), 305–324. https://doi.org/10.1093/forestry/75.3.305 Quinn, G.P., Keough, M.J. 2002. Experimental design and data analysis for biologists. Cambridge: Cambridge University Press. Ren, X., Yang, G., Qin, X., Bai, Y., Liu, Z., Zhao, S., Wang, D., 2012. Effects of environmental factors on species distribution and diversity in an Abies fargesii-Betula utilis mixed forest. Sheng Tai Xue Bao 32(2), 0605–0613. https://doi.org/ 10.5846/stxb201011291695 Reques, P. 1997. Población y territorio en Cantabria. Santander: Universidad de Cantabria. Rescia, J.A., Pons, A., Lomba, I., Esteban, C., Dover, J.W., 2008. Reformulating the social- ecological system in a cultural rural mountain landscape in the Picos de Europa region (northern Spain). Landsc Urban Plan 88, 23–33. https://doi.org/10.1016/j.landurbplan.2008.08.001 Root‐Bernstein, M., Svenning, J‐C., 2018. Human paths have positive impacts on plant richness and diversity: A meta‐analysis. Ecol Evol. 00, 1–11. https://doi.org/10.1002/ece3.4578 Rotherham, I.D. (Ed.), 2013. Cultural Severance and the Environment. The Ending of Traditional and Customary Practice on Commons and Landscapes Managed in Common. Springer, Dordrecht. Saeki, I., 2007. Effects of tree cutting and mowing on plant species composition and diversity of the wetland ecosystems dominated by the endangered maple, Acer pycnanthum. For. Ecol. Manage. 242, 733–746. https://doi.org/10.1016/j.foreco.2007.02.009 Sandom, C.J., Ejrnæs, R., Hansen, M.D.D., Svenning, J-C., 2014. High herbivore density associated with vegetation diversity in interglacial ecosystems. PNAS 111(11), 4162–4167. https://doi.org/10.1073/pnas.1311014111 Shultis, J., Heffner, S., 2016. Hegemonic and emerging concepts of conservation: a critical examination of barriers to incorporating Indigenous perspectives in protected area conservation policies and practice. J SUSTAIN TOUR 24(8-9), 1227–1242. http://dx.doi.org/10.1080/09669582.2016.1158827 Silva-Flores, R., Hernández-Díaz, J.C., Wehenkel, C., 2016. Does community-based forest ownership favour conservation of tree species diversity? A comparison of forest ownership regimes in the Sierra Madre Occidental, Mexico. For. Ecol. Manage. 363, 218–228. http://dx.doi.org/10.1016/j.foreco.2015.12.043 28 Spanish Government, 2005. Boletin Oficial del Estado 03/06/2005. https://www.boe.es/boe/dias/2005/06/03/pdfs/A18636-18637.pdf (accessed 5 March 2018) Spanish Government, 2015. Boletin Oficial del Estado 24/03/2015. https://www.boe.es/buscar/pdf/2015/BOE-A-2015-4458-consolidado.pdf (accessed 5 March 2018) StataCorp, 2013. Stata Statistical Software: Release 13. College Station, TX: StataCorp LP. Symonds, M.R.E., Moussalli, A., 2011. A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike’s information criterion. Behav Ecol Sociobiol 65, 13–21. https://doi.org/10.1007/s00265-010-1037-6 Tárrega, R., Calvo, L., Marcos, E., Taboada, A., 2006. Forest structure and understory diversity in Quercus pyrenaica communities with different human uses and disturbances. For. Ecol. Manage. 227, 50–58. https://doi.org/10.1016/j.foreco.2006.02.008 Viedma, O., Moity, N., Moreno, J.M., 2015. Changes in landscape fire-hazard during the second half of the 20th century: Agriculture abandonment and the changing role of driving factors. Agric Ecosyst Environ 207, 126–140. https://doi.org/10.1016/j.agee.2015.04.011 Vorčák, J., Merganič, J., Saniga, M., 2006. Structural diversity change and regeneration processes of the Norway spruce natural forest in Babia hora NNR in relation to altitude. Journal of Forest Science 52(9), 399–409. https://doi.org/10.17221/4520-JFS Voth, A., 2007. National Parks and Rural Development in Spain. In: Mose, I. (Ed.), Protected Areas and Regional Development in Europe. Ashgate, Aldershot, UK. West, P., Igoe, J., Brockington, D., 2006. Parks and Peoples: The social impact of Protected Areas. Annu Rev Anthropol 35, 251–277. https://doi.org/10.1146/annurev.anthro.35.081705.123308 Wilshusen, P.R., Brechin, S.R., Fortwangler, C.L., West, P.C., 2002. Reinventing a Square Wheel: Critique of a Resurgent ‘‘Protection Paradigm’’ in International Biodiversity Conservation. Soc Nat Resour 15, 17–40. https://doi.org/10.1080/089419202317174002 Wissel, C., 1992. Modelling the mosaic cycle of a Middle European beech forest. Ecol Modell 63, 29–43. https://doi.org/10.1016/0304-3800(92)90060-R Xu, J., Melick, D.R., 2007. Rethinking the Effectiveness of Public Protected Areas in Southwestern China. Conserv. Biol. 21(2), 318–328. https://doi.org/10.1111/j.1523- 1739.2006.00636.x Zuur, A.F., Ieno, E.N., Smith, G.M., 2007. Statistics for Biology and Health: Analysing Ecological Data. Springer, New York. 29