Encroachment drives facilitation at alpine shrublines

1.   Introduction

Widespread vegetation greening via changes in growth, recruitment and distribution range of woody plants in recent decades is considered a fingerprint of global climate warming (Piao et al., 2020; IPCC, 2022). In particular, the greening at high latitudes or elevations is among the world’s most important broad-scale biological responses to anthropogenic climate change (Gao et al., 2022). Heterogeneous vegetation browning and greening-to-browning reversal after 2000 have been also detected in some alpine and arctic regions, resulting from warming-induced soil moisture deficit and other factors (Sigdel et al., 2021; Liu et al., 2023). However, it is unknown how warming and consequent reductions in soil moisture have impacted productivity of alpine woody ecosystems. Such shifts in greening or unexpected browning could result from changes in the interactions between alpine shrubs or reduced recruitment due to drought stress (Liu et al., 2023). Variations of stem density and encroachment rates could alter the intensity or direction of plant-plant interactions and associated feedback processes on species range boundaries (e.g., alpine treelines and shrublines) (Liang et al., 2016; Wang et al., 2016; Pausas and Bond, 2022). However, field evidence remains scarce for alpine shrublines, which represent the uppermost distribution limit of woody plant existence (Wang et al., 2015; Myers-Smith and Hik, 2018).

Evidence is mounting that facilitative interactions play a critical role in driving spatial patterns of woody plants under harsh environmental conditions (Kikvidze et al., 2005; Grau et al., 2012; Wright et al., 2021). Positive interactions (facilitation) due to abiotic amelioration, nutrient enrichment and associational defense are strongly linked to aggregated distribution of plants in cold climate conditions (Choler et al., 2001; Vega-Álvarez et al., 2019). Facilitative interactions also contribute to foster reproduction, survival and growth rate of alpine plants (Callaway et al., 2002; Bruno et al., 2003; Fajardo and McIntire, 2011), thus expanding their realized niche (He and Bertness, 2014). In parallel, owing to the improved microenvironment, facilitators were found to have positive effects on their neighbors' functional traits (e.g. height, canopy cover, specific leaf area) and thus promote ecosystem functioning under stressful conditions (Wang et al., 2021a; Bashirzadeh et al., 2022). The aforementioned research focused mostly on the interspecific facilitative interactions at local and large scales, but the relationships between plant interactions and stem density changes were less studied. Indeed, very few large-scale studies have addressed whether conspecific positive interactions are closely linked to stem density variations and drive the spatial patterns and functioning of alpine woody plant population. Nevertheless, conspecific facilitation within alpine shrublands has been shown to stabilize shrubline dynamics at local scale on the central Tibetan Plateau (Wang et al., 2015).

The world’s highest alpine shrublines are found on the Tibetan Plateau (Lu et al., 2021), where rapid changes of montane vegetation greenness over recent decades provide the ideal context for linking species interactions, vegetation spatial organization and functioning. There, it was found that shrub facilitation was favorable for tree seedling establishment and regulated upward shift rate of alpine treelines (Lyu et al., 2016; Chen et al., 2020; Li et al., 2023). Shrub facilitation can also boost recruitment of accompanying herbs through soil enrichment (Li et al., 2021b). Facilitation from shrubs can increase the above-ground production of neighboring plants and improve ecosystem functioning (Chen et al., 2019). However, none of these studies assessed how conspecific facilitation drove the spatial patterns and functioning of alpine shrubline population at the regional scale, even though shrubs are often considered as the dominant nurse-life form structuring alpine woody plant communities (Filazzola and Lortie, 2014).

Salix oritrepha Schneid. (hereafter Salix) is a widespread alpine deciduous shrub species forming shrublines across the Tibetan Plateau (Fang et al., 2011; Han et al., 2022; Lu et al., 2023). Similar to other arctic or alpine shrublines, Salix shrublines are located below permanent snowlines and above treelines on the Tibetan Plateau (Lu et al., 2022). In contrast to the stability of evergreen Juniperus pingii var. wilsonii shrublines under climate change, however, upward shifts of deciduous Salix shrublines were observed on the Tibetan Plateau (Wang et al., 2021b). Based on a network of Salix shrubline plots, a previous study revealed that Salix shrub encroachment and infilling before 2010 promoted by climate warming was one of major drivers of vegetation greening on the eastern Tibetan Plateau (Wang et al., 2021b). By contrast, due to warming-related moisture stress, Salix shrub recruitment decline occurred in recent decade (Wang et al., 2021b). Nevertheless, the relationships among shrub recruitment dynamics, spatial patterns and conspecific interactions for Salix shrubline populations on the eastern Tibetan Plateau, remain poorly understood.

