Exploring Assembly Trajectories of Abandoned Grasslands in Response to 10 Years of Mowing in Sub-Mediterranean Context

: Abandoned semi-natural grasslands are characterized by lower plant diversity as a con-sequence of tall grasses spreading. Mowing is a widely used restoration practice, but its effects to maintain the restored diversity over time is poorly investigated in sub-Mediterranean grasslands. Since 2010 in the central Apennines, we fenced a grassland, invaded by Brachypodium rupestre , which was mowed twice a year. Before the experiment started, we recorded species cover in 30 random sampling units (0.5 m × 0.5 m). The sampling was repeated every two years for a total of ten years. We used linear mixed-effect models to investigate the trajectory of functional diversity and community weighted mean for traits related to space occupation, resource exploitation, temporal niche exploitation, and Grime’s CSR strategies. The reduction of the weaker competitor exclusion exerted by B. rupestre affected the functional plant community. In the short term (4–6 years), this fostered space occupation strategies, decreasing convergence of clonal strategies and horizontal space occupation types. In the longer term (8–10 years), mowing ﬁltered ruderal strategies, i.e., species with faster resource acquisition (lower leaf dry matter content, LDMC). LDMC and CSR strategies, initially convergent due to the dominance of B. rupestre , lowered convergence over time due to higher differentiation of strategies.


Introduction
Semi-natural grasslands are an essential part of the European cultural landscape, being hotspots of biodiversity resulting from centuries or millennia of land use [1]. Grassland biodiversity is endangered by different drivers mainly related to processes of intensification or abandonment [2,3], with the latter particularly relevant in Mediterranean mountains [2,4,5]. Here, after World War II, the cessation of traditional management practices led to the abandonment of large areas [6,7], which are left to spontaneous vegetation regrowth, such as the expansion of shrublands or forest ecosystems at large scales [5,6] and invasion of competitive tall grasses at community scales [8][9][10][11][12].
In the Italian peninsula, two perennial species of Brachypodium are fostered by land-use cessation: Brachypodium rupestre (Host) Roem. & Schult. in low mountains, and B. genuense (DC.) Roem. & Schult. over 1300-1400 m a.s.l., both having competitive and stress-tolerant strategies. The competitive success of these Brachypodium species is related to their high tiller density and branching frequency [13], as well as to their capacity for clonal growth and clonal integration strategy [14]. Their invasive behavior causes a change in community structure and species composition [8,10,15,16]. Overall, the increase of tall-grass species leads to a marked drop of species richness and species diversity [10,17,18], by a process called "weaker-competitor exclusion" [19]. This process suggests that the prevalence of Land 2021, 10,1158 3 of 20 the decrease of B. rupestre would lead to an increase of species diversity as evidenced by previous studies [8,10,27,31], by reducing the dominance of this species, we hypothesized that the reduction of the weaker competitor exclusion exerted by B. rupestre as a result of mowing twice a year would affect the functional plant community. In the short term, this would foster space occupation strategies to maximize the exploitation of the newly available resources, increasing divergence of clonal strategies, plant height, and horizontal space occupation types. We hypothesized also that in the longer term, mowing would shape a highly adapted plant community, positively filtering ruderal strategies, i.e., species with higher specific leaf area and lower leaf dry matter content and faster resource acquisition (a disturbance avoidance strategy), and also species that allocate resources in belowground clonal organs, which would ensure tolerance to mowing. We expected that an initial high convergence of leaf traits and CSR strategies due to the dominance of B. rupestre would have paved the way to higher differentiation of strategies in the community and, thus, a lowered convergence over time.
As far as leaf and flowering phenology are concerned, B. rupestre-dominated communities are characterized by temporal niche differentiation, which ensures species coexistence [44,45]. We expected that the decrease of B. rupestre led to the re-occupation of the most favorable period for flowering and photosynthesizing, i.e., in early summer at the maximum of phytomass production [46], reducing functional diversity (i.e., increased convergence). Specifically, our questions were: (I) Does mowing foster in the short term the divergence of spatial niche occupation strategies? (II) Does mowing reduce in the long term the convergence of temporal niche occupation strategies? (III) Does mowing favor the divergence of CSR strategies, leading to more ruderal communities, and how does this change reflect on leaf traits?

