Goffredo Filibeck
email: filibeck@unitus.it@tel.390761357215

Marta G Sperandii

Manuele Bazzichetto

Leone D Mancini

Francesco Rossini

Laura Cancellieri

Exploring the drivers of vascular plant richness at very fine spatial scale in sub-Mediterranean limestone grasslands (Central Apennines, Italy)

Keywords: Dry grasslands, Diversity metrics, Multi-scale approach, Soil properties, Species-area relationship

Secondary dry grasslands in Europe can host high levels of vascular plant richness at small spatial scales.

However, in Southern Europe their biodiversity patterns are largely unexplored. In this work, we aim at: i) estimating plant species richness patterns at very fine scales in montane dry grasslands, on limestone bedrock, in Abruzzo Lazio & Molise National Park (Central Apennines, Italy); ii) assessing the most important physical and edaphic drivers of biodiversity patterns at multiple plot sizes. We used randomly placed nested-plot series where we measured alpha-diversity at three different plot sizes (1 m 2 , 0.1 m 2 and 0.01 m 2 ) and within-plot beta-diversity (as expressed by the slope of the species-area curve across plot sizes). Variable selection was performed by means of Random Forests. Relationships between selected variables and diversity measures were then assessed using Regression Trees, Linear and Generalized Linear Models. Overall, results pointed to topographically-controlled edaphic factors (soil pH and silt fraction) as the main drivers positively influencing alpha-diversity at all spatial scales, with a positive effect of rock cover and slope inclination at smaller spatial grains. Beta-diversity was positively influenced by rock cover.

We suggest that high-pH, steep and/or rocky sites feature higher species richness because they lack A c c e p t e d m a n u s c r i p t competitive grass species. Our results are in agreement with previous works underlining the importance of less productive habitats for the conservation of secondary grassland biodiversity.

Introduction

Fine-scale species richness (species density or alpha-diversity) of vascular plants, i.e. the number of species occurring in fixed-area plots of small size (usually ≤ 100 m 2 ), and the quest for its environmental determinants have been the subject of important studies and considerable discussion by plant ecologists (reviewed e.g. by Grace 1999). A considerable number of such works focused on European secondary grasslands (e.g. Grime 1973;Al-Mufti et al. 1977;[START_REF] Van Der Maarel | Small-scale plant species turnover in a limestone grassland: the carousel model and some comments on the niche concept[END_REF]Pärtel and Zobel 1999;[START_REF] Bennie | Influence of slope and aspect on long-term vegetation change in British chalk grasslands[END_REF][START_REF] Löbel | Species richness of vascular plants, bryophytes and lichens in dry grasslands: The effects of environment, landscape structure and competition[END_REF][START_REF] Turtureanu | Scaleand taxon-dependent biodiversity patterns of dry grassland vegetation in Transylvania[END_REF][START_REF] Reitalu | Determinants of fine-scale plant diversity in dry calcareous grasslands within the Baltic Sea region[END_REF][START_REF] Palpurina | Patterns of fine-scale plant species richness in dry grasslands across the eastern Balkan Peninsula[END_REF]Güler et al. 2016). These are semi-natural habitats originated and maintained by disturbance such as livestock grazing, in areas that would be otherwise covered by forest vegetation; however, they are not man-made ecosystems, as livestock disturbance replaces the grazing and browsing by wild prehistoric megaherbivores and by extant ungulates (Pärtel et al. 2005;Feurdean et al. 2018). European secondary grasslands, and especially "dry" grasslands (i.e. those growing on free-draining substrates such as limestone), are of outstanding ecological value since they are among the most diverse plant communities in the world at plot sizes <100 m² (Pärtel et al. 2005;Wilson et al. 2012;[START_REF] Chytrý | The most species-rich plant communities in the Czech Republic and Slovakia (with new world records)[END_REF][START_REF] Roleček | Understanding the extreme species richness of semi-dry grasslands in east-central Europe: a comparative approach[END_REF]. Moreover, they are an essential part of the cultural landscape of the continent, often resulting from millennia of lowintensity land-use (Poschlod and Wallis DeVries, 2002;[START_REF] Dengler | Biodiversity of Palaearctic grasslands: a synthesis[END_REF]Feurdean et al. 2018); they host a high proportion of Europe's endemic plants [START_REF] Habel | European grassland ecosystems: threatened hotspots of biodiversity[END_REF]; and act as a network of "interglacial refugia" for Pleistocene (ice-age) steppic species [START_REF] Hájková | Prehistoric origin of the extremely species-rich semi-dry grasslands in the Bílé Karpaty Mts (Czech Republic and Slovakia)[END_REF][START_REF] Kajtoch | Phylogeographic patterns of steppe species in Eastern Central Europe: a review and the implications for conservation[END_REF]Cancellieri et al. 2017b). European secondary grasslands are currently threatened by grazing abandonment, conversion into cropland, afforestation and transformation into "improved" grasslands [START_REF] Dolek | Conserving biodiversity on calcareous grasslands in the Franconian Jura by grazing: a comprehensive approach[END_REF][START_REF] Habel | European grassland ecosystems: threatened hotspots of biodiversity[END_REF] A c c e p t e d m a n u s c r i p t landscapes, as they are less studied than their Central European counterparts (Apostolova et al. 2014).

Although richness patterns at very fine spatial scales (≤1 m 2 ) have to be collected and interpreted with caution because of many potential pitfalls (Gotelli and Colwell 2001;[START_REF] Palmer | Artifacts and artifictions in biodiversity research[END_REF]), these are the scales at which -in herbaceous vegetation -plant individuals interact with each other and build up ecosystem functioning and services (Grime 1998;[START_REF] Dengler | Biodiversity of Palaearctic grasslands: a synthesis[END_REF][START_REF] Borer | Herbivores and nutrients control grassland plant diversity via light limitation[END_REF]). In addition, using plots ≤1 m 2 ensures that the sampling units are small in relation to the spatial scale of variation in topography and soil heterogeneity (Grime 1973).

A c c e p t e d m a n u s c r i p t

In the Italian Peninsula, extensive surveys of dry grassland biodiversity patterns and their edaphic drivers at spatial grains ≤1 m 2 are lacking: the only paper we are aware of, that analyzed richness patterns and edaphic variables at fine scales [START_REF] Maccherini | Factors associated with species richness in a remnant calcareous grassland[END_REF], was focused on a very small study area (cf. also [START_REF] Dengler | Biodiversity of Palaearctic grasslands: a synthesis[END_REF][START_REF] Chelli | Climate change response of vegetation across climatic zones in Italy[END_REF]. Nevertheless, the large areas of secondary grasslands in the Apennines chain (the mountain range forming the backbone of the Italian peninsula) are extremely interesting, because of considerable climatic and edaphic gradients (Primi et al. 2016), biogeographical position (e.g. [START_REF] Rosati | The forest vegetation of the Alburni Mts. (Cilento N. P., Campania, Southern Italy)[END_REF][START_REF] Blasi | Submediterranean dry grasslands along the Tyrrhenian sector of central Italy: Synecology, syndynamics and syntaxonomy[END_REF], conservation values (e.g. [START_REF] Manzi | Storia dell'ambiente nell'Appennino centrale: la trasformazione della natura in Abruzzo dall'ultima glaciazione ai giorni nostri[END_REF]Primi et al. 2016;Cancellieri et al. 2017b) and a long history of sheep grazing and transhumance [START_REF] Barker | Ancient and modern pastoralism in central Italy: an interdisciplinary study in the Cicolano mountains[END_REF][START_REF] Brown | Multi-proxy study of Holocene environmental change and human activity in the Central Apennine Mountains, Italy[END_REF].

The present paper thus aims at analysing fine-scale diversity patterns and drivers from a network of secondary dry grasslands within the Central Apennine Mountains. We focused on the karst plateaux of the upper-montane bioclimatic belt: these habitats are characterized by a high topographical heterogeneity at micro-and meso-scale, thus featuring steep gradients of the potential controlling variables for species richness. We used randomized nested plots of three different sizes (0.01-0.1-1 m 2 ) to answer the following questions: i) what are the levels and patterns of very fine-scale species richness at multiple spatial grains (alpha-diversity)? ii) what are the patterns of spatial turnover of species across the considered spatial grains (within-plot beta-diversity)? iii) what are the main environmental drivers of these diversity metrics, and do they vary between spatial scales?

Methods

Study area

Data were gathered within the "Abruzzo Lazio e Molise" National Park, in the Central Apennine mountains, Italy (41°48'N, 13°50'E). This protected area was designated in 1922 because it was the last refuge for both the Marsican brown bear (Ursus arctos ssp. marsicanus) and the Apennine chamois (Rupicapra pyrenaica ssp. ornata); at the same time, this area was until the 1950's one of the most important centers in Italy for transhumant shepherding [START_REF] Sievert | Abruzzo National Park: Land of Dreams[END_REF]. The Park's flora comprises >2,000 species, including >30 taxa endemic to the Central Apennine mountains (Conti and Bartolucci 2015); most of the Park's grasslands A c c e p t e d m a n u s c r i p t correspond to ecosystem types protected under the EU Habitats Directive (Primi et al. 2016). Five patches of secondary dry grassland were chosen as study sites, for a total area of c. 500 ha (Fig. 1; additional information are given in Online Resource 1). Their elevation ranges between c. 1350-1750 m a.s.l. and they thus belong to the upper-montane phytoclimatic belt of the Apennines (Gerdol et al. 2008). All the studied grassland patches are very similar: they lie on Mesozoic limestone bedrock [START_REF] Bigi | Carta Geologica del Parco Nazionale d'Abruzzo[END_REF]) and are located in large bowl-shaped depressions, hosting in turn a system of different smaller landforms (see photos in Online Resource 1) originated by complex interactions between karst, tectonic and glacial processes [START_REF] Bigi | Carta Geologica del Parco Nazionale d'Abruzzo[END_REF]Galadini and Messina 1990;[START_REF] Cinque | Le tracce glaciali nel Parco Nazionale d'Abruzzo (Appennino centrale): nota preliminare[END_REF]). Such micro-and meso-scale landforms include shallow dolines, sinkholes, hillocks, rocky outcrops, karren, etc. Because of karst hydrology, no surface water bodies occur. The grassland patches are surrounded by beech (Fagus sylvatica) forest, which is considered the potential vegetation of the area (Filibeck et al. 2015). Within the considered elevation belt, dominant species in the Park's grasslands include Avenula praetutiana, Bromus erectus, Festuca circummediterranea, Koeleria lobata and Phleum ambiguum on shallow to medium-depth soils; while on deep soils (e.g. in karst hollows) Agrostis capillaris, Festuca microphylla and Nardus stricta often prevail [START_REF] Bruno | Note illustrative alla carta della vegetazione del Parco Nazionale d'Abruzzo[END_REF]. No phtyosociological studies targeted the study sites; however, relatively similar habitats in neighbouring areas were described by [START_REF] Lucchese | I prati a Bromus erectus Hudson dell'Appennino Laziale[END_REF], [START_REF] Blasi | On the vegetation series in the northern sector of the Simbruini mountains (Central Apennines)[END_REF] and Di [START_REF] Pietro | A geobotanical survey on acidophilous grasslands in the Abruzzo, Lazio and Molise National Park (Central Italy)[END_REF].

