Nurse species and indirect facilitation through grazing drive plant community functional traits in tropical alpine peatlands.

Facilitation among plants mediated by grazers occurs when an unpalatable plant extends its protection against grazing to another plant. This type of indirect facilitation impacts species coexistence and ecosystem functioning in a large array of ecosystems worldwide. It has nonetheless generally been understudied so far in comparison with the role played by direct facilitation among plants. We aimed at providing original data on indirect facilitation at the community scale to determine the extent to which indirect facilitation mediated by grazers can shape plant communities. Such experimental data are expected to contribute to refining the conceptual framework on plant-plant-herbivore interactions in stressful environments. We set up a 2-year grazing exclusion experiment in tropical alpine peatlands in Bolivia. Those ecosystems depend entirely on a few, structuring cushion-forming plants (hereafter referred to as "nurse" species), in which associated plant communities develop. Fences have been set over two nurse species with different strategies to cope with grazing (direct vs. indirect defenses), which are expected to lead to different intensities of indirect facilitation for the associated communities. We collected functional traits which are known to vary according to grazing pressure (LDMC, leaf thickness, and maximum height), on both the nurse and their associated plant communities in grazed (and therefore indirect facilitation as well) and ungrazed conditions. We found that the effect of indirectly facilitated on the associated plant communities depended on the functional trait considered. Indirect facilitation decreased the effects of grazing on species relative abundance, mean LDMC, and the convergence of the maximum height distribution of the associated communities, but did not affect mean height or cover. The identity of the nurse species and grazing jointly affected the structure of the associated plant community through indirect facilitation. Our results together with the existing literature suggest that the "grazer-nurse-beneficiary" interaction module can be more complex than expected when evaluated in the field.


| INTRODUCTION
Species interactions are known to drive species coexistence, their relative abundance, and ultimately ecosystem properties, such as productivity and resilience to perturbations (Agrawal et al., 2007). Improving our knowledge of the drivers of species interactions could help both our fundamental understanding and our predictive ability of ecosystem responses to global changes (HilleRisLambers, Harsch, Ettinger, Ford, & Theobald, 2013;Soliveres, Smit, & Maestre, 2015).
For example, Graff et al. (2007) showed that the risk of a palatable plant being eaten decreased when it was near an unpalatable plant, resulting in indirect facilitation; this effect was, however, reduced at relatively high pressures (see also Smit et al., 2007). Several studies found nonlinear, complex patterns of species interactions along grazing gradients (see Smit, Rietkerk, & Wassen, 2009 for a review, Le Bagousse-Pinguet et al., 2012). Those studies were, however, restricted to the measurement of sapling performances (e.g., survival, biomass). Thus, the extent to which indirect facilitation through grazing drives the dynamics and organization of plants at the community and ecosystem scales remains largely unknown ).
An increasing number of studies in community ecology have stressed the advantage of using functional traits and strategies rather than species to get a more detailed understanding of the mechanisms generating observed communities (Mcgill, Enquist, Weiher, & Westoby, 2006;Violle et al., 2012). Functional traits are measurable features of an organism, which are linked to their fitness (Violle et al., 2007). The existing trade-offs between them allow defining life strategies (Grime, 1977) and enable to compare the effects of ecological mechanisms across gradients and scales. The goal of the trait-based comparative framework is to find general patterns, which allow the prediction of species interactions according to their functional traits.
In this study, our objective was to contribute to the understanding of indirect interactions among plants mediated by herbivores using a functional trait approach at the community scale. We hypothesize that the effects of grazing on associated plant communities (1) depend on the identity (and therefore on the traits) of the nurse and (2) are visible on the functional traits measured in the associated plant communities.
Taking the current literature into account, we aim at integrating our results into the current conceptual framework of indirect interactions and contribute to refine it.
We set up a grazing exclusion experiment in the tropical alpine peatlands of Bolivia. These ecosystems highly depend on a few, structuring cushion-forming plants which host associated plant communities (Cooper, Kaczynski, Slayback, & Yager, 2015;Cooper et al., 2010).
Cushion species are recognized as obligatory nurse species in these ecosystems, operating a transition from a mineral and aquatic environment to an organic and terrestrial environment and being a refuge for a number of endemic plant species (Loza Herrera et al., 2015;Ruthsatz, 2012;Squeo, Warner, Aravena, & Espinoza, 2006).
We investigated the effect of indirect facilitation on the community structure of the associated plant communities through the grazing protection provided by the nurse species. We set fences over two nurse species chosen for their contrasted strategies to cope with grazing: one, Distichia muscoides, is a compact, short-leaved cushion, whose shape limits the removal of biomass by grazers (indirect defense, sensu Boege & Marquis, 2005) but whose defenses seem unlikely to extend to the associated community because the associated communities develop above its canopy; the other, Oxychloe andina Phil., is a loose cushion with long spiny leaves (direct defense, sensu Boege & Marquis, 2005), which can possibly benefit the species living in the cushions. Because of these contrasted strategies within a same life form, we assume that D. muscoides does not provide protection to its associated communities whereas O. andina does, thereby leading to indirect facilitation from O. andina but not from D. muscoides.
As functional traits are known to vary with grazing pressure, we asked (1) if, as hypothesized, D. muscoides was less defended against grazing than O. andina by studying variations in the nurse traits and (2) if indirect facilitation could decrease the effect of grazing on the trait distribution of the associated plant communities by comparing the species compositions and the trait values of the associated communities found in the two nurse species.