Positive interactions are generally reflected as clustered patterns of sessile organisms in stressful environments (Callaway et al., 2002; Elliott, 2011; Sigdel et al., 2020); thus, we hypothesized that conspecific facilitation mainly drove the clustered patterns of shrubs at alpine shrublines. Since encroachment and densification could alter plant-plant interactions which further feedback to plant functional traits (Wang et al., 2016), we also hypothesized that shrub density changes mediated conspecific interactions which further affected functional traits at the shrubline. The objectives of this study were to test these two hypotheses by evaluating the relationships between conspecific interactions, shrub spatial patterns and functional traits and changes of shrub density in 26 Salix shrubline plots located on the eastern Tibetan Plateau.

2.   Material and methods

2.1   Study area

We established a network of 26 Salix shrubline plots along wide latitudinal (30°–38° N) gradients in the eastern Tibetan Plateau to characterize shrub spatial patterns and interactions at regional scale (Fig. 1). Cold and continental climate conditions dominate the study region (Yao et al., 2012). January (mean temperature ca. −10 ​℃) and July (ca 10 ​℃) are the coldest and warmest months, respectively. A significant warming trend in summer was observed since 1950, while summer precipitation did not show significant changes (Wang et al., 2021b).

Figure  1.  Spatial locations of the studied shrubline plots (right) on the eastern Tibetan Plateau (upper left insets) (a), view of alpine shrubline in site HSX2 (taken by Y. Wang, July 2018) (b). Green squares represent Salix shrubline plots (a, see plot names in Table S1). The shrub plot (4579.9 m a.s.l.) in site HSX2 included the current shrubline and was located in topographically uniform areas (b). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

2.2   Shrubline plots and field sampling

We surveyed 1–5 shrubline plots in each subregion (n ​= ​10) of the eastern Tibetan Plateau to cover the major distribution range of alpine Salix shrubline (30°–38° N, 96°–103° E) (Fig. 1a, Table S1). Owing to the actual distribution conditions of shrubline, the number of shrubline plots varied within each subregion. The Salix shrubline ecotone in the mountain slope forms a narrow transition zone varying from 40 to 200 ​m which determines the longer side of plot size. A total of 26 rectangular plots across the Salix shrubline ecotone with 30 ​m wide and length changing from 50 to 210 ​m were established between July 2018 and August 2020 in topographically uniform areas (Fig. 1b, Table S1) (see also details in ​Wang et al., 2021b). The elevation of shrublines varied from 3645 to 4895 ​m a.s.l. These shrublines were mainly located on shady and semi-shady slopes (Table S1). These plots were not subjected to strong human disturbances. The longer side of each plot was placed parallel to the altitudinal gradient and included the current shrubline. We used a GPS to measure the lower and upper elevations of each plot. For each plot, the coordinate origin ((x, y) = (0, 0)) was situated in its lower-left corner (Camarero and Gutiérrez, 2004). We recorded the location (accuracy: ±0.1 ​m) of the main stem of each individual shrub within each plot. Three parameters for each shrub were measured: basal stem diameter of each individual, maximum shrub height and crown width (mean value of the two vertical projections along the x and y directions). Dead shrubs were not found within all the plots.

Willow (Salix spp.) shrubs have well-defined annual rings and have been investigated in alpine and arctic dendrochronological studies (Myers-Smith and Hik, 2018; Buchwal et al., 2019; Han et al., 2022). Cross sections at the stem base of all shrubs were collected in each plot. Then, these wood samples of Salix shrubs were processed (air dried, sanded, cross-dated and measured) following dendrochronological procedures (Cook and Kairiukstis, 1990).