The local climate can be defined as sub-Mediterranean (Filibeck et al., 2015): it features high annual precipitation (c. 1,200-1,500 mm/yr), with a rainfall maximum in autumn and a very sharp drop in summer.

This precipitation regime leads to significant drought stress, especially for grassland plants, during July and August [START_REF] Catorci | Effect of spatial and temporal patterns of stress and disturbance intensities in a sub-Mediterranean grassland[END_REF]Primi et al. 2016). Average annual temperature is between 6-8 °C; winters are characterized by marked frost and abundant snowfall (see Filibeck et al. 2015, for climate diagrams and details).

In the study area, secondary grasslands are maintained by grazing. Transhumant sheep grazing was widely practiced until the 1950's, when ovine farming started to decrease for socio-economic reasons; sheep stocking rate is nowadays drastically lower than in the early twentieth century [START_REF] Manzi | Storia dell'ambiente nell'Appennino centrale: la trasformazione della natura in Abruzzo dall'ultima glaciazione ai giorni nostri[END_REF]. The decrease of sheep husbandry has resulted in a steep increase in cattle and horse farming; such animals are usually A c c e p t e d m a n u s c r i p t raised for meat production only and left free-ranging night and day in the wild (Primi et al. 2016). Because of summer drought, mowing has never been a large-scale management system in the region [START_REF] Manzi | La gestione dei pascoli montani in Abruzzo e la Società delle Erbe Seconde di Pescasseroli ed Opi[END_REF], and is currently not practiced at any of our study sites.

Data collection