| Study area and target ecosystem
The study area was located in the Cordillera Real, a mountain range of the Bolivian Andes located between Amazonia and the Altiplano, close to Lake Titicaca (highest peak: Mt Illimani, 6,462 m; Figure 1). In this region, peatlands are found between 4,000 and more than 5,000 m a.s.l. and are surrounded by drylands (Squeo et al., 2006). Being located in tropical alpine regions, the high Andean peatlands experience a dry, windy climate, with daily frost and intense solar radiation (Squeo et al., 2006) but long growing season and the absence of persisting snow cover (Anthelme & Dangles, 2012). The vegetation of the tropical alpine peatlands is dominated by cushion-forming species, especially Juncaceae and Cyperaceae (Ruthsatz, 2012). This type of life form has been extensively recognized as nurse for other plants in alpine ecosystems worldwide, through direct facilitation . In particular, cushion plants increase richness , β-diversity (Kikvidze et al., 2015), and the intensity of facilitation increases with phylogenetic distance , particularly where local diversity is low. The cushion plants slowly accumulate up to 10-12 m of organic matter mostly in valley bottoms along watercourses (Buttolph & Coppock, 2004;Cooper et al., 2015), covering the terrestrial surface in its entirety (i.e., no bare soil). The local population uses traditionally those ecosystems to sustain livestock, mainly camelids (Llama glama L. and Llama pacos L.; Buttolph & F I G U R E 1 (a) The study site (red circle) is located near La Paz and Titicaca lake in the Cordillera real (gray area, elevation >3,000 m); (b) a picture of an alpine peatland of the study site. Note the presence of fences; (c) the two cushion species studied; (d) schema of the experimental design showing the four treatments which were replicated 10 times. From the left to the right: Oxychloe andina control, Oxychloe andina under fence, Distichia muscoides control, and Distichia muscoides under a fence. Coppock, 2004) because they are far more productive than the surrounding dry vegetation all year long (Squeo et al., 2006).
Our study site was located in the Palcoco Valley (16°08′50″S, 68°17′08″W) in the vicinity of La Paz and the Titicaca lake, at an elevation range between 4,300 and 4,500 m a.s.l. In our study site, the amount of rainfall is around 410 mm in the humid season (December to March) and 184 mm during the dry season (i.e., the remaining part of the year) and the temperature reaches 6.4°C during the humid season and 4.5°C during the dry season (Loza Herrera et al., 2015). Those peatlands are dominated by two cushions species: Oxychloe andina and Distichia muscoides. These two species have different strategies against grazing. O. andina has spiny cylindrical leaves repelling grazers, while D. muscoides has tiny leaves and a rounded shape making it a highly compact cushion impeding the removal of biomass. Like in other tropical alpine peatlands, these cushion species allow other plant species, the so-called associated plant communities, to grow inside them through direct facilitation (Loza Herrera et al., 2015;Ruthsatz, 2012).
Because almost no plant species is found to grow outside the cushions in these ecosystems, the two cushion species are considered obligate nurses for the majority of the species composing the associated communities.
It is noteworthy that it is therefore not possible to evaluate the amount of direct facilitation provided by the cushion species, for example by comparing the associated communities with and without the presence of cushion species (e.g., Cavieres et al., 2014). Neither

