# Cushion plants as soil engineers: the formation of islands of fertility differ between species in nival environments (feb,10.1007/s11104-025-07274-0,2025)

Keyvan Dumas, Nicolas Bonfanti, Aurélien Grange, Annie Millery-Vigues, The Orchamp Consortium, Christiane Gallet, Glenn Yannic, Jérome Poulenard, Sébastien Lavergne, Sébastien Ibanez

► To cite this version:

Keyvan Dumas, Nicolas Bonfanti, Aurélien Grange, Annie Millery-Vigues, The Orchamp Consortium, et al.. Cushion plants as soil engineers: the formation of islands of fertility differ between species in nival environments (feb,10.1007/s11104-025-07274-0,2025). Plant and Soil, 2025, 10.1007/s11104-025-07274-0. hal-04958206

HAL Id: hal-04958206

https://hal.science/hal-04958206v1

Submitted on 21 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Cushion plants as soil engineers: the formation of
islands of fertility differ between species in nival
environments

Keyvan Dumas*¹, Nicolas Bonfanti²,³, Aurélien Grange¹, Annie Millery-Vigues¹,², The Orchamp Consortium, Christiane Gallet¹, Glenn Yannic¹, Jérome Poulenard³, Sébastien Lavergne¹, and Sébastien Ibanez¹

¹Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA,
Laboratoire d'Écologie Alpine, 38000 Grenoble, France

²Univ. Savoie Mont Blanc, INRAE, CARRTEL, 73370 Le Bourget-du-Lac, France

³Univ. Savoie Mont Blanc, CNRS, Laboratoire Environnement Dynamique et Territoire de la Montagne (EDYTEM), 73370 Le Bourget-du-Lac, France

February 4, 2025

<footnote>*Corresponding author : dumaskvn@mailfence.com</footnote>
16

17

**Keywords**: Cushion plants, Plant-soil feedback, Ecosystem engineering, Alpine ecology, Alpine soils, Nival belt

18

# Summary

19

## Aims

20

21

22

23

24

In extreme environments such as the nival vegetation belt, isolated cushion plants create islands of soil fertility within a mostly mineral environment. It, however, remains unknown whether different co-occurring cushion plant species engineer soils with contrasted properties, and to what extent plant functional traits mediate this process.

25

## Methods

26

27

28

29

30

We sampled 117 individual cushions from seven species, together with 26 bare soil controls from 2600 to 3200m of elevation. Soil property (carbon and nitrogen content, soil C:N and available phosphorus) were compared between species, and the role of cushion plant traits was assessed using Structural Equation Modelling (SEM).

31

## Results

32

33

34

35

36

37

The soils beneath cushion plants showed similar characteristics to those found 1000m below, where vegetation is continuously distributed. Isotopic signatures showed that the soil organic matter beneath cushions originated from the plant itself. Soil characteristics (C, N, C:N, P) differed among plant species. SEM revealed that soil characteristics were not driven by leaf C/N or phenolic content, but rather by cushion size and bedrock type.

38

## Conclusions

39

40

41

Providing a unique overview of soil formation in nival ecosystems, our study provides novel insights into the mechanisms of ecosystem engineering by foundation plants and shows how interspecific variability contributes to this process.
## Introduction

While the soil environment is a highly complex system, its interaction with plants is pivotal in shaping terrestrial ecosystems. Relations between plants and soil are characterized by retroactive loops: on one hand, soil affects the composition of plant communities and functional traits (Ordoñez et al., 2009) while, on the other hand, plants shape soil characteristics. This feedback spans varied spatial (Bever et al., 1997) and temporal scales (Kardol et al., 2006). This interplay between plant and soil is particularly important in extreme environments where plant life reaches its limits of habitability (Ehrenfeld et al., 2005). In dry lands (Ochoa-Hueso et al., 2018; Schlesinger et al., 1990) and high alpine ecosystems (Zhao, 2022), dominant plant species cope with poor soil conditions by forming *islands of fertility*, i.e. patches of vegetation where the soil is locally richer in nutrient and water availability than surrounding bare ground.

In the European Alps, the nival belt extends above 2700m and is characterized by low annual temperatures and a short growing season. These harsh conditions prevent the formation of a continuous vegetation cover, and result in discrete patches of plants (Körner, 2021). These abiotic constraints are amplified by the shallow and coarse nature of high elevation soils (Baroni and Orombelli, 1996; Pintaldi et al., 2021). Pedogenesis in nival soils is described as slow and intermittent as it follows snow cover dynamics (Magnani et al., 2017b) and ice-cover fluctuations (Badino et al., 2018). In the harsh conditions of the nival belt, several plant species from different lineages have evolved toward a cushion life-form (Boucher et al., 2016). Cushion plants may create fertility islands in high mountain ecosystems by trapping dust and nutrients within the plant (Magaña Ugarte et al., 2024), and are the dominant plant life-form in the nival belt, where they increase organic matter (OM) and nutrient content in the soil (Mihoč et al., 2016; Roy et al., 2013; Shuo et al., 2024). Plant characteristics such as size or canopy structure are particularly important in shaping the soil characteristics (Mohebali et al., 2023; Zhao, 2022). Cushion plants are thus often described as engineer species (Badano and Cavieres, 2006; Chen et al., 2019; Wang et al., 2021).

While in multispecies plant communities, pedogenesis is a complex process involv-
ing multiple feedback between the soil and many individual plants belonging to different species (Fig. 1a), isolated cushion plants offer a more tractable study system, where a single individual plant exerts a strong influence on its own soil. This allows the comparison of the consequences of engineering between different plant species. Due to their discrete nature, the influence of cushion plants on soil is concentrated directly beneath their canopy, in the transition zone between the plant litter and the engineered soil (Fig. 1b).

According to (van Breemen and Finzi, 1998), ecosystem engineers correspond to "organisms that directly or indirectly modulate the availability of resources to other species, by causing physical state changes in biotic or abiotic materials". The presence of engineer plant species along an altitudinal gradient might modify the altitudinal distribution of soil characteristics. For example, in mountainous regions, soil organic matter and soil N and P content generally decrease with altitude (Körner, 2021). However, the engineering effect of cushion plant might mitigate this trend, due to the increase of soil organic matter and nutrient content beneath cushion plants, compared to soils not covered by cushion plants (Roy et al., 2013; Zhao, 2022).

Soils and habitats in the alpine and nival belts are heterogeneous due to micro-topographic variations (Choler, 2005; Li et al., 2024). It can therefore be difficult to decipher whether islands of fertility result from environmental modifications induced by plants, or from higher plant survival rates in richer patches (Schlesinger et al., 1990; Steinbauer et al., 2016). Similarly, in the context of plant facilitation it is unclear whether plants are found in association because of the facilitating effect of some species or because plants tend to be clumped in richer patches (Michalet et al., 2024). Although experiments described in Li et al. (2017) showed that abiotic processes alone do not produce island of fertility in shrub-like systems, evidence is still needed to test whether fertility islands result from habitat selection by plants in heterogeneous landscapes, or from ecosystem engineering by plants.

Soil characteristics in plant-engineered systems are strongly influenced by the functional traits of the dominant plant cover, either at the inter- or intra-specific level (Ochoa-Hueso et al., 2018; Walker et al., 2001; Whitham et al., 2006). For in-
stance, intraspecific variations in leaf tannin content within *Populus* spp. is known
to affect soil nutrient dynamics (Whitham et al., 2006). Polyphenolic compounds,
as C-based molecules are also known to influence the decomposability of litter
(Hättenschwiler and Vitousek, 2000). Plants with higher C:N or higher total phe-
nols should therefore show a slower decomposition process, resulting in a higher
soil carbon content. In addition, plant morphological features such as size may be
particularly important in shaping soil characteristics, with an expected increase
in soil organic carbon (SOC), nitrogen and phosphorus as cushions grow, due to
a faster input of dead OM compared to the decomposition rate (Zhao, 2022; Zhao
and An, 2021). Microhabitats, such as screes (Fig. 1b) or rock cracks (Fig. 1c),
might also play a role in ecosystem modification by cushions plants. Their engi-
neering effect is likely to be more spatially diffuse in screes than in cracks, where
pre-existing soil is extremely limited (Fig. 1). We therefore expect a stronger
engineering effect of cushions in cracks than in screes.

Eight cushion plant species were sampled from 2600m to 3200m in six sites through-
out the French Alps (117 individuals plus 26 bare soils). For each individual plant,
we recorded its species and measured several plant traits (leaf C:N, leaf phenolic
content, cushion size) and soil characteristics (SOC, total nitrogen (tN), soil C:N,
bioavailable phosphorus (Pavail) and pH). We first compared the soils formed
beneath cushions to soils found at lower elevations (montane to alpine ecosys-
tems), and hypothesized that soil characteristics beneath cushions correspond to
those usually found at lower elevations (H1). Then, we hypothesized that cushion
plants strongly affect their soil by increasing and concentrating dead OM (H2).
We then compared the characteristics of bare soils and soils beneath cushions,
and the characteristics of soils beneath different species of cushion plants. We
expected higher nutrient and OM content beneath cushions than in bare soils, and
notable inter-specific differences (H3). We then studied the role of plant traits
in the engineering effect of cushion plants, assuming that bigger cushions, and/or
ones with higher leaf C:N or total phenol content should produce soils enriched in
OM and nitrogen, with an increased engineering effect for crack-growing cushion
plants (H4).
(a) Multi-specific plant cover

<figure><img src="image_1.png" /></figure>

PLANT

COVER

SOIL

BEDROCK

<figure><img src="image_2.png" /></figure>

(b) In scree soil

<figure><img src="image_3.png" /></figure>

PLANT
COVER

SOIL

BEDROCK

(c) In rock cracks

<figure><img src="image_4.png" /></figure>

Cushion leaves

Cushion plant
litter

Engineered soil

