Traits of dominant plant species drive normalized difference vegetation index in grasslands globally

Aim: Theoretical, experimental and observational studies have shown that biodiversity– ecosystem functioning (BEF) relationships are influenced by functional community

The biodiversity-ecosystem functioning (BEF) debate has become even more relevant recently, given that the ongoing global biodiversity crisis (Ceballos et al., 2015;Díaz et al., 2019;McGill et al., 2015) requires us to understand how changes in the diversity and composition of biological communities (especially plants) will affect ecosystem processes and services (Pillar et al., 2018).
Over recent decades, the focus of BEF research has shifted from considering species richness alone to investigating how plant functional community structure (i.e., the distribution of trait values measured in a given community; Garnier et al., 2016) affects ecosystem processes (Cadotte et al., 2011(Cadotte et al., , 2013;;Dias et al., 2013;Fischer et al., 2016;Grigulis et al., 2013;Kröber et al., 2015).Such traitbased approaches have offered two main hypotheses for a mechanistic explanation of BEF relationships.The dominance hypothesis (also known as the mass ratio hypothesis) states that ecosystem processes are determined primarily by the traits of the dominant species in a community (Grime, 1998).It suggests that ecosystem functioning is linked to the community-weighted mean (CWM) of relevant traits (Garnier et al., 2004).The niche partitioning hypothesis, in contrast, refers to the variation of a trait value within a community.
It suggests that a higher range of trait values [i.e., functional diversity] reflects niche differences, allowing for enhanced use of biotope space, better resource use efficiency and, ultimately, enhanced ecosystem functioning (Petchey & Gaston, 2002).Dominance and niche partitioning effects are not mutually exclusive, and disentangling them has proved difficult.For example, it is challenging to manipulate functional community structure experimentally such that functional diversity and CWM vary independently of each other (Dias et al., 2013;Ricotta & Moretti, 2011).
Furthermore, in natural assemblages, both plant traits and ecosystem functioning are strongly influenced by environmental drivers, such as climate, and exhibit strong spatial autocorrelation.Traditional experimental approaches have tried to keep abiotic factors constant and have therefore been restricted to small spatial and temporal scales, whereas more recent approaches have specifically included or manipulated environmental variation, such as climate and nutrient supply (Chollet et al., 2014;Craven et al., 2016;Roscher et al., 2013;Zhou et al., 2017).Overall, experimental and observational studies tend to find dominance effects more frequently and with stronger evidence than niche partitioning effects (Chollet et al., 2014;Kröber et al., 2015;Lavorel, 2013;Mokany et al., 2008;Mouillot et al., 2011;van der Plas, 2019).In particular, plant traits related to the trade-off between resource acquisition and conservation show strong dominance effects on productivity (Díaz et al., 2004;Grigulis et al., 2013;Lavorel, 2013).Recent studies, however, have shown that both niche partitioning and dominance effects contribute to ecosystem multifunctionality (i.e., providing multiple functions; Hector & Bagchi, 2007;Mouillot et al., 2011;Valencia et al., 2015) and to biomass production outside the peak productivity season (Chollet et al., 2014).Despite these advances, there continues to be a lack of studies and syntheses focusing on large spatial scales (but see Cornwell et al., 2008), and we know little about how dominance and niche partitioning effects interplay with climatic factors to shape biomass production along larger biogeographical gradients.
Remote sensing observations might allow this knowledge gap to be closed (Franklin et al., 2017;Polley et al., 2020;Rocchini et al., 2021).Remote Earth observation products have become widely available and allow the assessment of ecosystem functionality and biodiversity (Lausch et al., 2016;Schmidtlein et al., 2012).

Satellite-derived normalized difference vegetation index (NDVI),
for instance, can be used to measure surface greenness over wide extents as a proxy of primary productivity (Paruelo et al., 1997;Running, 1990) and other ecosystem processes (Ustin et al., 2004), including the stability of productivity (White et al., 2020).As such, NDVI has been used in ecological studies to evaluate ecosystem responses to environmental changes (reviewed by Pettorelli et al., 2005), to measure ecosystem stability (De Keersmaecker et al., 2014) and to study the effect of plant functional community structure on resilience (Spasojevic et al., 2016).Yet, no study to date has linked NDVI to plant functional community structure over large extents to disentangle the relative contributions of the dominance and niche partitioning effects.
Here, we investigate the BEF relationship in grassland systems at the global scale.We rely on the world's most comprehensive the importance of community assembly processes for BEF relationships in natural communities.

Main conclusions:
Our analysis provides empirical evidence that plant functional community structure and global patterns in primary productivity are linked through the resource economics and size traits of the dominant species.This is an important test of the hypotheses underlying BEF relationships at the global scale.

K E Y W O R D S
biodiversity-ecosystem functioning, biodiversity, community-weighted mean, ecosystem functioning, functional diversity, sPlot, traits, vegetation repository of plant community data, sPlot (Bruelheide et al., 2019), coupled with the plant trait database TRY (Kattge et al., 2011), to quantify plant functional community structure in grassland communities around the globe.We used maximum annual NDVI as a proxy for primary productivity and modelled it using plant functional community data and climatic variables.Specifically, we were interested in the role of CWM (as a proxy for dominance) and functional diversity (as a proxy for niche partitioning) for predicting the biogeographical patterns of NDVI.In addition, we discerned the unique and shared contributions of our functional and climatic predictors using a deviance partitioning approach.

| Datasets
We based this study on the integration of data compiled from three global initiatives: (1) the vegetation plot database sPlot, v.2.1 (Bruelheide et al., 2019), combined with (2) the TRY database of plant functional traits (Kattge et al., 2011) and (3) MODIS satellite data providing NDVI time series (Didan, 2015).Additionally, we used the land cover map GlobCover 2009 (Arino et al., 2012) and temperature and precipitation data from the CHELSA project (Karger et al., 2016).

| Filtering of sPlot
The sPlot database, v.2.1, contains information on plant community composition in >1.1 million vegetation plots contributed by >100 vegetation plot databases.In line with a long tradition of BEF research in grassland ecosystems (Fraser et al., 2015;Hooper et al., 2005;Tilman et al., 1996), we chose grasslands as a model system for this analysis.We screened the database for vegetation plots in grassland communities.Any selected plot had to satisfy at least one of the following criteria: (1) the plot was marked as "grassland" by the vegetation survey that provided the sPlot entry; and (2) ≥90% of the plot was covered by species that were not trees or shrubs or that had a plant height >2 m (according to the TRY database).We excluded plots that were labelled as "forest", "shrubland" or "wetland", observations before the year 2000 (launch of the MODIS programme) and plots with high location uncertainty (>100 m) and imprecise GPS coordinates (less than four decimal places in decimal degrees).We used the land cover map GlobCover 2009 to identify and exclude plots on heterogeneous MODIS pixels that contained non-grassland land cover types.In total, 2941 grassland plots fulfilled the selection criteria and were included in the analysis (Figure 1).The median size plot size was 78.5 m 2 , with an interquartile range from 20 to 80 m 2 .

| Annual maximum NDVI
We used annual maximum NDVI as a proxy for primary productivity at the selected grassland plots.This value corresponds to the yearly maximum photosynthetic activity of the vegetation and is considered an indicator for ecosystem productivity (Pettorelli et al., 2005).
We obtained the annual maximum NDVI values from the MODIS product MOD13Q1, which has global coverage, a spatial resolution of 250 m and a temporal resolution of 16 days (Didan, 2015).year and ends in the next one.In cases where an NDVI pixel contained more than one vegetation plot, surveys from different years were treated as independent observations.Whenever there were multiple vegetation plots per year and NDVI pixel, one plot was selected at random.

| Functional community structure
We selected 18 ecologically relevant plant traits to characterize the functional community structure at all selected plots (Figure 2).The same set of traits has already been used by Bruelheide et al. (2018) and Testolin et al. (2021) to analyse global trait-environment associations.These traits represent different ecological trade-offs and plant strategies, such as the leaf economics spectrum (Wright et al., 2004), plant size and reproduction (Moles & Westoby, 2006).
Although our trait selection was ultimately driven by data availability, we note that most of the selected traits relate to the leaf economics spectrum of a plant (i.e., to its carbon gain strategy; Wright et al., 2004).We considered these traits as "effect traits" because of their role for biomass production (Suding et al., 2008).We acknowledge that other important aspects of plant physiology, particularly those related to C 3 /C 4 photosynthetic pathways, belowground traits and, among them, those associated with symbionts such as mycorrhiza or nitrogen-fixing bacteria, are underrepresented in our trait selection, although these traits have long been known to be related strongly to productivity in grasslands (Cadotte et al., 2009;Díaz & Cabido, 2001;Weigelt et al., 2021).For more information on the selected traits and their ecological relevance, see the Supporting Information (Appendix S1).Species mean traits were retrieved from TRY (Kattge et al., 2011) and gap-filled using Bayesian hierarchical probabilistic matrix factorization (Schrodt et al., 2015).The gapfilled trait data were available for 88.7% of all species occurrences in sPlot (Bruelheide et al., 2018).We ln-transformed all trait values for downstream analysis.
We calculated the CWM for all traits using the following equation: where n k is the number of species in plot k, p i,k is the relative cover of species i in plot k, and t i,l is the mean value of species i for trait l.The CWM is an abundance-weighted trait mean value, which is most affected by the dominant species in the community.Accordingly, CWMs are often linked to the mass ratio effect and dominance hypothesis (Garnier et al., 2004).
(1) To reduce the functional trait space described by the 18 CWM variables, we performed a principal components analysis (PCA) on the CWM values of the selected sPlot communities using the R package FactoMineR (Lê et al., 2008).The first two PCA axes of the CWM space captured 27.42% (CWM axis 1) and 20.88% (CWM axis 2) of the total variance in the data (Figure 2).Communities that had negative scores on CWM axis 1 were dominated by species with "acquisitive" traits, such as high specific leaf area and high leaf concentrations of nitrogen and phosphorus.Positive values on this axis were accompanied by a dominance of species with "conservative" traits, such as high leaf dry matter content and stem density (Supporting Information Appendix S2).This axis corresponds to the resource economics spectrum described by many studies (i.e., leaf economics spectrum, fast-slow spectrum; Díaz et al., 2016;Reich et al., 1997;Wright et al., 2004).The CWM axis 2 was related to overall plant size: communities with high scores on this axis were characterized by relatively tall species, with both large and heavy leaves and seeds.This corresponds to a plant size spectrum that is commonly found in plants (Díaz et al., 2016).We used the two axes as composite variables characterizing the community position along the resource economics and plant size spectrum.In the subsequent analysis, we interpreted associations between the CWM axes and NDVI as evidence for dominance effects.
To quantify functional diversity across all 18 traits, we calculated functional dispersion (hereafter, FD) using the R package FD (Laliberté et al., 2014) and taking Gower's distance as the dissimilarity measure (Podani, 1999).FD has the advantage of accommodating multiple traits in addition to species abundances.Furthermore, it is independent of species richness and relatively unaffected by outliers (Laliberté & Legendre, 2010).Here, we interpret associations between FD and NDVI as evidence for niche partitioning effects.

| Data analysis
We modelled annual maximum NDVI using a number of nested generalized additive models (GAMs; Wood, 2017) with subsets of predictors relating to traits (CWM axis 1, CWM axis 2 and FD), climate (temperature and precipitation) and space (latitude and longitude) (Table 1).We chose the GAM framework because it allowed us to include the spatial coordinates as smooth spherical splines, which accounts for the spatial structure of the NDVI response variable and spatial autocorrelation in the residuals.All other predictors were included as simple linear terms.As climatic variables, we included mean annual temperature (Temp) and the precipitation of the driest quarter (Prec).With the response variable NDVI ranging between zero and one, we specified the models using the beta-regression family and a logit link function, using the mgcv R package (Wood, 2017).All linear predictors were scaled and centred before entering the model.
After fitting all models, we used a deviance partitioning approach to disentangle the predictive effects of climate, traits and space on NDVI.To obtain the unique fraction explained by a set of predictors, we calculated the difference in model deviance of nested models with and without the given set of predictors and expressed it as a percentage of the total deviation of the full model.Shared fractions between two predictor sets were calculated by subtracting their unique contributions from the joint model deviance.Taking the same approach, we also partitioned the contribution of traits further

TA B L E 1 Fitted models predicting annual maximum normalized difference vegetation index (NDVI)
into the joint and unique effects of the two CWM axes and FD.All analyses were run in R v.3.5.3 (R Core Team, 2019).

| RE SULTS
The full model explained 86% of the deviance in the NDVI and had an adjusted R 2 of 0.87 (Table 1).All linear predictors except FD were significant (Figure 3a).Temperature and precipitation had positive effects on annual maximum NDVI.CWM axis 1 had a negative effect on annual maximum NDVI (i.e., communities dominated by species with acquisitive traits tended to have higher NDVI than those dominated by species with conservative traits), and CWM axis 2 (i.e., plant size) had a positive effect (Figure 3a).The model coefficients and significance were similar across the different nested models, and the effect of FD was significant only when no spatial term was included in the model (Supporting Information Appendix S3).All models including the spatial term showed by far the best fit (Table 1).Generally, the spatial term was highest in mid-to high latitudes of Eurasia and lowest in parts of Africa and Australia (Figure 3b).The importance of space was also reflected by the "space only" model, which explained nearly as much deviance as the full model.The largest fraction of the explained model deviance was shared between the three variable sets of climate, traits and space (38.16%), followed by the unique fraction of space (22.44%), the shared fractions of space and traits (13.52%) and space and climate (11.50%).The unique and shared fractions of climate and traits were relatively small, but traits had a larger unique importance than climate.Most of the unique fraction of traits was linked to the CWM (especially axis 1) and not to FD (Figure 4a).Notably, even the "only traits" and "only climate" models showed relatively high predictive power, with R 2 values of 0.54 and 0.55, respectively.Most of total deviance explained by traits was associated with the unique fraction of CWM axis 1, whereas FD and CWM axis 2 made smaller contributions to the trait effect.There also was a sizeable fraction shared between CWM axis 1 and FD (Figure 4b).For the pairwise correlations between all model variables, see the Supporting Information (Appendix S4).

| DISCUSS ION
A large body of research revolves around the question of how traits contribute to ecosystem functioning.Using data from grasslands all over the globe, we found strong support for a dominance effect for peak season productivity (i.e., dominant plant species contribute disproportionately to ecosystem productivity).At the same time, our findings failed to provide evidence for the niche partitioning hypothesis.For the large spatial extent and the traits covered by our study, we did not find that higher trait diversity resulted in enhanced biomass production.
Our results show that grasslands dominated by species with acquisitive traits have higher productivity than those dominated by species with conservative traits.This is in line with existing evidence from observational studies of smaller extents and experiments that point out the role of dominant species for determining ecosystemlevel primary productivity through their resource economics traits (Lavorel, 2013;Mokany et al., 2008).Acquisitive species with high specific leaf area and leaf nutrient concentrations generally exhibit rapid carbon and nitrogen turnover, which entails high ecosystemlevel rates of, for instance, biomass accumulation, decomposition and evapotranspiration (Reich, 2014).Our study suggests that this relationship also holds at a larger spatial scale, because we found this resource axis to be an important predictor of NDVI, which is a proxy for primary productivity.Traits related to overall plant, leaf and seed size were also positively associated with productivity.This shows that the NDVI signal is not captured fully by leaf economics traits alone but also requires a second trait dimension as an estimation of stand biomass (Chave et al., 2005).This axis is captured by plant height and seed size, which are allometrically linked (Moles & Westoby, 2006;Rees & Venable, 2007;Thompson & Rabinowitz, 1989).Although ecological theory and many experiments show that niche partitioning can enhance ecosystem functions (Hooper et al., 2005;Petchey & Gaston, 2002;van der Plas, 2019), we found no association between functional diversity and grassland productivity at the global scale, at least for the traits we used.Overall, the recent literature suggests that the effect of plant functional diversity on ecosystem processes is varied and subtle, especially in naturally assembled communities (van der Plas, 2019).Although a majority of studies report positive relationships (Cardinale et al., 2011;Mouillot et al., 2011;Petchey & Gaston, 2006;van der Plas, 2019), it is not uncommon to find negative (Thompson et al., 2005) or no associations (Chollet et al., 2014) between functional diversity and ecosystem functions.The effect of functional diversity is also known to vary among traits (Petchey et al., 2004).Although we covered the major plant trait axes identified by Díaz et al. (2016), we cannot exclude the possibility that some other traits would lead to a stronger association between functional diversity and NDVI.For example, we did not have data on nitrogen use (e.g., N fixation) strategies and photosynthetic pathways that are known to play a role for complementarity in grasslands (Cadotte et al., 2009).The selection of the functional traits affecting an ecosystem function remains difficult, because different traits might play different roles across communities or biomes.A new methodology has recently been published to select those traits that contribute to functional community assembly (Pillar et al., 2021).Whether this technique can prove useful also for selecting those traits related to a specific ecosystem function, however, remains to be tested.To test the sensitivity of our results to the choice of the functional diversity metric, we also performed a supplementary analysis using a functional diversity measure based on Rao's quadratic entropy (Rao, 1982) rather than FD.Overall, those results were qualitatively consistent with the findings presented here in the main text, but the effect of functional diversity was even smaller (Supporting Information Appendix S5).
Another possible explanation of our results relates to the use of annual maximum NDVI as a response variable.This value represents the high-productivity season, when resource turnover is expected to be at its maximum (Pettorelli et al., 2005).Previous work shows that niche partitioning effects in plant communities might be stronger outside the high-productivity season (Chollet et al., 2014;Doležal et al., 2019).Possibly, the role of functional diversity might be more substantial when quantifying productivity as an integrated measure over the whole vegetative season, rather than a snapshot at the peak of the vegetation season (Chi et al., 2017).
Although grassland productivity ultimately depends on plants, our models highlighted the existence of important confounding variables, namely climate and other unmeasured factors that are spatially structured.It is well known that precipitation (La Pierre et al., 2011;Sala et al., 1988) and temperature (Briggs & Knapp, 1995) are important determinants of grassland productivity and often outweigh the pure effect of biodiversity (van der Plas, 2019).Although substantial, the deviance explained by climate was mostly shared with plant traits, highlighting the importance of climate as an environmental filter shaping the functional structure of vegetation assemblages (Bruelheide et al., 2018).Indeed, the resource economics axis described by CWM axis 1 is known to be linked to droughttolerance strategies (Reich, 2014).Both climate and vegetation are spatially structured and therefore expected to be collinear (Currie et al., 2020).In fact, the largest fraction of the explained model deviance was shared between space, climate and traits, following Tobler's well-known first law of geography: "everything is related to everything else, but near things are more related than distant things" (Tobler, 1970).Nonetheless, even when we included the smooth spatial term in the model, the effect of the functional structure (i.e., CWM) of vegetation assemblages remained significant, which suggests that the relationship between plant traits and productivity goes beyond the mere spatial covariation in abiotic conditions and local species pools.
There is no doubt that the approach adopted in this study has some limitations.Not only does the use of trait databases lead to a very rough approximation of functional community structure (e.g., neglecting intraspecific and site-specific trait variation) but also the vegetation plot database had vast gaps in its global coverage, which was partly attributable to the strict filtering criteria applied here.South America and Africa are particularly underrepresented in our compilation.Nevertheless, our analysis has an unprecedented spatial coverage, and spans more than nine of the ten sPlot biomes (sensu Bruelheide et al., 2018).Furthermore, by explicitly including geographical coordinates in our models, we mitigated as much as The fact that we limited our analysis only to vegetation plots located within MODIS pixels representing uniform land cover classes, however, is expected to mitigate the negative consequences of this scale mismatch.Furthermore, NDVI, like other vegetation indices, is known to saturate in highly productive systems, which means that it might not correctly discern productivity differences in those areas (e.g., the tropics).However, this mostly concerns forests with high leaf area index values and should not be a major issue for our sites in open grasslands (Haboudane et al., 2004).Finally, although the selected functional traits account for two of the main trait dimensions in plants (leaf economics and size), data coverage was limited for other important aspects, such as photosynthetic pathways and association with nitrogen-fixing bacteria.We agree that it would have been very desirable to have more specific traits, example, on C 3 /C 4 photosynthetic pathways or bacterial associations.More importantly, we missed belowground traits in general.We cannot exclude, therefore, that such traits could increase the relative contribution of niche partitioning mechanisms to ecosystem productivity patterns.
In conclusion, our analysis provides empirical evidence that plant functional community structure and global patterns in primary productivity are mostly linked through the resource economics and size traits of the dominant species.Over large extents, the dominance hypothesis received substantially more support than the niche partitioning hypothesis, at least based on the functional traits included in the present study.This is consistent with the hypothesis that productivity mostly relates to the functional profile of dominant species, rather than to the functional diversity of a community.Our results also provide an empirical justification for focusing on dominant plant species when running Earth system models simulating the response of vegetation to changing conditions.Identifying the dominant traits and species and assessing their susceptibility to environmental change then becomes crucial to improve our ability to foresee, and possibly mitigate negative consequences on the provisioning of plant biomass by grasslands.

Furthermore
, owing to the high image frequency of MODIS, its NDVI product has a relatively low susceptibility to cloud cover, which is why we chose it over Landsat.For each of the selected vegetation plots, we retrieved NDVI time series covering the period 2000-2016 using Google Earth Engine (Google Earth Engine Team, 2015) and averaged the maximum annual NDVI values of the year when the vegetation plot was sampled and the following year.This averaging across two consecutive years was done to reduce the effect of anomalies and noise in the NDVI values and bearing in mind that the growing season in the Southern Hemisphere starts in one calendar F I G U R E 1 Locations of included vegetation plots coloured by annual maximum normalized difference vegetation index (NDVI).
Variables factor map on the first two principal components analysis (PCA) axes for community-weighted mean (CWM) traits.Arrows indicate the strength and direction of correlation between the trait CWM and the axes.We interpret CWM axis 1 as a resource economics spectrum (acquisitive vs. conservative) axis and CWM axis 2 as a plant size axis.Traits are as follows: specific leaf area (SLA), leaf P concentration (Leaf.P), leaf N concentration (Leaf.N), leaf area (Leaf.A), number of seeds of the reproductive unit (Rel.seed.num),leaf fresh mass (Leaf.fr.mass), stem conduit density (Stem.cond.dens),leaf δ 15 N (dN15), plant height (Height), dispersal unit length (Disp.u.l), wood vessel length (Vessel.l),seed length (Seed.l),leaf C per dry mass (Leaf.Cpmass), seed mass (Seed.m),leaf N/P ratio (NpP), leaf N per area (Leaf.NpA), leaf dry matter content (LDMC) and stem specific density (Stem.Dens).
E 3 Estimates of the full model.(a) Standardized effect sizes of linear predictors (mean and confidence intervals).Negative values are red; positive values are blue.(b) Estimated smooth spatial term at the plot locations.Abbreviations: CWM, community-weighted mean; FD, functional dispersion.

F
Venn diagram partitioning total explained model deviance (i.e., 85.99%) among sets of predictors.(a) Partitioning total deviance among climate (mean annual temperature and precipitation of the driest quarter), traits (i.e., CWM axis 1, CWM axis 2 and FD) and space (i.e., smooth term using spatial coordinates).(b) Partitioning the total fraction of traits [i.e., yellow circle in (a)] between CWM axis 1, CWM axis 2 and FD.Values are expressed as percentages.14668238, 2023, 5, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/geb.13644by JOSEP Penuelas -Readcube (Labtiva Inc.) , Wiley Online Library on [12/04/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License possible the bias deriving from the uneven distribution of plots, besides partly controlling for other spatially structured drivers of productivity, such as soil fertility and management.Moreover, there was a substantial scale mismatch between the vegetation surveys and the NDVI measure.Even though the MODIS images have a pixel size of 250 m, which is roughly two orders of magnitude coarser than a typical grassland plot, this was the best possible compromise in the trade-off between temporal coverage and spatial resolution.
sPlot has been initiated by sDiv, the Synthesis Centre of the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, funded by the German Research Foundation (FZT 118) and is now a platform of iDiv.We are grateful to thousands of vegetation scientists who sampled vegetation plots in the field or digitized them into regional, national or international databases.We appreciate the support of the German Research Foundation for funding sPlot as one of the iDiv research platforms (DFG FZT 118, 202548816).T.E., D.H. and F.M.S. were funded under this scheme.We thank Ludmilla F igu eir edo and the iDiv Data and Code Unit for assistance with the curation and archiving of the dataset.The study was supported by the TRY initiative on plant traits (http://www.try-db.org).T.E. was financially supported by the European Commission through the programme Erasmus Mundus Master Course, International Master in Applied Ecology (EMMC-IMAE; FPA 532524-1-FR-2012-ERA MUNDUS-EMMC).D.H. was s upp ort ed by th e C oor dination of Superior Level Staff Improvement (CAPES).V.P. was sup-