Local topographic and edaphic factors largely predict shrub encroachment in Mediterranean drylands

Local topographic and edaphic factors largely predict shrub encroachment in Mediterranean drylands Alice Nunes a ,MelanieKöbel a , Pedro Pinho a,b , Paula Matos a , Edoardo A.C. Costantini c , Cristina Soares a , Francesco de Bello d ,OtíliaCorreia a , Cristina Branquinho a, ⁎ a Centre for Ecology, Evolution and Environmental Changes, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, C2, Piso 5, 1749-016 Lisboa, Portugal b CERENA - Centro de Recursos Naturais e Ambiente, Universidade de Lisboa, Instituto Superior Técnico, Av. Rovisco Pais, 1049-001 Lisboa, Portugal c CREA-AA, Research Centre for Agriculture and Environment, Florence, Italy d Institute of Botany, Czech Academy of Sciences, CZ-37982 Trebon, Czech Republic HIGHLIGHTS •Shrub encroachment is a worldwide phenomenon that affects many ecosystem services. •Itscausesneedtobebetterunderstood and quantified in the Mediterranean Basin. •We studied its drivers, and impact on plant functional traits in Holm-oak woodlands. •We found that topo-edaphic factors, not climate, largely predict shrub encroachment. •We help to predict encroachment under climate change and to improve management actions. GRAPHICAL ABSTRACT abstractarticle info Article history: Received 15 July 2018 Received in revised form 29 November 2018 Accepted 30 November 2018 Available online 1 December 2018 Shrub encroachment influences several ecosystem services in drylands worldwide. Yet, commonly used strategies to reduce encroachment show a low medium-term success, calling for a better understanding of its causes. Previous works identified multiple drivers responsible for this phenomenon, including anthropogenic and environmental causes. However, the relative effect of climate, topography and edaphic factors on shrub encroachment is not fully understood nor has been properly quantified in Mediterranean Basin drylands. Also, understanding how these drivers lead to changes in plant communities' functional traits associated to shrub encroachment is crucial, considering traits influence ecosystem processes and associated ecosystem services. Here, we studied the understory of a Mediterranean dryland ecosystem composed of savanna-like Holm-oak woodlands, along a regional climatic gradient. We specifically assessed (i) how climatic, topographic and edaphic factors influence understory relative shrub cover (RSC) and (ii) their direct and indirect effects (via RSC) on plant functional traits. We studied the mean and diversity of 12 functional traits related to plant regeneration, establishment, and dispersal, at the community-level. We found that, under similar low-intensity land use, topographic and edaphic factors, namely slope variations and soil C:N ratio, were the most important predictors of shrub encroachment, determining communities' functional characteristics. Climate, namely summer precipitation, had a much lesser influence. Our model explained 52% of the variation in relative shrub cover. Climate had a stronger effect on a set of functional traits weakly involved in shrub encroachment, related to flowering Keywords: Aridity Climate change Holm-oak Functional traits Structural equation model Sustainable land management Science of the Total Environment 657 (2019) 310–318 ⁎Corresponding author. E-mail addresses: [email protected] (A. Nunes), [email protected] (M. Köbel), [email protected] (P. Pinho), [email protected] (P. Matos), [email protected] (E.A.C. Costantini), [email protected] (C. Soares), [email protected] (F. de Bello), [email protected] (O. Correia), [email protected] (C. Branquinho). https://doi.org/10.1016/j.scitotenv.2018.11.475 0048-9697/© 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Science of the Total Environment journal homepage: www.elsevier.com/locate/scitotenv and dispersal strategies. We show that shrub encroachment is largely predicted by topo-edaphic factors in Mediterranean drylands subject to conventional low-intensity land use. Hence, management strategies to reduce encroachment need to take these drivers into account for efficient forecasting and higher cost-effectiveness. Our results suggest that climate change might not greatly impact shrub encroachment in the Mediterranean Basin, but may affect functional structure and reduce functional diversity of plant communities, thus affecting ecosystem functioning. © 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http:// creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction Shrub encroachment, defined as the increase in density, cover and biomass of native woody species (Van Auken, 2009), has been intensively debated in recent literature (Eldridge et al., 2011;Sala and Maestre, 2014;Eldridge and Soliveres, 2015;Maestre et al., 2016). It has been reported mostly for dryland ecosystems, e.g. in North America (Van Auken, 2009), Australia (Eldridge and Soliveres, 2015), South Africa (Roques et al., 2001) and in the Mediterranean Basin (Castro and Freitas, 2009;Maestre et al., 2009;Caldeira et al., 2015). Shrub encroachment results from both natural and human-induced factors that are difficult to disentangle, leading to apparently contradictory results. Long-term overgrazing has been often considered a cause of shrub encroachment in grasslands in western United States or in southern Africa drylands (Schlesinger et al., 1990;Roques et al., 2001;MEA, 2005;Van Auken, 2009). Contrastingly, in the Mediterranean Basin, shrub encroachment is seen as a consequence of grazing exclusion (Castro and Freitas, 2009) or of the abandonment of the land, ceasing active shrub removal (Maestre and Cortina, 2005), i.e. as an expected outcome of a secondary succession. It has been also associated to past shrub clearing, reduced fire frequency, increases in CO 2 , N deposition and long-term climate change (Archer, 2010). At a global scale, mean annual precipitation sets an upper limit for woody cover (shrubs or trees) (Sankaran et al., 2005;Hirota et al., 2011). However, locally, topography and soil characteristics limit this potential (Colgan et al., 2012;Sala and Maestre, 2014), particularly in water-limited ecosystems (Gómez-Plaza et al., 2001), because they largely determine water availability to plants. Finally, biotic-abiotic interactions may also influence shrub cover through positive feed-backs. Shrub colonization may lead to a more heterogeneous distribution of soil water and nutrients, thus promoting further shrub encroachment (Schlesinger et al., 1990). In short, shrub encroachment depends on multiple factors, varying with the historical management (e.g. grazing pressure), environmental conditions (e.g. climate, soil properties, topography) and the scale of analysis (e.g. landscape or patch level) (Eldridge et al., 2011;Eldridge and Soliveres, 2015). Shrub encroachment has been described as a land degradation process (Schlesinger et al., 1990;MEA, 2005), mostly from the point of view of livestock production, a primary use of drylands. However, it may improve carbon sequestration (Daryanto et al., 2013), soil fertility and N mineralization rate (Maestre et al., 2009;Gómez-Rey et al., 2013), having a beneficial effect on ecosystem functioning at least until a certain threshold of relative woody cover (41–60%) is reached (Soliveres et al., 2014;López-Díaz et al., 2015). Several previous works support the idea that the consequences of shrub encroachment are contextdependent, varying e.g. with precipitation (Jackson et al., 2002;Knapp et al., 2008;Anadón et al., 2014). Recent works suggest that the effect of shrub encroachment on ecosystem functioning largely depends on the functional traits of the woody and herbaceous species involved (Rivest et al. 2011;Rolo et al., 2012;Soliveres et al., 2014;Maestre et al., 2016), which largely influence ecosystem processes (Díaz et al., 2007;de Bello et al., 2010). Shrub encroachment should involve changes in community functional traits that co-vary among shrub and herbaceous growth-forms, such as plant height, life-cycle, root depth, and specificleafarea(Valencia et al., 2015), reflecting trade-offs in resource allocation patterns (Díaz et al., 2016). Previous studies addressed the effects of multiple factors on shrub encroachment, mainly from outside of the Mediterranean Basin (Roques et al., 2001;D'odorico et al., 2012;Iannone et al., 2015). Some of them predict its increase under a climate change scenario of increased aridity (D'odorico et al., 2012), like the one forecasted for Mediterranean Basin drylands (Huang et al., 2016). However, to our knowledge, none has quantified jointly the relative effect of climate, topography and edaphic factors on shrub encroachment and associated changes in community functional traits in Mediterranean Basin drylands. We propose to overcome this knowledge gap by studying the plant community of savanna-like Holm-oak woodlands along a spatial climatic gradient, comprising sites with varying topo-edaphic characteristics. Hence, we adopted a ‘space-for-time substitution’, assessing changes along space to infer changes over time (Blois et al., 2013). Holm-oak woodlands consist of scattered trees with an understory of semi-natural grasslands and shrubland patches, and have a traditional low intensity silvo-pastoral use. These heterogeneous savanna-like ecosystems are called montados in Portugal or dehesas in Spain, and support a remarkably rich biodiversity, being very important for the economy of rural areas (Pereira and Da Fonseca, 2003). Specifically, we assessed (i) how climatic, topographic and edaphic factors influence relative shrub cover and (ii) their direct and indirect effects (via relative shrub cover) on plant functional traits. We studied the mean and diversity of 12 functional traits related to plant establishment, persistence, regeneration, and dispersal, at the community level. To test direct and indirect effects of environmental drivers on the plant community we built an a priori causal model based on our results and on observationsof previous studies, using structural equation modelling (SEM model) (Grace et al., 2012). Under similar low-intensity land use, climate and topo-edaphic variables are expected to exert a strong control on relative shrub cover (Sala and Maestre, 2014) at a regional and local scale, respectively. Our first hypothesis is that local topo-edaphic factors are as important as climate, as drivers of shrub encroachment. Grass and shrubdominated communities are expected to differ considerably in functional traits' means and diversity. Because functional traits are linked to ecosystem processes (e.g. primary productivity, decomposition) (de Bello et al., 2010), this information can provide clues about the consequences of shrub encroachment for ecosystem functioning. In addition, it is important to disentangle which functional consequences are mostly driven by local environmental factors, from those due to regional filters, i.e. climatic variables. Our second hypothesis is that the community response to local environmental filters via changes in relative shrub cover (i.e. indirectly) will affect a different set of functional traits than that due to regional climatic filters. 2. Methods 2.1. Study sites and environmental variables The study was developed in Mediterranean drylands dominated by Holm-oak woodlands (Quercus ilex L. subsp. rotundifolia) in southwestern Iberian Peninsula. Field sampling was carried out in 54 sites selected along a regional climatic gradient. A stratified random selection of sites along a macroclimatic gradient was made (Fig. S1), based on the United Nations aridity index (AI) (Middleton and Thomas, 1992). The AI represents the ratio of mean annual precipitation to annual potential 311A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318 evapotranspiration, and the lower the index value, the higher the aridity. The global aridity database (http://www.cgiar-csi.org/data/globalaridity-and-pet-database)(Trabucco and Zomer, 2009) was used to retrieve AI data for the period 1950–2000. Along the study area AI ranged from 0.42 to 0.56. The selection of the sampling sites was made in order to ensure that they had no agricultural activities over the last five years, and moderate to low grazing intensity (see Nunes et al., 2017 for more details about sampling site selection). We sampled the understory of Holm-oak woodlands with a sparse tree cover (b40 trees/ha, on average) (Amaral et al., 1997). Shrub cover ranged from 0% to 87% (average 16%) over all sampling sites. The most common species in shrub patches were Cistus ladanifer,Cistus salvifolius and Lavandula stoechas,and in grassland areas were Brachypodium distachyon,Vulpia geniculata, Gaudinia fragilis,Leontodon taraxacoides,Tolpis barbata,and Chamaemelum mixtum. Sampling sites were characterized by a set of 19 climate variables extracted from Worldclim database with a 1 km 2 resolution (Hijmans et al., 2005) (Table S1). Also, three topographic variables were computed for each site based on digital elevation models with 10 m resolution using ArcGIS 10.1 (ESRI, 2010): local slope values, the standard deviation of the slope within a 250 m buffer around the sampling site centroid, and a topographic wetness index (TWI) (Sørensen et al., 2005), to quantify topographic control on hydrological processes (Table S1). A higher TWI value indicates a higher potential for runoff. Soil samples were collected at each site from the upper 10 cm (composite sample of 5 subsamples). Soil texture (% of sand, clay and silt) and nutrient content (soil organic matter, soil N, and C:N ratio) were determined using standard procedures (Table S1). After exploring the correlation among the climatic, topographic and edaphic descriptors at each site (Tables S2 and S3), we latter summarized them into 8 variables showing a Pearson's r correlation b0.70, to be used in subsequent analysis (Table 1). 2.2. Vegetation sampling and community trait metrics Understory plant sampling was performed in April–June 2012, at the peak of standing biomass. Plant cover was assessed using the point– intercept method, avoiding drainage and flooding areas. Species cover was estimated along six 20 m transects established at each site, in points spaced every 50 cm (further detail on the sampling scheme may be found in Nunes et al., 2014). The year in which the sampling took place had a low average annual precipitation (370 ± 61 mm) when compared to the 50-year mean of the sampling sites altogether (561 ± 27 mm). We used the dominant species attaining ≥80% of the relative cover to characterize communities' functional traits (Pakeman and Quested, 2007) (95 species overall, and 30 ± 9 per site). Plant trait selection was based on plant strategies related to establishment, persistence, regeneration, and dispersal, totaling 12 functional traits (Table S4). Trait data was obtained through direct measurements in the field following standard protocols (Pérez-Harguindeguy et al., 2013), or retrieved from various bibliographic sources (Table S4). We calculated functional structure and diversity through the communityweighted-mean (CWM) (Garnier et al., 2007) and functional dispersion (FDis) (Laliberte and Legendre, 2010), respectively. The CWM is calculated as the average trait value in a community weighted by the relative abundance of the species carrying each value (Garnier et al., 2007). FDis corresponds to the weighted mean distance of individual species from the weighted centroid of all species in a multidimensional trait space, where weights are species relative abundances (Laliberte and Legendre, 2010). It reflects the degree of functional dissimilarity within the community (Laliberte and Legendre, 2010), and a higher FDis is expected to lead to increased complementarity in resource use between species and thus to increased ecosystem functioning (Tilman et al., 1997). Mean trait values per species were used to compute functional metrics. Prior to analysis, continuous traits were log transformed. To deal with continuous, ordinal, and categorical variables, as well as with missing values, the Gower distance was used in calculations. All calculations were done using the dbFD function of the FD package (Laliberté et al., 2014)inR(RCoreTeam,2015). 2.3. Statistical analyses We first analysed relative shrub cover through linear regression using the selected climatic (aridity index, summer precipitation, winter temperature), topographic (slope standard deviation, topographic wetness index) and edaphic (soil organic matter, C:N ratio, and sand content) predictors (Table 1), as well as two-way interactions among them. We also included a quadratic term for aridity, to account for nonlinear responses (Gross et al., 2013). Correlation among the predictors was below Pearson's r 0.7 in all cases (Tables S2 and S3). We used a multi-model inference approach and the best model was selected based on the Akaike information criterion (AIC), after ensuring no significant multicollinearity through variance inflation factors (VIF). Afterwards, we built regression models for the functional trait metrics using the same predictors and procedure. As we found strong correlations between some trait functional metrics (Tables S5 and S6), we conducted two principal component analyses (PCA), based on correlation matrices, using separately CWM and FD values of all traits measured. This allowed us to summarize trait variation, and identify the most informative axes of functional specialization (Díaz et al., 2016). We then used the two first PCA components of each analysis as a measure of the CWM and the FD of each community, i.e., as response variables (see Valencia et al., 2015 for a similar approach). Then, we used a variance decomposition analysis based on the best models selected to obtain the percentage of variance explained by climatic and topo-edaphic predictors and their interactions (Dubuis et al., 2013). Finally, to identify the mechanisms that control shrub encroachment, and to test direct and indirect effects (i.e., via changes in relative shrub cover) of climatic and topoedaphic factors on the functional traits of the plant community, we used structural equation modelling (SEM). This approach is well suited for studying hypotheses about processes with complex causal connections (Grace et al., 2012). It allows separating the direct and indirect effects of the predictors included in a model and estimating the strengths of multiple effects. We established an a priori model based on our current knowledge and on the best regression models previously built for each response variable. We hypothesized that (i) climatic and topoedaphic variables would directly affect relative shrub cover and the functional characteristics of the plant community and (ii) changes in Table 1 Environmental variables selected, their units and range in the study area. Environmental variables Description Range and units Climatic Aridity index Ratio of mean annual precipitation to annual potential evapotranspiration, for the period 1950–2000 0.42–0.56 (unitless) Summer precipitation Precipitation of driest quarter 17–34 mm Winter temperature Mean temperature of the coldest quarter 9.4–11.3 °C Topographic Slope standard deviation (slope SD) Standard deviation of the slope within a 250 m buffer around the sampling site (based on digital elevation models, 10 m resolution) 1.16–14.69° Topographic wetness index (TWI) Ln [(upslope area)/local slope], used to quantify topographic control on hydrological processes 6.32–14.44 (unitless) Edaphic SOM Soil organic matter content 2.3–15.6% C:N Soil carbon/nitrogen ratio 6.8–18.1 (unitless) Sand Soil sand content 17.9–71.4% 312 A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318 functional metrics would be inter-related and largely mediated by changes in relative shrub cover. Although we admit that other a priori model structures could be considered, and that assuming cause-effect relationships from observational studies has its pitfalls, the metrics of goodness-of-fit of our model, and ecological evidences drawn from ours and other studies, clearly support the plausibility of the relations we propose among variables. Before modelling, we examined the distributions of all our endogenous variables, and tested their normality. The first component of the PCA on CWM values was normalized (z-score) before analyses. To test the overall fit of SEM models we used the chisquare statistic and its significance (the model is rejected if p-value b 0.05), the comparative fit index (CFI N0.9), the root mean square error of approximation index (RMSEA b0.10) and the standardized root mean square residual (SRMR b0.10). Standardized path coefficients estimated by maximum likelihood were used to measure the direct and indirect effects of the predictors (Grace et al., 2012). These coefficients are interpreted as the size of an effect that one variable exerts upon another. SEM models were fit using the Lavaan package (Rosseel, 2012). All the analyses were performed under R statistical environment (R Core Team, 2015). All the predictors used in modelling were standardized and normalized (z-score) before analyses. 3. Results 3.1. Relative shrub cover Variation in relative shrub cover (RSC) was mostly explained by topo-edaphic variables (64% of the adj. r 2 = 0.52), increasing with slope standard deviation and soil C:N ratio, which showed additive effects (Table 2,Fig. 1).Loweraridity(highervaluesofthearidity index) and lower summer precipitation also contributed to higher RSC, although climatic predictors explained a small proportion of its variation along the climatic gradient either per se (15%), or together with topo-edaphic factors (21%) (Table 2). 3.2. Functional trait metrics The PCA of the CWMs of all traits segregated two main PCA components, which accounted for 57% of the total variance found in the data (Fig. 2A). The first component, hereafter called ‘CWM-PCA1’,explained 46% of the variance and described a gradient of RSC, being negatively correlated with it. Communities with lower RSC showed a higher proportion of annuals (graminoids and rosettes), more anemochory dispersal and higher SLA; communities with higher RSC showed higher plant height and maximum root depth, more perennial species and more barochory dispersal. This was shown by the strong correlations found between CWM values of these traits and the CWM-PCA1 (Fig. 2A, Table 3). The second component, hereafter called ‘CWMPCA2’, explained 10% of the variance. The traits more strongly correlated with CWM-PCA2 were flowering traits (onset and duration), growthform (erect and prostrate), dispersal strategy (ectozoochory) and Nfixing ability (Fig. 2A, Table 3). Community FDis values were explained by the two main PCA components, accounting for 62% of the total variance found in the data (Fig. 2B). The first component, hereafter ‘FDisPCA1’, accounted for 42% of the variance and separated communities based on life-cycle, root depth, height and SLA traits, which showed the strongest correlations with this axis (Fig. 2B, Table 3). The second component, hereafter ‘FDis-PCA2’, explained 20% of the variance and segregated communities according to flowering traits (onset and duration), seed mass and growth-form, all positively correlated with this axis (Fig. 2B, Table 3). The first components of both PCAs, i.e. CWMPCA1 and FDis-PCA1, were mostly explained by topo-edaphic factors (64% and 76% of the total variance explained by the models, respectively), while the second PCA components, i.e., CWM-PCA2 and FDisPCA2 were best predicted by climatic variables (100% and 61% of the total variance explained by the models, respectively) (Table 2). 3.3. Structural equation modelling (SEM) Our a priory SEMs were well fitted to our data, as indicated by goodness-of-fit statistics (Fig. 3). They explained 54% of the variation Table 2 Summary of the best models selected for each response variable (see Tables S7 and S8 for further detail). The sign of the coefficients of the selected predictors in each model is indicated. The proportion of variance explained by topo-edaphic and climatic predictors (and by both) was calculated using a variance decomposition analysis based on each model. Response variables Predictors Proportion of the explained variance (%) Topo-edaphic Climatic Topo-edaphic × climatic Slope SD Soil C/N Slope SD × soil C/N Sand Soil C/N × Sand Aridity Summer precip. Winter temp. Slope SD × summer precip. Ad. R 2 Topo-edaphic Both Climatic Relative shrub cover (+) (+) (−)(−) 0.52 64 21 15 CWM PCA 1 (−)(−) (+) (+) 0.55 64 14 22 CWM PCA 2 (−)(−) 0.20 0 0 100 FDis PCA 1 (+) (+) (−) (+) 0.42 76 17 7 FDis PCA 2 (−) (+) (−)(−) (+) 0.36 39 0 61 Aridity represents –the value of the aridity index, so that higher aridity values correspond to drier conditions, to facilitate results interpretation. The highest values of the proportion of the explained variance by different predictors for each response variable are highlighted in bold. Fig. 1. Relative shrub cover variation with slope standard deviation (slope SD) and soil C:N ratio. Planes represent the predicted values of a linear regression fitted to both variables. The colors of the predicted planes change from blue (low values of relative shrub cover) to red (high values of relative shrub cover). 313A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318 in RSC, which was mainly driven by topo-edaphic factors. Slope standard deviation and soil C:N ratio had a strong direct positive effect on RSC (Fig. 3). On the contrary, summer precipitation had a moderate negative direct effect on RSC. Changes in RSC largely determined CWMPCA1 variation, explaining 98% of its variance through a direct negative effect (Fig. 3), showing that the effects of topo-edaphic variables and of summer precipitation on CWM-PCA1 (Fig. S2) were mediated by changes in RSC. The CWM-PCA2 was related only with climatic factors, although they were only able to explain 21% of its variation; summer precipitation and winter temperature had a direct strong and moderate negative effect on CWM-PCA2, respectively (Fig. 3). Topo-edaphic factors were negatively related to FDis-PCA1, explaining 64% of its variation; about 3% of this effect was direct, through the interaction between slope standard deviation and soil C:N ratio (although individually these factors had a marginal positive effect), and 97% indirectly, via CWM-PCA1 (Fig. 3). FDis-PCA2 was related to climatic factors and also to CWM-PCA1, which jointly explained 38% of its variation; about 87% of this effects was direct, driven by a marginal positive influence of summer precipitation and by a negative effect of aridity, and 13% was due to an indirect positive effect of CWM-PCA1 (Fig. 3). 4. Discussion The causes of shrub encroachment in Mediterranean Basin drylands need to be better understood and quantified in order to improve encroachment forecasting and the success of removal programmes. Here, Relave shrub cover A Relave shrub cover B Fig. 2. Principal component analysis (PCA) of A) community-weighted mean (CWM) and B) functional dispersion (FDis) trait values. Vectors represent traits described in Table S4. Abbreviations: Ther = therophyte; Phan = phanerophyte; Cham = chamaephyte; Gram = graminoid; Roset = rosette; Erect = erect; SLA = specific leaf area; Anem = anemochory; Baro = barochory; Height = height; Rdepht = root depth; Smass = seed mass; Slong = seed persistence; Dispers = dispersal strategy; Lcycle = life cycle; Bflow = onset of flowering; Dflow = duration of flowering. Study sites are represented by points in a gray scale indicating the class of relative shrub cover (%) at each site (see legend). See Supporting information Table S5 for correlations among CWM values and Table S6 for correlations among FDis values. 314 A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318 we show that, under similar low-intensity land use, shrub encroachment is largely predicted by local topographic and edaphic factors, and influenced to a much lesser extent by climate in the studied communities. Hence, integrating topo-edaphic factors into management strategies to predict and revert shrub encroachment may contribute to improve their cost-effectiveness and sustainability in the long run. Also, according to our results, climate change predictions of increasing aridity for the Mediterranean Basin are not expected to greatly impact shrub encroachment, but may affect particular functional traits of the plant community, thereby affecting ecosystem functioning. 4.1. Main drivers of shrub encroachment and implications for management Our findings partially match our first hypothesis, as both topoedaphic and climatic variables influenced shrub encroachment, jointly accounting for 52% of the variation in relative shrub cover. Topoedaphic factors explained most of the variation in relative shrub cover (64% of the r 2 ), which increased towards areas with higher slope standard deviation, and areas with higher soil C:N ratio. In drylands, water limitations constrain vegetation cover, favoring runoff and erosion in slopes. Hence, topography largely controls water and nutrient flow paths from upslope to lowlands, where they tend to accumulate (Gómez-Plaza et al., 2001). This leads to coarser-textured soils on hillslopes, with less ability to retain water at upper soil depths, favoring deep-rooted species such as shrubs. Conversely, plain areas tend to have finer-textured soils, and shallower water profiles, favoring shallowrooted species such as most herbaceous species (Sala et al., 1997; McAuliffe, 2003). Soil C:N ratio, which showed no significant correlation with slope standard deviation, neither with aridity (Table S3), was also an important predictor of shrub encroachment. The soil C:N ratio reflects the dynamic interaction between slow-modifying soil features such as particle size distribution and mineralogy, and local site characteristics like hydrology, and vegetation. A higher soil C:N ratio in shrub-encroached areas may indicate a lower decomposition rate, probably due to lower water availability. Yet, this assumption would need other direct indicators of soil organic carbon turnover (e.g. microbial biomass) to be confirmed. Additionally, a high soil C:N in shrub-dominated areas may also be a result of higher inputs of above and belowground litter from shrubs, which have a usually higher C:N ratio and lignin content, that take longer to decompose, than that from herbaceous species (Bot and Benites, 2005). These conditions may create a positive feed-back loop promoting further shrub encroachment. One of the strategies most used to deal with encroachment is mechanical shrub removal, along with the use of herbicides or prescribed burning (Archer, 2010). Yet, these techniques often fail its goals in the long run (Rango et al., 2005;Archer and Predick, 2014), calling for a better understanding of the causes of shrub encroachment under different contexts. Considering the importance of topographic and edaphic factors as predictors of shrub encroachment in Mediterranean Basin drylands, we suggest that combining shrub clearing with other locally applied techniques that take into account topo-edaphic factors, involving the manipulation of water and nutrient flow pathways and accumulation, might contribute to their longer-term effectiveness. There are Table 3 Pearson correlation coefficients between community-weighted-means (CWM) and the two first components of the respective principal component analysis (CWM-PCA1 and CWM-PCA2); and between functional dispersion (FDis) and the two first components of the respective principal component analysis (FDis-PCA1 and FDis-PCA2). Traits Category CWM FDis PCA1 PCA2 PCA1 PCA2 Life cycle −0.99⁎⁎⁎ −0.02 0.95⁎⁎⁎ 0.04 Life-form Therophyte 0.96⁎⁎⁎ 0.10 0.78⁎⁎⁎ 0.16 Hemicryptoph 0.01 −0.35⁎⁎ Geophyte 0.04 0.05 Chamaephyte −0.34⁎0.34⁎ Phanerophyte −0.98⁎⁎⁎ −0.04 Growth-form Bulb 0.04 0.05 0.16 0.68⁎⁎⁎ Erect 0.44⁎⁎ −0.37⁎⁎ Graminoid 0.73⁎⁎⁎ 0.24 Prostrate 0.26 −0.74⁎⁎⁎ Rosette 0.62⁎⁎⁎ 0.19 Shrub −0.99⁎⁎⁎ 0.03 Max. height −0.94⁎⁎⁎ 0.10 0.87⁎⁎⁎ −0.05 SLA 0.94⁎⁎⁎ 0.06 0.84⁎⁎⁎ −0.28⁎ Onset flower. −0.37⁎⁎ 0.67⁎⁎⁎ −0.20 0.79⁎⁎⁎ Duration flower. 0.08 −0.52⁎⁎⁎ −0.08 0.81⁎⁎⁎ Dispers. strategy Anemochory 0.95⁎⁎⁎ 0.23 0.79⁎⁎⁎ 0.23 Barochory −0.95⁎⁎⁎ −0.15 Ectozoochory 0.16 −0.57⁎⁎⁎ Endozoochory −0.15 −0.08 Seed mass −0.24 0.04 −0.26 0.65⁎⁎⁎ Seed persistence −0.79⁎⁎⁎ 0.07 0.73⁎⁎⁎ 0.08 N-fixing ability −0.32⁎ −0.48⁎⁎⁎ 0.13 0.04 Max. root depth −0.95⁎⁎⁎ 0.00 0.89⁎⁎⁎ 0.17 Significant correlations are highlighted in bold. ⁎pb0.05. ⁎⁎ pb0.01. ⁎⁎⁎ pb0.001. Fig. 3. Structural equation model to explain relative shrub cover and the main axis of functional specialization of the plant community regarding functional structure (CWM-PCA1 and CWM-PCA2) and functional dispersion (FDis-PCA1 and FDis-PCA2). Overall goodness-of-fit statistics: χ 2(23) = 26.668, p = 0.270, root mean square error of approximation (RMSEA) = 0.054 (0.00–0.121), comparative fit index (CFI) = 0.989, standardized root mean square residual (SRMR) = 0.056. Arrow widths are proportional to the standardized path coefficients, which are presented. The R 2 next to response variables indicates the proportion of variance explained. *p b0.05; **p b0.01; ***p b0.001. Aridity = -Aridity index. 315A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318 many of such techniques (e.g. cross-slope barriers like terraces, stone bunds, water-harvesting structures, manure spreading) proposed for sustainable land management in drylands which show some evidence of success in increasing water and nutrient availability to plants (Schwilch et al., 2012), although its effectiveness in dealing with shrub encroachment remains to be tested. A sustainable land management in drylands, and improved strategies to promote the recovery of degraded land e.g. through soil and water conservation, is among the aims of the UN sustainable development goals, and our findings contribute to approach these goals. Topo-edaphic and climatic variables jointly explained 52% of the variation in relative shrub cover, suggesting that it may be influenced by other factors, such as biotic interactions or past land management. In Holm-oak woodlands, plain areas may be more grazed by livestock (particularly cows) than slopes, because they have in general more productive soils, providing higher quantity and quality of forage (Bailey et al., 1996). In addition, occasional mechanical shrub clearing that may have occurred in the past, is easier to perform in flatter areas than in slopes. Hence, although a similar low-intensity land-use was one of the criteria used to select sampling sites, these factors may have partially contributed to a higher shrub cover found in areas with higher slope variation. Climate was a weak predictor of shrub encroachment. In light of this result, climate change forecasts of increased aridity for Mediterranean Basin drylands (Huang et al., 2016) are not expected to greatly affect shrub encroachment. Still, relative shrub cover decreased with increasing summer precipitation across the studied communities. This may be because precipitation during the warmer season (summer) usually leads to low percolation of water into the soil, thus favoring shallower-rooted herbaceous species in detriment of shrubs (Sala et al., 1997). Previous studies based on manipulative experiments found increased woody plant encroachment with higher precipitation intensity (shorter events with higher precipitation), but not with the total precipitation amount (Kulmatiski and Beard, 2013). The authors argue that a higher precipitation intensity can push soil water deeper into the soil, in shrub rooting depth, thus favoring woody species over grasses. The intensity and frequency of extreme precipitation events and intra and inter-annual precipitation fluctuations are likely to increase under a climate change scenario (IPCC, 2014). However, as we did not assess variations in precipitation distribution in our observational field study, we cannot confirm nor discard its potential effect on shrub encroachment, based on our findings. Also, since our study was developed within a regional climatic gradient where the aridity index ranged from 0.42 to 0.56, further studies along broader climatic gradients comprising more extreme aridity levels would be needed to fully assess the potential effect of climate change on shrub encroachment. 4.2. Consequences of shrub encroachment for community functional traits Changes in relative shrub cover explained most of the variation in the functional structure and dispersion of the plant communities. Shrub encroachment was associated with changes in CWMs of traits differentiating perennial taller shrubs with low SLA and larger seeds, from herbaceous short-lived and mainly anemochorous species with high SLA. These CWM variations may be seen as axes of functional specialization from a resource-conservative to a resource-acquisitive strategy to cope e.g. with water limitations (Díaz et al., 2016). Interestingly, functional dispersion of the same traits, particularly of life-cycle, root depth, height and SLA, increased with shrub encroachment. This is probably due to the co-occurrence of herbaceous and shrub species in areas with higher relative shrub cover, increasing functional diversity, when compared to grasslands devoid of shrubs. In addition, other authors have found that at high levels of shrub cover, only herbaceous species with contrasting trait values are present (Rolo et al., 2016), which would contribute to a higher functional dispersion with increasing shrub cover. The effect of topo-edaphic variables and, to a lesser extent, of climatic variables, in explaining variation along CWM-PCA1 and FDisPCA1 axes was virtually all due to changes in relative shrub cover. This indicates that main changes in plant community functional traits were mostly driven indirectly by local topo-edaphic filters, through changes in relative shrub cover. 4.3. Effects of climate on community functional traits Climate affected mostly a set of traits not closely related to changes in relative shrub cover, described by the second components of CWM and FDis PCAs. Precipitation and temperature affect many plant traits (Moles et al., 2014), and can influence the timing of phenophases (Peñuelas et al., 2004). Higher summer precipitation and higher winter temperature led to an earlier flowering and longer flowering duration. This was probably due to an attenuation of water and temperature limitations, respectively, which may provide favorable conditions for longer growth and flowering periods (Crimmins et al., 2013;Ramos et al., 2015). Previous studies reported an earlier onset of flowering triggered by higher mean annual temperatures (Miller-Rushing and Primack, 2008), supporting our results. More favorable climatic conditions also increased the CWM of prostrate and erect growth-forms, leading to a higher FDis for growth-form, and increased the CWM of ectozoochorous and N-fixing species. Attenuation of water and temperature limitations in summer and winter, respectively, may allow a higher reproductive success of these species. N-fixing species, in particular, seem to be highly vulnerable to drought during their reproductive phase, which may shorten the duration of reproductive development, reducing seed number and weight (Daryanto et al., 2015). Higher aridity levels and lower summer precipitation were associated with a lower FDis of flowering traits and of seed mass. A lower FDis is expected to reflect a lower complementary in resource use between species, suggesting a reduction in ecosystem functioning (Mouillot et al., 2011;Valencia et al., 2015) and resilience to environmental change (Volaire et al., 2014). We should also notice that our study was conducted in a relatively dry year, which may have led to a lower functional dispersion of vegetative traits such as SLA and height, particularly of annual species (Carmona et al., 2015), possibly preventing the detection of changes on the FDis of those traits along the spatial aridity gradient. 5. Concluding remarks The causes of shrub encroachment need to be better understood and quantified in the Mediterranean Basin, to improve encroachment forecasting and the cost-effectiveness and sustainability of management strategies, and our work contributes to fill this knowledge gap. We show that, under similar low-intensity land use, and despite the likely influence of other factors not addressed in this study (e.g. past land management), topographic and edaphic factors were the best predictors of shrub encroachment, influenced to a much lesser extent by climate, jointly explaining 52% of the variation in relative shrub cover. Thus, we suggest that monitoring and management actions to reduce shrub encroachment that take into account topo-edaphic factors, involving, for instance, the manipulation of water and nutrient flow pathways and accumulation, are more likely to sustainably succeed, than shrub removal alone. Climate was a weak predictor of shrub encroachment, affecting mostly a set of functional traits not so directly involved in shrub colonization. In light of this result, climate change predictions of increased aridity in the Mediterranean Basin are not expected to greatly impact shrub encroachment. Yet, they may affect the functional structure and diversity of plant communities, compromising ecosystem functioning. This work provides important advances to better predict shrub encroachment now and under a climate change scenario, and to help define more cost-effective and sustainable strategies to deal with it in Mediterranean ecosystems. 316 A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318 Acknowledgements We kindly acknowledge the help of Susana Tápia during field and laboratory work and Andreia Ferreira during field work. This research was funded by EAA Grants - Programa AdaPT and by FCT-MEC through: project PTDC/AAG-GLO/0045/2014; A.N. by contract SFRH/BD/51407/ 2011, P.P. by contract SFRH/BPD/75425/2010. The preparation of the paper was supported by NitroPortugal project [European Union's Horizon 2020 research and innovation programme under grant agreement No 692331]. NitroPortugal is a co-financer for the project Towards INMS, implemented by UN Environment. Data accessibility Data available from Dryad Digital Repository. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.scitotenv.2018.11.475. References Amaral, M.R., Borges, C., Santos, G., 1997. Montados de sobro e azinho: estrategias para o séc. XXI. Congresos Forestales. Anadón, J.D., Sala, O.E., Turner, B.L., Bennett, E.M., 2014. Effect of woody-plant encroachment on livestock production in North and South America. Proc. Natl. Acad. Sci. 111 (35), 12948–12953. Archer, S.R., 2010. Rangeland conservation and shrub encroachment: new perspectives on an old problem. Wild Rangelands: Conserving Wildlife While Maintaining Livestock in Semi-arid Ecosystems, pp. 53–97. Archer, S.R., Predick, K.I., 2014. An ecosystem services perspective on brush management: research priorities for competing land-use objectives. J. Ecol. 102 (6), 1394–1407. Bailey, D.W., Gross, J.E., Laca, E.A., Rittenhouse, L.R., Coughenour, M.B., Swift, D.M., Sims, P.L., 1996. Mechanisms that result in large herbivore grazing distribution patterns. J. Range Manag. 49, 386–400. Blois, J.L., Williams, J.W., Fitzpatrick, M.C., Jackson, S.T., Ferrier, S., 2013. Space can substitute for time in predicting climate-change effects on biodiversity. Proc. Natl. Acad. Sci. 110, 9374–9379. Bot, A., Benites, J., 2005. The Importance of Soil Organic Matter: Key to Drought-resistant Soil and Sustained Food Production. Food & Agriculture Org. Caldeira, M.C., Lecomte, X., David, T.S., Pinto, J.G., Bugalho, M.N., Werner, C., 2015. Synergy of extreme drought and shrub invasion reduce ecosystem functioning and resilience in water-limited climates. Sci. Rep. 5, 15110. Carmona, C.P., Mason, N.W., Azcárate, F.M., Peco, B., 2015. Inter-annual fluctuations in rainfall shift the functional structure of Mediterranean grasslands across gradients of productivity and disturbance. J. Veg. Sci. 26 (3), 538–551. Castro, H., Freitas, H., 2009. Above-ground biomass and productivity in the Montado: from herbaceous to shrub dominated communities. J. Arid Environ. 73, 506–511. Colgan, M., Asner, G., Levick, S.R., Martin, R., Chadwick, O., 2012. Topo-edaphic controls over woody plant biomass in South African savannas. Biogeosciences 9, 1809–1821. Crimmins, T.M., Crimmins, M.A., Bertelsen, C.D., 2013. Spring and summer patterns in flowering onset, duration, and constancy across a water-limited gradient. Am. J. Bot. 100, 1137–1147. Daryanto, S., Eldridge, D.J., Throop, H.L., 2013. Managing semi-arid woodlands for carbon storage: grazing and shrub effects on above-and belowground carbon. Agric. Ecosyst. Environ. 169, 1–11. Daryanto, S., Wang, L., Jacinthe, P.-A., 2015. Global synthesis of drought effects on food legume production. PLoS One 10, e0127401. de Bello, F., Lavorel, S., Diaz, S., Harrington, R., Cornelissen, J.H.C., Bardgett, R.D., Berg, M.P., Cipriotti, P., Feld, C.K., Hering, D., da Silva, P.M., Potts, S.G., Sandin, L., Sousa, J.P., Storkey, J., Wardle, D.A., Harrison, P.A., 2010. Towards an assessment of multiple ecosystem processes and services via functional traits. Biodivers. Conserv. 19, 2873–2893. Díaz, S., Lavorel, S., de Bello, F., Quetier, F., Grigulis, K., Robson, M., 2007. Incorporating plant functional diversity effects in ecosystem service assessments. Proc. Natl. Acad. Sci. U. S. A. 104, 20684–20689. Díaz, S., Kattge, J., Cornelissen, J.H.C., Wright, I.J., Lavorel, S., Dray, S., Reu, B., Kleyer, M., Wirth, C., Colin Prentice, I., Garnier, E., Bönisch, G., Westoby, M., Poorter, H., Reich, P.B., Moles, A.T., Dickie, J., Gillison, A.N., Zanne, A.E., Chave, J., Joseph Wright, S., Sheremet'ev, S.N., Jactel, H., Baraloto, C., Cerabolini, B., Pierce, S., Shipley, B., Kirkup, D., Casanoves, F., Joswig, J.S., Günther, A., Falczuk, V., Rüger, N., Mahecha, M.D., Gorné, L.D., 2016. The global spectrum of plant form and function. Nature 529 (7585), 167. D'odorico, P., Okin, G.S., Bestelmeyer, B.T., 2012. A synthetic review of feedbacks and drivers of shrub encroachment in arid grasslands. Ecohydrology 5, 520–530. Dubuis, A., Rossier, L., Pottier, J., Pellissier, L., Vittoz, P., Guisan, A., 2013. Predicting current and future spatial community patterns of plant functional traits. Ecography 36, 1158–1168. Eldridge, D.J., Soliveres, S., 2015. Are shrubs really a sign of declining ecosystem function? Disentangling the myths and truths of woody encroachment in Australia. Aust. J. Bot. 62, 594–608. Eldridge, D.J., Bowker, M.A., Maestre, F.T., Roger, E., Reynolds, J.F., Whitford, W.G., 2011. Impacts of shrub encroachment on ecosystem structure and functioning: towards a global synthesis. Ecol. Lett. 14, 709–722. ESRI, 2010. ArcMap v. 10.1. ESRI. Garnier, E., Lavorel, S., Ansquer, P., Castro, H., Cruz, P., Dolezal, J., Eriksson, O., Fortunel, C., Freitas, H., Golodets, C., 2007. Assessing the effects of land-use change on plant traits, communities and ecosystem functioning in grasslands: a standardized methodology and lessons from an application to 11 European sites. Ann. Bot. 99, 967–985. Gómez-Plaza, A., Martınez-Mena, M., Albaladejo, J., Castillo, V., 2001. Factors regulating spatial distribution of soil water content in small semiarid catchments. J. Hydrol. 253, 211–226. Gómez-Rey, M.X., Madeira, M., Gonzalez-Prieto, S.J., Coutinho, J., 2013. Soil C and N dynamics in a Mediterranean oak woodland with shrub encroachment. Plant Soil 371, 339–354. Grace, J.B., Schoolmaster, D.R., Guntenspergen, G.R., Little, A.M., Mitchell, B.R., Miller, K.M., Schweiger, E.W., 2012. Guidelines for a graph-theoretic implementation of structural equation modeling. Ecosphere 3, 1–44. Gross, N., Boerger, L., Soriano-Morales, S.I., Le Bagousse-Pinguet, Y., Quero, J.L., GarciaGomez, M., Valencia-Gomez, E., Maestre, F.T., 2013. Uncovering multiscale effects of aridity and biotic interactions on the functional structure of Mediterranean shrublands. J. Ecol. 101, 637–649. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. Hirota, M., Holmgren, M., Van Nes, E.H., Scheffer, M., 2011. Global resilience of tropical forest and savanna to critical transitions. Science 334, 232–235. Huang, J., Yu, H., Guan, X., Wang, G., Guo, R., 2016. Accelerated dryland expansion under climate change. Nat. Clim. Chang. 6 (2), 166. Iannone, B.V., Heneghan, L., Rijal, D., Wise, D.H., 2015. Below-ground causes and consequences of woodland shrub invasions: a novel paired-point framework reveals new insights. J. Appl. Ecol. 52, 78–88. IPCC, 2014. In: Core Writing Team, Pachauri, R.K., Meyer, L.A. (Eds.), Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. IPCC, Geneva, Switzerland (151 pp.). Jackson, R.B., Banner, J.L., Jobbágy, E.G., Pockman, W.T., Wall, D.H., 2002. Ecosystem carbon loss with woody plant invasion of grasslands. Nature 418 (6898), 623. Knapp, A.K., Briggs, J.M., Collins, S.L., Archer, S.R., Bret-Harte, M.S., Ewers, B.E., ... Cleary, M.B., 2008. Shrub encroachment in North American grasslands: shifts in growth form dominance rapidly alters control of ecosystem carbon inputs. Glob. Chang. Biol. 14 (3), 615–623. Kulmatiski, A., Beard, K.H., 2013. Woody plant encroachment facilitated by increased precipitation intensity. Nat. Clim. Chang. 3 (9), 833. Laliberte, E., Legendre, P., 2010. A distance-based framework for measuring functional diversity from multiple traits. Ecology 91, 299–305. Laliberté, E., Legendre, P., Shipley, B., 2014. FD: Measuring Functional Diversity From Multiple Traits, and Other Tools for Functional Ecology. R Package Version 1.0-12. López-Díaz, M.L., Rolo, V., Benítez, R., Moreno, G., 2015. Shrub encroachment of Iberian dehesas: implications on total forage productivity. Agrofor. Syst. 89 (4), 587–598. Maestre, F.T., Cortina, J., 2005. Remnant shrubs in Mediterranean semi-arid steppes: effects of shrub size, abiotic factors and species identity on understorey richness and occurrence. Acta Oecol. 27 (3), 161–169. Maestre, F.T., Bowker, M.A., Puche, M.D., Belén Hinojosa, M., Martínez, I., García-Palacios, P., Castillo, A.P., Soliveres, S., Luzuriaga, A.L., Sánchez, A.M., 2009. Shrub encroachment can reverse desertification in semi-arid Mediterranean grasslands. Ecol. Lett. 12, 930–941. Maestre, F.T., Eldridge, D.J., Soliveres, S., 2016. A multifaceted view on the impacts of shrub encroachment. Appl. Veg. Sci. 19, 369–370. McAuliffe, J., 2003. The atmosphere–biosphere interface: the importance of soils in arid and semi-arid environments. Changing Precipitation Regimes in Terrestrial Ecosystems: A North American Perspective. University of Arizona Press, Tucson, Ariz, pp. 9–27. MEA, 2005. Ecosystems and Human Well-being: Desertification Synthesis. World Resources Institute. Middleton, N., Thomas, D., 1992. World Atlas of Desertification: United Nations Environmental Programme. Arnold. Miller-Rushing, A.J., Primack, R.B., 2008. Global warming and flowering times in Thoreau's Concord: a community perspective. Ecology 89, 332–341. Moles, A.T., Perkins, S.E., Laffan, S.W., Flores-Moreno, H., Awasthy, M., Tindall, M.L., Sack, L., Pitman, A., Kattge, J., Aarssen, L.W., 2014. Which is a better predictor of plant traits: temperature or precipitation? J. Veg. Sci. 25, 1167–1180. Mouillot, D., Villeger, S., Scherer-Lorenzen, M., Mason, N.W.H., 2011. Functional structure of biological communities predicts ecosystem multifunctionality. PLoS One 6. Nunes, A., Tápia, S., Pinho, P., Correia, O., Branquinho, C., 2014. Advantages of the pointintercept method for assessing functional diversity in semi-arid areas. iForest Biogeosci. For. 884. Nunes, A., Köbel, M., Pinho, P., Matos, P., de Bello, F., Correia, O., Branquinho, C., 2017. Which plant traits respond to aridity? A critical step to assess functional diversity in Mediterranean drylands. Agric. For. Meteorol. 239, 176–184. Pakeman, R.J., Quested, H.M., 2007. Sampling plant functional traits: what proportion of the species need to be measured? Appl. Veg. Sci. 10, 91–96. Peñuelas, J., Filella, I., Zhang, X., Llorens, L., Ogaya, R., Lloret, F., Comas, P., Estiarte, M., Terradas, J., 2004. Complex spatiotemporal phenological shifts as a response to rainfall changes. New Phytol. 161, 837–846. 317A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318 Pereira, P.M., Da Fonseca, M.P., 2003. Nature vs. nurture: the making of the montado ecosystem. Conserv. Ecol. 7 (3). Pérez-Harguindeguy, N., Díaz, S., Garnier, E., Lavorel, S., Poorter, H., Jaureguiberry, P., BretHarte, M., Cornwell, W., Craine, J., Gurvich, D., 2013. New handbook for standardised measurement of plant functional traits worldwide. Aust. J. Bot. 61, 167–234. R Core Team, 2015. R: A Language and Environment for Statistical Computing. URL. http:// www.R-project.org (Vienna, Austria). Ramos, A., Pereira, M.J., Soares, A., do Rosário, L., Matos, P., Nunes, A., Branquinho, C., Pinho, P., 2015. Seasonal patterns of Mediterranean evergreen woodlands (Montado) are explained by long-term precipitation. Agric. For. Meteorol. 202, 44–50. Rango, A., Huenneke, L., Buonopane, M., Herrick, J., Havstad, K., 2005. Using historic data to assess effectiveness of shrub removal in southern New Mexico. J. Arid Environ. 62, 75–91. Rivest, D., Rolo, V., López-Díaz, L., Moreno, G., 2011. Shrub encroachment in Mediterranean silvopastoral systems: Retama sphaerocarpa and Cistus ladanifer induce contrasting effects on pasture and Quercus ilex production. Agric. Ecosyst. Environ. 141 (3-4), 447–454. Rolo, V., López-Díaz, M.L., Moreno, G., 2012. Shrubs affect soil nutrients availability with contrasting consequences for pasture understory and tree overstory production and nutrient status in Mediterranean grazed open woodlands. Nutr. Cycl. Agroecosyst. 93 (1), 89–102. Rolo, V., Rivest, D., Lorente, M., Kattge, J., Moreno, G., 2016. Taxonomic and functional diversity in Mediterranean pastures: insights on the biodiversity–productivity tradeoff. J. Appl. Ecol. 53 (5), 1575–1584. Roques, K., O'connor, T., Watkinson, A., 2001. Dynamics of shrub encroachment in an African savanna: relative influences of fire, herbivory, rainfall and density dependence. J. Appl. Ecol. 38, 268–280. Rosseel, Y., 2012. lavaan: an R package for structural equation modeling. J. Stat. Softw. 48, 1–36. Sala, O.E., Maestre, F.T., 2014. Grass–woodland transitions: determinants and consequences for ecosystem functioning and provisioning of services. J. Ecol. 102, 1357–1362. Sala, O., Lauenroth, W., Golluscio, R., 1997. Plant functional types in temperate semi-arid regions. Plant Functional Types: Their Relevance to Ecosystem Properties and Global Change, pp. 217–233. Sankaran, M., Hanan, N.P., Scholes, R.J., Ratnam, J., Augustine, D.J., Cade, B.S., ... Ardo, J., 2005. Determinants of woody cover in African savannas. Nature 438 (7069), 846. Schlesinger, W.H., Reynolds, J.F., Cunningham, G.L., Huenneke, L., Jarrell, W., Virginia, R., Whitford, W., 1990. Biological feedbacks in global desertification. Science 247, 1043–1048. Schwilch, G., Hessel, R., Verzandvoort, S. (Eds.), 2012. Desire for Greener Land. Options for Sustainable Land Management in Drylands. University of Bern - CDE, Alterra - Wageningen UR, ISRIC - World Soil Information and CTA - Technical Centre for Agricultural and Rural Cooperation, Bern, Switzerland, and Wageningen, The Netherlands. Soliveres, S., Maestre, F.T., Eldridge, D.J., Delgado-Baquerizo, M., Quero, J.L., Bowker, M.A., Gallardo, A., 2014. Plant diversity and ecosystem multifunctionality peak at intermediate levels of woody cover in global drylands. Glob. Ecol. Biogeogr. 23 (12), 1408–1416. Sørensen, R., Zinko, U., Seibert, J., 2005. On the calculation of the topographic wetness index: evaluation of different methods based on field observations. Hydrol. Earth Syst. Sci. Discuss. 2, 1807–1834. Tilman, D., Knops, J., Wedin, D., Reich, P., Ritchie, M., Siemann, E., 1997. The influence of functional diversity and composition on ecosystem processes. Science 277, 1300–1302. Trabucco, A., Zomer, R., 2009. Global aridity index (global-aridity) and global potential evapo-transpiration (global-PET) geospatial database. CGIAR Consortium for Spatial Information Published online, available from the CGIAR-CSI GeoPortal at:. http:// www.csi.cgiar.org/. Valencia, E., Maestre, F.T., Le Bagousse-Pinguet, Y., Luis Quero, J., Tamme, R., Boerger, L., Garcia-Gomez, M., Gross, N., 2015. Functional diversity enhances the resistance of ecosystem multifunctionality to aridity in Mediterranean drylands. New Phytol. 206, 660–671. Van Auken, O., 2009. Causes and consequences of woody plant encroachment into western North American grasslands. J. Environ. Manag. 90, 2931–2942. Volaire, F., Barkaoui, K., Norton, M., 2014. Designing resilient and sustainable grasslands for a drier future: adaptive strategies, functional traits and biotic interactions. Eur. J. Agron. 52, 81–89. 318 A. Nunes et al. / Science of the Total Environment 657 (2019) 310–318