Fine dust

PLANT
COVER

SOIL

BEDROCK

<figure><img src="image_5.png" /></figure>

Galibier

<figure><img src="image_6.png" /></figure>

(e)

<figure><img src="image_7.png" /></figure>

(f)

**Fig. 1** Conceptual scheme of the different effects of distinct engineer species on ecosystems, with a special emphasis on the contribution of plants to soil creation. (a) In alpine meadows, the soil is influenced by many individuals belonging to sometimes several dozens of different plant species. In the nival zone, each individual cushion plant influences its own soil, in a more diffuse way in scree than in rock cracks (b and c). (d) Example of sampling design of our study at a given site : 2-3 plots (red triangles) sampled per site, with ca. 200 meters between each plot. Illustration is taken from the Galibier site, with lower and higher sites elevation of 2780m a.s.l. and 3190m a.s.l., respectively. (e) Typical habitat of alpine cushion plants, Archeboc site, 3000m. (f) Bare soil sampled in Galibier
135 **Material and Methods**

136 **Sampling**

137 In July 2022, 117 cushion plants belonging to eight different species, along with 138 23 plant-free soils ("bare soils" hereafter; Fig. 1f), were sampled in six altitudinal 139 gradients (hereafter "site") ranging from 2600 to 3200m throughout the French 140 Alps (Table 1). These sites represent a range of climatic and geological conditions, 141 representative of the diversity of the French western Alps (Table 2). On each site, 142 two to three plots were selected every ca. 200m of elevation, in order to cover 143 the altitudinal variation within the nival zone (Fig. 1d). At each plot, three 144 individuals of each cushion plant species were sampled. To avoid confusion due to 145 pluri-specific cushion or facilitated species, we focused on mono-specific cushions 146 only. Some samples lacking material to conduct all analyses were discarded, leading 147 to the final sampling design presented in Table 1. On each plot, two bare soils, 148 uncovered by cushion plants, were also sampled.

149 The following species were sampled : *Silene acaulis* subsp *bryoides* (Jord.) Nyman, 150 1878 (*Silene acaulis* hereafter), *Cherleria sedoides* L., 1753, *Eritrichium nanum* 151 (L.) Schrad. ex Gaudin, 1828, *Saxifraga exarata* Vill., 1779, *Saxifraga bryoides* 152 L., 1753, *Saxifraga oppositifolia* L., 1753, *Androsace saussurei* Dentant, Lavergne, 153 Boucher & Ibanez, 2021, and *Androsace alpina* (L.) Lam., 1779. Data gathered on 154 Androsace alpina and A. saussurei were pooled and referred as *Androsace spp.* to 155 increase the statistical power in later analyses. These species are indeed morpho- 156 logically cryptic species, described as separate species in 2021 only, on the basis of 157 molecular data (Boucher et al., 2021). The following analyses are therefore based 158 on seven distinct species or taxa.