Study Area
The study area is located in the central Apennines (337398.00 m E, 4757200.00 m N, UTM-WGS84), encompassing north-facing slopes on limestone bedrocks. The area ranges from 1000 to 1300 m a.s.l. The climate is sub-Mediterranean [47], characterized by winter cold stress and summer drought stress. Mean annual rainfall is 1100 mm and the average annual temperature is 10-11 • C. Soils are neutral-sub acidic and their depth ranges from 20 to 40 cm. Its texture is mainly characterized by sand (8-10%), silt (32-35%), and clay (50-60%). Skeleton ranges from 25 to 40%, whereas total nitrogen averages from 6 to 7 g Kg −1 [48].
Grasslands The area has had a century-old history of grazing by sheep and cattle, mainly during the period between late spring and late autumn, but in the last decades the number of grazing animals strongly decreased. Besides grazing, mowing has been carried out during the period of early summer on semi-flat slopes [48].

Sampling Design
In 2010, we selected the following elements as stratifying criteria to reduce macroenvironmental variables among sampling units: topography features (altitude, ranging from 1000 to 1300 m a.s.l.; aspect, from west-north-west to east-north-east; slope, 20 • to 50 • ). Using remote sensing images and the QGIS software, the areas matching the stratifying criteria were overlaid with a grid of 100 m × 100 m and we selected those cells with a homogeneous aspect, slope angle, and sward cover. Following these criteria, we chose two adjacent cells of 1 ha each, with altitude equal to 1220-1280 m a.s.l., aspect north-westfacing and slope angle of 25 • . In order to avoid disturbances by domestic and wild animals the cells were fenced.
As management practice, we selected two mowing events a year as previous research found the highest vascular plant species richness when the grassland was mown twice a year and hay was removed [49]. Moreover, experimental research suggests the adequacy of summer mowing to avoid reallocation of nutrients to storage organs [50].
One of the two selected cells was mown twice a year (i.e., end of June, before B. rupestre flowering, and at the end of October) starting from 2010 until 2019. Aboveground phytomass was removed and the sward was cut 2 cm above soil level using a grass trimmer and the mown phytomass and litter were hand removed every mowing event (i.e., twice a year). The other adjacent cell was left unmown and used as a control. Inside the fenced cells, using a QGIS generator of random points, a 5 m x 5 m grid was overlaid, and 30 points present at the crosses of the grid were randomly selected in each cell. Points close to the fence (<5 m) were excluded to avoid edge effects. Considering the findings of previous studies on community assembly rules [51], we assumed that fine-scale spatial resolution (0.5 m × 0.5 m) offered the most appropriate scale of study to detect the effect of biotic interactions. The 30 selected points were the lower left-hand corner of 30 sampling units (0.5 m × 0.5 m) per cell. For each year of data collection (2010,2012,2014,2016,2018, and 2020), 30 new sampling units were randomly placed to avoid temporal autocorrelation. In the unmown area, data were collected in 2010, before the start of the experiment, and in 2020, i.e., ten years later.

