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

Research

# Phenotypic effects of the nurse *Thylacospermum caespitosum* on dependent plant species along regional climate stress gradients

Xingpei Jiang, Richard Michalet, Shuyan Chen, Liang Zhao, Xiangtai Wang, Chenyue Wang, Lizhe An and Sa Xiao

X. Jiang (http://orcid.org/0000-0003-0228-1275), S. Chen, X. Wang, C. Wang, L. An (http://orcid.org/0000-0003-0228-1275) (jianghenchi@hotmail.com) and Sa Xiao (http://orcid.org/0000-0002-8427-752X), State Key Laboratory of Grassland and Agro-ecosystems, School of Life Sciences, Lanzhou Univ., 222 Tianshuinan Road, CN-730000 Lanzhou, PR China. – R. Michalet, Univ. of Bordeaux, U.M.R. CNRS, Allée Geoffroy Saint-Hilaire Pessac, France. – L. Zhao, Ministry of Education Key Laboratory of Three Gorges Reservoir Region's Eco-Environment, Chongqing Univ., Chongqing, PR China.

Oikos

127: 252–263, 2018

doi: 10.1111/oik.04512

Subject Editor: Miguel Verdú

Editor-in-Chief: Dries Bonte

Accepted 25 July 2017

Contrasting phenotypes of alpine cushion species have been recurrently described in several mountain ranges along small-scale topography gradients, with tight competitive phenotypes in stressful convex topography and loose facilitative phenotypes in sheltered concave topography. The consistency of phenotypic effects along large-scale climate stress gradients have been proposed as a test of the likely genetic bases of the differences observed at small-scale. Inversely, plastic phenotypic effects are more likely to vanish at some points along climate stress gradients. We tested this hypothesis for two phenotypes of the alpine cushion species *Thylacospermum caespitosum* at four points along regional gradients of cold and drought stress in northwest China. We measured the traits of the two cushion phenotypes and quantified their associated plant communities and environmental variables along the regional temperature and aridity gradients. Cushion height, convexity and stem density overall showed significant effect of phenotypes. Difference in tightness of cushions between phenotypes was consistent across climate conditions, whereas differences in cushion convexity and height between phenotypes increased with increasing cold stress. Phenotypic effects on species richness and abundance were consistent along both climate gradients but not effects on species composition, while there were no phenotypic effects on environmental variables. Additionally, RII (relative interaction index) curves were linear along the drought gradient but unimodal along the temperature gradient, likely due to the occurrence of contrasting species pools at the different sites. We conclude that the consistency of phenotypic effects of *T. caespitosum* was high for species richness and abundance and mainly explained by differences in interference mediated by likely heritable differences in cushion tightness. Additionally, our study shows that the shapes of the relationship between plant responses to neighbours and environmental stresses are not necessarily driven by niche-based deterministic factors.

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

Pioneer facilitation study has shown that phenotypic dif-
ferentiation within dominant species may alter plant–plant
interactions and drive community composition and diver-
sity (Callaway et al. 1991). More recently there has been a
growing interest in linking community ecology and evolu-
tionary sciences (Valiente-Banuet et al. 2006, Bronstein
2009, Thorpe et al. 2011, Ikeda et al. 2017), together with
an emergence of community genetics as a research discipline
(Whitham et al. 2003, 2006, Pakeman et al. 2006, Johnson
and Stinchcombe 2007, Lau et al. 2016). Community genet-
ics assesses the community consequences of genetic differen-
tiations within foundation species (sensu Ellison et al. 2005),
i.e. community phenotypes, their heritability and feed-back
effects on the genetics of foundation species, with an over-
all goal to improve our understanding of community evo-
lution (Whitham et al. 2006). Pioneer community genetics
studies have focused on multi-trophic interactions (reviewed
by Whitham et al. 2006, Bailey et al. 2009), but there is an
increasing interest in assessing in plant communities the con-
sequences of phenotypic differentiation within dominant
foundation species for dependent plant species (Pakeman
et al. 2006, Lankau and Strauss 2007, Crutsinger et al. 2010,
Michalet et al. 2011, Aschehoug and Callaway 2014, Ehlers
et al. 2016).

Alpine cushion plant systems have particularly stimulated
community genetics studies (Michalet et al. 2011, Al Hayek
et al. 2014, 2015a, b, Michalet et al. 2017), likely due to
the steepness of environmental gradients enhancing phe-
notypic differentiation (Choler and Michalet 2002, Körner
2003, Schöb et al. 2013, Bonanomi et al. 2016, Michalet
et al. 2016) and eventually ecological speciation (Choler et al.
2004). Alpine community genetics studies have shown how
small-scale variation in topography may induce heritable phe-
notypic differentiation within dominant cushion plants and
in turn contrasting effects on dependent plant species with
feedbacks on the fitness of the foundation cushion species
(Michalet et al. 2011, Al Hayek et al. 2014).

Interestingly, similar patterns and processes have been
described worldwide in different alpine climate conditions:
in northern Arizona (USA) with a subtropical continental
climate (Michalet et al. 2011), in western Pyrenees
(France) with a temperate oceanic climate (Al Hayek et al.
2014, 2015b) and in Mount Lebanon (Lebanon) with a
Mediterranean climate (Al Hayek et al. 2015a, Michalet
et al. 2017). In all these alpine or subalpine systems two
main phenotypes were observed in contrasting topogra-
phies, tight phenotypes with compact morphology and
high competitive effect on dependent species in dry and
unfertile convex topography and loose phenotypes with
loose morphology and high facilitative effect in wetter and
more fertile concave topography. In these different studies
the genetic basis of the phenotypic differentiation within
the dominant cushion plant and the contrasting pheno-
typic effects on other species have been assessed directly

in common-gardens controlling for environmental plastic-
ity effects only in the case of fast-growing cushion species
(Michalet et al. 2011, Al Hayek et al. 2015b). However,
most alpine cushion plants are very slow-growing species
and the genetic basis of their differences in competitive and
facilitative effects can be assessed in common-gardens only
beyond a threshold of size which might need several decades
to be reached. Al Hayek et al. (2015a) proposed an alter-
native indirect method for these particular species hardly
to be studied in common-gardens; this method assumes
that phenotypic effects are more likely to be inconsistent
along gradients of environmental stress when they are due
to environmental plasticity. In contrast, phenotypic effects
should be more consistent along climate stress gradients if
they have a genetic basis, as heritable phenotypes from con-
trasting habitats have consistent traits when transplanted in
common gardens (Clausen et al. 1940). Schöb et al. (2013)
and Bonanomi et al. (2016) found that phenotypic traits
and effects changed along elevation stress gradients for *Are-
naria tetraquetra* in the Sierra Nevada (Spain) and *Silene
acaulis* in the Apennines (Italy), respectively. Al Hayek et al.
(2015a) did not find a decrease in facilitation for the loose
phenotype of the cushion shrub *Onobrychis cornuta* along
a local gradient of climate drought stress due to exposure
and concluded that small-scale phenotypic effects of this
cushion shrub were very likely due to heritable differences
in traits between phenotypes. However, phenotypic effects
of the shrub were not fully consistent and rather increased
with increasing drought stress since the loose phenotype
was even more facilitative and the tight phenotype more
competitive, which suggests that phenotypic effects may
inversely vanish in low stress conditions.

Thus, our knowledge of the consistency of phenotypic
effects along stress gradients is still limited. In particular, there
has been no study assessing the consistency of phenotypic
effects both along regional climate stress gradients and with
several points along the gradients. Variation in competition
and facilitation have been proposed to follow either a linear
(e.g. Stress gradient hypothesis, Bertness and Callaway 1994),
or unimodal curve (e.g. the collapse of interactions, Michalet
et al. 2006, Holmgren and Scheffer 2010) and assessing the
consistency of phenotypic effects with likely nonlinear varia-
tions in competition and facilitation along environmental
gradients needs a design including at least three points along
an environmental gradient (Michalet et al. 2014a).

Our main goal in this study is to assess the consistency
of the phenotypic effects of the alpine cushion species
*Thylacospermum caespitosum* on dependent plant species
along regional climate stress gradients. Although loose flat
and tight dome-shaped phenotypes of the species have been
described in Ladakh (Dvorský et al. 2013) and northern
China (Michalet et al. 2016) suggesting the occurrence of
contrasting effects on other species, these phenotypic effects
have not been specifically quantified. An assessment of the
consistency of phenotypic effects along regional gradients of
climate stress should allow us to progress in our knowledge of
the consequences of phenotypic variation within community
dominants for plant–plant interactions and species richness
in alpine environments.

# Material and methods

## Species and sites

Our model species, *Thylacospermum caespitosum* (Caryophyllaceae), is one of the most drought- and cold-tolerant alpine cushion species of the world (de Bello et al. 2011). It is widely distributed, between 3600 m and 5900 m of elevation, from the south (dry western Himalayan Mountains in Ladakh, India and wetter southeast Tibetan ranges (de Bello et al. 2011, Chen et al. 2015) to the north (Qilian Mountains and Tianshan Mountains, in northwestern China, Michalet et al. 2016) of Central Asia.

Two phenotypes of *T. caespitosum* have been described, first in the southwest part of its range in Ladakh, India (Dvorský et al. 2013) and then north of the Qinghai-Tibet Plateau, in the Qilian Mountains, China (Michalet et al. 2016). In both mountain ranges the two phenotypes have similar contrasting morphologies, with either flat loose or dome-shaped tight cushions. Additionally, and similarly to what have been described for other alpine cushion species from different mountain ranges (Michalet et al. 2011, Al Hayek et al. 2014, 2015a, b), the loose cushions seem to be more facilitative for other alpine species and the tight cushions more competitive (Michalet et al. 2016). In our study sites, common subordinated species are *Elymus nutans*, *Draba altaica*, *Rhodiola juparensis*, *Poa sibirica* and *Festuca coelestis*.

In order to assess variation in phenotypic effects of *T. caespitosum* along climate stress gradients, we conducted late July 2015 a field study at four sites with contrasting climate conditions (Fig. 1). The sites were chosen in order to vary both temperature and rainfall, because several experiments

<figure><img src="image_3.png" /><figcaption>Figure 1. Location of the four sampling sites in northwest China. From west to east, Tianshan (TS, 43°7′N, 86°49′E, 3772 m a.s.l.); Qilian low (QL, 39°35′N, 96°25′E, 3696 m a.s.l.); Qilian high (QH, 39°30′N, 96°30′E, 4260 m a.s.l.); Maya (MY, 37°4′N, 102°40′E, 3800 m a.s.l.).</figcaption></figure>

conducted in dry alpine climates have shown that cold and drought stresses can produce contrasting variations in facilitation and competition along elevation gradients (Cavieres et al. 2006, Michalet et al. 2014b). China is subjected to subtropical monsoonal oceanic influences coming from the southeast and which strength decreases with increasing rainshadow towards the northwest, where the climate becomes much colder and drier. Thus, the warmest and wettest site Maya (hereafter MY, 37°4′N, 102°40′E, 3800 m a.s.l.), was located at low alpine elevation close to the eastern border of the Qilian Mountains (Qilianshan national nature reserve, Gansu Province). The driest site, Qilian low (hereafter QL, 39°35′N, 96°25′E, 3696 m a.s.l., Yanchiwan national nature reserve, Gansu Province), was also located at low alpine elevation but more west in the most continental part of the Qilian mountain range. The two coldest sites were located both at high alpine elevation. Qilian high (hereafter QH, 39°30′N, 96°30′E, 4260 m a.s.l.) was located in the same mountain range than QL but 560 m higher in elevation and Tianshan (hereafter TS, 43°7′N, 86°49′E, 3772 m a.s.l., Xinjiang Province) at the northwest end of our study zone in a slightly wetter climate than QH due to northwestern wet influences coming from Europe and western Asia. The climate of the sampling sites is overall semiarid continental (Wang 1983) with 167–528 mm (TS, 312 mm; QL, 167 mm; QH, 230 mm; MY, 528 mm) of annual rainfall and –8.3 – –2.3°C (TS, –7.2°C; QL, –3.9°C; QH, –8.3°C; MY, –2.3°C) of mean annual temperature (Hijmans et al. 2005). Among the four sites, only QL and MY suffer very slight grazing from yaks *Bos grunniens* and sheep (small tailed han sheep, *Ovis aries*).

## Data collection

To quantify variation in the strength of associations between the cushions of the two phenotypes of *T. caespitosum* and the dependent species along climate stress gradients, we randomly selected, at each of the four sites and for each phenotype, 40 cushion replicates and 40 paired adjacent open areas. Cushion individuals were considered as tight or loose at each site depending on their shape, dome-shaped or flat. Phenotype was strongly associated to microtopography since tight phenotypes were mostly present in convex topographies and loose phenotypes in concave topographies (Michalet et al. 2011).

For all cushion individuals, we measured height (H, using two perpendicular rulers), maximum (a) and minimum (b, perpendicular to maximum axis) axes of their area and calculated the approximate cushion area ($A = a \times b$). In order to assess the tightness of each cushion, we quantified the cushion stem density, by counting the number of stems in a $2 \times 2$ cm² wire frame, which was randomly located on the cushion surface. We also calculated an index of cushion convexity (i.e. ICC):

$$ICC = H/A$$

Where H is the height of the cushion in cm and A is the area of the cushion in dm².
Meanwhile, two phenotypes cushion leaves were randomly collected (stored in an incubator with ice bags) and took back to measure physiological traits: leaf fresh weight (LFW), leaf area (LA), leaf dry weight (LDW), leaf total organic carbon (LTOC), leaf total nitrogen (LTN). Leaf area was measured with open source software ImageJ (Schindelin et al. 2015), leaf total organic carbon was determined according to an improved K₂Cr₂O₄–H₂SO₄ method (Wu and Tao 1993) and leaf total nitrogen was digested by concentrated H₂SO₄ at 120°C for 10 min, followed by semi micro-Kjeldahl method, using an auto chemistry analyzer. Then, we calculated leaf dry matter content (LDMC, g g⁻¹, the ratio LDW/LFW), specific leaf area (SLA, g cm⁻², the ratio LA/LDW) and the ratio between leaf total carbon and nitrogen (i.e. LCN, the ratio LTOC/LTN).

We then haphazardly selected in the close vicinity of each cushion (less than a meter) and in similar micro-topographical conditions (i.e. concave open plots paired with loose cushions and convex open plots paired with tight cushions) an open plot of the same area. We counted all individuals of each vascular species growing within the cushions and within the open paired plots.

We also sampled soils in five replicates of each combination of the site by phenotype by cushion treatments. Cushions and open plots were randomly selected and soil samples were collected 5–10 cm deep at least three days after a rain event. Soil samples were fully mixed in plastic bags, stored in an incubator with ice bags and transferred to the lab for chemical analysis. Soil water content (SWC) was quantified by gravimetric method, by drying 5 g of fresh soil at 105°C for 24 h. The remaining soil samples were air-dried and then sieved with 100 meshes (0.15 mm). Soil pH was measured with a pH meter using a 1:2 ratio of soil to deionized water. Soil total organic carbon (STOC) was quantified following the Walkley–Black chromic acid wet oxidation method (Agricultural Ministry of P R China 2006). Soil total nitrogen (STN) and phosphorus (STP) were both digested by concentrated H₂SO₄ at 375°C for 3 h and 45 min, respectively, followed by semi micro-Kjeldahl and Mo–Sb antispetrophotography method, using an auto chemistry analyzer.

## Climatic data and plant–plant interaction index

We calculated mean temperature (T, °C) and mean precipitation (P, mm) during the growing season (June–August) based on monthly climatic data obtained from the Worldclim database (<www.worldclim.org>; Hijmans et al. 2005, ArcGIS ESRI 2014) with a spatial resolution of 30 seconds (ca 1 km²).

To quantify the effect of cushions on abundance and richness of dependent species we used the relative interaction index (RII) (Armas et al. 2004):

$$RII_{richness/abundance} = (N_{cushion} - N_{open})/(N_{cushion} + N_{open})$$

where N_cushion and N_open are the richness or abundance of non-cushion species in cushion and paired open areas, respectively;

RII values range from –1 to 1, positive values indicating facilitation, zero indicating neutral effect and negative values indicating competition.

## Statistical analyses

We used a two-way analysis of variance to assess the effects of the site, phenotype treatments and their interactions on cushion traits (area, height, ICC, stem density, SLA, LDMC, LCN). To explore variation in phenotypic effects on richness and abundance of understory species along climate stress gradient (temperature and drought), we used two-way analysis of covariance (ANCOVA) models, with the mean summer temperature/ precipitation (T/P, covariate) and phenotype (factor) as factors and the RII_richness and RII_abundance as response variables (Maalouf et al. 2012). Quadratic terms of covariates (T², P²) were also adopted for potential curved responses of RII_richness and RII_abundance along the temperature and drought gradients. One sample t-tests was used to evaluate if RII values were significantly different from zero and univariate analyses were followed by Tukey tests when necessary.

We used principal coordinates analysis (PCoA, Gower 1966) to assess at each site differences in species composition between loose and tight habitats and between cushions and open plots. Species with occurrence lower than 3 at each site (160 plots) were excluded from species matrix before analysis. At last, MY, QL, QH and TS sites covered 23, 14, 15, 9 species, respectively (Supplementary material Appendix 1 Table A5). To further evaluate phenotypic effects on understory community composition we also conducted a npMANOVA with permutation of 9999 times to assess effects of the two treatments (cushion and habitat) and their interaction at each site. Because both phenotypes occur in different environmental conditions we had to eliminate the effects caused by the difference in habitats between the two phenotypes to really assess phenotypic effects on species composition. Therefore, we defined the difference between the two habitats as habitat effect, the difference between cushion and open plots as cushion effect, phenotypic effects being revealed by the occurrence of a significant interaction between the habitat and cushion treatments.

We assessed differences in environmental conditions between sites, between habitats and between cushion and open plots using principal component analysis (PCA) on soil (SWC, STOC, STP, C/N ratio (the ratio of STOC/STN) and pH, Supplementary material Appendix 1 Table A1) and climate (T and P) variables. We also used npMANOVA with permutation of 9999 times to assess effects of the three treatments (site, cushion and habitat) and their interactions.

All dependent variables (except RII) were log-transformed before univariate analyses. All statistical analyses were conducted with R ver. 3.3.0 (<www.r-project.org>), with the module vegan (Oksanen et al. 2016), stats (<www.r-project.org>) for multivariate analyses, and ggplot2 (Wickham 2009) for drawings.
## Data deposition

Data available from the Dryad Digital Repository: <http://dx.doi.org/10.5061/dryad.nm13c> (Jiang et al. 2017).

## Results

Cushion area only differed across sites with the largest cushions observed at the warmest and wettest site (MY), intermediate cushions at the warmest and driest site (QL) and the smallest cushions at the two coldest sites (QH and TS) (highly significant site effect, Fig. 2a). The tallest cushions were also observed at MY, with no difference in cushion height among the other three sites (highly significant site effect, Fig. 2b). Additionally, tight phenotypes were taller than loose ones at all sites, but differences between the two phenotypes were stronger at QH, the coldest site,

than at the two warmest sites (MY and QL), with intermediate difference at TS (highly significant phenotype effect and significant site by phenotype interaction, respectively, Fig. 2b). Similar interaction among treatments were observed for the convexity index (ICC) with the highest difference in convexity between the two phenotypes at QH, intermediate difference at TS and the lowest difference at MY and QL (Fig. 2c). However, overall phenotypes the most convex cushions were observed at TS with no difference in cushion convexity among the other three sites (highly significant site effect, Fig. 2c). In contrast, stem density was only affected by the phenotype and site treatments, with overall higher stem density for tight than loose cushions and the highest stem density at QH, intermediate values at QL and TS and the lowest stem density at MY (Fig. 2d). Overall these results show that the difference in tightness of cushions between phenotypes was rather consistent across sites and thus climate conditions but that phenotypic effects on cushion shape were

<figure><img src="image_4.png" /><figcaption>Figure 2. Traits (Mean ± SE, n = 40) of the two cushion phenotypes of *T. caespitosum* at the four sites: (a) cushion area (dm²); (b) cushion height (cm); (c) index of cushion convexity; (d) stem density (no. stems cm⁻²). Significant effects of the site, phenotype treatments and their interactions in the two-way ANOVAs are indicated above each panel. Capital letters give the results of the Tukey tests on the effect of the site treatment and small letters give the results of the Tukey tests on the effect of the site by phenotype interaction treatment when significant. *, p < 0.05; **, p < 0.01; ***, p < 0.001.</figcaption></figure>
slightly inconsistent, with differences in cushion convexity and height between phenotypes increasing with increasing cold stress. There was a highly significant site effect on LDMC with higher LDMC at the driest site (QL) (Supplementary material Appendix 1 Table A6). For SLA there was only a significant site effect with higher SLA at the two coldest sites (QH and TS) (Supplementary material Appendix 1 Table A6). QL had higher LCN values than the other three sites (site effect in Supplementary material Appendix 1 Table A6) while loose cushions had higher LCN than tight cushions (phenotype effect in Supplementary material Appendix 1 Table A6).

The main result of the ANCOVAs on variation in RII along climate stress gradients was the occurrence of a highly significant phenotype effect for both abundance and richness and along both drought and cold gradients (Supplementary material Appendix 1 Table A2). This was explained by the consistency of the higher positive effect or lower negative effect

(depending on the sites) of the loose phenotype as compared to the tight phenotype along the stress gradients (Fig. 3). The second result was the contrasting shapes of the curves of increasing facilitation with increasing climate stress depending on the considered stress. Along a gradient of increasing cold stress the curves were overall unimodal, due to the occurrence of a collapse of facilitation in the coldest site (QH), as shown by the more significant quadratic than linear temperature effect (Supplementary material Appendix 1 Table A2, Fig. 3a, c). In contrast, along a gradient of increasing drought stress the curves were overall linear (except RII<sub>abundance</sub> for tight cushion), due to the highest facilitation found at the driest site (QL), as shown by the more significant linear than quadratic P effect (Supplementary material Appendix 1 Table A2, Fig. 3b, d). There was only a weakly significant T × phenotype interaction observed for RII<sub>abundance</sub>, due a slightly more linear curve observed for the tight than for the loose phenotype (Supplementary material Appendix 1 Table A2, Fig. 3c).

<figure><img src="image_5.png" /><figcaption>Figure 3. Variation of mean RII<sub>richness</sub> (upper panels) and mean RII<sub>abundance</sub> (lower panels) along the temperature (left panels) and drought stress (right panels) gradients. Grey quadratic regression curves are fitted by the raw data of loose cushions. Asterisks above or below dots give the significance of the one sample t-test on RII means at each site. Significances of effects of treatments in the analysis of covariance (ANCOVA) are indicated above each panel: (*), p < 0.1; *, p < 0.05; **, p < 0.01; ***, p < 0.001. Detailed results of ANCOVA are given in Supplementary material Appendix 1 Table A2.</figcaption></figure>
The results of PCoA showed that there was a significant interaction between the cushion and habitat treatments, i.e. a phenotypic effect on species composition, only at MY and QL (Fig. 4a–b, respectively). At MY, the wettest and warmest site, this phenotypic effect was observed because cushion effects of the two phenotypes occurred on different axes (i.e. cushion effects of loose cushions occurred on the first axis whereas cushion effects of tight cushions occurred on the second axis). RII results showed that the effect of the tight phenotype for abundance and richness was negative at MY (Fig. 2a, c) and this negative cushion effect on species composition was likely due to a high frequency and abundance of pioneer species between tight cushions at this mild site (e.g. *Potentilla fragarioides*, *Cerastium fontanum*, *Viola biflora*, *Corydalis hamata*, *Stellaria arenaria* at the positive end of axis 2, Fig. 4a). In contrast, in QL, the driest site, this phenotypic

effect was due to a stronger cushion effect in the loose than in the tight habitat observed on axis 1. RII results showed that the effect of the loose phenotype was positive at QL and this stronger cushion effect of the loose phenotype was likely due to a high frequency and abundance of stress-intolerant and competitive alpine species, among them many graminoids, (e.g. *Leontopodium nanum*, *Kobresia pygmaea*, *Kobresia royleana*, *Festuca coelestis* at the negative side of axis 1, Fig. 4b) within loose cushions at this stressed site. At QH and TS, there was only a significant cushion effect (Fig. 4c–d, respectively).

The first two PCA axes explained 38% and 21% of the total variance in environmental variables among sites (see the details of soil properties in Supplementary material Appendix 1 Table A1), habitats and cushions (Fig. 5). PC1 was mainly explained by cold stress and soil phosphorus,

<figure><img src="image_6.png" /><figcaption>Figure 4. Principal coordinates analysis (PCoA) ordination of plant communities in *T. caespitosum* cushions and open areas in each habitat (tight and loose) at each of the four sites: (a) MY; (b) QL; (c) QH; (d) TS. Different symbols indicate different habitats while solid and empty symbols indicate cushion and open, respectively. Results of treatment effects in the npMANOVA are given in the lower left of each panel: (*), *p* < 0.1; *, *p* < 0.05; **, *p* < 0.01; ***, *p* < 0.001. The results of npMANOVA (permutation = 9999) are shown in Supplementary material Appendix 1 Table A3. Four-italic letters are the abbreviation of species listed in appendix (Supplementary material Appendix 1 Table A5).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5. Principal component analysis (PCA) on environmental variables (SWC (%), soil water content; STOC (mg g⁻¹), soil total organic carbon content; STP (mg g⁻¹), soil total phosphorus content; C/N, ratio of total organic carbon content and total nitrogen content; pH; P (mm), mean summer precipitation; T (°C), mean summer temperature) at the four sites: MY; QL; QH and TS. Different symbols indicate different habitats, whereas solid and empty symbols represent cushion and open area conditions. Dashed ellipses delimit sites (and cushion and open plots at QH only) and arrows give the contributions of environmental variables. PCA axis 1 accounts for 38.25% of the multivariate variation and axis 2 for 21.32%. Asterisks refer to the significance of non-parametric multivariate analysis of variance (npMANOVA, permutation = 9999): ***, *p* < 0.001, **, 0.001 < *p* < 0.01, *, 0.01 < *p* < 0.05, (*), 0.05 < *p* < 0.1. The results of npMANOVA are showed in Supplementary material Appendix 1 Table A4.</figcaption></figure>

opposing the two warmest sites (MY and QL) in its positive part to the two coldest ones (TS and QH) in its negative part. PC2 was mainly explained by drought stress and C/N ratio, mainly opposing the driest site (QL) in its negative part to the three others in its positive part. Thus, there was a highly significant site effect in the npMANOVA (Supplementary material Appendix 1 Table A6). There was a highly significant cushion effect because the open plots were overall colder and drier than cushion plots, as suggested by their more negative scores than cushion plots on axes 1 and 2, respectively. There was also a weakly significant habitat effect because tight plots were also colder and drier than loose ones, as suggested by their also slightly more negative scores than loose plots on axes 1 and 2, respectively. Finally, there was a highly significant site by cushion interaction in the npMANOVA because the overall wetter and warmer environmental conditions of the cushion plots as compared to the open ones were not observed at TS (Fig. 5). There were no significant phenotypic effects on environmental variables, neither variation in phenotypic effects across sites (no significant habitat by cushions interactions, Fig. 5), which show that there were, both overall and at each site, no stronger positive effects of the loose

phenotypes than of the tight phenotypes on environmental conditions.

## Discussion

### The consistency of phenotypic effects

All measured cushion traits of *Thylacospermum caespitosum* were significantly affected by climate stress gradients and in particular cold stress, with cushions being larger and taller at the mildest low alpine and wet site, MY, and more convex and tighter at the two coldest high alpine sites, QH and TS. Intraspecific changes in plant functional traits along environmental gradients are widespread (Albert et al. 2010) and can be due either to phenotypic plasticity or adaptation (Nicotra et al. 2010). Several alpine cushion plants have been shown to become more convex and tighter as cold stress increases with increasing elevation (Schöb et al. 2013, Bonanomi et al. 2016). These changes are consistent with the hypothesis that the cushion architecture allows plants to better adapt to the adverse alpine abiotic conditions, in particular through mitigating extreme temperatures inside the cushion and increasing soil moisture beneath it (Körner 2003, Cavieres et al. 2006, Anthelme et al. 2014). At the same time and at each site, cushions of *T. caespitosum* become more convex and tighter from the more sheltered and mild loose habitat to the more exposed and stressful tight habitat, as observed for other alpine cushion species displaying similar phenotypic differentiation between loose and tight phenotypes (Michalet et al. 2011, Al Hayek et al. 2014, 2015a). Cross-transplantation experiments and reciprocal common-gardens allow disentangling environmental plasticity effects from heritable effects, which are both likely to be involved in such phenotypic differentiation along small- or large-scale environmental gradients (Al Hayek et al. 2015b). However, because most alpine cushion plants are too slow-growing to be tested in such design, Al Hayek et al. (2015a) proposed that a consistency of phenotypic effects observed along climate stress gradients suggests a likely heritability of phenotypic differentiation occurring at small scale. Following this argumentation, in our experiment only the difference in tightness between loose and tight cushions but not the difference is convexity was consistent across the four sites with contrasting climate conditions, which suggests that only the former trait was likely heritable.

At the same time, we found that phenotypic effects on dependent species were consistent along both climate stress gradients only for species richness and abundance but not for species composition. For species composition, phenotypic effects occurred only at the two most extreme sites along the drought gradient, MY and QL, but were not observed at the two other sites, QH and TS. Phenotypic effects on species composition, i.e. community phenotypes, have rarely been observed in community genetics studies on plants. Al Hayek et al. (2015a) found in Lebanon very weak phenotypic effects of the alpine cushion species *Onobrychis cornuta* on species composition of the dependent community but strong
effects on species abundance and richness. Most evidence of
community phenotypes has come from interactions across
trophic levels, and in particular for plant–insect herbivore
interactions. In the case of plant–plant interactions, the
occurrence of true community phenotypes has been observed
when phenotypic effects are due to contrasting leaf chem-
istry (Pakeman et al. 2006, Lankau and Strauss 2007), but
not to contrasting plant architectures (Callaway et al. 1991,
Crutsinger et al. 2010, Michalet et al. 2011).

Interestingly, we did not find any phenotypic effects on environmental conditions. Our PCA result showed that both cushion types mitigated both cold and drought stresses at three of the four sites (not in TS, likely because of the young and poor soils of this site located close to a glacier), but these positive effects were not different among phenotypes. There was a significant habitat effect because the loose habitat was milder than the tight one at all sites, but these milder conditions of the loose habitat were not due to contrasting cushion effects since the same difference in environmental conditions between habitats was observed both in the open plots and within cushions. These differences in environmental conditions between habitats were consistent with the contrasting topographies of the two phenotypes, with the tight phenotype more exposed to drought stress in convex topography than the loose one occupying more concave slopes (Michalet et al. 2016).

We did not measure the likely effects of cushions on microclimate, although contrasting effects of alpine cushion species have often been shown to be explained by subtle differences in air temperature and humidity at cushion surfaces (Cavieres et al. 2006, Molenda et al. 2012, Schöb et al. 2013). Bonanomi et al. (2016) have shown that cushions of *Silene acaulis* had varying effects on microclimate related to changes in phenotypes along an elevation gradient in central Italy. However, these changes in cushions effects on microclimate were not related to biotic effects since in their study facilitation peaked at medium elevation, whereas the most convex cushions with the most important positive effects on temperature occurred at the highest elevation. They concluded that the collapse of facilitation observed at the highest elevation was due to the importance of interference effects (due the tightness of the cushions) occurring at this particularly stressful site and balancing the positive cushion effects on temperature. The net effects of neighbours measured in plant-plant interaction studies are known to be the balance of hidden positive and negative effects that are not necessarily quantified in the field (Callaway and Walker 1997).

Interestingly, Callaway et al. (1991) also found in
California (USA) that, although not varying in direct posi-
tive effects on environmental conditions, two phenotypes of
*Quercus douglasii* had contrasting net effects on understorey
species, because negative trees had strong interference effects
due to shallow roots counterbalancing their positive effects
on microclimate and soil resources. Michalet et al. (2017)
have shown in Lebanon that the negative effects of the tight
phenotype of the alpine cushion shrub *Onobrychis cornuta*
was mainly explained by a high interference litter effect due

to slower litter decomposition rate on convex topography
than in concave topography occupied by the positive loose
phenotype. They argued that such interference effects are also
likely involved in the competitive ability of other described
tight phenotypes of alpine species growing on dry and nutri-
ent-poor soils from rocky outcrops (Michalet et al. 2011, Al
Hayek et al. 2015a). In the case of typical very tight cushions
like *T. caespitosum* in our study, *Silene acaulis* in the study
of Bonanomi et al. (2016) or Azorella species in the Andes
(Arroyo et al. 2003), the interference effect is likely due to the
difficulty of seeds to germinate and seedlings to root between
the stems of the cushions, an effect named the “cushion qual-
ity effect” by Dvorský et al. (2013) in their study on *T. caespis-
tosum*. Cushion quality has also been proposed to be a main
driver of collapse of facilitation by cushion species (Pugnaire
et al. 1996, Schöb et al. 2013) and the lack of facilitation of
tight cushions might also be explained by a decrease in the
vigor of benefactor cushions. However, our measurements of
cushion leaf traits showed that tight cushions had lower C/N
than loose cushions. This suggests that tight cushions were
more healthy and fast-growing than loose cushions, a func-
tional evidence rather supporting the interference hypothesis
of Bonanomi et al. (2016).

The shapes of the variation in effects of
*T. caespitosum* along climate stress gradients

Interestingly, for both species richness and abundance we
found a unimodal relationship between RII and tempera-
ture stress and a linear curve along the drought gradient, the
former supporting the collapse of facilitation model (Michalet
et al. 2006, 2014a, Holmgren and Scheffer 2010) and the latter
supporting the Stress gradient hypothesis (SGH, Bertness and
Callaway 1994). Contrasting response curves along environ-
mental gradients have been found in plant–plant interaction
experiments and these discrepancies have induced an impor-
tant and still not solved debate in the literature (Maestre et al.
2009, Callaway 2013, He and Bertness 2014, Soliveres and
Maestre 2014, Michalet et al. 2014a, Michalet and Pugnaire
2016, Liancourt et al. 2017). The most accepted explana-
tion is that unimodal curves should be more frequent along
resource or disturbance gradients, due to the switchback of
competition in resource-poor conditions or to the collapse
of facilitation in highly disturbed conditions, respectively,
whereas linear increases in facilitation should be more frequent
along non-resource gradients such as temperature gradients
(Michalet 2007, Maestre et al. 2009, Michalet et al. 2014a).

However, these arguments cannot explain the different curves
observed in our study, because we found linear curves along the
drought gradient (a resource gradient) and unimodal curves
along the temperature gradient (a direct factor). The linear curves
observed along the drought gradient were obviously due to the
two extreme cushion effects found in our study 1) the highly
competitive effect of the tight phenotype at the wettest site, MY
and 2) the highly facilitative effect of the loose phenotype at the
driest site, QL. In contrast, when QL was positioned at interme-
diate position along the temperature gradient and QH and TS at
the cold extreme, the curves turned to unimodal, because facilitation was lower at QH and TS than at QL. Interestingly, results on species composition have shown that phenotypic effects were only observed at these two extreme sites, mostly because of the abundance of both pioneer ruderal species in the open plots exacerbating negative effects of the tight phenotypes at MY and stress-intolerant competitive graminoid species in cushions plots exacerbating positive effects of the loose phenotypes at QL. In contrast species richness was lower in the two coldest sites, in particularly at TS, which might explain lower RII values than at the two other sites. Al Hayek et al. (2015a) also showed in Lebanon that the increase in competitive effects of the tight cushions of *Onobrychis cornuta* with increasing drought stress from north to south exposure was rather explained by a higher species pool of ruderal species preferring open plots in south than in north exposure. Species pools effects have been also proposed in other systems for explaining similar collapses of facilitation, such as the US salt marshes (Pennings et al. 2003). This suggests that variation in plant-plant interactions along environmental gradients are not necessarily driven by deterministic niche-based mechanisms but might also be influenced by stochastic regional factors. In other words, even in the case of strong potential effects of neighbours due to particular favorable environmental conditions (such as high cold stress enhancing positive effects) and nurse traits, interactions may collapse because of a too low species pool and the failure of the recipient community to recording this effect (Michalet et al. 2015, Liancourt et al. 2017).