159 At each sampling point, we noted the GPS coordinates using a Garmin© GPSMAP 160 65s portable device. The micro-site of each sampled plant was classified as scree 161 or crack (Fig.1b, c). GPS coordinates measured on field were used to retrieve 162 elevation from IGN (French National Geographic Institute) LiDAR HD digital 163 elevation model at 0.2 m resolution (https://geoservices.ign.fr/lidarhd). 164 Finally, plant diameter was estimated from the average of the longest and shortest
<table><caption>Table 1: Number of samples per site and sample type, i.e. plant, by species, or bare soils. The family of each plant species is also given.</caption><thead><tr><th>Species</th><th>Family</th><th>Viso</th><th>Galibier</th><th>Puy-gris</th><th>Bellecote</th><th>Passon</th><th>Archeboc</th><th>Total</th></tr></thead><tbody><tr><td>Eritrichium nanum</td><td>Boraginaceae</td><td>-</td><td>2</td><td>2</td><td>-</td><td>-</td><td>-</td><td>4</td></tr><tr><td>Cherleria sedoides</td><td>Caryophyllaceae</td><td>4</td><td>2</td><td>-</td><td>5</td><td>6</td><td>5</td><td>22</td></tr><tr><td>Silene acaulis</td><td>Caryophyllaceae</td><td>2</td><td>5</td><td>3</td><td>4</td><td>4</td><td>3</td><td>21</td></tr><tr><td>Androsace spp</td><td>Primulaceae</td><td>-</td><td>-</td><td>-</td><td>1</td><td>3</td><td>4</td><td>8</td></tr><tr><td>Saxifraga exarata</td><td>Saxifragaceae</td><td>3</td><td>7</td><td>-</td><td>4</td><td>2</td><td>1</td><td>17</td></tr><tr><td>Saxifraga bryoides</td><td>Saxifragaceae</td><td>4</td><td>-</td><td>4</td><td>1</td><td>6</td><td>4</td><td>19</td></tr><tr><td>Saxifraga oppositifolia</td><td>Saxifragaceae</td><td>3</td><td>9</td><td>4</td><td>5</td><td>1</td><td>4</td><td>26</td></tr><tr><td>Bare soil</td><td>-</td><td>3</td><td>5</td><td>4</td><td>5</td><td>3</td><td>3</td><td>23</td></tr><tr><td>Total</td><td>-</td><td>19</td><td>30</td><td>17</td><td>25</td><td>25</td><td>24</td><td>140</td></tr></tbody></table>

perpendicular dimension of each individual plant, as in Zhao and An (2021).

<table><caption>Table 2: Main characteristics of sampling site and samples distribution among sites. Minimal and maximal elevation of samples for each site are reported. MAT : Mean Annual Temperature, Prec : Precipitation, SCD : Snow Cover Duration</caption><thead><tr><th>Site</th><th>Elevation (m)</th><th>Lon.</th><th>Lat.</th><th>MAT (°C)</th><th>Prec (mm)</th><th>SCD (days)</th><th>Bedrock</th></tr></thead><tbody><tr><td>Viso</td><td>2745-3094</td><td>7.06</td><td>44.70</td><td>-1.10</td><td>39.74</td><td>193</td><td>Ophiolite</td></tr><tr><td>Galibier</td><td>2784-3192</td><td>6.44</td><td>45.06</td><td>-0.70</td><td>39.35</td><td>246</td><td>Dolomite</td></tr><tr><td>Puy-gris</td><td>2597-2772</td><td>6.14</td><td>45.29</td><td>1.10</td><td>52.13</td><td>209</td><td>Gneiss</td></tr><tr><td>Bellecote</td><td>2899-3121</td><td>6.75</td><td>45.27</td><td>-1.70</td><td>36.72</td><td>215</td><td>Dolomite &amp; Quartzite</td></tr><tr><td>Passon</td><td>2847-3079</td><td>6.98</td><td>45.98</td><td>-4.20</td><td>48.67</td><td>239</td><td>Gneiss</td></tr><tr><td>Archeboc</td><td>2658-2985</td><td>6.98</td><td>45.59</td><td>-0.70</td><td>41.58</td><td>255</td><td>Schist sandstone coal</td></tr></tbody></table>

167 Soil analysis

A variable quantity of bulk soil (10 to 100 grams, depending on available material)

was sampled in the first 10 cm layer of the soil beneath each cushion plant or in bare
soil using a clean tool, and sealed immediately in a plastic bag. For cushion plants growing in cracks, soil was undifferentiated from the cushion litter, and samples were taken indifferently. Soil samples were let to air-dry at lab temperature within 24h after sampling, until no mass loss was observed. Water pH measurement was performed on raw bulk soil following NF ISO 10390-2005 (F), as in Martinez-Almoyna et al. (2020). For all other analyses, soil samples were sieved at 2-mm before grinding (NF ISO 11464 standard). After decarbonation of the samples, C, N, $\delta^{13}C$ and $\delta^{15}N$ contents were measured at the Silvatech facility (Université de Lorraine, AgroParisTech, INRAE, SILVA, F-54000 Nancy, France) on an EA/GA-IRMS analyser. Bioavailable phosphorus was determined using the Olsen method (Olsen, 1954), involving phosphorus extraction with 0.5 M $NaHCO_3$ (NF ISO-11263). Sample with enough soil material were processed for X-Ray fluorescence (XRF) elemental analysis (S8 Tiger Series 2, Bruker) (Towett et al., 2013).

## 183 Leaf trait measurements

On each cushion, a few stems were sampled and air-dried at room temperature until mass stabilized. Leaves were ground to fine powder for C and N content and isotopic analysis at the Silvatech facility on an EA/GA-IRMS analyser. Phenols and flavonoids were extracted with 1.5 mL of 70% ethanol from 20 mg of leaf material, using a Precellys-Evolution extractor-grinder, in 2 cycles of 15 sec at 5800 rpm. Following Lefebvre et al. (2016), for polyphenolic contents, 100 µL of extract were mixed with 3.4 mL of distilled water and 0.5 mL of Folin-Ciocalteu reagent. After vortexing, 1 mL $Na_2CO_3$ was added, and the mixture was heated to 40°C. The absorbance peak was measured at 760 nm with a Jenway UV-Vis 6800 Spectrophotometer. Due to high phenol concentrations in Saxifragaceae, samples were diluted four times before measurement. Absorbance results were calibrated using gallic acid, hence results are given in Gallic Acid Equivalent (GAE) mass per dry matter mass (mgGAE g$^{-1}$). For flavonoids, 1 mL of the ethanolic extract was diluted in 4 mL of distilled water and 0.3 mL of 5% $NaNO_2$, with 0.3 mL of $AlCl_3$ at 10% concentration being added after 5 minutes still. After 1 minute, 2 mL of 1 M NaOH and 2.4 mL of distilled water were added, and absorbance measured at 510 nm. Final values were expressed as Catechin Equivalent (EC) mass per unit
of dry matter (mgCE g⁻¹).

## ORCHAMP dataset

Comparative analysis along an elevational sequence utilized data from the Or-champ long-term monitoring programm (www.orchamp.osug.fr). It consists of multiple altitudinal gradients throughout the western Alps, along which multiple soil parameters and biodiversity, including plants, are measured every 200 m of elevation. Soil description and analysis protocols of top soil layer use the same methods as decribed above. A detailed description of this program can be found in Thuiller et al. (2024), Calderón-Sanou et al. (2023) and Leclerc et al. (2023). Among the several sites sampled in the Orchamps program, we selected the five sites (Archeboc, Bellecote, Passon, Puy-Gris and Viso) that were located just be-low our sampling sites.

## Statistical analyses

All the following analysis were performed using the R software version 4.4.0 (R Core Team, 2024). SOC, soil total nitrogen and bioavailable phosphorus were log-scaled, as were leaf phenolic and flavonoid content, and plant diameter. Elevation, leaf total carbon and nitrogen, and isotopic values of δ¹³C and δ¹³N were centered around zero and their standard deviation scaled to one. We did not apply correction for unequal sample size since linear mixed models (LMMs) are robust to unbalanced sampling design (Pinheiro, 2014).

## Soil characteristics and cushion plant traits variability

Characteristics of soils sampled in this study, as well cushion plant traits, were first investigated using Principal component analysis (PCA). One PCA focused on the geochemical characteristics of soils, based on the five main elements (SiO₂, CaO, Al₂O₃, MgO and Fe₂O₃). The relative nature of XRF measurement was taken into account into the PCA using the *composition* package (van den Boogaart et al., 2024). Mean coordinates of each sampling plot along the first axis of the resulting PCA were used as a "bedrock" composite variable in following analyses. A second
229

230

PCA using SOC, total nitrogen, soil C:N and bioavailable Phosphorus was used to investigate variation between our samples.

231

232

233

234

235

236

237

Inter-specific variation of cushion species traits was analyzed using a PCA on all plant related variables : plant diameter, leaf total phenols and flavonoids content, phenol: flavonoid ratio, leaf carbon, nitrogen and C:N ratio, δ¹³C and δ¹⁵N isotopic signature. The multivariate variation in cushion plant traits was analyzed using Permutational analysis of variance (PERMANOVA), as a function of plant identity, bedrock composition, elevation and geographic position, whose effects were assessed using 9999 permutations.

238

## Influence of engineer plants on an elevational soil sequence

239

240

241

242

243

244

245

246

247

248

249

250

251

252

Combining Orchamp and cushion plants datasets, we obtained the soil characteristics of 151 samples ranging from 1400 m to 3090 m on five distinct gradients. Samples from the Orchamp dataset and nival bare soils (referred as "soil") represent elevation sequences with limited or diffuse engineering effect. Soil under cushion plant (labelled as "cushion") showed the elevational evolution of soils under the concentrated influence of an ecosystem plant engineer. Altitudinal variation in soil characteristics (SOC, tN, C:N or Pavail) was analyzed using LMMs, with an asymmetric quadratic predictor modeling the effect of elevation and another fixed term for the type of soil sequence considered, i.e. soil or cushion. Elevation was centered around zero and scaled to a standard-deviation of one. Sampling site was also included in models as a random effect. The best model was selected after a step-wise variable selection procedure based on corrected Aikake Index Criteria (AICc) and parsimony. Model assumption were checked using DHARMa (Hartig and Lohse, 2022) and MuMln (Bartoń, 2023) packages.

253

## Isotopic signature of plant and soils and enrichment signal

254

255

256

257

258

To test the origin of OM in soils (H2), we focused on δ¹⁵N isotope which reflect sources of nitrogen acquisition, in both plants and soil (Spasojevic and Weber, 2021). We fitted a LMMs to model the soil isotopic signature against a two-way interaction between soil total nitrogen content and cover type, i.e. bare soils or cushion plant. Plot was considered as a random variable in all models.
259 **Effect of the plant species identity on soil characteristics**

260 As different plant species are expected to have different effects on their environ-
261 ment (H3), we investigated how SOC, total nitrogen, C:N ratio and bioavailable
262 phosphorus contents of soil beneath different cushions varied according to plant
263 species, and in comparison to bare soil control samples. Each soil variable was
264 fitted against plant identity using LMMs with site as a random effect. Estimated
265 marginal means (EMM) were computed to obtain a mean value per cushion species,
266 which were compared pairwise using Tukey's HSD test. Resulting p-values were
267 corrected using Benjamini and Hochberg false detection rate (FDR) procedure for
268 multiple comparisons.

269 **Structural Equation Modelling**

270 In the framework of ecosystem engineering, biotic effect from ecosystem engi-
271 neer and abiotic effect from the environment interact in a complex way to shape
272 ecosystems properties. Structural Equation Model (SEM) allows modeling com-
273 plex causal relationships between multiple variables, including direct and indirect
274 effects, based on relationships that are assumed or inferred from previous knowl-
275 edge. Here, we used the piecewiseSEM package (Lefcheck, 2016).
276

276 We built a hypothetical causality structure (Appendix A: Table A1) taking into
277 account Elevation, Bedrock and micro-habitat (μ-habitat in the figures, and rep-
278 resenting screes vs. crack) as abiotic predictors, plant size, leaf C:N and total
279 phenolic content as plant traits predictors, and soil C:N, total nitrogen content
280 (variable total N) and pH as soil response variables. We considered error corre-
281 lation between phenolic content and leaf C:N, as these are related through plant
282 carbon balance. To account for the mutual influence between total nitrogen and
283 soil C:N, a correlated error term was set between these variables.

284 All relationships between variables were fitted using LMMs, with plot and species
285 identiy as random variables. Individual model assumptions were checked using
286 the DHARMa (Hartig and Lohse, 2022) and MuMln packages (Bartoń, 2023).
287 Standardized effect size were estimated by bootstrapping (999 iterations) using
288 the semEff package (Murphy, 2022), and standardized to express change in the
response variable per unit of its standard deviation. Models goodness of fit was assessed using directed separation test, using a χ² test and the Root Mean Square Error of Approximation (RMSEA) index of fit. In a SEM framework, a χ² p-value higher than 0.05 is a sign of good model fit. RMSEA is a [0,1] bounded index, decreasing with the goodness of fit (a value below 0.1 is usually considered as good). The results of the SEM were plotted using the tidySEM package (van Lissa, 2023).

## 296 Results
<table><caption>Table 3: Soil and leaf parameters measured in each cushion plant species and in bare soils in the present study. Values are given as mean values ± standard deviation. Abbreviations: SOC = Soil Organic Carbon, total N = soil total nitrogen, bioavail P = bioavailable phosphorus, Phenols = total polyphenols content, Flavo = total flavonoids content.</caption><thead><tr><th colspan="8">Soil parameters</th></tr><tr><th>Species</th><th>SOC<br>(g.kg<sup>-1</sup>)</th><th>total<br>(g.kg<sup>-1</sup>)</th><th>N</th><th>C:N</th><th>bioavail P<br>(mg.kg<sup>-1</sup>)</th><th>δ<sup>13</sup>C‰</th><th>δ<sup>15</sup>N‰</th></tr></thead><tbody><tr><td>Bare soil</td><td>12.7±11.5</td><td>1.1±1</td><td></td><td>11.9±9.5</td><td>4.3±3.5</td><td>-25.4±1.3</td><td>1.1±2</td></tr><tr><td>Androsace spp</td><td>19.8±13.4</td><td>0.8±0.4</td><td></td><td>28.6±24.2</td><td>4.1±3.6</td><td>-24.4±1.2</td><td>0.1±1.3</td></tr><tr><td>Eritrichium nanum</td><td>132.1±62.7</td><td>10.6±3.3</td><td></td><td>12.1±3.2</td><td>14.2±3.9</td><td>-25.8±0.8</td><td>-3±2.3</td></tr><tr><td>Cherleria sedoides</td><td>42.1±74.6</td><td>2.5±3.4</td><td></td><td>14±4.1</td><td>5.6±5.8</td><td>-25.2±0.9</td><td>-1.1±1.8</td></tr><tr><td>Silene acaulis</td><td>133.7±124.4</td><td>5.5±4.1</td><td></td><td>21±8.5</td><td>13.9±13.2</td><td>-25.1±1</td><td>-2.2±2.4</td></tr><tr><td>Saxifraga bryoides</td><td>48.9±53.1</td><td>3.3±3.4</td><td></td><td>14.9±5</td><td>7.5±6.3</td><td>-25.5±0.8</td><td>-1.5±1.5</td></tr><tr><td>Saxifraga exarata</td><td>45.5±61</td><td>3±2.9</td><td></td><td>13.5±3.7</td><td>7.6±5.7</td><td>-26.4±1.2</td><td>-1.2±1.8</td></tr><tr><td>Saxifraga oppositifolia</td><td>76.5±87.7</td><td>4.3±4.5</td><td></td><td>16.8±6.6</td><td>9.3±7.8</td><td>-26±1.2</td><td>-1.5±2</td></tr></tbody></table><table><thead><tr><th colspan="8">Leaf parameters</th></tr><tr><th>Species</th><th>C (g.kg<sup>-1</sup>)</th><th>N (g.kg<sup>-1</sup>)</th><th>C:N</th><th>Phenols<br>(mgGAE.g<sup>-1</sup>)(mgEC.g<sup>-1</sup>)</th><th>Flavo</th><th>δ<sup>13</sup>C‰</th><th>δ<sup>15</sup>N‰</th></tr></thead><tbody><tr><td>Androsace spp</td><td>430.8±44.3</td><td>16.3±2.7</td><td>26.8±3.6</td><td>17.2±8</td><td>9.2±6.7</td><td>-27.2±1.3</td><td>-5.8±1.1</td></tr><tr><td>Eritrichium nanum</td><td>419±13.9</td><td>14.7±0.9</td><td>28.7±1.9</td><td>9.7±7.1</td><td>7±6.2</td><td>-27.8±0.1</td><td>-6±1.2</td></tr><tr><td>Cherleria sedoides</td><td>438±35.2</td><td>12.5±3.4</td><td>36.9±8.2</td><td>9.5±3.4</td><td>1.6±0.6</td><td>-26.2±0.8</td><td>-5.2±0.9</td></tr><tr><td>Silene acaulis</td><td>410.5±42.9</td><td>14.2±5.4</td><td>32.3±10.2</td><td>14.4±5.9</td><td>6.5±3.1</td><td>-25.3±0.8</td><td>-5±1.5</td></tr><tr><td>Saxifraga bryoides</td><td>416.1±25</td><td>9.3±2.7</td><td>48±12.7</td><td>26.3±9.6</td><td>6±2.2</td><td>-27.8±0.9</td><td>-6.5±1.2</td></tr><tr><td>Saxifraga exarata</td><td>454.5±17.7</td><td>15.7±4.1</td><td>30.6±7.2</td><td>57.1±17.3</td><td>19±7.5</td><td>-28.4±1.3</td><td>-4±2.1</td></tr><tr><td>Saxifraga oppositifolia</td><td>440.1±21.3</td><td>13.2±3.4</td><td>35.2±8.5</td><td>87.9±25.4</td><td>14.1±4.3</td><td>-28±1.3</td><td>-5.3±1.8</td></tr></tbody></table>
297 Soil characteristics

The first two PCA axis explained 93% of the variance in soil geochemical compo-
sition. The first axis separated soils from carbonated to siliceous bedrocks, with
much more dispersion observed in carbonated sites (Appendix A, Fig. A1a), while
the second axis discriminated metamorphic bedrocks based on their content in
Al₂O₃, MgO and Fe₂O₃. Notably, within the Bellecote site, plots had very differ-
ent characteristics, one being of carbonated type while the other was almost pure
silicate (Appendix A, Fig. A1a). The variations of the five soil variables were
accounted for ca. 88% by the first two PCA axis (Appendix A, Fig. A1b). For
bare soils only, SOC ranged from 1.4 g.kg⁻¹ to 43.3 g.kg⁻¹, but reached 398.50
g.kg⁻¹ under cushion plants. Total nitrogen ranged from 0.2 g.kg⁻¹ to 4.0 g.kg⁻¹ in
bare soils, and reached 16.5 g.kg⁻¹ beneath cushions. All soils had low bioavailable
phosphorus content, ranging from 0.3 mg.kg⁻¹ to 45.5 mg.kg⁻¹. Soil C:N varied con-
siderably as well, as detailed in Table 3. As expected from the different bedrock
nature of study sites, soil pH varied greatly (Appendix A, Fig. A1b), ranging from
5.0 to 9.1. These results highlight the high heterogeneity of nival soils in the Alps.

313 Cushion plant traits

The Permanova analysis showed that plant traits were mainly driven by species
identity which accounted for over 40% of explained variance while elevation and site
explained less than 10% altogether (Table A2). The strong effect of species identity
on traits was also revealed by the PCA (Fig. 2. Inter-specific differences mainly
concerned secondary compound content and carbon isotope fractioning (Fig. 2a).
Saxifragaceae had much higher total phenolic content than other cushion plants
(Fig. 2b). Conversely, Saxifragaceae displayed a low Flavonoid:Phenol ratio, and
a low δ¹³C signature. Leaf C:N varied between and within species, ranging from
15.7 to 73.7 (mean leaf C:N = 35.61 ± 10.8). The highest C:N values were reached
by the nitrogen-depleted *Saxifraga bryoides* (mean leaf C:N = 48.0 ± 12.7, Fig.
2c). Mean plant traits are given in Table 3.
<figure><img src="image_8.png" /><figcaption>**Fig. 2** Cushion plant traits. (a) First two axis of a Principal Component Analysis including leaf carbon content, leaf nitrogen content, leaf C:N ratio, total phenolic and flavonoids content, plant diameter, leaf δ¹³C and leaf δ¹⁵N, showing a strong separation between Saxifragaceae and other study plant families. (b) Leaf C:N ratio and (c) total Phenols content (mg Gallic Acid Equivalent per gram of dry matter) values for all cushion plant species. Total Phenols x-axis is log scaled. Species abbreviation: And spp = *Androsace* spp., Eri nan = *Eritrichium nanum*, Che sed = *Cherleria sedoides*, Sil aca = *Silene acaulis*, Sax bry = *Saxifraga bryoides*, Sax exa = *Saxifraga exarata*, Sax opp = *Saxifraga oppositifolia*</figcaption></figure>
## Soils beneath cushions compared to soils along an altitudinal gradient

Since elevation is likely to be a strong driver of local environments in alpine ecosystems, we expected that cushion plants buffer these changes (H1). SOC, total nitrogen (total N), soil C:N and bioavailable phosphorus (bioavail P) all responded to elevation, and all but soil phosphorus increased under the influence of cushion plants (Appendix A, Table A3). At the exception of soil C:N, elevation and cover type accounted for more than 20% of the variance, while site explained only 10% (Table A3). Soil C:N and bioavailable phosphorus steadily decreased with altitude (Fig. 3c, d), while SOC and total nitrogen showed a bell-shaped curve with a maximum in the alpine vegetation belt, at around 2000 m (Fig. 3a, b). In the presence of cushion plants, SOC, total nitrogen and soil C:N were much increased compared to bare soils (Fig 3a, b, c). Cushion plants thus modify the altitudinal trend of SOC, total nitrogen and soil C:N on the upper part of the sequence and considerably increase the OM content of the soil by their presence, soil reaching levels similar to those observed 1000m lower.

## Nitrogen isotopic signature under the influence of cushion plants

Using nitrogen isotopic signatures of soil and leaves, we aimed to demonstrate how plants contribute to soil OM. The greater this contribution, the closer the soil $\delta^{15}$N should be to the leaves' $\delta^{15}$N (H2). Leaves had a lower $\delta^{15}$N values than soils (Leaf $\delta^{15}$N = -5.47 ± 1.63 ‰; Soil $\delta^{15}$N = -1.03 ± 2.19 ‰; P-value < 0.001; Fig. 4A), as well as lower $\delta^{13}$C values (Leaf $\delta^{13}$C = -27.1 ± 1.53 ‰, Soil $\delta^{13}$C = -25.5 ± 1.17 ‰; P-value < 0.001; Fig. 4A). Bare soils and soils beneath cushions had similar $\delta^{13}$C signature (Bare soil $\delta^{13}$C = -25.4 ± 1.24 ‰, Cushion soil $\delta^{13}$C = -25.5 ± 1.16 ‰; P-value = 0.32), but $\delta^{15}$N was higher in bare soils than in soils beneath cushions (Bare soil $\delta^{15}$N = 1.07 ± 1.93 ‰, Cushion soil $\delta^{15}$N = -1.45 ± 1.98 ‰; P-value < 0.01). Linear regression showed contrasted dynamics between bare soils and cushion plant soils (P-value = 0.001; Fig 4B). The evolution of soil isotopic signature was well described by the interaction between cover type and
<figure><img src="image_9.png" /><figcaption>**Fig. 3** Observed and modeled values of (a) Soil Organic Carbon (SOC), (b) total nitrogen (total N), (c) soil C:N ratio and (d) bioavailable Phosphorus (bioavail P) along an altitudinal sequence and under the effect of a high alpine engineer species. Regression lines from Linear Mixed Models fitting the response variable to an asymmetric quadratic term of elevation and the type of elevational sequence, under typical alpine vegetation or under cushion plant influence. Light grey points represent classical sequence from the montane to the nival belt, dark grey points alpine and nival soils under the influence of cushion plants. Dashed lines represent 95% confidence interval. All y-axis are log₁₀ scaled</figcaption></figure>
<figure><img src="image_10.png" /></figure>

**Fig. 4** a) Carbon and nitrogen isotope biplot of soils and leaf samples, for δ¹³C and δ¹⁵N. Top and right plots show the distribution of isotopic signature for δ¹³C and δ¹⁵N respectively, for leaf and soil samples. Among soil samples, the ones corresponding to bare soil are highlighted as triangles. b) Effect of soil organic matter increase and cover type (Bare soil vs cushion plant) on soil δ¹⁵N signature. Solid lines indicate results from Linear Mixed Model regression including site as random effect, dashed lines indicate 95% confidence interval. Horizontal dashed lines represent the isotopic signature of atmospheric nitrogen and the mean value of cushion plant leaves, as measured in the present study