The mean series correlations for the cross sections of mature shrubs varied between 0.34 and 0.60, indicating a reliable cross-dating of the shrub ring-width sequences. Finally, the germination age of all the shrub individuals was identified within each plot. The static age structure within each plot represented the balance between the mortality and survival of shrubs; thus, it can be used to reconstruct shrub recruitment dynamics (Wang et al., 2015; Myers-Smith and Hik, 2018; Lu et al., 2022).

2.3   Data analysis

2.3.1   Spatial point pattern analysis

Spatial point-pattern analysis was employed to calculate the spatial patterns of Salix shrubline populations (Wiegand and Moloney, 2014). The inhomogeneous version of the O-ring statistics (O(r)) was used because of the heterogeneous environment across the shrubline ecotone. We used the univariate O₁₁(r) statistic to explore the spatial distribution patterns of two Salix size classes (juveniles: basal stem diameter <10 ​mm, adults: basal stem diameter ≥10 ​mm), while the bivariate O₁₂(r) statistic was used to investigate the spatial associations of these two size classes (adults vs. juveniles). In addition, based on the age structure of each shrubline plot (see also Wang et al., 2021b), mean O₁₁(r) values for shrub individuals after 2010 (age ≤10 ​yr) and before 2010 (age ≤10 ​yr) was calculated to explore the spatial distribution patterns of juveniles over these two time periods. In the univariate and bivariate cases, the observed values were compared to the lower or upper limit of the 95% confidence envelopes based on 999 Monte Carlo simulations (Wiegand and Moloney, 2004). The observed values of O₁₁(r) situated below or above the lower and upper limit of the envelopes indicate significant regularity or clustering at the studied spatial distance, respectively, whereas the observed values of O₁₂(r) situated below or above the lower and upper limit of the envelopes signify significant repulsion or attraction, respectively. The spatial point pattern analyses were implemented in the Programita software (Wiegand and Moloney, 2014). Additionally, one-way ANOVA tests were used to test if significant differences existed between the mean intensity of spatial patterns of juveniles and adults within the study plots.

2.3.2   Relating spatial patterns to abiotic and biotic variables

Following previous studies, the mean intensity of univariate and bivariate spatial point patterns can be related to abiotic and biotic variables (López et al., 2016). For instance, the intensity of woody plant interactions was linked with biotic and climatic factors (Muhamed et al., 2013). In this study, mean O₁₁(r) values for the two considered size classes at scales of 1–10 ​m and mean O₁₂(r) values between them in each plots were used to characterize the mean intensity of univariate spatial patterns and conspecific interactions. Given that Salix shrubline dynamics are driven by summer warming but mediated by summer soil moisture level on the eastern Tibetan Plateau (Wang et al., 2021b), the mean intensity of univariate spatial patterns and conspecific interactions were related to these temperature and moisture variables (spring and summer warming rate, spring and summer moisture condition, spring, summer and annual precipitation, spring, summer and annual VPD). Seasonal warming rates and precipitation in the 26 study plots were derived from meteorological data and the CRU 0.5°-gridded database, whilst soil-moisture content at depth of 10 ​cm was extracted from the NASA GLDAS-2 (see more details in ​Wang et al., 2021b). Atmospheric vapor pressure deficit (VPD, an indicator of evaporative demand) since 1954 ​at near plot was calculated from CRU TS4.03 gridded data (0.5° ​× ​0.5° grids). Several site factors (aspect of each shrubline ecotone, northness (sin(aspect) ​× ​cos(slope)) and rainfall continentality) influencing local moisture regimes were also related to the mean intensity of univariate spatial patterns and conspecific interactions. Note that rainfall continentality, which depends on site elevation and annual precipitation, usually increased (decreased) with increasing (decreasing) distance from the Pacific Ocean leading to rainshadowed continental (oceanic) climates and it is an important driver of alpine plant communities in the Tibetan Plateau and other ranges (Li et al., 2021a; Michalet et al., 2021). The linear regression equation for shrub recruitment series in each plot was calculated and its slope was used as a recruitment rate (Table S5, Fig. S1). Considering that spatial patterns of plants are partly density dependent (Cameron et al., 2019), we calculated the Pearson correlations between rate of Salix shrub recruitment in recent decades and mean intensity of O₁₁(r) and O₁₂(r) values for the two Salix size classes. Pearson correlations between summer warming rate and mean shrub recruitment rate over the past decades were calculated to explore their relationships. Additionally, the mean intensity of O₁₁(r) values for shrub juveniles before and after 2010 (age ≤10 ​yr) was related to shrub recruitment before and after 2010, respectively.