## Conclusion

We conclude that phenotypic effects of *T. caespitosum* were overall consistent along both regional stress gradients and influencing mainly the abundance and richness of dependent alpine plants, but weakly the community composition because of contrasting species pools across sites. Additionally, these biotic phenotypic effects appeared not to be related to contrasting phenotypic effects on environmental conditions but rather explained by differences in tightness between phenotypes and thus contrasting interference effects. Our study increases our knowledge on the consistency of phenotypic effects with changing environmental conditions along regional stress gradients, and highlights after previous studies, the importance of interference effects in stressful habitats of alpine cushion species. Our work also shed new light on the ecological factors shaping variation in plant responses to neighbours along environmental gradients, and in particular on the importance of species pools in regional studies.

*Acknowledgements* – We thank Peng Guo, Zilong Chen, Meiqin Xie, Yingfeng Ding, Letao Bo and Aifeng Guo for helping during field work.

*Funding* – The project was supported by the State Key Program of the National Natural Science of China (31230014, 41430749), the National Natural Science Foundation of China (31470492, 31670435, 41671038, 31670437), the Central University Special Fund (Izujbky-2016-87, Izujbky-2014-k09, Izujbky-2014-200,

Izujbky-2016-81, 106112016CDJXY210006), the Program for New Century Excellent Talents in University (NCET-13-0265) and the Gansu Province Natural Science Foundation (1606RJZA071).

