Changes in taxonomic and functional diversity of plants in a chronosequence of Eucalyptus grandis plantations

Tree plantations have become one of the fastest-growing land uses and their impact on biodiversity was evaluated mainly at the taxonomic level. The aim of this study was to analyze environmental changes after the Eucalyptus plantation in an area originally covered by natural grasslands, taking into account the alpha and beta (taxonomic and functional) diversity of plant communities. We selected nine plantation ages, along a 12 years chronosequence, with three replicates per age and three protected grasslands as the original situation. At each replicate, we established three plots to measure plant species cover, diversity and environmental variables. Results showed that species richness, and all diversity indices, significantly declined with increasing plantation age. Canopy cover, soil pH, and leaf litter were the environmental drivers that drove the decrease in taxonomic and functional diversity of plants through the forest chronosequence. Based on the path analyses results, canopy cover had an indirect effect on plant functional diversity, mediated by leaf litter depth, soil pH, and plant species richness. The high dispersal potential, annual, barochorous, and zoochorous plant species were the functional traits more affected by the eucalypt plantations. We recommend two management practices: reducing forest densities to allow higher light input to the understory and, due to the fact that leaf litter was negatively associated with all diversity facets, we recommend reducing their accumulation or generate heterogeneity in its distribution to enhance biodiversity.

The human population is growing more rapidly than ever before and consequently increasing the demand for goods and services provided by ecosystems processes 1,2 . Among human activities that imply intensive land use, tree plantations with non-native species are expanding worldwide because of the global demand for timber, pulp and biofuel products 3 . The effects of land conversion on ecosystem services will largely depend on how native biodiversity responds to the resulting environmental changes. Some studies revealed that human land uses that preserve vegetation structure and composition of natural biomes are more used by native species than those anthropogenic habitats that cause drastic environmental changes 4,5 . Specifically, tree plantations have more structural similarities to natural forests than crop fields, then promoting suitability for forest biodiversity but less than native forest 6 . In addition, studies on biodiversity changes among different plantation ages showed that the biome where the tree plantations are developed is also relevant 7 . Mature plantations contributed to maintain biodiversity better than young plantations in forest biome 8 . The opposite pattern was observed for bird and ant communities in tree plantations developed in grassland sites, where the richness and abundance were the lowest in mature plantations and highest in the grasslands 4,5,8,9 . Several studies have already underlined the negative effects on native species and their ecosystem functions due to the expansion of forest plantations 10 . Conversely, others have pointed out that non-native tree species can restore degraded ecosystems faster than native tree plantations and have the potential to enhance above-ground biomass accumulation 11 . Plant richness, species diversity, and heterogeneity increased significantly during succession in abandoned croplands of the Loess Plateau in China 12 . However, there is evidence that that non-native Traditional beta diversity was defined as the variation in species composition 24 , and their study attempt to reveal the assembly mechanisms that drive the variation. More recently, functional beta diversity has been defined as the change in ecological functions or species traits between assemblages 25 . The inclusion of both facets of beta diversity is useful for understanding how human activities impact on biodiversity. Using this concept, Vaccaro et al. 26 have shown that bird composition and functional diversity in cattle pastures, were the most similar to native pampean grasslands than other human land uses. However, the effects of tree plantations on the functional diversity of understory vegetation during the forest cycle developing in grasslands are still unclear, especially considering that they are among the most threatened and less conserved habitats 27 . Plants are often used as biodiversity indicators because they are sensitive to abiotic conditions and land-use changes, are at the base of food webs, and provide habitat for animals influencing animal diversity and distribution 28 .
The study aims to analyze the effects of the Eucalyptus plantation chronosequence, on alpha, beta, taxonomic and functional diversity of understory plant communities. Recognizing the structural and microclimatic changes that occur during the eucalypt forest cycle, we hypothesized that the plant taxonomic and functional similarity with the predecessor natural habitat (Pampean grasslands), will decline with decreasing environmental similarity along the plantation cycle chronosequence. Therefore, we expect that (1) as plantation age increases, environmental conditions typical of grasslands will change through the forest cycle (e.g. increase in vegetation stratification, leaf litter depth, canopy cover, and decrease in soil pH); (2) species richness and alpha functional diversity will be higher in young plantations than mature plantations; and (3) plant taxonomic and functional similarity (the inverse of beta diversity) between forest plantation and native grassland, will be higher in young plantations due to the higher environmental similarity. Given that the capability of surviving in disturbed habitats is related with life-history traits (dispersion, establishment, and persistence), our predictions about functional traits are that (1) plant communities in young plantations will have similar functional characteristics to pioneer species e.g. high potential dispersal and shorter lifespan (2) loss of functional traits typical of grassland plant species will occur as plantation age increases and (3) late-successional species may persist at the end of the forest cycle 29,30 .