Based on the correlations between spatial patterns or associations of shrub population and different environmental variables, the most important temperature or moisture variables were used in the principal component analysis (PCA). The scores of PCA with high explained variance were used as indicators of climatic stress and were related to vegetation patterns or interactions.

2.3.3   Bivariate point pattern analysis with one quantitative mark

Bivariate point pattern analysis considering a quantitative mark (e.g., a plant trait) is frequently used to examine the influence of the pattern of large individuals on the pattern of small individuals (Wiegand and Moloney, 2014). For instance, the normalized bivariate mark correlation functions K_m1m2(r) with a quantitative mark can be applied to investigate the spatial correlation of marks associated with pairs of points at the given spatial scale. Values of K_m1m2(r) ​< ​1 or K_m1m2(r) ​> ​1 represent mutual inhibition or attraction for the quantitative marks of individuals, respectively (Wiegand and Moloney, 2014). We calculated the bivariate patterns (adults vs. juvenile) of the following variables or quantitative marks: shrub height, basal stem diameter, and crown width. Then, we also calculated a goodness-of-fit (GoF) test to assess the relative intensity of each analyzed pattern (Wiegand and Moloney, 2014). Lastly, we calculated the percentage value using the number of plots with attraction divided by the total number of plots (n ​= ​26).

2.3.4   Linking conspecific interactions, patterns, functioning and shrub recruitment