## References

Agricultural Ministry of PR China. 2006. Method for determination of soil organic matter. – In: Soil testing, part 6 (in Chinese). Agricultural Ministry of PR China, Beijing, China.

Al Hayek, P. et al. 2014. Phenotypic differentiation within a foundation grass species correlates with species richness in a subalpine community. – Oecologia 176: 533–544.

Al Hayek, P. et al. 2015a. Differential effects of contrasting phenotypes of a foundation legume shrub drive plant-plant interactions in a Mediterranean mountain. – J. Veg. Sci. 26: 373–384.

Al Hayek, P. et al. 2015b. Disentangling the heritable and plastic components of the competitive and facilitative effects of an alpine foundation species. – J. Ecol. 103: 1172–1182.

Albert, C. H. et al. 2010. Intraspecific functional variability: extent, structure and sources of variation. – J. Ecol. 98: 604–613.

Anthelme, F. et al. 2014. Facilitation among plants in alpine environments in the face of climate change. – Front. Plant Sci. 5: 387.

ArcGIS ESRI 2014. Release 10.3 Environmental Systems Research Institute. – Redlands, CA.

Armas, C. et al. 2004. Measuring plant interactions: a new comparative index. – Ecology 85: 2682–2686.

Arroyo, M. T. K. et al. 2003. Positive associations between the cushion plant *Azorella monantha* (Apiaceae) and alpine plant species in the Chilean Patagonian Andes. – Plant Ecol. 169: 121–129.