Data Collection
Relevés were all carried out at the end of June, before mowing. In each sampling point, we recorded the cover of each vascular species (visually estimated in percentage scale).
To understand the functional trajectory of the plant community induced by mowing twice a year, we selected a core of plant traits capturing different functional dimensions: competition for light and space occupation (plant vegetative height, H; horizontal architecture, HA) [52,53]; resource exploitation (specific leaf area, SLA; leaf dry matter content, LDMC) [53] and regulation of leaf temperature and water-use efficiency during photosynthesis (leaf area, LA) [53]; temporal niche exploitation (flowering phenology, FP; leaf phenology, LP) [27,52]; clonal strategy (clonal position, CP; clonal spread, CS) (adapted from [54]). Moreover, leaf traits (LA, SLA, and LDMC) were also used to calculate plant CSR (competitive, stress tolerant, and ruderal) strategies [39], by integrating both leaf economics (typical of S and R strategy selection) and leaf size (typical of the C strategy) variation spectra in a three-way trade-off [36]. These selected traits are related to both competitive ability and response to mowing [52][53][54][55]. A description of each trait and their respective list of states and data sources are reported in Table 1. Data on SLA, LA, LDMC, and H were drawn from field measurements at the end of June 2019 in the surrounding areas sharing similar topographic conditions, following the protocols described by Perez-Harguindeguy et al. [56], i.e., taking measurements from 10 fully expanded leaves without signs of damage of different individuals per species. In total, we measured 860 individuals of 86 species. Even if we did not consider intraspecific variation, the approach adopted, which considered a mean trait value fixed for each species, is more advantageous compared to the use of a regional or global trait database since they rarely provide climatic or habitat information on trait collection sites [57]. Since measuring these quantitative traits for all the species was not feasible, to perform the following statistical analyses we considered the species whose relative cumulative plot cover (in decreasing order of abundance) reached at least 80% of the total plot cover [58]. This choice relies on the assumption that the ecosystem functioning and services depend strongly on the dominant species (biomass-ratio hypothesis) [59]. These trait values were then log 10transformed before running the analysis [58].
On the contrary, other traits were categorized by means of fuzzy coding in which species show multiple states of a trait (e.g., flowering time, Table 1). The percentage value of a species (0-100%) for each state was estimated dividing 1 by the number of the states. For example, species that can flower in both spring and summer would have a score of 0.5 for spring flowering and 0.5 as summer flowering [60]. Furthermore, Grime CSR strategies, calculated with the StrateFy CSR classification tool [36], have been treated as a fuzzy coded trait in which the C-, S-, and R-scores are expressed as a percentage (0-100%). Categoric traits were obtained through databases for all the species: LEDA for flowering time, leaf phenology, horizontal architecture [61]; CLO-PLA for clonal position and clonal spread [54] and supplemented by observations in the Herbarium Universitatis Camerinensis (CAME).