soil total nitrogen ($R^2 = 0.47$; P-value = 0.001, Appendix A, Table A4). As OM built-up in soils beneath cushion plants, we found that the δ¹⁵N isotopic signature of soils decreased, tending toward strongly negative values (min δ¹⁵N<sub>soil</sub> = -6.17 ‰), similar to the isotopic signature of plant tissues (mean δ¹⁵N<sub>plant</sub> = -5.47 ± 1.63 ‰), thus confirming H2. No such relationship between soil δ¹⁵N and total nitrogen was observed in bare soils (P-value ~ 0.9; Appendix A, Table A4; Fig. 4B).
362
363

Interspecific differences in soil characteristics beneath cushion plants, and comparison with bare soils

364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381

Considering the importance of traits in ecosystem functioning, we hypothesized that soil should differ between cushion plant species, and be poorer in bare soils (H3). Cushion species identity influenced all considered response variables (Fig. 5), explaining 14%-28% of the variance, while site explained 17%-41%. Bare soils had low SOC (12.7 ± 11.5 g.kg⁻¹) and total nitrogen content (1.10 ± 3.81 g.kg⁻¹) while soil beneath cushions were enriched in SOC (Effect size = 4.70 ± 2.41, SOC = 69.3 ± 88.1 g.kg⁻¹; Appendix A, Table A5) and total nitrogen (Effect size = 3.19 ± 1.54, total nitrogen = 3.81 ± 3.98; Appendix A, Table A5). Measures of soil C:N ratio varied substantially important between soils sampled beneath cushions (16.9 ± 9.14) and bare ground (11.9 ± 9.52). This result was also observed for SOC and total nitrogen, but to a lesser extent. SOC and total nitrogen beneath *Androsace* species did not differ from bare soils. *Eritrichium nanum* soils were rich in SOC (132.0 ± 62.7 g.kg⁻¹) and total nitrogen content (10.6 ± 3.27 g.kg⁻¹), but had similar C:N ratio than bare soils (12.1 ± 3.20). *Silene acaulis* had the highest values of SOC (134.0 ± 124.0 g.kg⁻¹), total nitrogen (5.52 ± 4.07 g.kg⁻¹), soil C:N (21.0 ± 8.45) and bioavailable phosphorus (13.9 ± 13.2 mg.kg⁻¹), although with an important variability between samples. Other species had a lower effect on soil bioavailable phosphorus (Appendix A, Table A5).