Data were collected in early summer in 2013 and 2014. Sampling design was based on a grid of 300 x 300 m cells, overlaid on the National Park area, with one random point selected within each cell (tessellation stratified sampling: see e.g. Barabesi et al. 2012;[START_REF] Chiarucci | Scale dependence of plant species richness in a network of protected areas[END_REF]. The five grassland clearings to be investigated were chosen a priori, following field reconnaissance and inspection of detailed geological maps that showed they had comparable features. All the random points falling within the five selected grassland patches were retained, resulting in a total of 50 points (Fig. 1). At each GPS-located point (hereafter "sampling point"), vascular plants were surveyed using a nested-plot series based on the design proposed by [START_REF] Dengler | A flexible multi-scale approach for standardised recording of plant species richness patterns[END_REF], as modified by Cancellieri et al. (2017a): one 1 m 2 squared plot was centred on the random point, and oriented N-S; two 0.1 m 2 subplots were placed at the NE and SW corners of the 1 m 2 plot, and two 0.01 m 2 subplots were placed within each of the 0.1 m 2 plots (i.e., there are four 0.01 m 2 subplots at each sampling point) (Fig. 2 and Online Resource 1). The rationale for sampling a higher number of subplots of smaller areas within each nested-plot series is twofold [START_REF] Dengler | A flexible multi-scale approach for standardised recording of plant species richness patterns[END_REF]): i) it allows to better account for the micro-environmental heterogeneity contained within the largest plot; ii) it allows the use of mean species richness values within each nested series, thus liberating the resulting species-area curve from stochasticity. Moreover, while species counts often do not follow a normal distribution, averaging alpha-diversity values among sub-plots leads to normally distributed data (e.g. [START_REF] Quinn | Experimental Design and data Analysis for Biologists[END_REF].

Each nested-plot series was precisely delimited in the field by means of a plastic frame, and was surveyed by at least two people to reduce inter-observer variability (e.g. [START_REF] Morrison | Observer error in vegetation surveys: a review[END_REF]. For each plot and subplot, we listed all vascular plant species (presence/absence data only), following the "shoot presence" criterion: that is, a species is recorded as present if the vertical projection of any above-ground part of an individual falls within the plot (Cancellieri et al. 2017a). Plant specimens that could not be identified in the field were A c c e p t e d m a n u s c r i p t collected and determined in the laboratory. Identification was mostly based on Pignatti (1982); plant names were subsequently standardized following [START_REF] Conti | An Annotated Checklist of the Italian Vascular Flora[END_REF]. Given the still poorly-known taxonomy and distribution of the genus Festuca in Italy, and considering that this taxon is one of the main builders of the local plant communities, all Festuca "species" within each nested-plot series were collected for identification in the laboratory through leaf section, following Foggi andRossi (1996) andFoggi et al. (1999).

For each nested-plot series, the following environmental variables were measured at 1-m 2 spatial scale (Tab. 1). Elevation was measured with a GPS receiver. Slope inclination was measured by a clinometer.

Slope aspect was measured with a compass and used to calculate heat-load index (based on latitude, aspect and inclination) following eq. 2 in [START_REF] Mccune | Equations for potential annual direct incident radiation and heat load[END_REF]. Soil depth was measured at 5 random points with a 65-cm long iron rod (soils deeper than rod length were arbitrarily given a value of 65 cm); the median of the 5 measurements was used in the analyses. Percentage of stone and rock cover was visually estimated. Land-form position was visually assessed as a categorical variable as follows: "hollow" (indicating the flat or gently sloping position in concave morphologies such as the bottom of dolines, other karst depressions or gullies), "hillside" (for the sloping sides of both meso-and macro-scale landforms), and "hump" (for convex landforms such as hillocks, ridges or rocky outcrops). Soil samples were collected with a hand auger (Eijkelkamp, NL) from the upper 10-cm layer as a mixed sample from 3 randomly placed points within each 1 m 2 plot, and subsequently air-dried and crushed to pass a 2-mm sieve. Soil pH was measured in water (soil-solution suspension of 1:2.5). Soil texture (clay, silt and sand fraction) was determined through the standard sedimentation techniques (Gee and Bauder 1986) adopting the ISSS dimensional limits, without pre-treatment with H2O2. Soil organic matter was determined using Walkley-Black method.

The plot data are available in GrassPlot, the Database of Scale-Dependent Phytodiversity Patterns in Palaearctic Grasslands (GIVD ID EU-00-003) [START_REF] Dengler | GrassPlot-a database of multi-scale plant diversity in Palaearctic grasslands[END_REF].

Data analysis

For each nested-plot series, the following diversity metrics were calculated from the collected data, to be used as response variables in the statistical analysis: species richness of the whole 1 m 2 plot (hereafter, "1 A c c e p t e d m a n u s c r i p t m 2 alpha-diversity"); mean species richness of the two 0.1 m 2 subplots (hereafter, "0.1 m 2 alpha-diversity"); mean species richness of the four 0.01m 2 subplots (hereafter, "0.01m 2 alpha-diversity"); and "z-value", i.e. the slope of the species-area relationship calculated from the three spatial scales at each nested-plot series.

We assumed for the species-area curve the Arrhenius (1921) power function: S = c A z and we calculated the z-values with linear regression in double-log space (e.g. [START_REF] Chiarucci | Effects of productivity on species-area curves in herbaceous vegetation: evidence from experimental and observational data[END_REF], i.e. using the equation: log S = log c + z log A (where S = number of species, A = plot size; c and z are the estimated parameters, i.e. the intercept at A=1 and the slope of the regression line, respectively).

Fitting of the species-area function was performed for each nested-plot series, using for the 0.1 m 2 size the mean richness values of the two subplots, and for the 0.01 m 2 size the mean richness values of the four subplots. The z-values of nested-plot series are useful for comparing the degree of spatial species turnover between different ecological situations (e.g. [START_REF] Chiarucci | Effects of productivity on species-area curves in herbaceous vegetation: evidence from experimental and observational data[END_REF][START_REF] De Bello | Grazing effects on the species-area relationship: variation along a climatic gradient in NE Spain[END_REF][START_REF] Dengler | Pitfalls in small-scale species-area sampling and analysis[END_REF][START_REF] Turtureanu | Scaleand taxon-dependent biodiversity patterns of dry grassland vegetation in Transylvania[END_REF]Polyakova et al. 2016). They are a measure of beta-diversity (i.e. of the variation in the identities of species between plots) at fine spatial scale, as the shape of the species area-curve reflects the level of spatial homogeneity of species distribution and community composition [START_REF] Storch | The theory of the nested species-area relationship: geometric foundations of biodiversity scaling[END_REF].

Collinearity between the predictor variables, as well as the strength of the univariate relationships between each predictor variable and the diversity metrics, was preliminarily checked using pairwise Spearman correlations. We then excluded the percentage of sand from further analyses because by definition the three soil texture fractions sum up to 100 and, of the three texture components, sand exhibited the highest degree of correlation with other predictor variables. Since the study area had been deliberately defined within a specific elevation belt, and the range of altitudinal variation among plots was thus small, elevation was not included in the analyses. One sampling point presented outlying values from soil analysis, leading to suspect an error occurred in laboratory processing, and was excluded from statistical analysis on diversity-environment relationships, that are thus based on n=49.

A c c e p t e d m a n u s c r i p t

We performed our main statistical analysis through Recursive Partitioning approaches ("Regression Trees" and "Random Forests": Strobl et al. 2009a). Recursive Partitioning methods have many desirable properties that make them attractive for the ecologist [START_REF] De'ath | Classification and regression trees: a powerful yet simple technique for ecological data analysis[END_REF]: they can successfully handle very high number of predictors, are robust to multicollinearity, independent from any assumptions on the distribution of errors and can successfully handle and reveal complex interactions and non-linear relationships [START_REF] Sandri | Variable selection using random forests[END_REF][START_REF] Hapfelmeier | A new variable selection approach using random forests[END_REF].

As a first step, in order to identify the most relevant explanatory variables we used Random Forests (i.e. an ensemble method combining several Classification or Regression Trees), which several studies highlighted as an effective variable selection approach [START_REF] Sandri | Variable selection using random forests[END_REF][START_REF] Hapfelmeier | A new variable selection approach using random forests[END_REF][START_REF] Svitok | Environmental thresholds and predictors of macrophyte species richness in aquatic habitats in central Europe[END_REF]. Random Forests were run using R package "party" (function "cforest": Hothorn et al., 2006b;[START_REF] Strobl | Bias in Random Forest Variable Importance Measures: Illustrations, Sources and a Solution[END_REF][START_REF] Strobl | Conditional variable importance for random forests[END_REF], as it implements conditional inference trees which are designed to reduce the bias that such tools otherwise have towards highly correlated variables [START_REF] Strobl | Conditional variable importance for random forests[END_REF](Strobl et al. , 2009b)). The number of "trees" in the "forests" was set to 1000, while the number of randomly preselected predictor variables for each split was set to 3, following Strobl et al. (2009a). As Random Forests are, by definition, "random", repeating the process (i.e. re-running the algorithm) is recommended in order to ensure the robustness and stability of the model (Strobl et al. 2009a;Shih 2011). For this reason, for each of the four response variables a Random Forest was run 10 times. For each Random Forest, a conditional measure of variable importance was computed for all predictors using function "varimp". We defined as "important" and retained for further analysis only those environmental variables which, in at least six out of the 10 Random Forest tries, featured positive values exceeding the absolute value of the most negative score (cf. Strobl et al. 2009a).

As a second step of analysis, individual Regression Trees were built for each response variable, including as predictors those environmental variables identified by the Random Forests as important for each biodiversity measure. Regression Trees have recently proven to be very effective in studying plant diversity drivers (e.g. [START_REF] Chytrý | Plant species richness in continental southern Siberia: effects of pH and climate in the context of the species pool hypothesis[END_REF][START_REF] Palpurina | Patterns of fine-scale plant species richness in dry grasslands across the eastern Balkan Peninsula[END_REF][START_REF] Svitok | Environmental thresholds and predictors of macrophyte species richness in aquatic habitats in central Europe[END_REF]. We used them to identify the environmental characteristics most strongly associated with each response variable, assess the sign of the relationship, detect interactions among predictors and estimate the environmental thresholds for species A c c e p t e d m a n u s c r i p t richness [START_REF] De'ath | Classification and regression trees: a powerful yet simple technique for ecological data analysis[END_REF]Strobl et al. 2009a). We performed Regression Tree analysis with package "party" in R (function "ctree", Hothorn et al. 2006a), implementing "Conditional Inference" Trees: these were preferred to the more classical approach, based on the "CART" algorithm, as the latter may be biased towards continuous predictors measured over large scales or categorical predictors with many categories (Hothorn et al. 2006a;Strobl et al. 2009a). Besides, "Conditional Inference" Trees, implementing a statistically robust stopping criterion, do not suffer from overfitting (Strobl et al. 2009a). The algorithm tests the global null hypothesis of independence between any of the input variables and the response; it stops if this hypothesis cannot be rejected, otherwise selects the input variable with strongest association to the response. This association is measured by a p-value corresponding to a test for the partial null hypothesis of a single input variable and the response. The stop criterion is by default based on multiplicity adjusted p-values (Bonferroni correction), but univariate p-values can be alternatively used. This approach ensures that the right sized "tree" is grown and no "pruning" or cross-validation is needed (Hothorn et al. 2006a).

Additionally, we assessed the strength of the relationship between the response variables and the selected predictors through Generalized Linear Models (GLM) and multiple Linear Models. Whenever two predictors (among those selected through the Random Forest output) were highly correlated with each other (i.e. had a Spearman correlation coefficient > 0.6), we retained only the most biologically informative one. Since the explanatory variables are measured in different units, standardized regression coefficients (obtained by subtracting, from each observation, the variable mean and dividing by the standard deviation: e.g. Quinn and Keough 2002) were computed to allow comparison of the effect sizes on the response variable ; for Linear Models we standardized both explanatory and response variables, while for Generalized Linear Models we only transformed predictors. For 1 m 2 alpha-diversity, a Poisson Generalized Linear Model was fitted, setting "log" as link function. Overdispersion was tested using R package "AER" (function "dispersiontest": Kleiber and Zeileis 2008). We used an adjusted Pseudo-R 2 (Guisan and Zimmermann 2000)

to quantify the improvement in the deviance explained by the model including the selected predictors over the intercept-only model (R package "modEvA", function "Dsquared": [START_REF] Barbosa | modEvA: Model Evaluation and Analysis[END_REF]). For 0.1 m 2 and 0.01 m 2 alpha-diversity, as they were obtained averaging the species counts among subplots of the same A c c e p t e d m a n u s c r i p t size within each nested series, we assumed a Gaussian distribution of errors: thus multiple Linear Models were fitted and adjusted-R 2 was used to quantify explained variability. Z-values were log-transformed before model fitting in order to meet linear-model assumptions.

Finally, to help ecological interpretation of the richness patterns and elucidate their relationship with floristic composition, a PCA (based on a correlation matrix) of the plot × species matrix was run, adding the 1 m 2 alpha-diversity as a covariate, in order to underline which taxa were associated with different levels of species richness.

Note that in the results section we shall use p-values as a continuous measure of fit of the null hypothesis to the data, instead of using them for a "significant/non-significant" dichotomy (Wassertein and Lazar 2016;Greenland et al. 2016).

Results

Overall diversity

Across all 50 sampling points, a total of 198 vascular plant species were recorded; 64 species (32% of the total recorded flora) were found at only one nested series (singletons). Mean species richness per plot was 10.3, 18.5 and 28.0 for 0.01 m 2 , 0.1 m 2 and 1 m 2 , respectively (Tab.2). The mean slope (z) of the speciesarea relationship was 0.238 (Tab.2); i.e., richness increased by c. 1.7 times, on average, with a ten-fold increase in plot area. The most frequent species (in terms of proportion of sampling plots where the taxon was recorded) are listed in Tab. 3. The full 1m²-plot × species matrix (including the frequency values of all species) is provided as Online Resource 3.

Diversity-environment relationships

Descriptive statistics of the environmental variables are shown in Tab.1. Pairwise Spearman's correlation matrix (Tab.4) highlighted that many predictor variables were correlated with each other: most notably, soil pH was positively correlated with stone and rock cover (rs=0.64) and with slope inclination (rs=0.48), and negatively correlated with soil depth (rs= -0.61); stone and rock cover was positively related to slope inclination (rs=0.54). As for the relationships between predictor and dependent variables, species richness A c c e p t e d m a n u s c r i p t at 1 m 2 spatial scale (1 m 2 alpha-diversity) was positively correlated with silt % (rs=0.38), soil pH (rs=0.37) and stone and rock cover (rs=0.30). However, all these relationships became weaker moving to 0.1 and 0.01 m 2 spatial scales. The z-values were positively correlated with stone and rock cover (rs=0.58); in addition, they showed a positive correlation with slope inclination (rs=0.37) and a negative relationship with soil depth (rs=-0.48).

In Tab. 5 it is reported, for each of the four response variables, how many times, out of the 10 Random Forest tries, each predictor variable turned out to be significant. Both soil pH and silt % were found to be significant predictors across all 10 runs for all three alpha-diversity spatial scales. Stone and rock cover and slope inclination were stable outcomes for 1 m 2 and 0.1 m 2 ; soil depth and land-form position were constantly significant only at 1 m 2 spatial scale. As for z-values, soil pH was never a significant predictor, while land form, stone and rock cover, inclination and soil depth were significant in all runs. Heat-load index was never significant for any of the dependent variables.

The Regression Tree for 1 m 2 alpha-diversity identified a significant threshold value for pH at 6.4, i.e. the plots with soil pH smaller than this value hosted, on average, significantly lower richness. The model further partitioned the group of high-pH sites based on soil texture, showing that the maximum species richness is attained at sites featuring simultaneously high pH values and high silt % (fig. 3). However, as shown above, within the study area soil pH is positively correlated with slope inclination and stone and rock cover, and negatively correlated with soil depth; thus, high pH values can also indicate steep slopes or rocky humps with shallow soil. For 0.1 m 2 and 0.01 m 2 alpha-diversity, Regression Tree analysis performed with "ctree" default rules didn't yield any statistically significant partition. If univariate p-values were used (see Methods) instead of the stop criterion based on Bonferroni correction (cf. Nakagawa 2004), then the resulting tree for the 0.1-m 2 sub-plots highlighted the role of slope inclination and soil pH (and their interaction). Steeper slopes (indicating also stony soils, as explained above) showed greater species richness at 0.1 m 2 spatial scale. The lowest richness was found when a low slope inclination was combined with a pH value lower than 6.2, indicating flat bottoms of concave landforms with deep, sub-acid soils (fig. 4). As for the 0.01 m 2 alphadiversity, no significant partitions were detected even after the Bonferroni correction was removed.

A c c e p t e d m a n u s c r i p t

In the Regression Tree for z-values, stone and rock cover was identified as a strongly significant predictor, with two subsequent partitions based on the same variable (fig. 5), thus underlining a regular increase of the slope of the species-area curve with increasing soil rockiness, and showing that low z-values were consistently found when stone and rock cover was less than 3%.

For the 1 m 2 alpha-diversity GLM, the predictor variables selected on the basis of Random Forest output were soil pH, land-form, silt %, slope inclination, soil depth and stone and rock cover. However, as both soil depth and stone and rock cover were highly correlated with soil pH (Spearman's rs > 0.6), they were not included in the model. The inclusion of pH as quadratic term (based on a preliminary graphical exploration, that showed a slightly hump-shaped trend) did not improve significantly the fit of the models as checked via likelihood ratio test for the GLM and ANOVA for the Linear Models; therefore we retained only the firstorder term for pH in the final models. GLM results (Tab. 6) showed that 1 m 2 alpha-diversity was positively associated with silt % (standardized coefficient=0.13, p<0.001), and was on average higher in "hump" landforms than in hollows and hillsides. Slope inclination and pH were positively, but very weakly, associated with 1 m 2 alpha-diversity (stand. coeff.=0.05, p=0.1, for both predictors); however, as soil pH was markedly correlated with two omitted variables, its parameter estimate may be biased. The adjusted Pseudo-R 2 of the model was 0.53 (i.e. the predictors explained 53% of the deviance in 1 m 2 alpha-diversity).

As for 0.1 m 2 alpha-diversity, the predictor variables selected for fitting the Linear Model were: soil pH, slope inclination, silt %, clay % and stone and rock cover. As stone and rock cover was highly correlated with soil pH (rs>0.6), it was removed from the model. Silt fraction showed a positive relationship with 0.1 m 2 alpha-diversity (stand. coeff.=0.47, p<0.001); a smaller positive effect, with a high error, was showed by slope inclination (stand. coeff.=0.23, p=0.1) and soil pH (stand. coeff.=0.22, p=0.1); for this latter, the same caveats hold as in the 1m 2 model (Tab. 6). The model explained 31% of variability.

For 0.01 m 2 alpha-diversity, variable selection based on Random Forests led to include only silt % and soil pH. Species richness at this spatial scale showed a positive relationship with silt % (stand. coeff.=0.37; p<0.01); a smaller positive effect was showed also by soil pH, but with a high standard error and a large pvalue (stand. coeff.=0.18; p=0.2) (Tab. 6). The model explained only 13% of variability. In the PCA of the 1-m 2 plots × species matrix, species richness was correlated with PC2 (loading = -0.57).

The correlation coefficients of the species on this factor, allowing to highlight the taxa that were consistently associated with either the richest or the poorest plots, are listed in Online Resource 2.

Discussion

Overall diversity

Especially at the finest spatial scales, mean and maximum species richness per plot ranked quite high within the context of Palaearctic grasslands, supporting earlier hypothesis that S-European grasslands can host high levels of alpha-diversity (Apostolova et al. 2014): for instance, our maximum of 20 vascular plant species in 0.01 m 2 quadrats is the second highest among all the Eurasian datasets reported by Dengler et al.

(2016) (cf. also [START_REF] Chytrý | The most species-rich plant communities in the Czech Republic and Slovakia (with new world records)[END_REF]. At 1 m 2 scale, the maximum richness value (41 species) was not exceptionally high (cf. [START_REF] Merunková | White Carpathian grasslands: can local ecological factors explain their extraordinary species richness?[END_REF]Wilson et al. 2012), although still within the "upper range" at European scale [START_REF] Chytrý | The most species-rich plant communities in the Czech Republic and Slovakia (with new world records)[END_REF][START_REF] Dengler | Scale-dependent plant diversity in Palaearctic grasslands: a comparative overview[END_REF]). However, our data stem from a randomized sampling design, while most datasets that contain or approach the European "richness records" originate from plots that were placed subjectively with the aim to record the most species-rich sites (e.g. [START_REF] Merunková | White Carpathian grasslands: can local ecological factors explain their extraordinary species richness?[END_REF]). Beta-diversity as expressed by z-values was within the range typically reported for dry grasslands in Europe, and our mean z-value (0.24) was very close to what found in many previous studies (see e.g. Baumann et al. 2016, and references therein).

The probabilistic sampling used in this study allowed to quantitatively define [START_REF] Chiarucci | To sample or not to sample? That is the question... for the vegetation scientist[END_REF]) the most common species within the landscapes studied (Tab. 3). Among grasses, the most frequent taxon was Festuca circummediterranea, followed by Poa alpina and Koeleria lobata, pointing to an interesting codominance of Mediterranean and Boreal geographic elements. Similarly, among forbs, the Apennine endemic Viola eugeniae turned out to be present in 2/3 of the samples, but the most common species was Hieracium pilosella (found in 80% of the quadrats), a typical example of the European geographic element.

A c c e p t e d m a n u s c r i p t

Diversity-environment relationship: alpha-diversity

Soil pH and silt fraction were positively correlated with species richness at all spatial scales. Slope inclination, stone and rock cover and "hump" land-form were positively influencing species richness as well, but only at 1 m 2 and/or 0.1 m 2 spatial grains, and with weaker effect sizes.

There is a huge body of works on the positive (or unimodal) relationship between soil pH and plot-scale diversity in Eurasian grasslands (see references in: [START_REF] Chytrý | Plant species richness in continental southern Siberia: effects of pH and climate in the context of the species pool hypothesis[END_REF]Palpurina et al. 2017), although in most studies the quadrat size was larger (10-100 m 2 ) than in the present work (but see Grime 1973). This relationship has been explained assuming that during the glacial phases the acidophilous pool of species underwent a much higher extinction rate than the acidifuge contingent because of the unbalanced distribution of soil types across the ice-age biomes), thus leading to a smaller pool of acidophilous than calcicolous species in the present-day regional floras (Grime 1973;Grace 1999;Ewald 2003;[START_REF] Tyler | Some ecophysiological and historical approaches to species richness and calcicole/calcifuge behaviour -contribution to a debate[END_REF][START_REF] Chytrý | Local and regional patterns of species richness in central European vegetation types along the pH/calcium gradient[END_REF]).

However, since many of our environmental variables showed correlated gradients in the field (cf. [START_REF] Bennie | Influence of slope and aspect on long-term vegetation change in British chalk grasslands[END_REF][START_REF] Baumann | Diversity and functional composition of alpine grasslands along an elevational transect in the Gran Paradiso National Park (NW Italy)[END_REF], it is difficult to separate the role of pH per se from that of other environmental variables. Sites with steep and rocky or shallow soils have necessarily high pH values because of the surfacing limestone material (e.g. Palpurina et al. 2017); on the other hand, because of the widespread calcareous bedrocks, in our study area low pH values are usually confined to concave landforms (hollows, dolines, bases of slopes), in turn typically associated with deep soils, low rock cover and low slope inclination. Such land-forms can host acid soils, in spite of the base-rich bedrock, mainly because they accumulate clayey and decarbonated soil through colluvial deposition and karstic processes (e.g. Calandra 1999); in addition, they may trap huge amounts of wind-blown siliceous materials such as Saharan quartz silt and volcanic ashes (Giraudi 1995; Sauro 2004), and/or host a thick snow cover, favoring podzolization [START_REF] Hiller | Snowbeds on silicate rocks in the Upper Engadine (Central Alps, Switzerland)-pedogenesis and interactions among soil, vegetation, and snow cover[END_REF]. At the same time, dolines and hollows will receive more water and nutrients than convex land-forms or steep slopes (Araya et al. 2011;Batori et al. 2017), host soil types with higher clay fraction and higher available water content (Di [START_REF] Pietro | A geobotanical survey on acidophilous grasslands in the Abruzzo, Lazio and Molise National Park (Central Italy)[END_REF], and undergo lower levels of natural disturbance at micro-scale via summer drought, water run-off and erosion [START_REF] Bennie | Influence of slope and aspect on long-term vegetation change in British chalk grasslands[END_REF].

A c c e p t e d m a n u s c r i p t

A large number of studies found that grassland species richness in Europe and N-America was reduced on soils with high moisture and nutrients, while attained maximum values at sites with relatively dry, nutrientpoor and shallow soils (reviews by Grace 1999 and by Keddy 2005; see also [START_REF] Bennie | Influence of slope and aspect on long-term vegetation change in British chalk grasslands[END_REF][START_REF] Bonanomi | Short-term effects of nitrogen enrichment, litter removal and cutting on a Mediterranean grassland[END_REF][START_REF] Maccherini | Factors associated with species richness in a remnant calcareous grassland[END_REF]Moeslund et al. 2013b;[START_REF] Baumann | Diversity and functional composition of alpine grasslands along an elevational transect in the Gran Paradiso National Park (NW Italy)[END_REF]Filibeck et al. 2016;[START_REF] Kleinebecker | And the winner is….! A test of simple predictors of plant species richness in agricultural grasslands[END_REF]Palpurina et al. 2019). Even in a different floristic kingdom, [START_REF] Cantero | Species richness, alien species and plant traits in Central Argentine mountain grasslands[END_REF] found that in mountain grasslands in Argentina, species richness of 1-m 2 quadrats had a negative relationship with soil depth and available water, and a positive correlation with stone cover. This is consistent with the classical model of humped-back relationship between productivity and richness, usually explained assuming that in very productive sites a few strong competitors will create closed canopies and accumulate litter, reducing the amount of light reaching soil surface and outcompeting stress-tolerant species that have lower growth rates (e.g. Al-Mufti et al. 1977;[START_REF] Janssens | Relationship between soil chemical factors and grassland diversity[END_REF]Grytnes 2000;[START_REF] Stevens | Resource quantity, not resource heterogeneity, maintains plant diversity[END_REF]Pärtel et al. 2005;[START_REF] Kleinebecker | And the winner is….! A test of simple predictors of plant species richness in agricultural grasslands[END_REF]; but see Šímová et al. 2013 for a review and empirical testing of some alternative or concurrent explanations). Thus, we surmise that, in our study area, the positive effect by soil pH on richness may partly be explained with the different levels of stress controlled by the same topographical factors that influence soil pH (Moeslund et al. 2013a): dry and steep land-forms may promote a higher species density because they favor the coexistence of many small-sized species that colonize "gaps" and are not able to grow in more productive sites dominated by large, mat-forming and litter-accumulating species (Grubb 1976;[START_REF] Bennie | Influence of slope and aspect on long-term vegetation change in British chalk grasslands[END_REF][START_REF] Habel | European grassland ecosystems: threatened hotspots of biodiversity[END_REF][START_REF] Dengler | Biodiversity of Palaearctic grasslands: a synthesis[END_REF]. This hypothesis is supported by the PCA of the plot × species matrix (Online Resource 2), revealing that the poorer plots were characterized by perennial grasses known for their typical "competitive" characteristics (Grime 1973) such as Agrostis capillaris and Festuca microphylla. The former is a patch-forming species, capable of establishing a dense leaf canopy and assuming dominance due to its ability for rapid lateral spread (Grime et al. 2007); the latter belongs to the group of fescues with extra-vaginal branching strategy (Pignatti 1982; for a definition of this trait see e.g. [START_REF] Stace | The distinction between the Festuca ovina L. and Festuca rubra L. aggregates in the British Isles[END_REF], a character leading to the development of expanded tussocks and mats. Notably, both A. capillaris and F. microphylla are also well known for being bound to acidic to sub-acid soils throughout their geographic range (e.g. Pignatti 1982;

A c c e p t e d m a n u s c r i p t Grime et al. 2007) and, within the pedo-climatic context of the Apennines, to fine-earth deep soils (e.g. [START_REF] Bruno | Note illustrative alla carta della vegetazione del Parco Nazionale d'Abruzzo[END_REF]Pignatti 1982;[START_REF] Blasi | On the vegetation series in the northern sector of the Simbruini mountains (Central Apennines)[END_REF]. Thus, in our study area, the observed effect of pH on species richness might actually be caused by vegetation structure patterns led by these two grasses, unable to establish at high-pH and dry sites. Consistently, among the species significantly associated (in the PCA) with the poorer plots, there were some herbs characterized by creeping growth-forms (Potentilla rigoana, Veronica serpyllifolia), a strategy that allows to exploit the shaded micro-habitat under the grass canopy or to colonize the surface of the thick litter (Grime 1998). On the other hand, the perennial grass most significantly associated in the PCA with the high-richness plots was Festuca circummediterranea, a fescue that can only develop into isolated tufts with a reduced lateral spread, because of intra-vaginal shoot system (e.g. Pignatti 1982). Indeed, the other species correlated with the richer plots included either opportunistic annuals (e.g. Arenaria serpyllifolia, Cynosurus echinatus, Geranium molle, Medicago lupulina, Trifolium scabrum) that may fill the "gaps" between the perennial grasses, or typically stress-resistant, slow-growing plants (e.g. the CAM species Sedum acre).

Somewhat similar conclusions were drawn by a study on plant traits and assembly rules in another area of the Central Apennines [START_REF] Habel | European grassland ecosystems: threatened hotspots of biodiversity[END_REF][START_REF] Dengler | Biodiversity of Palaearctic grasslands: a synthesis[END_REF], showing that Sesleria nitida, a grass capable of forming high-density carpets, was the main factor suppressing species coexistence at extremely fine spatial scale (0.0025 m 2 ) -but only on N-facing slopes, because on the dryer S-facing slopes it built single, isolated tussocks leaving large spaces available to other species. However, since in the present study we didn't collect data on species cover or traits, further work is necessary to confirm our hypotheses.

Finally, an interesting problem concerns silt fraction: this predictor showed a positive effect on species richness, even at the 0.01 m 2 grain, and it was the only significant environmental variable to be uncorrelated with the pH gradient. The strong role of this predictor was quite unexpected and is not easy to explain. Only a small number of studies on grassland diversity analyzed soil texture, and even these often failed to find a correlation with species richness (e.g. [START_REF] Turtureanu | Scaleand taxon-dependent biodiversity patterns of dry grassland vegetation in Transylvania[END_REF]. The results by Stohlgreen et al. (1999) and [START_REF] Cousins | The influence of management history and habitat on plant species richness in a rural hemiboreal landscape, Sweden[END_REF], showing a positive effect of silty soils on species richness, can't be compared with our observations because of their plot sizes being much larger and their pedo-climatic context very different. Moreover, we acknowledge that the soil texture data used in our study were A c c e p t e d m a n u s c r i p t obtained without the prior destruction of soil organic matter (SOM) with H2O2, leading to possible overestimation of silt and underestimation of clay fraction because of silt-sized micro-aggregates made up of SOM-clay complexes [START_REF] Jensen | Soil texture analysis revisited: Removal of organic matter matters more than ever[END_REF]). However, silty soils are generally better drained than clayey soils, and extremely silt-rich soils may be scarcely fertile as they do not show any structural arrangement because of lack of the binding properties provided by clay (e.g. [START_REF] Thompson | Soils And Soil Fertility[END_REF]). In a phytosociological study performed in a neighboring area, silt fraction was found to be higher in the xerophytic, sparse grasslands of steep slopes, than in the mesic, dense communities of karst hollows (Di [START_REF] Pietro | A geobotanical survey on acidophilous grasslands in the Abruzzo, Lazio and Molise National Park (Central Italy)[END_REF].

Thus, the effect of silt might point once again to a negative correlation -within the study region -between species richness and biomass or competition, in turn controlled by topographic-hydrological gradients.

Diversity-environment relationship: spatial scales

We found that the relative importance of some predictors varied with plot size. This was already shown by previous studies dealing with multiple-scale sampling in grasslands (e.g. Auestad et al. 2008;Polyakova et al. 2016). It has been suggested that at extremely fine scales the availability of nutrients will influence the identity of the dominant species, in turn influencing richness; while at larger plot sizes, factors controlling richness via habitat heterogeneity will prevail (Auestad et al. 2008). Although we only had a small range of plot sizes, this may explain while stone and rock cover (a proxy of habitat heterogeneity) had a relevant effect only at 1 m 2 .

The overall explanatory power of the predictor variables became weaker with decreasing plot size. This is not unexpected, because very small spatial scales are unavoidably associated with the increasing effects of factors other than the physical variables, such as plant individual number and size [START_REF] Palmer | Artifacts and artifictions in biodiversity research[END_REF][START_REF] Storch | The theory of the nested species-area relationship: geometric foundations of biodiversity scaling[END_REF], stochastic disturbance and diaspore availability (Zobel et al. 2000), and fine-scale plant interactions [START_REF] Dengler | Biodiversity of Palaearctic grasslands: a synthesis[END_REF]. Moreover, environmental variables were measured at the scale of the largest plot, so they might not be representative of the micro-habitat variability between the subplots.

Diversity-environment relationship: z-values

A c c e p t e d m a n u s c r i p t

A high z-value (slope) of the species-area curve implies that there is a low mean number of individuals per species [START_REF] Storch | The theory of the nested species-area relationship: geometric foundations of biodiversity scaling[END_REF], i.e. each species occurs with scattered individuals leading to a high spatial heterogeneity of the vegetation. This can be in some cases a result of physical habitat heterogeneity. Indeed, we found that z-values were positively correlated with stone and rock cover and with "humps", and negatively correlated with soil depth: this can be explained as rocky and shallow soils give rise to a "rough" microtopography; while the flat, deep soils of the doline bottoms are rather homogeneous. Although stone and rock cover can also influence z-values via a trivial mechanism (i.e., on very rocky soils there is a higher chance that the subplots of the smallest size will fall on a completely bare surface, thus inflating the steepness of the species-area curve), in our study area rock-free soils are usually found only in concave land-forms, i.e. at sites with high pH and deep soils: hence, there is some overlapping between sites that curtail species richness (alpha-diversity) and sites featuring a low turnover of species between spatial scales (beta-diversity). Our results are thus consistent with [START_REF] Chiarucci | Effects of productivity on species-area curves in herbaceous vegetation: evidence from experimental and observational data[END_REF], who found a negative relationship between z-values and productivity in Tuscany and in Germany. Polyakova et al. (2016), in primary steppe grasslands in S-Siberia, found a strong positive effect on z-values of soil roughness at microscale. In contrast, both [START_REF] Turtureanu | Scaleand taxon-dependent biodiversity patterns of dry grassland vegetation in Transylvania[END_REF] (in secondary grasslands in Romania) and [START_REF] Baumann | Diversity and functional composition of alpine grasslands along an elevational transect in the Gran Paradiso National Park (NW Italy)[END_REF] (in high-elevation grasslands in the Italian Alps) found that their measured environmental predictors explained little of the variability in z-values. However, comparisons of z-values between different studies have to be interpreted very cautiously, as the results will depend on the size of the plots used [START_REF] Chiarucci | Effects of productivity on species-area curves in herbaceous vegetation: evidence from experimental and observational data[END_REF].

Conclusions

The importance ranking and effect size of the environmental variables were partly different across alphadiversity spatial scales or between alpha-and beta-diversity; nevertheless, the sites with acid and deep soils had on average lower levels of all the diversity measures, while maximum alpha-and beta-diversity were usually attained at sites with neutro-basic and shallow soils, often on steep or rocky land-forms. These relationships seem to be mediated by competition: at the low-diversity sites, species coexistence is probably suppressed by a small number of carpet-forming species, in turn bound to topographically-A c c e p t e d m a n u s c r i p t controlled "islands" of acidic or deep soils. Further studies will be necessary in the Apennine grasslands to confirm this explanation. However, the observed patterns point to the importance of the less productive and dryer sites for the conservation of secondary grassland plant diversity, as underlined by many previous works (e.g. Grime 1973; [START_REF] Moore | Conservation of wetlands: do infertile wetlands deserve a higher priority?[END_REF]Wassen et al. 2005;[START_REF] Bennie | Influence of slope and aspect on long-term vegetation change in British chalk grasslands[END_REF][START_REF] Cingolani | Predicting alpha, beta and gamma plant diversity from physiognomic and physical indicators as a tool for ecosystem monitoring[END_REF]Filibeck et al. 2016;[START_REF] Kleinebecker | And the winner is….! A test of simple predictors of plant species richness in agricultural grasslands[END_REF]Palpurina et al. 2019). In the Apennine mountains, these oligotrophic habitats seem to be not (yet) significantly threatened by atmospheric nitrogen deposition [START_REF] Bonanomi | Short-term effects of nitrogen enrichment, litter removal and cutting on a Mediterranean grassland[END_REF][START_REF] Nogueira | Extended autumn drought, but not nitrogen deposition, affects the diversity and productivity of a Mediterranean grassland[END_REF]; but see Ochoa-Hueso et al. 2017), a factor heavily impacting some C-and N-European grasslands (e.g. [START_REF] Stevens | Nitrogen deposition threatens species richness of grasslands across Europe[END_REF]). They may be nevertheless at risk because of grazing abandonment [START_REF] Bonanomi | Short-term effects of nitrogen enrichment, litter removal and cutting on a Mediterranean grassland[END_REF]): many sectors of the Apennines have witnessed the development of a closed woodland in a very short time [START_REF] Bracchetti | Land-cover changes in a remote area of central Apennines (Italy) and management directions[END_REF]Amici et al. 2013). Moreover, in many European and Mediterranean landscapes, shrub and tree colonization seems to be hindered in land-forms with deep clayey soils and to take place preferentially on shallow or rocky soils, probably because of reduced competition by grasses and larger availability of suitable microsites (Grove and Rackham 2001: 57;Rackham 2006: 201-202;Pärtel and Helm 2007). This might lead to a proportionally larger loss of those grassland types that host the highest levels of alpha-and beta-diversity.

  s c r i p t Finally, z-values were modelled as a function of land-form position, stone and rock cover, slope inclination, and soil depth. Model results (Tab. 6) showed a strong positive effect (stand. coeff.=0.63, p<0.001) of stone and rock cover. The model explained 47% of variability.

  . and L.C. conceived the study and the sampling methodology; L.C. coordinated field work; L.D.M., M.G.S. and G.F. participated in field sampling; L.C. revised critical vascular plants; F.R. supervised soil analysis; M.G.S., M.B. and G.F. performed data analysis; G.F. supervised the overall work and led the writing; all authors discussed the results and contributed to text editing.

AFIGURESFig. 1

 1 FIGURES

Fig. 2

 2 Fig. 2 Layout of the nested-plot series surveyed at each random sampling point.

Fig. 3

 3 Fig. 3 Regression Tree for 1 m 2 alpha-diversity as response variable. Variables responsible for tree splitting

Fig. 4

 4 Fig. 4 Regression Tree for 0. 1 m 2 alpha-diversity as response variable. See caption in Fig. 3 for explanations.Unlike the tree for 1 m 2 , here "univariate p-values" (instead of Bonferroni-corrected) was set as stopping

Fig. 5

 5 Fig. 5 Regression Tree for z-value as response variable. See caption in Fig. 3 for explanations

  Number of times, out of 10 Random Forest tries, that each predictor variable turned out to be important for each dependent variable, i.e. featured positive values exceeding the absolute value of the most negative score. Variables are listed in decreasing order of number of significant results A c c e p t e d m a n u s c r i p t

  Summary statistics of the environmental variables measured at each nested-plot series and used as predictor variables in the statistical analysis Summary statistics of the vascular plant richness encountered in all plots at different scales, and of the derived biodiversity metrics used as dependent variables in the statistical analysis Pair-wise correlation table between all predictor and dependent variables. Below the diagonal: Spearman's correlation coefficients (rs). Above the diagonal (in italic): p-values. Note that these are not corrected for multiple testing. Correlation coefficients corresponding to a p-value <0.05 are in bold

	TABLES
			Species richness values
	Area (m 2 )	n	Mean ±SD	Median Min	Max
	Continuous variables Elevation (m) Slope inclination (°) Heat-load index Soil depth (median) (cm) Stone and rock cover (%) Soil pH Silt (%) Clay (%) Sand (%) Soil organic matter (%) Categorical variables Land-form position (number of cases in each category) Tab. 1 A c c e p t e d Mean 1523 10.3 0.9 20.7 10.3 6.5 20.2 14 65.9 12.08 Hollow (16), Hillside (19), Hump (15) SD Min. Max. 138 1331 1767 8.5 0 35 0.1 0.57 1.04 13.4 4 60 11.9 0 40 0.7 4.9 7.6 12.1 0 68.5 7.7 0.06 41 16.5 3.3 89.8 5.4 0.21 21.5 m a n u s c r i p t 1 50 28.0 ±6.38 27 17 41 0.1 100 18.46 ±5.03 17 9 31 0.01 200 10.26±3.48 10 0 20 Diversity metrics Name n Mean ±SD Median Min Max 1 m 2 alpha-diversity 50 28.0 ±6.38 27 17 Elevation Slope inclination Heat-Load index Stone and rock cover Soil depth Clay Silt Sand Soil organic matter Soil pH 1 m 2 alpha-diversity 0.1 m 2 alpha-diversity 0.01 m 2 alpha-diversity 41 0.1 m 2 average alpha-diversity 50 18.46 ±4.79 17 10.5 29.5 0.01 m 2 average alpha-diversity 50 10.26±2.87 10 4.5 16.25 z-value 50 0.238± 0.063 0.224 0.150 0.443 m Elevation 0.002 0.0406 0.0001 0.3229 0.5566 0.8269 0.7875 0.1089 <0.0001 0.1496 0.1043 0.6268 0.4198 z-values Slope inclination 0.425 0.0064 <0.0001 0.0838 0.1252 0.0553 0.0424 0.0700 0.0004 0.0987 0.0600 0.7375 0.0231 Heat-Load index -0.291 -0.380 0.0888 0.7618 0.2853 0.4740 0.4275 0.1519 0.4418 0.3290 0.5415 0.9119 0.1512 Stone and rock cover 0.537 0.562 -0.243 0.0007 0.0375 0.8169 0.4132 0.0325 <0.0001 0.0167 0.0947 0.5773 <0.0001 Soil depth -0.143 -0.247 0.044 -0.464 0.8117 0.6150 0.6368 0.3589 <0.0001 0.1114 0.7386 0.3567 0.0005 Clay -0.085 -0.220 0.154 -0.295 -0.035 0.4055 0.0001 0.0849 0.4342 0.2178 0.3512 0.4836 0.1790 Silt 0.032 -0.273 0.104 -0.034 0.073 0.120 <0.0001 0.0310 0.5116 0.0329 0.0729 0.0775 0.9096 a n u s c r i p t Sand 0.039 0.288 -0.115 0.118 -0.068 -0.530 -0.868 0.0207 0.1899 0.2058 0.2712 0.2725 0.6742 Soil organic matter 0.229 0.258 -0.206 0.303 -0.133 -0.246 -0.305 0.326 0.0016 0.3513 0.2896 0.5206 0.6732 Soil pH 0.590 0.481 -0.111 0.642 -0.611 -0.113 -0.095 0.189 0.435 0.0063 0.0309 0.4272 0.0314 1 m 2 alpha-diversity 0.207 0.236 -0.141 0.337 -0.228 -0.177 0.302 -0.182 0.135 0.381 <0.0001 <0.0001 0.1323 0.1 m 2 alpha-diversity 0.232 0.268 -0.088 0.239 -0.048 -0.135 0.256 -0.159 0.153 0.306 0.877 <0.0001 0.6193 0.01 m 2 alpha-diversity 0.070 0.049 0.016 -0.081 0.133 -0.101 0.252 -0.158 0.093 0.115 0.702 0.844 0.0002 z-values 0.117 0.321 -0.206 0.579 -0.477 -0.193 -0.016 0.061 0.061 0.305 0.216 -0.072 -0.502 Tab.2 A c c e p t e d Tab.4

Acknowledgements

We are grateful to Abruzzo Lazio & Molise National Park agency for having authorized field work (special thanks to C. Gentile and C. Sulli for scientific advice and field support), G. Ubertini for laboratory analysis, S. Soliveres and M. Brandi for advice on recursive partitioning, S. Maccherini and M. Marignani for useful