Aschehoug, E. T. and Callaway, R. M. 2014. Morphological variability in tree root architecture indirectly affects coexistence among competitors in the understory. – Ecology 95: 1731–1736.

Bailey, J. K. et al. 2009. From genes to ecosystems: a synthesis of the effects of plant genetic factors across levels of organization. – Phil. Trans. R. Soc. B. 364: 1607–1616.

Bertness, M. D. and Callaway, R. 1994. Positive interactions in communities. – Trends Ecol. Evol. 9: 191–193.

Bonanomi, G. et al. 2016. Cushion plant morphology controls biogenic capability and facilitation effects of *Silene acaulis* along an elevation gradient. – Funct. Ecol. 30: 1216–1226.

Bronstein, J. L. 2009. The evolution of facilitation and mutualism. – J. Ecol. 97: 1160–1170.

Callaway, R. M. 2013. Life at the edge, cooperation in Antarctica. – J. Veg. Sci. 24: 417–418.

Callaway, R. M. and Walker, L.R. 1997. Competition and facilitation: a synthetic approach to interactions in plant communities. – Ecology 78: 1958–1965.

Callaway, R. M. et al. 1991. Facilitation and interference of *Quercus douglasii* on understory productivity in central California. – Ecology 72: 1484–1499.

Cavieres, L. A. et al. 2006. Positive interactions between alpine plant species and the nurse cushion plant *Laretia acaulis* do not increase with elevation in the Andes of central Chile. – New Phytol. 169: 59–69.