Results
We recorded a total of 247 plant species/morphospecies along the eucalypt forest cycle and grassland sites in the Mesopotamic Pampa (30 sites in total). We identified 47 species/morphospecies in grassland sites, 32 species/ morphospecies in both grassland and eucalypt plantations and 168 species/morphospecies only in eucalypt plantations (Appendix A). Moreover, we identified 49 families and the most dominant were Poacea (47 spp.), Compositae (47 spp.) and Leguminosae (16 spp.).
Regarding environmental changes, the first two axes of PCA analyses explained 87% of the variation (Fig. 1). PCA axis 1 arranged sites according to plantation age from negative to positive values. Grassland sites, 0-year-old and 1-year-old plantations have similar environmental conditions: high levels of soil pH, phosphorus content and temperature. Moreover, the rest of the plantations were high canopy cover, DBH and leaf litter depth.
Our results showed that plant species richness was negatively affected by soil pH, leaf litter depth, and nitrogen content (Table 1). In addition, we found that overall taxonomic similarity was also explained by similarity in leaf litter depth and canopy cover between grassland sites and the eucalypt forest plantations. The decrease of alpha functional diversity was affected by the increase of canopy cover and leaf litter depth. Finally, the decrease in functional similarity was explained by the increase of similarity in leaf litter depth, soil pH, and nitrogen content. www.nature.com/scientificreports/ Piecewise SEM analysis revealed that plant species richness affected positive and directly on alpha functional diversity (Fig. 3). In addition, the model showed that the effect of the canopy cover on plant species richness is mediated by leaf litter depth and soil pH. Leaf litter depth had a direct and indirect influence on plant species richness. The analysis reproduced the data well based on a comparison of the Fisher's C statistic to a chi-square distribution (C 8 = 5.56, p = 0.7).
The fourth corner analysis revealed a positive and strong association between canopy cover and anemozoochory as dispersion syndrome and, a negative and strong association with barochorous species (Fig. 4). In addition, leaf litter depth had a negative and strong association with plants with high dispersal potential, perennial plants and zoochory and barochory as dispersal syndromes. Soil pH had a negative association with perennial plants and herbs whereas N had a negative and strong association with annual plants and anemozoochory. Finally, nitrogen content in soil had a negative strong association with annual plants and anemozoochory as dispersal syndrome.