Structural equation models (SEMs) are used to explore complex multivariate linkages among interrelated variables through path analysis (Lefcheck, 2016). We used SEMs to explore the relationships among population interactions, patterns, traits and shrub recruitment rate before 2010 and after 2010. Specifically, climatic (summer warming rate, summer VPD changes) and non-climatic variables (mean intensity of conspecific facilitation and clustering of juveniles, shrub recruitment rate before and after 2010, mean values for the mark correlation function with one shrub trait (height, basal diameter, crown width) were used in establishing the SEMs. The goodness of fit of the SEMs was assessed by a Chi-square test (p) and Akaike Information Criterion (AIC). The model with p ​> ​0.05 and lower AIC value was considered as the best-fitted model.

3.   Results

3.1   Spatial patterns of shrub size classes

Juveniles and adults usually formed dense patches at the shrubline (Fig. 2). Shrub recruitment showed an increasing trend over the past decades in most (88%) plots (Fig. S1a–k, Table S5). However, an increasing recruitment trend before 2010 but a declining trend after 2010 was also detected in the studied plots (Fig. S1a–k, Table S5).

Figure  2.  Point patterns of all shrub individuals mapped in 26 shrubline plots situated on the eastern Tibetan Plateau. Shrubs are grouped in two size classes: blue circles denote juveniles (basal stem diameter <10 mm), and red rhombuses represent adults (basal stem diameter ≥10 mm). The x axis was set along the contour, and increasing altitude corresponds to higher values in the y axis. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

Juveniles in 84.6% of the studied plots were spatially aggregated at scales of 1–10 ​m, but were randomly or regularly distributed up to 11–15 ​m (Table S2). Adults were spatially aggregated at scales of 1–10 ​m in 65% of the plots, but were randomly or regularly distributed from 11 to 15 ​m in 88% of the plots (Table S2). Adults and juveniles showed significant positive associations at scales of 1–10 ​m in 73% of the plots, but they were not spatially related from 11 to 15 ​m in 81% of the plots (Table S3). The intensity of spatial patterns of juveniles or adults and their spatial associations varied among the study plots (Fig. 3). According to one-way ANOVA analysis, the mean intensity of spatial patterns of juveniles were significantly higher than that of adults (F ​= ​6.34, p ​= ​0.015).

Figure  3.  Bivariate spatial patterns based on the O₁₂(r) statistic and calculated for the 26 shrubline plots. We considered two size classes of shrub individuals (adults, juveniles). Solid lines with symbols indicate the O₁₂(r) statistics, and the dashed lines represent the upper and lower 95% confidence envelopes. Values of the O₁₂(r) statistic situated below or above the and lower or upper envelopes represent significant repulsion or attraction, respectively (see also Table S3).

DownLoad: Full-Size Img PowerPoint

3.2   Relationships between spatial patterns and abiotic and biotic variables

The mean intensity for the spatial patterns of juveniles and spatial associations between juveniles and adults differed among the study plots (Figs. 2 and 3). Then the mean intensities of the patterns were negatively related to summer soil moisture content (Table S6). Non-linear functions with two inflection points described the relationships between mean intensity of the spatial patterns of juveniles and spatial associations between juveniles and adults and summer soil moisture content (Fig. 4a). The mean intensity for the spatial patterns of juveniles and the spatial associations between juveniles and adults showed positive and significant relationships with shrub recruitment rate in recent decades (Fig. 4b). Summer warming rate was positively associated with mean shrub recruitment rate in recent decades (r = 0.52, p < 0.01). However, the mean intensity for the spatial patterns of adults was not significantly related to summer warming rate or soil moisture content (Fig. S2, Table S6). Spring and summer VPD positively affected the mean intensity for the spatial patterns of juveniles and the spatial associations between juveniles and adults (spring VDP: p ​< ​0.05; Summer VPD: p ​> ​0.05) (Table S6). A significant impact of seasonal precipitation on vegetation patterns or interactions was not detected. Three site factors (aspect, northness, rainfall continentality) did not have significant effect on the spatial patterns of juveniles and the spatial associations between juveniles and adults (Table S6). In addition, the mean intensity for the spatial patterns of juveniles before and after 2010 was positively and negatively associated with shrub recruitment before 2010 and after 2010, respectively (Fig. S3).

Figure  4.  Significant relationships between mean intensity for spatial patterns of juveniles and spatial associations between juveniles and adults and summer soil moisture content (a) and shrub recruitment rate (b) in recent decades. The fitted lines with different colours are shown for juveniles (JV, blue line), and juveniles-adults (JV-AD, black line). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

Results of PCA including three important climatic factors (mean summer soil moisture content, changes of mean spring VPD, summer warming rate) showed that the first three components explained 77.48% (PCA1), 17.16% (PCA2), 5.36% (PCA3) of the variance in the climatic matrix, respectively (Table S7). The loading of PCA1 was mainly contributed by mean summer soil moisture. PCA1 score as the indicator of summer soil moisture was negatively associated with mean intensity for spatial patterns of juveniles and spatial associations between juveniles and adults (p ​< ​0.05 in both cases), but it did not have significant effects on the mean intensity for spatial patterns of adults (Table S7).

3.3   Bivariate mark correlation functions of shrub size traits

Juvenile shrubs were higher than expected at scales of 1–15 ​m from adults in the 77% of the studied plots excluding 6 cases (QL3, DL2-3, CDUO1, LWQ2, BD1) (Fig. 5, Table S4). Juvenile shrub individuals showed thicker basal stems than expected when located 1–15 ​m from adults in the 69% of the studied plots excluding 8 cases (WL1-3, DL2-3, BZ1-2, LWQ2) (Fig. 5, Table S4). Juveniles were larger in crown width than expected when located 1–15 ​m from adults in the 46% of the plots excluding 14 cases (QL3, WL2-3, HZ1-2, DL1-3, CDUO1-2, BZ2, LWQ1-2 and BD1) (Fig. 5, Table S4).

Figure  5.  Mark correlation functions (K_m1m2(r)) based on the spatial patterns of shrub height within 26 shrubline plots. Solid lines with symbols indicate the K_m1m2(r) statistics, and the dashed lines represent the upper and lower 95% limit of the simulated confidence envelopes. Values of the K_m1m2(r) statistic located below or above the and lower or upper envelopes represent significant inhibition or attraction, respectively. Values of K_m1m2(r) approaching 1 represent that the mark-pair products at a given distance only depend on the spatial location of Salix shrub individuals. The goodness-of-fit (GoF) tests and their significance levels (p) are shown in Table S4.

DownLoad: Full-Size Img PowerPoint

3.4   Linking spatial patterns and shrub recruitment

We used SEMs to reveal the linkages among spatial patterns and shrub recruitment dynamics. The first SEM (p ​= ​0.076, AIC ​= ​50.93) indicated that conspecific facilitation positively affected clustering of juveniles/adults and shrub size-related traits (shrub height, basal stem diameter, and crown width) (Fig. 6a). Note that conspecific facilitation and clustering of juveniles mutually promoted each other (Fig. 6a). The second SEM (p ​= ​0.065, AIC ​= ​29.21) indicated that shrub recruitment before 2010 promoted by summer warming but inhibited by summer VPD had a positive effect on clustering of juveniles before 2010 (Fig. 6b). The third SEM (p ​= ​0.843, AIC ​= ​22.83) indicated that shrub recruitment after 2010 limited by summer moisture level had a negative effect on clustering of juveniles after 2010 (Fig. 6c).

Figure  6.  Relationships among population interactions, patterns and traits, shrub recruitment and climate change. Structural equation models show the relationships among conspecific facilitation, clustering of juveniles/adults, and shrub traits (a), and the relationships among climate change, shrub recruitment before 2010 (b) and after 2010 (c) and conspecific clustering patterns. Arrow thickness denotes the path coefficient value.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

The stress gradient hypothesis states that facilitation (positive interactions) will increase as increasing environmental stress does (Bertness and Callaway, 1994; He et al., 2013). However, the refined stress gradient hypothesis suggested that maximum of facilitation usually occurred at intermediate stress levels (Michalet et al., 2006, 2014; Holmgren and Scheffer, 2010). Alpine shrublines are subject to very harsh conditions (e.g., low temperature, short growing season, seasonal drought, thaw-freeze cycles) (Körner, 2012a); thus, it is reasonable to consider that facilitation may play a key role in structuring alpine shrubline populations. In this study, the mean intensity of clustering of juveniles and inferred conspecific facilitation between juveniles and adults were negatively associated with summer soil moisture content but positively associated with spring and summer VPD, illustrating that facilitation intensity increased with increasing moisture stress. However, non-linear curves with two inflections well fitted the changes of mean intensity for clusters of juveniles and conspecific interactions along the summer soil moisture gradient, thus partly supporting the refined stress gradient hypothesis (Michalet et al., 2006, 2014; Holmgren and Scheffer, 2010). Juveniles may benefit from intraspecific clustering and facilitation through the amelioration of soil moisture levels. It was reported that soil moisture increased more under or near nurse shrubs than in open areas in semi-arid and in alpine habitats due to reduced thermal amplitude and decreased evaporation rate, thus facilitating seed germination and growth of juveniles (Armas and Pugnaire, 2005; Mihoč et al., 2016). Structural modifications in the soil induced by neighbors also lead to a greater water and nutrient availability under or near them than in bare areas in Mediterranean sites (Verdú and García-Fayos, 1996). However, we did not find significant relationships between mean intensity for clustering of adults and summer soil moisture, likely because adults could have more access to deep soil moisture reserves and better tolerate water shortage (Wright et al., 2014). In addition, the mean intensity for spatial patterns of juveniles and adults and conspecific facilitation were not related to summer warming rate, indicating that climate warming did not have direct impacts on spatial patterns and related interactions. Given the significant linkages between summer warming and shrub recruitment rates found in this study and before (Wang et al., 2021b), we cannot exclude the possibility that summer warming rate may indirectly drive vegetation spatial patterns or associations through altering shrub recruitment rates (Büntgen et al., 2015; Lu et al., 2022).

Plant clustering and facilitation are widely reported in harsh environments such as alpine or arctic regions (Choler et al., 2001; He et al., 2013; Verdú et al., 2021). In particular, several studies emphasized the key role of facilitation in broadening the niche of plants at their distribution limits (He and Bertness, 2014). In this study, clustered distribution and conspecific facilitation at scales of 1–10 ​m dominated the studied Salix shrubline populations. It has been suggested that seed dispersal distance, habitat heterogeneity and facilitation could lead to the formation of clustered patterns (Lara-Romero et al., 2016). Field data showed that Salix shrub seed could be dispersed upslope by as much as 65 ​m in elevation from parent shrubs based on the colonization of few Salix seedlings located above the current shrubline (Wang et al., 2021b). Therefore, seed dispersal distance might not be the primary factor for clustering pattern formation at the shrublines. Both inhomogeneous and homogeneous versions of the point pattern model used in this study showed the similar distribution patterns (Law et al., 2009; Wiegand and Moloney, 2014), therefore, habitat heterogeneity may not be the main driver of clustered distributions at the shrubline. Many studies have confirmed that clustering as a form of spatial self-organization in severe environment usually involves positive interactions which contribute to ameliorate abiotic stress and enhance individual performance (Trautz et al., 2017). Indeed, the fitted models suggested that conspecific facilitation positively affected clustering of juveniles and adults. In such cases, clustering patterns of alpine Salix shrubs could be largely attributed to conspecific facilitation, thus supporting our first hypothesis. Nevertheless, the random distribution of the two compared size classes and the lack of spatial relationships between them at scales of 11–15 ​m may indicate the combined effects of seed dispersal and low post-dispersal survival (Nathan et al., 2002; Rietkerk and Van de Koppel, 2008; Shen et al., 2018). Furthermore, clustering intensity significantly differed between juveniles and adults, suggesting that vegetation patterns would vary with life stage (young vs. mature shrubs) which involves different stress-tolerant ability and environmental sensitivity during shrub ontogeny (He et al., 2013; Wright et al., 2014). Our results agree with studies of alpine treelines suggesting that tree clustering or conspecific facilitation usually dominates tree population and tree spatial patterns are size or age-dependent (Elliott, 2011; Harper et al., 2011; Sigdel et al., 2020).

Positive plant-plant interactions might be modulated by changes in plant recruitment through time (Wieczorek et al., 2017). Many studies have documented that encroachment and densification of woody plants at high altitudes and high latitudes are occurring in the past decades (Camarero and Gutiérrez, 2004; Liang et al., 2016; Devi et al., 2020; Shi et al., 2022), which is one of major drivers of vegetation greening at such marginal sites. Considering that plant-plant interactions are highly density-dependent (Cameron et al., 2019), variations of vegetation greenness and productivity due to plant recruitment shifts could alter the intensity or direction of plant-plant interactions. It was found that conspecific competition between mature and juvenile cohorts intensified with increasing tree recruitment at treelines since 1950s on the southeastern Tibetan Plateau (Wang et al., 2016). In contrast, conspecific facilitation increased with climate warming at treelines across Canada (Harper et al., 2011). Although shrub range shifts and changes of shrub recruitment have been widely documented in alpine and arctic regions across the Earth (Myers-Smith and Hik, 2018; Lu et al., 2022), there are very few studies linking shrub recruitment with plant-plant interactions. In this study, the fitted models indicated that summer warming positively related to conspecific clustering before 2010 through increased shrub recruitment. Conversely, summer warming and related evaporative demand after 2010 negatively drove conspecific clustering through decreased shrub recruitment. Given that conspecific facilitation was closely coupled to conspecific clustering, the fitted models revealed that shrub recruitment dynamics under climate change largely modulated population interactions and patterns, thus partly confirming our second hypothesis. Conceivably, facilitation changes induced by shrub recruitment shifts may amplify the effects of low-frequency climatic variation, potentially resulting in slow responses of shrublines to climate (Butterfield et al., 2010). Our results are consistent with other studies highlighting the benefits of living in high-density aggregated patterns in some harsh environments including arctic and alpine communities (Choler et al., 2001; Bruno et al., 2003; He and Bertness, 2014; Silliman et al., 2015).

Positive plant-plant interactions could have considerable effects on species functional traits (Verdú et al., 2021). Field evidence showed that facilitation from nurse shrubs substantially improved leaf and stem traits of forb species at or above the treeline (García-Cervigón et al., 2015; Wang et al., 2021a). In parallel, transplant experiments demonstrated that intraspecific facilitation markedly increased survival or height growth of tree seedlings in comparison with open areas at and beyond the treeline (Fajardo and McIntire, 2011; Ettinger and HilleRisLambers, 2017). Moreover, nurse plants in alpine habitats enhanced survival and reproduction of protected plants (Callaway et al., 2002). In this study, juvenile shrubs growing around adults tend to be higher and thicker in the majority of the studied plots, while juveniles located around adults tend to have larger crowns in nearly half of the studied plots, indicating that conspecific facilitation would increase the range of functional traits related to plant size. Generally, plant facilitation may create favorable microenvironments, thus leading to suitable niche spaces for co-existing individuals and act as selective force for trait evolution (Verdú et al., 2021; Pausas and Bond, 2022). Equally importantly, size-related traits such as height and stem diameter are important indicators for aboveground carbon storage, stand structure and resource-use strategy (Enquist et al., 2009; Bjorkman et al., 2018; Conti et al., 2019), whereas plant crown width is closely related to carbon assimilation, respiration load, nutrient cycling, surface albedo and heat budget (Poorter et al., 2006; Fichtner et al., 2017; Li et al., 2017). Therefore, the increased values in size-related traits due to conspecific facilitation may have substantial impacts on the functioning of alpine shrublines. As mentioned above, shrub recruitment changes mediated conspecific facilitation, which was positively related to shrub size-related traits, thus supporting our second hypothesis.

Warming-recruitment relationships shifted from positive before 2010 to negative coupling after 2010, while recruitment-clustering relationships changed from positive before 2010 to negative coupling after 2010, suggesting that shrub recruitment shifts under rapid climate change might modify the direction of shrubline population responses. Recruitment-clustering correlations were higher after 2010 than before 2010, implying that the degree of shrubline population responses could increase with continuous climate warming. In particular, a weakening of facilitation due to climatic warming and drying may enhance the impacts of regional climate trends on alpine shrub ecosystems which are usually shaped by microclimatic conditions (Körner, 2012b). Given the key role of facilitation in maintaining vegetation structure and functioning in harsh environments (Cavieres et al., 2016; Wright et al., 2021), the negative coupling of climate warming, recruitment dynamics and vegetation pattern could alter the future development of alpine shrub ecosystems.

5.   Conclusions

Using a network of 26 Salix shrubline plots located along a broad geographical gradient spanning ca. 900 ​km in the Tibetan Plateau, we revealed the linkages among shrub recruitment, shrub spatial patterns, conspecific interactions and vegetation functioning under climate change. Clustering patterns and conspecific facilitative interactions dominate Salix shrublines across the Tibetan Plateau, thus providing further support to the stress gradient hypothesis. Previous studies found that plant functional traits acted as determinants of population stability (Májeková et al., 2014). This study elucidated that conspecific facilitation tended to enhance size-related functional traits such as height, stem diameter, and crown width, thereby driving the stability of shrubline population. It has also been reported that conspecific facilitation contributes to stabilize shrubline populations on the central Tibetan Plateau (Wang et al., 2015). Therefore, plant facilitation acting as an insurance policy plays an important role in sustaining the functioning and stability of alpine ecosystem. Further, shrub recruitment dynamics under environmental change (climate warming, increase evapotranspiration rates) would alter population spatial patterns and interactions. It should be noted that the shrub recruitment decline detected in recent decades can lead to negative cascade effects on population patterns and functioning with implications for the stability of alpine ecosystems. Further work is needed to clarify the mechanisms for the presented patterns (e.g., controlled experiments on facilitation, impact of facilitation on leaf and wood traits, microenvironment monitoring).

Data availability

Shrubline datasets analyzed in the current study are deposited on the website of the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/en/disallow/5bd43be4-a31e-46be-be5c-67b09e6f98ae/)

CRediT authorship contribution statement

Yafeng Wang: Writing – original draft, Investigation, Funding acquisition, Formal analysis, Conceptualization. Eryuan Liang: Writing – review & editing, Formal analysis, Conceptualization. J. Julio Camarero: Writing – review & editing, Writing – original draft, Methodology, Formal analysis, Data curation.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.fecs.2024.100168.