Chen, J. et al. 2015. Cushion plants can have a positive effect on diversity at high elevations in the Himalayan Hengduan Mountains. – J. Veg. Sci. 26: 768–777.
Choler, P. and Michalet, R. 2002. Niche differentiation and distribution of *Carex curvula* along a bioclimatic gradient in the southwestern Alps. - J. Veg. Sci. 13: 851-858.

Choler, P. et al. 2004. Genetic introgression as a potential to widen a species' niche: insights from alpine *Carex curvula*. - Proc. Natl Acad. Sci. USA 101: 171-176.

Clausen, J. et al. 1940. Experimental studies on the nature of species. I. Effect of varied environments on western North American plants. - Carnegie Inst. of Washington Publ. 520, Washington DC, USA.

Crutsinger, G. M. et al. 2010. Genetic variation within a dominant shrub species determines plant species colonization in a coastal dune ecosystem. - Ecology 91: 1237-1243.

de Bello, F. et al. 2011. Cushions of *Thylacospermum caespitosum* (Caryophyllaceae) do not facilitate other plants under extreme altitude and dry conditions in the northwest Himalayas. - Ann. Bot. 108: 567-573.

Dvorský, M. et al. 2013. Testing the stress-gradient hypothesis at the roof of the world: effects of the cushion plant *Thylacospermum caespitosum* on species assemblages. - PLoS One 8: e53514.