Discussion
Our findings reveal that changes in soil pH, canopy cover, nitrogen content in soil, and leaf litter from native grassland to mature plantations were the major drivers of the decreasing taxonomic and functional diversity of plants. To our knowledge, this is the first time to simultaneously compare plant taxonomic and functional diversities and their relationship with environmental variables providing new insights into the effects of plantation cycle on grassland plant communities. Furthermore, we propose a model that includes soil (pH) and structural variables (litter depth, canopy cover) to explain the variation in plant taxonomic and functional diversity along the eucalypt forest cycle. Eucalypt tree growth imposed a novel abiotic conditions to pampean grasslands and could modify resource availability to plant species (light, water, and nutrients) 31 .
As expected, species richness and alpha functional diversity decreased over plantation age, and it was explained by soil acidification, increasing leaf litter depth, nitrogen content in soil and increasing canopy cover along the eucalypt forest cycle. Decreasing in understory plant richness in the first four years after eucalypt plantation establishment was found within the Danling region in China where the landscape mosaic is composed of commercial plantations, cultivated land and unmanaged forests 32 . Contrary to our results, the richness and abundance of plant communities was significantly greater in mature plantations of Robinia pseudoacacia (nonnative tree) than in the native tree plantations in the Loess Plateau in China 13 . The management practices that characterized the tree plantations caused different levels of disturbance affecting the plant community composition during the secondary succession. In eucalypt plantations, the accumulation of leaf litter has been considered a factor of allelopathy as a response to competition, through a reduction in seed germination and light availability for understory species seedling 33 . Certain phenolic acids and volatile oils are released from the leaves, bark and roots of some Eucalyptus spp. act as allelopathic agents and are harmful to other plant species 34 . Jobbágy et al. 35 found that eucalypt plantations in grasslands acidify the soil mainly as a result of the release of phenolic acids. Soil acidification can also change the bioavailability of nutrients, being beneficial for certain types of species,  www.nature.com/scientificreports/  www.nature.com/scientificreports/ as well as detrimental to others 36 . On the other hand, canopy cover was a relevant factor in decreasing plant functional diversity and consequently may impact on functions that plant species performed. Furthermore, light availability on the forest floor, microclimatic conditions (i.e. temperature and humidity) and rainfall patterns are closely dependent on canopy cover 37,38 . The forest canopy generate a rainfall partitioning into interception, throughfall and stemflow and then affect soil moisture patterns, infiltration, groundwater recharge and water yield 39 . Considering that high nitrogen availability in the soil could stimulate microbial activity and then an increase of decomposition rate 40 , the positive association between plant richness and nitrogen content in soil would be related to the high microbial activity allowing that certain soil nutrients are available. Therefore, the eucalypt forest cycle modified habitat characteristics limiting the establishment of some plant species and their functional traits as their ecological niches were not more compatible with the new environmental conditions. Regarding taxonomic and functional similarity results, canopy cover, leaf litter depth, soil pH, and nitrogen content in soil were the most important factors which negatively affect grassland plant species. As we have reported previously, leaf litter could produce soil acidification and then allelopathic effects on plant species and would avoid that typical grassland species could remain. The presence of some plant species in eucalypt plantation and their absence in grassland sites (Appendix A) could be explained by the dispersion of some species from riparian forests of the Uruguay River 41 to mature plantations. For example, we have recorded the presence of zoochory species like Blepharocalyx salicifolius, Butia yatay, Teucrium vesicariumand Eugenia myrcianthes in mature plantations, whose dispersion would be favored by high similarity in habitat structure and microclimatic conditions with riparian forest 8 . Besides, studies have reported that tree plantations are hot-spots of plant invasion and threaten the remnants of semi-natural vegetation 42,43 . In addition, studies that evaluated changes in plant species composition developed in landscapes dominated by forest have shown that, the species composition changes considerably due to the increasing of invasive species in commercial plantation [44][45][46][47][48] . On the other hand, other important factors that could influence the variability in functional diversity are management practices carried out in eucalypt plantations of pampean grasslands. Particularly forest management that invokes increased light availability through the opening spaces, for instance thinning, may effectively increase overall plant species and so biodiversity at different scales 42,[49][50][51] . Nevertheless, the authors have remarked that this increase would be consistent with the increase of invasive plant species within tree plantations, then the original native ecosystem functions may have been damaged. Our findings suggested that thinning (in 2-year and 7-year eucalypt plantation) not increase plant diversity due to all diversity indexes showed a marked decrease in their values.
The proposed model explaining the effect of the eucalypt forest cycle on alpha functional diversity indicated that canopy cover had an indirect relationship mediated by leaf litter, soil pH, and species richness. Furthermore, the model has shown that the species richness decline was linked with the loss of unique functional traits suggesting low functional redundancy in grassland plant communities of Mesopotamic Pampa. The increase of leaf litter depth had a direct and indirect relationship with species richness would be due likely to the soil acidification and physical restrictions for plant growth. Based on the coefficients of the model, the impact caused by the physical restriction on species richness (and functional diversity) would be higher than the soil acidification and the allelopathic effects. Future studies may incorporate other abiotic factors such as nitrogen in the soil available to plants and humidity, to disentangle the effects of environmental changes on functional diversity.
The associations between life-history traits and environmental variables showed that high dispersal potential, annual, barochorous and zoochorous species were the functional traits further affected by the growth of eucalypt plantations. Furthermore, leaf litter, canopy cover, pH and nitrogen content in soil were the environmental variables with a large impact on life-history traits that could act as ecological filters and promoting that, only species with certain combinations of life traits, may pass through the filters. In addition, species that occurred in young plantations had similar characteristics to fast-growing pioneer species common in early-successional stages: high amount of seed production, annual lifeform that support our hypothesis. This concurs with the general pattern that greater dispersal ability allows quicker response following a disturbance 52 . Moreover, the increase of leaf litter depth could be a filter for grassland species which are characterized by high dispersal potential. On the other hand, zoochory and barozoochory were the dispersal syndromes affected by canopy cover and leaf litter and may impact on ecosystem services. However, we expected that anemochorous and anemozoochorus species would have been presented in young plantations, not only for their high dispersal abilities, but also as they are characteristics of grassland and can tolerate their environmental conditions. Our results are consistent with other studies that found plant traits positively associated with disturbance intensity 53 . Assemblages after disturbance comprised few plants with wind-dispersed seed, consistent with selection for species with better dispersal ability. The observed trait composition could be the response to high disturbance intensity (i.e. formicide and glyphosate application) in the early years of eucalypt plantation. This set of traits was more likely to colonize new environments and have better abilities to find the necessary resources (Duncan et al. 2003). Changes in trait composition towards generalist species in response to land-use intensification have also been reported for other traits and arthropod taxa in grasslands of Germany 54,55 . At the end of the eucalypt forest cycle, the plant species presented traits related to late-successional species (less dispersal abilities, perennial lifeform) and different dispersal syndrome with respect to younger plantations. It is important to notice that successional processes are driving plant community changes, not only the environmental changes caused by the growth of eucalypt plantations, and our study design cannot separate both effects.
Assessing the changes in taxonomic and functional diversity of plant species along the eucalypt forest cycle highlighted the importance of integrating approaches. Our results demonstrated that facets and components of native grassland plant diversity were negatively affected by novel environmental conditions provided by eucalypt plantation. Leaf litter depth was the most important driver in plant diversity decline may cause loss of native species in pampean grasslands and consequently could affect resource availability for herbivores, species interactions, trophic networks and ecosystem functioning. Our findings could be useful in places where eucalypt plantation is expanding within open habitats, especially in South America and Africa. We recommend to regional stakeholders www.nature.com/scientificreports/ that conservation efforts should be addressed to reduce litter accumulation and/or creating heterogeneity in its distribution within eucalypt plantations. Finally, our approaches and results may contribute to develop practical management guidelines aimed at enhancing the value of plantations for biodiversity conservation.