Statistical Analysis
Preliminarily, we calculated the mean and standard error of the mean for the number of species, the total vegetation cover (sum of species cover values), and the B. rupestre absolute and relative cover, for each year.
Then, to investigate the functional trajectory of the plant community along the experimental trail, we computed the functional diversity (hereafter, FD) with Rao's Quadratic Entropy [62]. We selected Rao's Q because it is the most widely used index for functional investigation [62,63] and it expresses the expected dissimilarity between two individuals Land 2021, 10, 1158 6 of 20 of a given species assemblage selected at random with replacement. We calculated the FD for each plant trait as follows: where S is the number of species, p i and p j are the proportions of i-th or j-th species in the sampling unit, and d ij is the functional distance or dissimilarity between the i-th and j-th species ranging from 0 (trait value of species i and j is exactly the same) to 1 (when two individuals have completely different trait values). As a measure of distance (d ij ) we used the Gower distance [64], implemented by de Bello et al. [60], since this is the appropriate distance for functional analysis being able to handle fuzzy coded traits [60]. Since detecting assembly rules needs to compare the observed values with those expected by chance [65], we calculated for each FD a standardized effect size (SES), as (Iobs − Isim)/σ sim , where I obs is the observed value of FD, I sim is the mean of the expected FD, and σ is the standard deviation of expected FD. A distribution of SES values lower than 0 indicates functional convergence, while a distribution of SES above 0 indicates functional divergence, according to a two-tailed t-test [66]. Expected FD values were generated, shuffling 999 times trait values among the species in the entire data set ("between-plot randomization strategy") [65]. This algorithm, assuming a null model in which each species may occur in each plot with any abundances, and in combination with Rao's Quadratic Entropy, is sensitive to habitat filtering and limiting similarities, and is less affected by Type I error rates [65]. Since we did not expect any clear relationship (e.g., linear or non-linear) between response variables (i.e., SES-FD) and the predictor [67] (year), we used generalized additive mixed models (GAMMs; with a Gaussian distribution) to analyze the trajectory of the standardized effect size of functional diversity (SES-FD) for each plant trait in relation to the number of years of mowing (fixed effect). In the model, we set the knot parameter (k) to 4 to prevent overfitting. To account for the random distribution of the sampling units, we included as random intercept the identities of samples (each sample being composed of 30 sampling units, as one different random sample was taken every two years).
In order to have a better understanding of the meaning of FD variations, we computed the community weighted mean values (hereafter, CWM). This index is widely used to quantify the shifting in the functional strategy along environmental gradients and provides complementary information to FD [63]. We calculated CWMs as follows: where CWM is the community weighted mean value of a given plant trait (for quantitative traits) or relative contribution of each state (for categorical traits), S is the number of species, p i is the relative abundance of species I (i = 1,2, . . . , S), and x i is the trait value for species i.
In the case of the CWMs, using observed values may lead to biased biological information due to the spurious correlation with the species composition matrix [68] resulting in higher inflated error Type I rates [69]. For this reason, we integrated generalized additive mixed models with the "Max Test" procedure which involved two independent permutation tests: one testing the species attributes-species composition link ("column-based permutation") in which 999 times the species trait across species are shuffled randomly, and the other testing the plot attributes-species composition link ("row-based permutation") in which 999 plot attributes are shuffled randomly across all plots [69]. A p-value for each permutation procedure (p i , row-based permutation; p c , column-based permutation) is estimated separately as the probability that the marginal adjusted R 2 value of the observed model is significantly higher than a distribution of 999 expected marginal adjusted R 2 values, according to a one-tailed t-test (α < 0.05). By picking the largest of the two p-values as results (p max = max (p i , p c )), the type I error is controlled.
To exclude that the results that were related to factors unrelated to the treatment (e.g., climate variability), we performed a BACI (Before/After Control/Impact) analysis [70], where sampling is conducted at simultaneous time periods in treatment and control sites Land 2021, 10, 1158 7 of 20 before and after treatment (2010 and 2020) [71]. We used mixed-effects models, modeling the two-level categorical variables Control/Impact (unmown vs. mown area) and Before (2010)/After (2020) treatment and their interaction term as fixed effects. A significant interaction term indicates a significant effect of the treatment. As sampling units were spatially nested in two locations and changed every year, we used the identities of samples (each sample being composed of 30 sampling units), as one different random sample was taken every year, before and after treatment nested inside the two sampled areas as random intercept. We ran the mixed-effects models for each of the SES-FD indices using the lme function of nlme R-package. Significance values were obtained using the Type II Wald Chi-square test and the Anova function (car R-package, version 2.1-5). For CWMs, we used the above mentioned approach (i.e., Max Test) (lme function, nlme R-package).
All statistical analyses were carried out using the R software, version 4.0.3 (R Foundation for Statistical Computing, Vienna, Austria, http://www.R-project.org). For the calculation of the FD index, we used the Rao function provided by de Bello et al. [58] and for the CWM, the functcomp function (FD package). To calculate Gower distance handling fuzzy coded traits, we used the gawsid function (gawdis package). To run the generalized additive mixed models, we used the gamm4 function (gamm4 package). In Tables S1 and S2 we reported the values of SES-FD and observed CWM calculated at plot level.

Results
In total, we found 138 species across all the years of the treatment in the mown area with mean values of 17 ± 0.54 (mean ± standard error mean) species per plot before the treatment, a peak after six years of 27 ± 0.59 and a slight decrease at the end of the experimental trial with 25 ± 0.79 species (Figure 1a). On the contrary, the total plot cover decreased linearly from 151 ± 5.6% at the beginning of the experiments to 113 ± 3.9% at the end of the manipulation (Figure 1b). Brachypodium rupestre cover decreased consistently for both absolute and relative cover (Figure 1c,d), with a proportional cover ranging from 34 ± 2.6% before the experiment to 4 ± 0.8% at the end of the experiment.

Functional Diversity Variation
Regarding the analysis of functional trajectories, we did not find any significant change, for the SES-FD of plant height (H), specific leaf area (SLA), leaf dry matter content (LDMC), leaf phenology (LP), and clonal spread (CS) ( Table 2). Leaf area (LA), horizontal architecture (HA), flowering period (FP), and clonal propagation (CP) changed significantly during the experiment (Table 2). In detail, for space occupation (Figure 2), HA saw a cubic trend (edf = 2.6; Adj. R 2 = 0.12, AIC = 522.0, Table 2 Figure 5a). Finally, also the Grime's CSR strategies ( Figure 6) showed significant changes in a quadratic trend (edf = 1.8; Adj. R 2 = 0.01, AIC = 352.5, Table 2, Figure 6a). Detailed results for all plant traits (edf, Adj. R 2 , and model AIC values) are reported in Table 2. Finally, the analysis on the distribution of SES-FD pointed out mainly a functional convergence pattern for all traits across all years (SES-FD < 0). A shift from a convergence pattern to a random pattern occurred for HA and CP from 4 to 10 years (Figures 2c and 5a), and in the last year of LDMC (Figure 3c). BACI analysis confirmed a significant effect of the treatment on SES-FD, for CP (Adj. R 2 = 0.35; p < 0.05), HA (Adj. R 2 = 0.58; p < 0.05), and Grime's CSR strategy (Adj. R 2 = 0.35; p < 0.05). with mean values of 17 ± 0.54 (mean ± standard error mean) species per plot before the treatment, a peak after six years of 27 ± 0.59 and a slight decrease at the end of the experimental trial with 25 ± 0.79 species (Figure 1a). On the contrary, the total plot cover decreased linearly from 151 ± 5.6% at the beginning of the experiments to 113 ± 3.9% at the end of the manipulation (Figure 1b). Brachypodium rupestre cover decreased consistently for both absolute and relative cover (Figure 1c,d), with a proportional cover ranging from 34 ± 2.6% before the experiment to 4 ± 0.8% at the end of the experiment.
In Appendix A we report the graphs with non-log10-transformed CWM values of H, SLA, LA, and LDMC.

Discussion
The functional structure of the plant community was influenced differently by mowing over time (Figures 2-6). In the first 4-6 years, mowing influenced the plant community probably through the weakening of the dominance of Brachypodium rupestre, highlighted by the decrease of its average cover value (Figure 1c,d) and by the fading of its weaker competitor exclusion effect, as suggested by the increase of species richness (Figure 1a). Contrariwise, after 6 years, the species richness value showed a fluctuation rather than a further increase (Figure 1a). In this second phase, the plant community was no longer dominated by B. rupestre, whose cover values tended to decrease to lower values (Figure 1c,d), and several SES-FD curves showed a variation in their trend (Figures 2a, 3c, 4c and 5c). This variation in some cases was more marked in the last two years, giving rise to a clear inversion in the trend. This is consistent with Tardella et al. [27], who stated that the execution of mowing twice a year influences the functional structure of the plant community, leading to a condition strongly adapted to this type of disturbance.

Shorter-Term Variations
In the first phase of the experiment (up to 4-6 years), the plant community showed a variation of SES-FD for most of the considered traits. In particular, the decreasing deviations from the null model for plant height, horizontal space occupation, and leaf phenology types (i.e., decreasing convergence) seem to underline a rapid increase in the available spatial and temporal niches. Actually, the reduction of B. rupestre coverage likely exerted a positive effect on the appearance of gaps within the turf (spatial niche implementation) [10], while the aboveground phytomass and litter removal reduced the space occupied by the long-lasting (living or dead) leaves of the B. rupestre individuals, making available new temporal niches for exploitation of resources [52]. In particular, plant height was linked to competition for light, which is size-asymmetric, meaning that taller species (i.e., B. rupestre) can remove the shorter ones, leading to a plant community characterized by species functionally similar to the dominant ones [17] (Figure 2a,b). Our results underlined indeed a variation from a grassland characterized by a taller and less stratified canopy (year 0) to a grassland with a shorter canopy but characterized by more variable values, indicating probably a vertically stratified community (Figure 2a,b).
The invasion of B. rupestre forces the spring-flowering species in shelter niches where few individuals of each species remain [15], reducing the ecological spatial gaps necessary for species establishment [72,73]. The persistence of high amounts of aboveground phytomass of dominant species and of a litter layer fosters the coexistence of species that use different temporal niches for flowering [74,75]. The increase of convergence of flowering phenology (FP) soon after the beginning of phytomass and litter removal is consistent with this observation. We found indeed an increasing trend of mid-spring and early summer flowering species (Figure 4c). Such findings are in line with the "mid-domain hypothesis" [76], which suggests that a peak of flowering occurs when the site conditions are the most favorable.

Longer-Term Effects
In the second part of the experiment (4-6 to 8-10 years), the plant community showed different trends. Functional variation could be ascribed directly to the effect of mowing, with B. rupestre dramatically reduced (Figure 1c,d) and no longer able to exclude weaker competitors [17]. Interestingly, plant community adaptation to mowing seems to be driven mainly by specific traits. In this phase, indeed, some functional strategies such as HA and CP, whose SES-FD values were not significantly different from those expected by chance, seem to be driven by stochastic events (i.e., random assemblage) rather than by deterministic processes. Contrariwise, some traits linked to disturbance avoidance and tolerance strategies were largely affected [52,77]. In particular, the effect of mowing in shaping the plant community is well distinguished by the pattern of clonality. In fact, we observed a CWM increase for species having buried clonal organs and a decrease for species having above-soil clonal organs (Figure 5b). Belowground clonal organs (mostly epigeogenous rhizomes) may allow plants to minimize the loss of biomass. Such a trend is consistent with the results of investigations performed in the Mediterranean [52] and in the temperate climatic contexts [78], which asserted that short clonal organs are fundamental for plant resprouting ("tolerance mechanism"), as they allow a considerable storage of carbon [79]. Finally, the filtering effect of mowing is underlined by the more negative SES-FD values of plant height. Indeed, recurrent mowing tends to select a community characterized by species functionally similar in their size (i.e., short-sized species) (Figure 2a,b). This is consistent with the adaptation used by the plant to reduce the amount of aboveground phytomass removed [17,29].

Variation of Grime's CSR Strategies
During the experimental trial we observed significant trends of the Grime's strategies, with the decrease of stress-tolerant strategies and the increase of ruderal ones ( Figure 6). This trend is in line with the general theory that disturbance repeated over time favors species with an overall ruderal strategy, i.e., species with avoidance and/or tolerance mechanisms [38,39]. Specifically, ruderal species do not invest a large proportion of resources in the individual growth but in propagules (i.e., seeds or clonal organs) by which the population can regenerate in case of disturbance [39]. Contrariwise, stress-tolerant species invest mainly in the capacity to retain resources, being present in resource-poor environments [38,39]. The higher CWM of stress tolerators at the beginning of the experimental trial results mainly in the higher proportion of B. rupestre, characterized mainly by a stress-tolerant performance (54%) in terms of resource exploitation and, secondly, by ruderal (26%) and competitive (20%) ones. The higher leaf plasticity allows it to spread in environments that experience different levels of stress intensity [80]. The combination of stress-tolerant resource exploitation strategies and tall stature has been seen as a functional portfolio for invasive species [81]. The increase in abundance of species having ruderal strategies is consistent with the pattern of leaf traits. Indeed, we observed the increase over time of species having lower LDCM values (Figure 3b,d,f), indicating a tendency toward fast-growth strategies [56].

Conclusions
The ongoing loss of biodiversity due to land-use variation is a worldwide concern that makes the recovery of abandoned grasslands a compelling issue. Since reintroducing the traditional land use in the sub-Mediterranean region is no longer economically sustainable for farmers, alternative practices should be found. In our study, we investigated the effect of recurrent mowing as an alternative practice to support the recovery of functional structure of abandoned grasslands. Overall, we found that the functional dimension of the studied plant community was influenced differently over time. In the first 4-6 years, mowing influenced indirectly the plant community, mostly through the weakening of the dominant ability of Brachypodium rupestre. This phase is characterized by the increase of species having fast-growing strategies. These species may suffer the competition with B. rupestre in abandoned grasslands, and also because they share the same flowering niche. In the second phase of the experiments (6-8 years), longer-term effects of mowing select directly the functional dimension of the plant community in that it tended to acquire functional adaptations based on avoidance strategies (e.g., lower LDMC) or tolerance strategies (e.g., belowground clonal organs). These trends are in line with Grime's theory [39], as we observed that the plant community shifted towards ruderal strategies over time.
Finally, we have to note that concerning the latest years of the experimental trial we detected the start of a different trend for most of the considered traits. This is a hint that, after 10 years of recurrent mowing, the plant community might start to shift toward other different functional and ecological adaptations. This observation highlighted the importance of monitoring the plant community over the next years.