Ehlers, B. K. et al. 2016. Intraspecific genetic variation and species coexistence in plant communities. - Biol. Lett. 12: 20150853.

Ellison, A. M. et al. 2005. Loss of foundation species: consequences for the structure and dynamics of forested ecosystems. - Front. Ecol. Environ. 3: 479-486.

Gower, J. C. 1966. Some distance properties of latent root and vector methods used in multivariate analysis. - Biometrika 53: 325-338.

He, Q. and Bertness, M.D. 2014. Extreme stresses, niches, and positive species interactions along stress gradients. - Ecology 95: 1437-1443.

Hijmans, R. J. et al. 2005. Very high resolution interpolated climate surfaces for global land areas. - Int. J. Climatol. 25: 1965-1978.

Holmgren, M. and Scheffer, M. 2010. Strong facilitation in mild environments: the stress gradient hypothesis revisited. - J. Ecol. 98: 1269-1275.

Ikeda, D. H. et al. 2017. Genetically informed ecological niche models improve climate change predictions. - Global Change Biol. 23: 164-176.

Jiang, X. et al. 2017. Data from: Phenotypic effects of the nurse *Thylacospermum caespitosum* on dependent plant species along regional climate stress gradients. - Dryad Digital Repository, <http://dx.doi.org/10.5061/dryad.nm13c>.