Material and methods
Study area and forest management. The study was conducted in the Entre Ríos Province that includes the northern portion of the Pampa region (the Mesopotamic Pampa subregion), Argentina (Fig. 5). Climate is temperate humid, with an average annual temperature of 18 °C and annual average precipitation ranging from 1000 to 1200 mm 56 . The study area presents a rolling relief, well-defined streams bordered by riparian forests (Cabrera 1971;Rodriguez et al. 2017Rodriguez et al. , 2018. The natural dominant vegetation type was a grassland composed of Axonopus sp., Paspalum sp., Digitaria sp., Schizachyrium sp., and Bothriochloa sp. (Bilenca and Miñarro 2004). The general history of anthropogenic land use in the region began with an historic extensive cattle ranching since the nineteenth century, followed by soybean cultivation and citrus crops in the twentieth century 60 . In the last 50 years, the expansion of eucalypt plantations in the Mesopotamic Pampa is replacing mainly cattle pastures, and native grassland fragments. The Pampas region in South America is one of the greatest temperate grassland of the world and harbored one of the fastest growth planted forests in the American continent 59 . Eucalypt grandis W.Hill ex Maiden is the primary species planted in the region, promoted by high biomass yield, public incentives, and an emerging market of carbon sequestration 9,35 .
Specific Eucalytus plantations were selected with the objective of reducing the error associated with environmental variables and management practices that do not help to focus on answering the questions posed in the objectives and hypotheses. Stand selection: only plantations whose origin were natural grasslands were selected, avoiding those coming from abandoned citrus orchard or agricultural cultivation. From those, stands located in sandy loam soils and with similar site quality were filtered and finally selected. About management practices, the main criteria used in this study were, even-aged plantation; 4 × 2.5 m plantation frame which is equivalent to stocking densities averaging 1000 trees/ha 61 . The use of herbicide before planting was carried out in the planting line to reduce root competition between eucalyptus and natural grasses and herbs. In addition, formicide www.nature.com/scientificreports/ application was implemented to avoid herbivory from leaf-cutting ants 62,63 . Both practices were applied at 0-1 year after planting. Two pruning practice: the first one during the 3rd year, clearing all lateral branches up to 2 m height, and the second during the 6th year up to 4 m height 61 . Two thinning practices are usually conducted at 3 and 7-8 years after plantation; the final forest density is usually around 400-600 trees/ha in plantations prior to harvest 64 . Eucalypt plantations, in this area, were harvested for pulp and timber between 10 and 12 years, after that, harvest residues (leaves, bark, twigs, and branches) are not removed from stands.
Study design, plant sampling and environmental variables. Eucalypt plantations of nine ages (0, 1, 2, 3, 4, 5, 7, 9, 12 years) were selected to represent the chronosequence; and natural grassland sites to use as reference habitat. The zero-year-old plantations referred to the period after planting but less than one year. There were three replicates per plantation age and grasslands for a total of 30 sites. At each site, we established three 16-m 2 plots (total 300 plots) at least 50 m apart from plantation borders to avoid edge effects. During the blooming season from the late spring to early summer (December to January), plant species were identified at each plot, assigned percentage cover and measured environmental and dasometric variables related with the forest cycle changes. The collected specimens were identified to species or morphospecies level consulting the main territorial floras 65,66 . The study complies with relevant institutional, national, and international legislation.
Moreover, canopy cover of eucalypt trees was estimated using five digital photos taken from 1.5 m above the ground toward the canopy at each site. The percent canopy on each photo was analyzed with ImageJ as percentage of pixels with vegetation 67 . The average of leaf litter depth (mostly composed by eucalypt leaves) was estimated by measuring the 10 random points per site. Similarly, average of tree height and diameter at breast height (DBH) were measured using 10 eucalypt trees. Temperature was measured with miniature temperature logger (i-button) within one replicate per age and grassland site every five min for the entire 15 days collection period. All temperature loggers were in direct contact with soil placed 10 cm above the ground level. Temperature is one of the primary factors explaining both species richness and abundances of plant likely through physiological constraints and resource availability 68 . Finally, 10 soil subsamples at 0-20 cm depth were taken by walking a 'zig-zag' pattern from each site and mixed for homogenization. From each soil sample, total soil carbon, nitrogen, phosphorus (Bray& Kurtz), and soil pH were determined in the laboratory.
Taxonomic diversity. To estimate alpha taxonomic diversity, we counted plant species richness per site.
Moreover, to estimate beta taxonomic diversity, we calculated the Sorensen dissimilarity index (β sor ) between each plantation age, and the native grassland pooled the three grassland sites to better represent the regional pool of species. The results of beta diversity are presented in terms of similarity in composition between communities (1-β SOR ).
Selection of functional traits. Plant species that occurred only once were removed from total species recorded. The main reason for this is that, dominant plant traits, are known to strongly influence ecosystem processes, indicating how the community in general responds to the environment 69 . In general, plant communities have a typical structure with a relatively small number of dominant species and a large number of subordinate or minor species that account for a low proportion of the biomass 70,71 . Plant species recorded on the plots were scored for four functional traits associated with community responses to disturbances: growth form, dispersal syndromes, dispersal potential and life history (see Table 2) 72 .

Functional diversity indices.
We calculated functional dispersion (FDis) based on previous functional traits selected using the FD package 73 . FDis is defined as the mean distance in a multidimensional trait space of individual species to the centroid of all species, weighted by their relative abundances 73 . Thus, functional dispersion describes one of the facets of functional diversity, showing how far the most abundant species move away from the center of the functional space. FDis is unaffected by species richness by construction 73 .
To estimate beta functional diversity along eucalypt chronosequences (considering grassland sites as the reference habitat), we used trait matrix to build a dendrogram tree using the Jaccard index because functional traits are categorical variables and to avoid that not sharing traits increase the similarity between two species 74 . Then, both dendrogram tree and presence/absence data per site were used to calculate the Fsor index using the PICANTE package 75 , which is analogous to a traditional Sorensen's Index. Fsor represents the proportion of dendrogram branch lengths shared by two assemblages 76 . To test changes in species richness (alpha taxonomic diversity) along the eucalypt forest cycle, we performed Generalized Linear Models (GLM) with Negative Binomial distribution because we had a count variable (species richness) and to reduce overdispersion. To assess changes in alpha functional diversity along the eucalypt forest cycle, we modeled FDis with Normal Distribution.
We modeled taxonomic (1 − β SOR ) and functional similarity (Fsor) with beta distribution and performed beta regression models to explain the variability along the eucalypt forest cycle. Beta regression models are based on the assumption that the response variable is beta-distributed and it was developed for situations where, the dependent variable, can take values between 0 and 1 as the case of the Sorensen index 77 .
To identify the main environmental variables that explain changes of taxonomic diversity along the eucalypt forest cycle, we used a model-selection approach with Akaike Criterion (AICc), as we handled the corrected version for small sample sizes 78 . Models with ∆AICc < 2 were considered as equivalent to the minimum AICc model and hence more robust to explain the observed variability 79 . Goodness-of-fit was evaluated by examining plots of standardized residuals versus predicted and checked normal distribution. Environmental similarity indices (1-Gower index), between each environmental variable respect to grassland sites, were calculated for beta taxonomic analyses. The models with ΔAICc > 2 have been present in Appendix B, C, D, and E.
As noted previously for taxonomic diversity analyses, the model selection approach was implemented to identify the environmental variables that explained the variability in alpha functional diversity and functional similarity along the eucalypt plantation cycle.
Structural equation models (SEMs) are a statistical technique that unite multiple variables in a single causal network, thereby allowing simultaneous tests of multiples hypotheses 80 . SEMs test the direct and indirect effects on pre-assumed causal relationships 81 . Piecewise SEM expands upon traditional SEMs by introducing a flexible mathematical framework that can incorporate a non-normal distributions, hierarchical structures and different estimation procedures 82 . We then performed a piecewise SEM to analyze the direct and indirect effects on the alfa functional diversity (FDis), including the environmental variables with the lowest AICc from the model selection and the plant species richness. Some researchers recommend a minimum sample size of five cases per free parameter in the model resulting in six possible relationships to be tested 81 . The direct relationship between species richness and alpha functional diversity was modeled with Poisson distribution. Fisher's C statistic was calculated to evaluate the global goodness-of-fit of the model (P > 0.05) using the 'psem' function from the 'piecewiseSEM' package 83 .
In order to identify associations between functional traits and environmental variables, we used fourth-corner analysis, i.e. measured through the interaction terms between each environmental variable and trait variables 84 . The sign and magnitude of the interaction coefficients denoted the nature and strength of the association of the trait-environment relationship, respectively. We fitted a linear regression model, using all functional traits and environmental variables and the LASSO penalty was estimated via cross-validation. LASSO is a method of penalized likelihood which imposes a constraint on estimates of model parameters 85 .
All statistical analyses were conducted in R 3.2.0. Beta regression was performed using functions in the R package 'betareg' 77 , and candidate models were selected with the R package 'MuMIn' 86 .

Data availability
The presence of plant species in eucalypt plantations and grasslands sites are in the Supplementary Information. More detail information of this study is available from the corresponding author on reasonable request.