382
383

Causal relationships between abiotic constraints, plant traits
and soil characteristics

384
385
386
387
388
389
390
391

We assumed that plant traits variation drove inter-specific variations in the pro-
cess of soil engineering (H4), and tested this using a SEM. The SEM model was
consistent with the data (χ² = 7.41, p-value = 0.49; RMSEA = 0; Fig 6). No
important path was omitted according to the directed separation test (Fisher's C
= 15.51, p-value= 0.48). Furthermore, the three output variables (Soil C:N, soil
total nitrogen and pH) were well explained by the model (total N R²m = 0.38, R²c
= 0.52; C:N ratio R²m = 0.28, R²c = 0.56; pH R²m = 0.42, R²c = 0.86; Fig.6a). In
the following, effect size values are expressed as standardized bootstrapped effects
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

**Fig. 5** Values of a) Soil Organic Carbon (SOC), b) total nitrogen (total N), c) C:N ratio d) and bioavailable Phosphorus (bioavail P) in cushion plants and bare soils. Squares and error bars indicate the estimated marginal means, and standard deviation for each trait, and for each species. Points indicate measured values. All x-axis are log-scaled. Compact letters indicate means that were significantly different following Tukey's HSD tests, once corrected for false discovery rate. P-value < 0.001 = ***
(± standard error), and are hence unit-less.

Cushion plant traits were weakly influenced by abiotic conditions. Total phenolic content was slightly higher in crack-growing cushions (0.09 ± 0.06; Fig. 6a) while leaf C:N decreased with elevation (-0.23 ± 0.11; Fig. 6a), although this relationship was weakly significant (P-value = 0.07). These plant traits were negatively correlated (Errors correlation = -0.19; Fig. 6a). However, cushion plant size significantly increased with both soil C:N (0.34 ± 0.08; Fig. 6b) and soil total nitrogen (0.23 ± 0.1; Fig. 6c). While soil C:N was independent of pH, total nitrogen decreased in acidic conditions (-0.17 ± 0.10; P-value = 0.08; Fig. 6c). Soil characteristics were strongly influenced by bedrock (Appendix A, Fig. A1a). Total nitrogen content was lower on siliceous than on calcareous bedrocks (-0.25 ± 0.11; Fig. 6c). Soil pH was strongly linked to the geological nature of the bedrock (-0.68 ± 0.13; Fig. 6d). Micro-habitat of the cushion affected soil total nitrogen directly (0.24 ± 0.08; Fig. 6c), with cushions growing in cracks being richer. However, soil C:N remained constant, indicating that SOC followed the increase in total nitrogen.

## Discussion