Johnson, M. T. J. and Stinchcombe, J. R. 2007. An emerging synthesis between community ecology and evolutionary biology. - Trends. Ecol. Evol. 22: 250-257.

Körner, C. 2003. Alpine plant life: functional plant ecology of high mountain ecosystems. - Springer.

Lankau, R. A. and Strauss, S.Y. 2007. Mutual feedbacks maintain both genetic and species diversity in a plant community. - Science 317: 1561-1563.

Lau, M. K. et al. 2016. Genotypic variation in foundation species generates network structure that may drive community dynamics and evolution. - Ecology 97: 733-742.

Liancourt, P. et al. 2017. SGH: stress or Strain gradient hypothesis? Insights from an elevation gradient on the roof of the world. - Ann. Bot. 120: 29-38.

Maalouf, J.-P. et al. 2012. The interplay of stress and mowing disturbance for the intensit and importance of plant interactions in dry calcareous grasslands. - Ann. Bot. 110: 821-828.

Maestre, F. T. et al. 2009. Refining the stress-gradient hypothesis for competition and facilitation in plant communities. - J. Ecol. 97: 199-205.

Michalet, R. 2007. Highlighting the multiple drivers of change in interactions along stress gradients. - New Phytol. 173: 3-6.

Michalet, R. and Pugnaire, F. I. 2016. Facilitation in communities: underlying mechanisms, community and ecosystem implications. - Funct. Ecol. 30: 3-9.