| Sampling design
Our sampling design was implemented in two peatlands of the Palcoco Valley, separated by 1.5 km and taken as replicates. The spatial location of the plots was chosen randomly so that spatial variation in grazing intensity is unlikely to be a confounding effect. The sites were visited by livestock all year round (llama Llama glama L., alpaca Llama pacos L., and a few sheep), reaching approximately 300 individuals.
A grazing exclusion experiment was set up in February 2014 using metal fences (1.5 × 1.5 × 0.5 m) on top of 10 O. andina and 10 D. muscoides cushions (i.e., five cushions in each of the two peatland sites for each cushion species; Figure 1b). The fence structure was made of aluminum bars, and the five faces above the soil were covered with galvanized hexagonal mesh netting (5 cm long). The fences allowed excluding large and medium herbivores present in the sites, including wild ones such as Lagidium sp. Each vegetation plot had a size of 1 square meter (smaller than the diameter of a cushion); this allowed eliminating border effects within exclusion fences. The cushions covered all the surface of the plots. Each "grazing exclusion plot" (i.e., fence) was paired with a grazed plot of the same size located 1-2 m away from the fence border. In total, we monitored 20 ungrazed and 20 grazed (control) cushions, half of which from each of the two cushion species. Our protocol included a variable "cushion" nested within a variable "herbivory", each of which being represented by two treatments (O. andina vs. D. muscoides and herbivory vs. no herbivory) and each treatment having 10 replicates (i.e., five in each of the two sites).
A power analysis confirmed that our protocol, despite the relatively low number of replicates per treatment, was able to detect a reasonable significant effect of grazing exclusion (66 mg/g difference of Leaf Dry Matter Content (LDMC) with a power of 0.80 for example, Table   S1).
Because of its long spiny leaves, our hypothesis was that the protection of O. andina against grazing would extend to its associated communities, while the one of the D. muscoides would provide low or no protection. This hypothesis was tested (and confirmed, see

| Data collection
The field sampling took place 22 months after the experimental setup, that is, in December 2015. The paired plots had a similar species composition at the beginning of the experiment (time of the fence setting). Within each plot, we identified all the species of the associated communities and estimated their relative cover using a 10*10 cm mesh (see Garcia, Meneses, Naoki, & Anthelme, 2014, for details). We measured plant functional traits on the nurses themselves and on all the species of the associated community: LDMC, leaf thickness, and maximum height.
The LDMC characterizes plant resource acquisition strategies and is negatively correlated with maximum relative growth rate, because of a trade-off between resource acquisition speed and nutrient conservation (Diaz et al., 2004;Westoby et al., 2002). High LDMC plants are slow-growing species typical of stressed environments and abandoned pastoral lands (Cruz et al., 2010;Navarro et al., 2006;Peco et al., 2005) and high LDMC values are linked with the stress-tolerator syndrome. On the contrary, low LDMC is characteristic of competitors and ruderals with fast resource acquisition strategies. LDMC has also been linked to palatability (Louault et al., 2005). To get a reliable measure of LDMC, we used the 12 leaves of a given species collected in each plot. There is strong evidence that leaf mechanical properties play a role in deterring herbivores (Read & Stokes, 2006), notably the toughness. Leaf thickness accounts for a large part of the physical resistance of the leaves (Pérez-Harguindeguy et al., 2013). So, we used leaf thickness to assess the presence of a resistance shift in the grazing exclusion treatment. Among the 12 leaves used to measure LDMC, we measured leaf thickness on four. We performed a transverse section in the middle of the leaf length and measured leaf thickness at the middle distance between the midrib and the border of the leaf. The measurement was made using a microscope. Finally, we measured the maximum height. Although often described as a trait reflecting competition for light (Westoby, 1998), we took it also as the most direct indicator of grazing effects (see Díaz, Noy-Meir, & Cabido, 2001). We measured the height of four individuals for each species in each plot.
These individuals were the same as those used for the leave collection.
The same protocol was applied to the nurse species and the associated communities (31 species).
In each plot, we collected three leaves per individual on four individuals for each species (associated and nurse, 357 individuals sampled in total). For species with very small leaves, we collected the entire individual. Each sample was stored in hermetic plastic bags with a humid paper. Following field sampling, the samples were stored in a fridge until further analysis (within 1 week maximum). Our trait collection and measurements followed the guidelines of Pérez- Harguindeguy et al. (2013). The collected samples belonged to fully developed healthy individuals (without grazing damage or apparent disease). In particular, in the case of maximum height, we did not measure individuals that had been visibly grazed; this means that the observed differences in maximum height reflect a strategy modification rather than the direct impact of grazing.

| Data analysis
To confirm our hypothesis that O. andina is a better nurse than D. muscoides, we first analyzed the effect of grazing on the mean trait variation of the nurse itself, as well as on the percentage of overall cover of associated communities in the plots and on the relative cover of associated species in the cushions. The analysis of the species relative cover in the associated communities was performed using a Principal Component Analysis (PCA).
To investigate the effects of the structural defense of O. andina on the functional structure of the associated communities, we computed the community-weighted mean (CWM, Equation 1) and the community-weighted variance (CWV, Equation 2) for each of the three measured traits in each plot. The CWM is computed as the sum of the species mean trait in each of the plots weighted by the relative cover of the species in the given plots (Equation 1). This is a parameter of position of the trait distribution, which is known to vary along environmental gradients (e.g., Sonnier et al., 2010;Wright et al., 2004) and has been suggested to be an indicator of the local optimal trait value (Muscarella & Uriarte, 2016). The CWV is computed as the trait dispersion around the CWM weighted by the relative cover of the species (Equation 2). It is parameter of dispersion which is also known to vary along environmental gradients (e.g., Sonnier et al., 2010) and has been suggested to reflect trait convergence or divergence when it respectively decreases or increases (Sonnier et al., 2010). The second index, the CWV, is an indicator of the trait convergence or divergence in the community. To make the different plots comparable, we divided all the species relative covers by the sum of the relative cover of the associated communities in a given plot (i.e., the sum of the covers of all associated species in a plot is equal to 1). Note that the computation of the index uses the mean trait value by plot and by species (t ijk , i.e., taking into account intraspecific variability), unlike the typical method which takes one trait value per species (i.e., no intraspecific variability, Kattge et al., 2011). In other words, the method used here takes into account the intraspecific trait variation (Violle et al., 2012). The two indices were computed as follows (Sonnier et al., 2010;Violle et al., 2007): where a ik is the relative cover of species i in plot k, t ijk is the value of trait j of species i in plot k, and S is the number of species.
All the statistical analyses were performed using R (R Core Team, 2016) v.3.3.0. To compare the effect of the different treatments, we performed an analysis of variance (ANOVA). The variables grazing, cushion, and site were the independent variables. The grazing variable was nested in the cushion one because the fenced and the control plots were paired. Because the plots were paired, we used paired t tests as post hoc analysis instead of the more classical Tukey HSD. To compare the treatments between the two cushions, we used unpaired t tests. We corrected the obtained p-values with the Benjamini and Hochberg method, which controls for the false discovery rates (Benjamini & Hochberg, 1995), as do other similar post hoc procedures such Tukey HSD. Because the distribution of errors did not verify the normality assumptions, post hoc analyses were performed with Mann-Whitney nonparametric tests for CWM Height and for CWV (variance) of all the traits. (1)

| RESULTS
The treatments (grazed/not grazed, identity of the nurse) explained an important part of the variation in mean height of the nurse and in CWM Height of their associated communities (r 2 adj >70% for each, Tables S2 and S3), in mean LDMC of the nurse and in CWM LDMC of the associated communities (r 2 adj 30% for both), showing that the treatments captured well those trait variations. However, the treatments failed to explain variations in mean thickness (nurse), CWM Thickness , CWV Thickness, and CWV Height (r 2 adj = 0). Site had a significant effect on the CWM LDMC (F = 11.5, p = .002, Table S5). The cushion species had a significant effect on CWM Height (F = 89.2, p < .001), the CWV Height (F = 35.9, p < .001), and the CWM LDMC (F = 12.0, p = .001). The interaction between grazing status and cushion species had a significant effect on the CWM Height (F = 6.3, p = .005) and a marginal one on the CWM LDMC (F = 2.6, p = .089).

| Species richness and composition
The species richness of the associated communities was not significantly affected by the identity of the cushion species or the grazing status (from D. muscoides to O. andina: −1 species, p = .2; from grazed to ungrazed plots: +0 species, p > .7; Table S7). The difference of total cover of the associated communities between the ungrazed and the grazed plots was of the same order of magnitude in both cushions (O. andina: +13%, p = .004; D. muscoides: +15%, p = .002), showing that the associated communities of the two cushion species were affected by the grazing (Figure 2).
In total, 30 associated species were recorded in the study. A PCA performed on the species composition (normalized per plot) showed that the associated communities differed between the two cushions along the two first axes of the PCA (relative inertia: 13.70% and 10.82%, respectively, for the axes 1 and 2; Figure 3). Apart from a core of common species, including Oritrophium limnophilum, Werneria apiculata, and Aciachne pulvinata, most species were rep-

| Nurse traits
The nurse trait analysis showed significant differences between the two cushions species reflecting the different strategies to cope with herbivory ( Figure 4; Table S6). First, the leaves of O. andina were significantly thicker (+0.7 mm, p < .0001) and longer (+1.4 square root transformed, p < .0001) than those of D. muscoides. The leaves of O. andina also had a significantly higher mean LDMC than those of D. muscoides, but only in the presence of grazers (grazed plots: +96 mg/g, p < .0001; ungrazed plots: −1 mg/g, p = .9).
Grazing had a significant effect on the measured traits of D. muscoides but not on those of O. andina, suggesting that the protection of D. muscoides against grazing is less efficient than the one of O. andina.
The mean LDMC was lower, and the mean individual vegetative height was higher in D. muscoides in ungrazed than in grazed plots (mean LDMC: −61 mg/g, p < .01; mean height: +0.2 square root transformed, p = .04). However, grazing exclusion did not significantly affect mean leaf thickness in both cushions (D. muscoides: +0.020 mm, p = .80; O. andina: +0.025 mm, p = .80). Note that the LDMC variations observed here are of the same order of magnitude or greater than those reported in other studies on grazing effects on plant communities (Cruz et al., 2010;Whitworth-Hulse, Cingolani, Zeballos, Poca, & Gurvich, 2016).

| Community traits
The associated communities growing in O. andina cushions were clearly less affected by grazing than those growing in D. muscoides ( Figure 5). However, the responses were not homogeneous across traits.
Grazing affected the CWM LDMC of the associated communities only in D. muscoides. In grazed conditions, the CWM LDMC of the associated communities of D. muscoides was lower than in ungrazed condition (−53 mg/g, p = .03), whereas there was no significant effect of grazing on the CWM LDMC of the associated communities of O. andina (+11 mg/g, p = .53). In ungrazed plots, the CWM LDMC of the associated communities was not significantly different between the two cushion species (ungrazed plots: 25 mg/g, p = .47). However, the CWV LDMC F I G U R E 2 Percentage of cover of the associated communities inside the nurse cushion species for control (dark gray) and grazing exclusion (light gray). Error bars represent the 95% confidence interval. The bullet points are outliers of 95% distribution. Letters represent the significantly different groups according to the post hoc contrasts (Table S4)  Regarding leaf thickness, CWM Thickness and CWV Thickness were not significantly affected by grazing or the identity of the cushion (Figure 5b,d).
F I G U R E 3 PCA biplot of the community-plot matrix containing the relative abundances of the associated species. Circles and triangles represent the plots of, respectively, D. muscoides and O. andina cushion species. Control plots are in dark gray, and grazing exclusion plots are in light gray.

| Evidence for indirect facilitation at the community level
The analysis of the nurse functional traits showed that only D. muscoides was affected by grazing (mean LDMC and height), suggesting that D. muscoides is not fully protected against grazing. Together with prior information on the physical structure of D. muscoides (very compact shape with tiny leaves, Figure 4c) and the observed changes in the composition of the associated communities in the presence of grazing (Figure 3), our results suggest that the protection offered by D. muscoides to its associated community is weak or null, confirming our hypothesis. Overall, we showed that (1) the nurse species (and therefore its traits) affects the capacity of the nurse to provide indirect facilitation to its associated communities and that (2)  Previous studies about cushion plant effects (e.g., Cavieres et al., 2014;Kikvidze et al., 2015) have found that the effect of different cushion species on their associated communities was positive and relatively homogeneous (e.g., Cavieres et al., 2014;Kikvidze et al., 2015). An important result of our study is that the identity of the nurse species, even when the species belong to the same life form, can generate different outcomes in terms of plant-plant interactions in grazed conditions. This corroborates the results of a recent study in the dry Central Andes with two different cushion species (Anthelme et al., 2017), in which the authors showed a different outcome in direct facilitation.
We also found that the effect of indirect facilitation through grazers on the associated communities seems to depend on the trait of the associated plants considered. While O. andina significantly reduced grazing effects on CWM LDMC and CWV Height , it did not for F I G U R E 5 Community-weighted means (CWM) (a-c) and community-weighted variances (CWV) (d-f) of the associated communities growing inside the nurse cushion species for control (dark gray) and grazing exclusion (light gray). From the left to the right: Leaf Dry Matter Content (LDMC), leaf thickness, and maximum height. Error bars represent the 95% confidence interval. The bullet points are outliers of 95% distribution. Letters represent the significantly different groups according to the post hoc contrasts (Tables S8 and S9) CWM Height . Grazing was found to decrease height dispersion in agreement with previous literature (Díaz et al., 2001;Sonnier et al., 2010 (Tables S4 and   S8) increased of the same magnitude for each cushion. As cover and maximum height are positively correlated with aboveground biomass (Catchpole & Wheeler, 1992;Muukkonen et al., 2006), it is likely that, if we would have measured the aboveground biomass, we would not have detected any indirect facilitation effect.
Another result of this study is that the associated communities could have a negative impact on the nurse in the absence of grazing pressure in our study system. Maximum height has been linked to the ability to compete for light, as a main axis of the plant life strategy (Westoby, 1998;Westoby et al., 2002). The fact that CWM Height and cover of the associated communities increased in ungrazed compared to grazed plots in both cushions suggests that facilitation could also have a cost for the nurse, for example in terms of competition for light by maintaining potential competitors nearby, in agreement with previous studies (Michalet et al., 2016;Schöb et al., 2014). For example, Michalet et al. (2016) showed that the negative feedback of the beneficiaries on their benefactor was proportional to the cover of the beneficiaries.
Interestingly, the LDMC of the associated communities responded in the opposite way than the LDMC of the nurse in D. muscoides cushions. Indeed, D. muscoides had a higher LDMC in grazed treatments and its associated communities had a lower one. One possible explanation could be the difference of strategy of the nurse and the associated communities in face of grazing. D. muscoides seemingly developed a defense strategy (cf Section 2.1), while the associated communities could have a tolerance strategy as found by previous studies in grassland communities (Louault et al., 2005). Plants with a tolerance strategy tend to have faster regrowth traits (i.e., lower LDMC) in grazed conditions (Louault et al., 2005). On the contrary, plants with a defense strategy tend to have slower regrowth traits (higher LDMC) in grazing conditions compensated by a lower digestibility (Louault et al., 2005;Pontes et al., 2007). Descombes et al. (2016) also found a negative correlation between palatability and LDMC.
Our study is therefore in agreement with a link between tolerance and defense strategies against grazing. The fact that the leaf thickness of the nurses and their associated communities was not affected by grazing is surprising because of evidence from previous literature that physical resistance plays a role in deterring herbivores (Read & Stokes, 2006) and that leaf thickness represents an important part of the physical resistance of the leaves (Pérez-Harguindeguy et al., 2013). Our measurements were possibly not precise enough to capture intraspecific variations in leaf thickness, pinpointing the difficulty to infer functional traits related to the shape of leaves (see also SLA, Loza Herrera et al., 2015). Another explanation could be that LDMC is a more plastic trait than thickness.

| Integration of our results in the current framework on indirect facilitation
A main objective of the present study was to revisit the current framework on indirect interactions at the community level using a functional trait approach. In our study, the configuration we are interested in, with only three components-two plants and a domestic herbivoreactually results in a rather complex set of interactions. Indeed, grazers had a negative effect on the associated communities but this effect was reduced by the nursing effects provided by the nurse plants, nursing (Stachowicz, 2001). However, there is increasing evidence in the literature of nurses that do pay a cost due to nursing (García, Bader, & Cavieres, 2016;Schöb et al., 2014). Our results showed that, without grazers, the associated communities increased in height and cover, thereby possibly capturing light at the expense of the nurse, which would be in agreement with these previous studies.
Our study belongs to the trait-based facilitation framework , and more broadly to trait-based comparative ecology. Indirect facilitation from a nurse to an associated community was found to be trait dependent in our study. Facilitation ecologists have argued that indirect facilitation through grazing depended on the nurse, notably on its palatability or architecture (Catorci et al., 2016). More empirical studies measuring functional traits of both nurses and associated communities are needed to reach a more subtle understanding of the effects of indirect facilitation on associated plant communities, notably on the different life strategy axes (sensu Westoby et al., 2002). In parallel, measurements of functional traits should also be realized along grazing gradients to test if there are linear (Smit et al., 2009), nonlinear (Verwijmeren, Rietkerk, Wassen, & Smit, 2013), or idiosyncratic relationships between indirect facilitation and grazing intensity , and how those relationships vary depending on the functional traits considered.

| CONCLUSION
This study is a contribution to the emerging conceptual framework on trait-based approaches to study positive interactions Gross et al., 2009;Schöb et al., 2012Schöb et al., , 2013. We also contribute to bridge a research gap by showing how indirect facilitation through grazing can affect the structure of associated communities (Filazzola & Lortie, 2014), and we argue for the inclusion of a functional approach at the community level in the more general framework of indirect positive interactions (Aschehoug et al., 2016;Aschehoug & Callaway, 2015). This study notably shows that indirect interactions including a grazer can generate a complex set of interactions within an ecosystem dominated by an apparently homogeneous cover of cushionforming species. However, the ecological consequences of the variations found in ecosystem functions remain to be explored. Finally, this study is in line with the suggestion to integrate trophic and nontrophic, positive and negative interactions in a single framework (Kéfi et al., 2012), a necessity to better understand community and ecosystem dynamics in a global change context.

CONFLICT OF INTEREST
The authors declare no conflict of interest.

AUTHOR CONTRIBUTION
FA, RM, AD, and SK designed the study; FA and RM set up the experimental design; AD collected the data and performed the analysis; AD wrote the first full draft of the manuscript. All authors contributed to the different versions of the drafts and gave final approval for publication.