This study provides novel insights into the mechanisms of ecosystem engineering by foundation plants. In particular, it shows how inter-specific trait variability among ecosystem engineers leads to the formation of contrasted soils. In the nival vegetation belt, cushion plants create fertility islands enriched in SOC, total nitrogen and C:N, which mimicked soil conditions found at lower elevations, according to our expectations (H1, Fig. 3). In line with H2, nitrogen isotopic signatures showed that the OM found in cushion plant soils mainly originates from the cushion itself. Indeed, the isotopic signature of soils tend to be similar to that of leaves when total nitrogen content in soil increased beneath cushion plants, whereas the isotopic signature and nitrogen content were uncorrelated in bare soils (Fig. 4). Almost all cushion species increased SOC, total nitrogen and C:N ratio in soils (Fig. 5). Bioavailable phosphorus was increased by some species (*Silene acaulis*, but generally did not differ from those found in bare soils. The diversity of ecosystems
<figure><img src="image_15.png" /><figcaption>**Fig. 6** Path Analysis of direct effects of cushion plant traits and abiotic factors on soil parameters. (a) Causal relationship between predictors. Edge sizes are proportional to their effect. Direct effect of predictors on (b) soil C:N, (c) total nitrogen (total N) and (d) soil pH</figcaption></figure>
tem engineers thus contributed to the heterogeneity of fertility islands (in line with H3) through their species-specific effects (Fig. 5), but this was not explained by leaf C:N and leaf phenolic content (contrary to H4, Fig. 6). Nevertheless, plant size increased soil C:N and total N content (Fig. 6), the latter being also higher in rock crack microhabitats, compared to scree slopes.

Several works previously studied soil characteristics beneath cushion plants. The range of values we found for SOC under *Silene acaulis* cushions (Table 3, Fig. 5) was in line with the results obtained by Roy et al. (2013) (SOC ∈ [150,250] g.kg⁻¹) or Bonanomi et al. (2016) (389.5 ± 99 g.kg⁻¹), but were depleted in total nitrogen when compared to these two studies (6 to 9 g.kg⁻¹ and 9.7 ± 3.1 g.kg⁻¹, respectively). *Saxifraga oppositifolia* soil's in Svalbard had total nitrogen content ranging from 0.15 to 10 g.kg⁻¹, in line with our findings (Kume et al., 2003). Soils under *Thylacospermum caespitosum* cushions in the Himalaya region displayed contrasted characteristics, with lower SOC and total nitrogen (SOC ∈ [0.5,45] g.kg⁻¹, tN ∈ [0.5,2] g.kg⁻¹) described in Zhao (2022), but values similar to ours were found in Zhao and An (2021). *Potentilla articulata* and *Arenaria polytrichoides* had lower total nitrogen (0.3 ± 0.03 g.kg⁻¹ and 0.2 ± 0.01 g.kg⁻¹, respectively) (Chen et al., 2015a), while soils under *Androsace tapete* (He et al., 2014) had values similar to ours (SOC = 68.1 ± 7.7 g.kg⁻¹, tN = 4.4 ± 0.2 g.kg⁻¹, C:N = 15.1 ± 1.0). Data regarding growth and metabolic rate of cushion plants are missing. However, given the very conservative nature of these plants, very low rates of resource acquisition by cushion plants could contribute to explain their potential to generate fertility islands.

Although the magnitude of soil amelioration by cushion plants is generally variable (Momberg and le Roux, 2020), several previous studies observed soil amelioration by cushion plants when compared to bare soils (Bonanomi et al., 2016; Roy et al., 2013; Zhao, 2022; Zhao and An, 2021), as supported by our findings. The positive influence of plant size on SOC, total nitrogen and phosphorus was also assessed in *Thylacospermum caespitosum* (Zhao, 2022; Zhao and An, 2021), although these studies also detected an effect of elevation on soil characteristics, and an effect of plant size on pH, which we did not observe (Fig. 6), the latter being probably due to the diversity of pH across sites in the present study. However, Roy et al.
454

455

(2013) described a change in soil characteristics between siliceous and calcareous bedrocks, a finding in line with our results (Fig. 6a)

456

457

458

459

460

461

It should be acknowledged that sampling constraints might have led to an over-
estimation of the engineering effect of some cushion plants. Since larger cushions
or those located in cracks increase the proportion of litter in their volume, it is
likely that this litter was sampled in bigger proportion than for smaller cushion or
cushions located in screes. This could artificially inflate the OM content of soils
beneath larger cushions and those in cracks.

462

Soil characteristics in the nival vegetation belt

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

Bare soils sampled in this study showed similar SOC (around 10 g.kg⁻¹), total
nitrogen (around 1 g.kg⁻¹), C:N or available Phosphorus than previous studies
on upper horizons of soils sampled in the nival belt in the Alps (Bockheim and
Munroe, 2014; Haselwandter et al., 1983; Huber et al., 2007; Magnani et al., 2017a;
Pintaldi et al., 2021; Praeg et al., 2019; Zumsteg et al., 2013). The organic and
nutrient characteristics of these soils result from the decomposition of past plant
cover, together with atmospheric deposition of carbon and nitrogen (Baroni and
Orombelli, 1996; Haselwandter et al., 1983; Huber et al., 2007). Though, the SOC
and nutrient content of the present samples (Fig. 5), was more variable than pre-
viously reported. In most studies a single site was sampled (Huber et al., 2007;
Pintaldi et al., 2021), while the present study compared several sites having con-
trasted geologies and elevations. Conversely, soil pH was relatively stable within
siliceous and carbonated soils, and was never more than two units from seven, re-
flecting the raw nature of the nival soils (Huber et al., 2007; Munroe and Bockheim,
2001).

478

479

Evidence of plant engineering of soil through isotopic anal-
ysis

480

481

482

Our measurements of bare soil δ¹⁵N ranged within previous findings at such el-
evations (Huber et al., 2007). The nitrogen isotopic signature of soils and leaves
were independent from elevation, unlike previous results by Huber et al. (2007).
This might be explained by wider elevational range sampled in the present study (2600-3200m), while Huber et al. (2007) sampled a smaller gradient (2900-3100m) in a single location. Values from cushion plant leaves fell within the lower range of what is expected for alpine plants $\delta^{15}$N = -10‰ to 5‰ (Huber et al., 2007; Körner et al., 2016). In the study of facilitation and engineering processes, it is often difficult to distinguish the actual amelioration of the ecosystem by a particular engineer species from habitat selection within an heterogeneous landscape (Schlesinger et al., 1990; Steinbauer et al., 2016). By comparing soil $\delta^{15}$N signature and soil nitrogen, we demonstrated that cushion plants influence soil nitrogen dynamics in the nival belt, likely through multi-factorial mechanisms. First, the tendency for soils $\delta^{15}$N to converge toward low values could result from plant supplying $^{15}$N depleted OM to bulk soil. But since cushion plant size and total nitrogen are correlated, it is possible that cushion plants reduce soil $\delta^{15}$N by trapping atmospheric deposition within their structure more so than bare soils. Both cases highlight how cushion plants modify their environment by their physical structure, i.e. by being autogenic engineers (Jones et al., 1994)

## Plants characteristics structure their ecosystem

Considering the importance of leaf functional traits in litter decomposition and soil formation processes (Chen et al., 2015b; Chomel et al., 2016), we hypothesized that these traits would also drive the soil engineering process (H4). We found that leaf C:N ratio and phenolic content did not impact soil characteristics, but slightly responded to elevation, with negative and positive correlations respectively (Appendix A, Table A2). This is consistent with expected plant adaptations to abiotic constraints and habitat heterogeneity in the alpine zone (Choler, 2005; Magaña Ugarte et al., 2019; Stark et al., 2017).

Conversely, we found that cushion diameter, a trait more related to the physical influence of the engineer species on its environment, positively affected soil properties. Larger plant size increase the spatial scale of the engineering effect of cushion plants, but in the case of cushion plants, size is also related to cushion age. Although this relation may be loose, as many parameters such as exposition, soil type (scree or crack) or plant species influence growth rate (Rai et al., 2024),
bigger plants may generally reflect a longer period of ecosystem engineering when compared to smaller ones. This is likely to result in larger organic matter accumulation, as observed in *Thylacospermum caespitosum* cushions (Zhao, 2022), since decomposition rates are extremely slow in the nival belt (Ohtsuka et al., 2008). In addition, the dynamic of organic matter deposition by cushion plants in relation to their age is unclear. Cushion plant growth rate decrease with age (Rai et al., 2024), but the quantity of deposited litter per year is likely increased in larger and older plants, because they have larger canopy and leaves surface. Bigger cushion plants might thus deposit more litter each year, and since a longer time, than smaller ones. It is noteworthy that *Silene acaulis*, known for its cushions lasting several hundred of years (Doak and Morris, 2010; Morris and Doak, 1998) contributed the most to soil OM in our study system.

Morphological and architectural features of cushion plant likely play a key role in OM accumulation by influencing both its deposition and loss rate. Deposition is mostly driven by the plant growth characteristics, which is expected to be slow due to the conservative nature of plant life at such elevations. Plant architecture can control several processes, from the composition of soil seed bank (Erfanzadeh et al., 2020) to litter loss rate, while cushion plant shape can evolve under abiotic constraint such as wind (Combrinck et al., 2020). Looser, and hence more rugged canopy, could enhance the trapping of wind-borne organic matter and soil particles. However, below cushion plants from the nival belt, most of the SOC and soil nutrient origin from the cushion itself, as shown using isotopic analysis. Closed canopy are hence likely to be more efficient in trapping dead litter from the cushion and in limiting the leaching of nutrients by water runoffs and snowmelt or aeolian erosion. Species with a denser or more sealed canopy, such as *Silene acaulis* and *Cherleria sedoides*, might thus retain more OM (Bonanomi et al., 2016) thereby increasing their engineering effect and the formation of fertility islands. Lastly, since the plant traits measured here only explained part of the inter-specific variation observed, other unmeasured traits, particularly those related to plant architecture such as canopy density or stem structure, might also drive the engineering process.
544

Insights into the evolutionary consequences of plant-soil feed-back
back

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

Plants and soils are intricately linked by feedback loops, since plants contribute to drive soil characteristics (as shown here) and, reciprocally, soil influence plant growth. On the long term, such feedbacks might have impacted the evolutionary dynamics of plants and soil organisms (Schweitzer et al., 2014; van Breemen and Finzi, 1998). Such recursive modifications of the selective pressures experienced by organisms due to the alteration of the environmental conditions has been termed "niche construction" (Odling-Smee et al., 2003). Theoretical studies have shown that niche construction is stronger in constrained environments (Arnoldi et al., 2020; Kylafis and Loreau, 2008). While the dense architecture of cushion plants is more likely an adaptation for the resistance to harsh climatic conditions, soil engineering may have been an evolutionary by-product. This is supported by the finding that the leaf functional traits studied here do not influence the quality of the soil built by the cushions. Subsequently, soil engineering could have influenced the selective pressures experienced by plants, potentially triggering the evolution of traits that are yet to be discovered.

561

562

563

564

565

566

567

Given the results of this study, it is highly likely that ecosystem engineering has been, and continues to be, one of the main drivers of biodiversity maintenance in the harsh and hostile conditions of high elevation ecosystems. To better understand this process, it is essential to accurately characterize the communities hosted by cushion plants at multiple taxonomic levels. This would enable us to decipher the respective contribution of plants and soil organisms to the engineering process of the belowground environment.

568

Acknowledgments

569

570

571

572

The authors thank the SILVATECH platform (Silvatech, INRAE, 2018. Structural and functional analysis of tree and wood Facility, doi: 10.15454/1.55724 00113627854E12) at INRAE Grand-Est - Nancy for providing carbon and nitrogen content and isotopes analysis. SILVATECH is supported by the French
573 National Research Agency through the Cluster of Excellence ARBRE (ANR-11-LABX-0002-01). We would also like to thank Anne-Lise Develle and Alice Molard (Edytem) for the XRF measurements on the soil samples, Ilan Bourgeois (Cartel) for discussion on nitrogen isotopy, Sylvain Coq (CEFE) for his encouraging comments, Alexy Rosa and Bastien Féaud who helped during field work and sample management, and Marine Allimant who prepared samples before analysis by SILVATECH. Finally, we are thankful to the three anonymous reviewers whose comments greatly helped to improve this manuscript.
574
575
576
577
578
579
580

References

582 Arnoldi, J.-F., Coq, S., Kéfi, S., and Ibanez, S. (2020). Positive plant-soil feedback
583 trigger tannin evolution by niche construction: A spatial stoichiometric model.
584 *Journal of Ecology*, 108(1):378-391.

585 Badano, E. I. and Cavieres, L. A. (2006). Impacts of ecosystem engineers on
586 community attributes: Effects of cushion plants at different elevations of the
587 Chilean Andes. *Diversity and Distributions*, 12(4):388–396.

588 Badino, F., Ravazzi, C., Vallè, F., Pini, R., Aceti, A., Brunetti, M., Champvillair, E., Maggi, V., Maspero, F., Perego, R., and Orombelli, G. (2018). 8800 years of high-altitude vegetation and climate history at the Rutor Glacier forefield, Italian Alps. Evidence of middle Holocene timberline rise and glacier contraction. *Quaternary Science Reviews*, 185:41–68.

593 Baroni, C. and Orombelli, G. (1996). The Alpine “Iceman” and Holocene Climatic
594 Change. *Quaternary Research*, 46(1):78–83.

595 Bartoń, K. (2023). MuMIn: Multi-Model Inference.

596 Bever, J. D., Westover, K. M., and Antonovics, J. (1997). Incorporating the
597 Soil Community into Plant Population Dynamics: The Utility of the Feedback
598 Approach. *Journal of Ecology*, 85(5):561–573.

599 Bockheim, J. G. and Munroe, J. S. (2014). Organic Carbon Pools and Genesis of
600 Alpine Soils with Permafrost: A Review. *Arctic, Antarctic, and Alpine Research*,
601 46(4):987–1006.

602 Bonanomi, G., Stinca, A., Chirico, G. B., Ciaschetti, G., Saracino, A., and Incerti,
603 G. (2016). Cushion plant morphology controls biogenic capability and facilita-
604

605

tion effects of Silene acaulis along an elevation gradient. *Functional Ecology*, 30(7):1216-1226.

606

607

608

Boucher, F. C., Dentant, C., Ibanez, S., Capblancq, T., Boleda, M., Boulangeat, L., Smyčka, J., Roquet, C., and Lavergne, S. (2021). Discovery of cryptic plant diversity on the rooftops of the Alps. *Scientific Reports*, 11(1):11128.

609

610

611

Boucher, F. C., Lavergne, S., Basile, M., Choler, P., and Aubert, S. (2016). Evolution and biogeography of the cushion life form in angiosperms. *Perspectives in Plant Ecology, Evolution and Systematics*, 20:22-31.

612

613

614

615

616

Calderón-Sanou, I., Ohlmann, M., Münkemüller, T., Zinger, L., Hedde, M., Ligonnet, C., Martinez-Almoyna, C., Saillard, A., Renaud, J., Le Guillarme, N., Gielly, L., Consortium, O., and Thuiller, W. (2023). Mountain soil multitrophic networks shaped by the interplay between habitat and pedoclimatic conditions. *Soil Biology and Biochemistry*, page 109282.

617

618

619

620

Chen, J., Yang, Y., Stöcklin, J., Cavieres, L. A., Peng, D., Li, Z., and Sun, H. (2015a). Soil nutrient availability determines the facilitative effects of cushion plants on other plant species at high elevations in the south-eastern Himalayas. *Plant Ecology & Diversity*, 8(2):199-210.

621

622

623

Chen, J.-G., He, X.-F., Wang, S.-W., Yang, Y., and Sun, H. (2019). Cushion and shrub ecosystem engineers contribute differently to diversity and functions in alpine ecosystems. *Journal of Vegetation Science*, 30(2):362-374.

624

625

626

627

Chen, Y., Sun, J., Xie, F., Wang, X., Cheng, G., and Lu, X. (2015b). Litter chemical structure is more important than species richness in affecting soil carbon and nitrogen dynamics including gas emissions from an alpine soil. *Biology and Fertility of Soils*, 51(7):791-800.

628

629

Choler, P. (2005). Consistent Shifts in Alpine Plant Traits along a Mesotopographical Gradient. *Arctic, Antarctic, and Alpine Research*, 37(4):444-453.

630

631

632

633

Chomel, M., Guittonny-Larchevêque, M., Fernandez, C., Gallet, C., DesRochers, A., Paré, D., Jackson, B. G., and Baldy, V. (2016). Plant secondary metabolites: A key driver of litter decomposition and soil nutrient cycling. *Journal of Ecology*, 104(6):1527-1541.

634

635

636

Combrinck, M. L., Harms, T. M., McGeoch, M. A., Schoombie, J., and le Roux, P. C. (2020). Wind and seed: A conceptual model of shape-formation in the cushion plant Azorella Selago. *Plant and Soil*, 455(1):339-366.
637 Doak, D. F. and Morris, W. F. (2010). Demographic compensation and tipping points in climate-induced range shifts. *Nature*, 467(7318):959–962.
638

639 Ehrenfeld, J. G., Ravit, B., and Elgersma, K. (2005). Feedback in the plant-soils system. *Annual Review of Environment and Resources*, 30(1):75–115.
640

641 Erfanzadeh, R., Hazhir, S., and Jafari, M. (2020). Effect of cushion plants on the soil seed bank in overgrazed semiarid regions. *Land Degradation & Development*, 31(8):990–1000.
642
643

644 Hartig, F. and Lohse, L. (2022). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models.
645

646 Haselwandter, K., Hofmann, A., Holzmann, H. P., and Read, D. J. (1983). Availability of nitrogen and phosphorus in the nival zone of the Alps. *Oecologia*, 57(1-2):266–269.
647
648

649 Hättenschwiler, S. and Vitousek, P. M. (2000). The role of polyphenols in terrestrial ecosystem nutrient cycling. *Trends in Ecology & Evolution*, 15(6):238–243.
650

651 He, Y., Xu, X., Kueffer, C., Zhang, X., and Shi, P. (2014). Leaf litter of a dominant cushion plant shifts nitrogen mineralization to immobilization at high but not low temperature in an alpine meadow. *Plant and Soil*, 383(1):415–426.
652
653

654 Huber, E., Wanek, W., Gottfried, M., Pauli, H., Schweiger, P., Arndt, S. K., Reiter, K., and Richter, A. (2007). Shift in soil–plant nitrogen dynamics of an alpine–nival ecotone. *Plant and Soil*, 301(1):65–76.
655
656

657 Jones, C. G., Lawton, J. H., and Shachak, M. (1994). Organisms as ecosystem engineers. In *Ecosystem Management*, pages 130–147. Springer.
658

659 Kardol, P., Martijn Bezemer, T., and van der Putten, W. H. (2006). Temporal variation in plant–soil feedback controls succession. *Ecology Letters*, 9(9):1080–1088.
660
661

662 Körner, C. (2021). *Alpine Plant Life: Functional Plant Ecology of High Mountain Ecosystems*. Springer Nature.
663

664 Körner, C., Leuzinger, S., Riedl, S., Siegwolf, R. T., and Streule, L. (2016). Carbon and nitrogen stable isotope signals for an entire alpine flora, based on herbarium samples. *Alpine Botany*, 126(2):153–166.
665
666

667 Kume, A., Bekku, Y. S., Hanba, Y. T., and Kanda, H. (2003). Carbon Isotope Discrimination in Diverging Growth Forms of Saxifraga oppositifolia in Different Successional Stages in a High Arctic Glacier Foreland. *Arctic, Antarctic, and Alpine Research*, 35(3):377–383.
668
669
670
671

672

Kylafis, G. and Loreau, M. (2008). Ecological and evolutionary consequences of niche construction for its agent. *Ecology Letters*, 11(10):1072–1081.

673

674

675

676

Leclerc, L., Calderón-Sanou, I., Martinez-Almoyna, C., Paillet, Y., Thuiller, W., Vincenot, L., and Kunstler, G. (2023). Beyond the role of climate and soil conditions: Living and dead trees matter for soil biodiversity in mountain forests. *Soil Biology and Biochemistry*, 187:109194.

677

678

679

Lefcheck, J. S. (2016). piecewiseseem: Piecewise structural equation modeling in r for ecology, evolution, and systematics. *Methods in Ecology and Evolution*, 7(5):573–579.

680

681

682

Lefebvre, T., Millery-Vigues, A., and Gallet, C. (2016). Does leaf optical absorbance reflect the polyphenol content of alpine plants along an elevational gradient? *Alpine Botany*, 126(2):177–185.

683

684

685

Li, J., Gilhooly III, W. P., Okin, G. S., and Blackwell III, J. (2017). Abiotic processes are insufficient for fertile island development: A 10-year artificial shrub experiment in a desert grassland. *Geophysical Research Letters*, 44(5):2245–2253.

686

687

688

Li, X., Li, X., Shi, Y., Zhao, S., Liu, J., Lin, Y., Li, C., and Zhang, C. (2024). Effects of microtopography on soil microbial communities in alpine meadows on the Qinghai-Tibetan Plateau. *CATENA*, 239:107945.

689

690

691

Magaña Ugarte, R., Escudero, A., and Gavilán, R. G. (2019). Metabolic and physiological responses of Mediterranean high-mountain and alpine plants to combined abiotic stresses. *Physiologia Plantarum*, 165(2):403–412.

692

693

694

Magaña Ugarte, R., Gavilán, R. G., and Rubio, A. (2024). Assessing the occurrence of soil improvement and its relationship to the dominant life form in the high mountains of Central Spain. *Geoderma Regional*, 36:e00744.

695

696

697

Magnani, A., Viglietti, D., Balestrini, R., Williams, M. W., and Freppaz, M. (2017a). Contribution of deeper soil horizons to N and C cycling during the snow-free season in alpine tundra, NW Italy. *CATENA*, 155:75–85.

698

699

700

701

Magnani, A., Viglietti, D., Godone, D., Williams, M. W., Balestrini, R., and Freppaz, M. (2017b). Interannual Variability of Soil N and C Forms in Response to Snow—Cover duration and Pedoclimatic Conditions in Alpine Tundra, Northwest Italy. *Arctic, Antarctic, and Alpine Research*, 49(2):227–242.

702

703

704

705

Martinez-Almoyna, C., Piton, G., Abdulhak, S., Boulangeat, L., Choler, P., Delahaye, T., Dentant, C., Foulquier, A., Poulenard, J., and Noble, V. (2020). Climate, soil resources and microbial activity shape the distributions of mountain plants based on their functional traits. *Ecography*, 43(10):1550–1559.
706 Michalet, R., Losapio, G., Kikvidze, Z., Brooker, R. W., Butterfield, B. J., Callaway, R. M., Cavieres, L. A., Lortie, C. J., Pugnaire, F. I., and Schöb, C. (2024). Assessing the accuracy of paired and random sampling for quantifying plant-plant interactions in natural communities. *Population Ecology*, 66(1):39-52.

710 Mihoč, M. A. K., Giménez-Benavides, L., Pescador, D. S., Sánchez, A. M., Cavieres, L. A., and Escudero, A. (2016). Soil under nurse plants is always better than outside: A survey on soil amelioration by a complete guild of nurse plants across a long environmental gradient. *Plant and Soil*, 408(1):31-41.

714 Mohebali, A., Erfanzadeh, R., and Jafari, M. (2023). Comparison of sub-canopy topsoil properties between three woody plant species, a case study of the Baladeh Watershed, Mazandaran Province, Iran. *Caspian Journal of Environmental Sciences*, 21(5):1017-1025.

718 Momberg, M. and le Roux, P. C. (2020). Testing for consistency in ecosystem engineering: Do cushion plants always turn up the heat? *Acta Oecologica*, 104:103532.

721 Morris, W. F. and Doak, D. F. (1998). Life History of the Long-Lived Gynodioecious Cushion Plant Silene acaulis (Caryophyllaceae), Inferred from Size-Based Population Projection Matrices. *American Journal of Botany*, 85(6):784-793.

724 Munroe, J. S. and Bockheim, J. G. (2001). Soil Development in Low-Arctic Tundra of the Northern Brooks Range, Alaska, U.S.A. *Arctic, Antarctic, and Alpine Research*, 33(1):78-87.

727 Murphy, M. V. (2022). *semEff: Automatic Calculation of Effects for Piecewise Structural Equation Models*. R package version 0.6.1.

729 Ochoa-Hueso, R., Eldridge, D. J., Delgado-Baquerizo, M., Soliveres, S., Bowker, M. A., Gross, N., Le Bagousse-Pinguet, Y., Quero, J. L., García-Gómez, M., Valencia, E., Arredondo, T., Beinticinco, L., Bran, D., Cea, A., Coaguila, D., Dougill, A. J., Espinosa, C. I., Gaitán, J., Guuroh, R. T., Guzman, E., Gutiérrez, J. R., Hernández, R. M., Huber-Sannwald, E., Jeffries, T., Linstädter, A., Mau, R. L., Monerris, J., Prina, A., Pucheta, E., Stavi, I., Thomas, A. D., Zaady, E., Singh, B. K., and Maestre, F. T. (2018). Soil fungal abundance and plant functional traits drive fertile island formation in global drylands. *Journal of Ecology*, 106(1):242-253.

738 Odling-Smee, F. J., Laland, K. N., and Feldman, M. W. (2003). *Niche Construction: The Neglected Process in Evolution*. Number 37. Princeton University Press.
741 Ohtsuka, T., Hirota, M., Zhang, X., Shimono, A., Senga, Y., Du, M., Yonemura, S., Kawashima, S., and Tang, Y. (2008). Soil organic carbon pools in alpine to nival zones along an altitudinal gradient (4400–5300m) on the Tibetan Plateau. *Polar Science*, 2(4):277–285.

745 Olsen, S. R. S. R. (1954). *Estimation of Available Phosphorus in Soils by Extraction with Sodium Bicarbonate*. Washington, D.C.: U.S. Dept. of Agriculture.

747 Ordoñez, J. C., Van Bodegom, P. M., Witte, J.-P. M., Wright, I. J., Reich, P. B., and Aerts, R. (2009). A global study of relationships between leaf traits, climate and soil measures of nutrient fertility. *Global Ecology and Biogeography*, 18(2):137–149.

751 Pinheiro, J. C. (2014). Linear Mixed Effects Models for Longitudinal Data. In Kenett, R. S., Longford, N. T., Piegorsch, W. W., and Ruggeri, F., editors, *Wiley StatsRef: Statistics Reference Online*. Wiley, 1 edition.

754 Pintaldi, E., D'Amico, M. E., Colombo, N., Colombero, C., Sambuelli, L., De Regibus, C., Franco, D., Perotti, L., Paro, L., and Freppaz, M. (2021). Hidden soils and their carbon stocks at high-elevation in the European Alps (North-West Italy). *CATENA*, 198:105044.

758 Praeg, N., Pauli, H., and Illmer, P. (2019). Microbial Diversity in Bulk and Rhizosphere Soil of Ranunculus glaucis Along a High-Alpine Altitudinal Gradient. *Frontiers in Microbiology*, 10.

761 R Core Team (2024). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria.

763 Rai, S., Breme, N., Jandova, V., Lanta, V., Altman, J., Ruka, A. T., Rixen, C., and Dolezal, J. (2024). Growth dynamics and climate sensitivities in alpine cushion plants: Insights from Silene acaulis in the Swiss Alps. *Alpine Botany*.

766 Roy, J., Albert, C. H., Choler, P., Clément, J.-C., Ibanez, S., Lavergne, S., Saccone, P., Zinger, L., and Geremia, R. A. (2013). Microbes on the Cliff: Alpine Cushion Plants Structure Bacterial and Fungal Communities. *Frontiers in Microbiology*, 4.

770 Schlesinger, W. H., Reynolds, J. F., Cunningham, G. L., Huenneke, L. F., Jarrell, W. M., Virginia, R. A., and Whitford, W. G. (1990). Biological feedbacks in global desertification. *Science (New York, N.Y.)*, 247(4946):1043–1048.

773 Schweitzer, J. A., Juric, I., van de Voorde, T. F. J., Clay, K., van der Putten, W. H., and Bailey, J. K. (2014). Are there evolutionary consequences of plant–soil feedbacks along soil gradients? *Functional Ecology*, 28(1):55–64.
Shuo, X., Yongtao, H., Ben, N., Xingliang, X., Qian, S., and Yingfan, W. (2024). Effect of a Cushion Plant, Androsace tapete, on Soil Net Nitrogen Mineralization and Enzyme Activities during the Growing Season. *Journal of Resources and Ecology*, 15(2).

Spasojevic, M. J. and Weber, S. (2021). Variation in δ13C and δ15N within and among plant species in the alpine tundra. *Arctic, Antarctic, and Alpine Research*, 53(1):340–351.

Stark, J., Lehman, R., Crawford, L., Enquist, B. J., and Blonder, B. (2017). Does environmental heterogeneity drive functional trait variation? A test in montane and alpine meadows. *Oikos*, 126(11):1650–1659.

Steinbauer, M. J., Beierkuhnlein, C., Arfin Khan, M. A., Harter, D. E., Irl, S. D., Jentsch, A., Schweiger, A. H., Svenning, J.-C., and Dengler, J. (2016). How to differentiate facilitation and environmentally driven co-existence. *Journal of Vegetation Science*, 27(5):1071–1079.

Thuiller, W., Saillard, A., Abdulhak, S., Augé, V., Birck, C., Bonet, R., Choler, P., Delestrade, A., Kunstler, G., Leccia, M.-F., Lienard, B., Poulenard, J., Valay, J.-G., Bayle, A., Bonfanti, N., Brousset, L., Bizard, L., Calderón-Sanou, I., Dentant, C., Desjonquères, C., Gielly, L., Guéguen, M., Guiter, F., Hedde, M., Hustache, E., Kedhim, N., Lapenu, P., Guillarme, N. L., Marchal, L., Mahieu, C., Martin, G., Martinez-Almoyna, C., Miele, V., Murienne, J., Paillet, Y., Rome, M., and Renaud, J. (2024). ORCHAMP: an observation network for monitoring Biodiversity and Ecosystem Functioning Across Space and Time in Mountainous Regions. *Comptes Rendus Biologies - Académie des sciences*.

Towett, E. K., Shepherd, K. D., and Cadisch, G. (2013). Quantification of total element concentrations in soils using total X-ray fluorescence spectroscopy (TXRF). *Science of The Total Environment*, 463–464:374–388.

van Breemen, N. and Finzi, A. C. (1998). Plant-soil Interactions: Ecological Aspects and Evolutionary Implications. *Biogeochemistry*, 42(1):1–19.

van den Boogaart, K. G., Tolosana-Delgado, R., and Bren, M. (2024). *compositions: Compositional Data Analysis*. R package version 2.0-8.

van Lissa, C. J. (2023). *tidySEM: Tidy Structural Equation Modeling*. R package version 0.2.6.

Walker, L. R., Thompson, D. B., and Landau, F. H. (2001). Experimental Manipulations of Fertile Islands and Nurse Plant Effects in the Mojave Desert, Usa. *Western North American Naturalist*, 61(1):25–35.
811

Wang, Y., Sun, J., Liu, B., Wang, J., and Zeng, T. (2021). Cushion plants as critical pioneers and engineers in alpine ecosystems across the Tibetan Plateau.

812

813

*Ecology and Evolution*, 11(17):11554-11558.

814

Whitham, T. G., Bailey, J. K., Schweitzer, J. A., Shuster, S. M., Bangert, R. K., LeRoy, C. J., Lonsdorf, E. V., Allan, G. J., DiFazio, S. P., Potts, B. M., Fischer, D. G., Gehring, C. A., Lindroth, R. L., Marks, J. C., Hart, S. C., Wimp, G. M., and Wooley, S. C. (2006). A framework for community and ecosystem genetics:

815

816

817

818

From genes to ecosystems. *Nature Reviews Genetics*, 7(7):510-523.

819

Zhao, R. (2022). The influence of cushion plants on soil fertility islands under their canopies with ontogenesis in an extreme alpine ecosystem. *Plant and Soil*.

820

821

Zhao, R. and An, L. (2021). Plant size of the alpine cushion Thylacospermum caespitosum affects soil amelioration at different elevations. *Plant Ecology*, 222(3):323-335.

822

823

824

Zumsteg, A., Bååth, E., Stierli, B., Zeyer, J., and Frey, B. (2013). Bacterial and fungal community responses to reciprocal soil transfer along a temperature and soil moisture gradient in a glacier forefield. *Soil Biology and Biochemistry*, 61:121-132.

825

826

827

# Statements and Declarations

828

## Fundings

829

830

This work was funded by a grant from the Universe Science Observatory of Grenoble (OSUG). KD was supported by a PhD grant from the University Savoie Mont Blanc. Part of the analysis for the Orchamps monitoring program came from the French Agence Nationale de la Recherche (ANR) through the GlobNets (ANR-16-CE02-0009) and from two ‘Investissement d’Avenir’ grants (Trajectories: ANR-15-IDEX-02; Montane: OSUG@2020: ANR-10-LAB-56).

831

832

833

834

835

## Competing interests

836

The authors declare no conflicts of interest.

837
## Author contributions

838 **Keyvan Dumas** : Conceptualization, Formal analysis, Investigation, Data Curation, Writing - Original Draft, Visualization, **Nicolas Bonfanti** : Conceptualization, Investigation (ORCHAMPS dataset), Writing - Review & Editing, **Aurélien Grange** : Investigation, **Annie Millery-Vigues** : Investigation, **The Orchamp Consortium** : Investigation, **Jérome Poulenard** : Conceptualization, Supervision, Writing - Review & Editing, Funding Acquisition, **Glenn Yannic** : Supervision, Writing - Review & Editing, **Christiane Gallet** : Supervision, Writing - Review & Editing, **Sébastien Lavergne** : Conceptualization, Supervision, Writing - Review & Editing, **Sébastien Ibanez** : Conceptualization, Supervision, Writing - Review & Editing, Funding Acquisition

## Data availability

850 Data associated with this study will be deposited on Zenodo during the reviewing process. For now, data about cushion plant and their soils, as well as R script for all analyses, figures and tables, are available at https://gitlab.com/dumaskv n/script_plso. Data concerning the soils studied in the Orchamps project are available at https://orchamp.osug.fr/

## Supporting Information

856 Fig.A1 : Principal Component Analysis describing soil geochemistry, and nutrient and organic matter composition of soils.

858 Fig.A2 : Description of the causal relationship expected *a priori* in the SEM model

859 Table A1 : Justification for each path hypothesized in Fig. A2

860 Table A2 : PERMANOVA analysis of cushion plant traits

861 Table A3 : Linear Mixed Model parameters for the altitudinal variation of soil characteristics

863 Table A4 : Linear Mixed Model parameters describing the variation of δ¹⁵N as
function of soil total nitrogen.

Table A5 : Summary statistics of linear mixed models exploring the inter-specific

effect of cushion plants on soil characteristics