Michalet, R. et al. 2006. Do biotic interactions shape both sides of the humped-back model of species richness in plant communities? - Ecol. Lett. 9: 767-773.

Michalet, R. et al. 2011. Phenotypic variation in nurse traits and community feedbacks define an alpine community. - Ecol. Lett. 14: 433-443.

Michalet, R. et al. 2014a. Two alternatives to the stress-gradient hypothesis at the edge of life: the collapse of facilitation and the switch from facilitation to competition. - J. Veg. Sci. 25: 609-613.

Michalet, R. et al. 2014b. Partitioning net interactions among plants along altitudinal gradients to study community responses to climate change. - Funct. Ecol. 28: 75-86.

Michalet, R. et al. 2015. Competition, facilitation and environmental severity shape the relationship between local and regional species richness in plant communities. - Ecography 38: 335-345.

Michalet, R. et al. 2016. Beneficiary feedback effects on alpine cushion benefactors become more negative with increasing cover of graminoids and in dry conditions. - Funct. Ecol. 30: 79-87.

Michalet, R. et al. 2017. Direct litter interference and indirect soil competitive effects of two contrasting phenotypes of a spiny legume shrub drive the forb composition of an oromediterranean community. - Oikos 126: 1090-1100.

Molenda, O. et al. 2012. The alpine cushion plant *Silene acaulis* as foundation species: a bug's-eye view to facilitation and microclimate. - PLoS One 7: e37223.

Nicotra, A. B. et al. 2010. Plant phenotypic plasticity in a changing climate. - Trends Plant Sci. 15: 684-692

Oksanen, J. et al. 2016. vegan: community ecology package. R package. - <https://CRAN.R-project.org/package=vegan>.

Pakeman, R. J. et al. 2006. The extended phenotype of Scots pine *Pinus sylvestris* structures the understorey assemblage. - Ecography 29: 451-457.

Pennings, S. C. et al. 2003. Geographic variation in positive and negative interactions among salt marsh plants. - Ecology 84: 1527-1538.

Pugnaire, F. I. et al. 1996. Facilitation and succession under the canopy of a leguminous shrub, *Retama sphaerocarpa*, in a semi-arid environment in southeast Spain. - Oikos 76: 455-464.

Schindelin, J. et al. 2015. The imageJ ecosystem: an open platform for biomedical image analysis. - Mol. Reprod. Dev. 82: 518-529.

Schöb, C. et al. 2013. Variability in functional traits mediates plant interactions along stress gradients. - J. Ecol. 101: 753-762.

Soliveres, S. and Maestre, F.T. 2014. Plant–plant interactions, environmental gradients and plant diversity: a global synthesis of community-level studies. - Perspect. Plant. Ecol. 16: 154-163.

Thorpe, A. S. et al. 2011. Interactions among plants and evolution. - J. Ecol. 99: 729-740.
Valiente-Banuet, A. et al. 2006. Modern quaternary plant lineages promote diversity through facilitation of ancient tertiary lineages. - Proc. Natl Acad. Sci. USA 103: 16812-16817.

Wang, Q. 1983. Distribution of arid and semiarid areas in China and their major climatic characteristics (in Chinese). - Agric. Res. Arid Areas 1: 11-24.

Whitham, T. G. et al. 2003. Community and ecosystem genetics: a consequence of the extended phenotype. - Ecology 84: 559-573.

Whitham, T. G. et al. 2006. A framework for community and eco-system genetics: from genes to ecosystems. - Nat. Rev. Genet. 7: 510-523.

Wickham, H. 2009. ggplot2: elegant graphics for data analysis. - Springer.

Wu, L. H and Tao, Q. N. 1993. Study on an improved method for the determination of organic carbon in plants (in Chinese). - Chin. J. Soil Sci. 6: 286-287.

Supplementary material (available online as Appendix oik-04512 at <www.oikosjournal.org/appendix/oik-04512>). Appendix 1.