# Changes in the Strength of Associations Between Tree Seedlings and Understory Shrubs Along a Regional Drought Gradient in Lebanese Coniferous Forests

Andrea Maamary, Florian Delerue, Richard Michalet

► To cite this version:

Andrea Maamary, Florian Delerue, Richard Michalet. Changes in the Strength of Associations Between Tree Seedlings and Understory Shrubs Along a Regional Drought Gradient in Lebanese Coniferous Forests. Journal of Vegetation Science, 2025, 36 (4), 10.1111/jvs.70054. hal-05479066

HAL Id: hal-05479066

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

Submitted on 20 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

OPEN ACCESS

# Changes in the Strength of Associations Between Tree Seedlings and Understory Shrubs Along a Regional Drought Gradient in Lebanese Coniferous Forests

Andrea Maamary ID |

Florian Delerue ID |

Richard Michalet ID |

Université Bordeaux, CNRS, Bordeaux INP, EPOC, UMR 5805, Pessac, France

**Correspondence:** Richard Michalet (richard.michalet@u-bordeaux.fr)

**Received:** 2 January 2025 |

**Revised:** 23 May 2025 |

**Accepted:** 14 July 2025

**Funding:** This work was supported by the UMR CNRS 5805 EPOC, France.

**Keywords:** competition |

facilitation |

Juniperus oxycedrus |

Mediterranean forests |

regeneration |

VPD

## ABSTRACT

**Aims:** Do the effects of an understory shrub on woody seedlings explain the distributions of relict conifer trees at a regional scale?

**Location:** Three natural *Cedrus libani* reserves from Lebanon.

**Methods:** Using an observational approach, we quantified associations of recruits of woody species at the center and periphery of the spiny shrub *Juniperus oxycedrus* in forests and gaps along a north-south climatic gradient of both winter rainfall continent- tality and summer drought. We quantified the associations of woody seedlings with shrubs using the Relative Interactions Index (RII) calculated for recruits of all species together (community scale) and for four groups of functionally different species. We then analyzed the variation of RIIs along different treatments. Vapor Pressure Deficit (VPD) was quantified in each treatment combination, and its relationship with species-group RIIs and abundance was evaluated.

**Results:** At the community-scale, shrubs showed likely higher negative effects on woody seedlings at their center than at their periphery. In contrast, there were contrasting responses at the species-group level, with overall negative associations for *Cedrus*, positive associations for a group including *Abies cilicica*, and intermediate responses for two other groups. Positive associations strongly varied along the drought gradient and with the canopy treatment, while negative associations only slightly increased from forests to gaps. Positive associations varied together with VPD, peaking at an intermediate position along the drought gradient coinciding with the southernmost limit of *Abies* and declining at the dry extreme of the gradient where *Cedrus* is still present.

**Main Conclusions:** VPD appeared as a crucial driver of tree species occurrence and likely facilitation. Although manipulative experiments are needed to confirm the real occurrence of facilitative and competitive effects of the shrub, our study highlights the potential role of facilitation in explaining tree species distribution at regional scales.

## 1 | Introduction

Facilitation and competition are important drivers of commu- nity composition and diversity (Chesson 2000; Lortie et al. 2004;

Michalet, Maalouf, et al. 2015; Soliveres et al. 2015). Plant-plant interactions vary along environmental gradients, and pioneer models have proposed that competition should be dominant when stress and disturbance are low (Grime 1973) and shift

**Co-ordinating Editor:** Stephen Roxburgh

This is an open access article under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2025 The Author(s). *Journal of Vegetation Science* published by John Wiley & Sons Ltd on behalf of International Association for Vegetation Science.
to facilitation with increasing environmental severity (i.e., the Stress-Gradient-Hypothesis, Bertness and Callaway 1994; Bertness et al. 2024). However, facilitation has been shown to decrease in extreme conditions of stress or disturbance (e.g., Kitzberger et al. 2000; Tielbörger and Kadmon 2000). Therefore, facilitation may peak at mid position along environmental gradients where diversity is the highest (Michalet et al. 2006; Holmgren and Scheffer 2010). In addition, Michalet et al. (2014) emphasized the importance of distinguishing two scenarios and types of gradients: (i) the decline in facilitation, which is the occurrence of null interactions in extreme conditions of non-resource stress or disturbance (Chaieb et al. 2021) (ii) and the switch back to competition, which is the return to competition in extreme conditions of resource stress (e.g., nutrient, water) (Tielbörger and Kadmon 2000; Maestre and Cortina 2004). In the former scenario, facilitation decreases to null interactions because nurse species are too weak to ameliorate environmental conditions for the beneficiary species (e.g., Chaieb et al. 2021) or because all remaining species are strongly stress-tolerant and do not require positive effects provided by other species (Liancourt et al. 2017); whereas in the latter, competition increases because of decreasing resource availability, consistent with the MacArthur and Wilson (1967) r/K strategy model (e.g., Davis et al. 1998).

Plant communities include species of contrasting functional strategies that may have contrasting responses to neighbors (Chapin III and Shaver 1985; Michalet, Chen, et al. 2015). Although interactions may be null or very weak when quantified at the community level, strong positive and negative interactions may exist at the species or functional group level within a single community (Bai et al. 2021). Liancourt et al. (2005) have shown that species tolerating drought but not shade have generally negative responses to neighbors because, for these species, the benefit of ameliorating physical stress is low while the cost of decreasing light is high. For opposite reasons, species tolerating shade but not drought have generally positive responses to neighbors (Liancourt et al. 2005). Therefore, to better understand changes in plant-plant interactions along spatial and temporal climate stress gradients, it is crucial to assess the responses of different functional groups of species while identifying the direct mechanism of interactions for each species. The quantification of spatial associations is a crucial preliminary method to assess the likely interactions between understory and nurse species since they include the response of all community members, in opposition to experimental removal studies that only include a few target species (Michalet, Chen, et al. 2015). However, spatial associational measurements can be confounded by habitat-sharing effects if not conducted with careful verifications of environmental homogeneity within the communities between vegetated and unvegetated patches (Michalet 2006; Steinbauer et al. 2016; Michalet et al. 2024). Therefore, the combination of both methods is certainly the most appropriate way of assessing plant interactions in plant communities.

The assessment of variation in plant-plant interactions along drought gradients needs to distinguish two different drought stresses: the soil and atmospheric stresses (Saccone et al. 2009; Aguirre et al. 2021). Focusing on the atmospheric stresses, high Vapor Pressure Deficit (VPD) can lead to water (and temperature) stress in plants, affecting their growth, photosynthesis,

and overall health (Novick et al. 2016). Nurse plants, through decreasing temperature with shade and increasing air humidity with transpiration, can decrease VPD beneath their canopy, thereby mitigating atmospheric stress for understory species (De Frenne et al. 2019). However, atmospheric drought stress has been overlooked in plant-plant interaction studies because many researchers have instead focused on soil drought (but see Muhamed et al. 2013; Guignabert et al. 2020). Aguirre et al. (2021) investigated the role of microclimate in biodiversity experiments. Simulating natural atmospheric drought scenarios, they found that, during periods of high VPD, facilitative interactions become more prominent in plant communities with greater diversity. In contrast, soil water availability stress is rather involved in competitive interactions, since species may compete more when resource (water, nutrients and light) availability decreases (MacArthur and Wilson 1967). Because VPD is a non-resource ecological factor, species are less likely to compete in conditions of high VPD (Michalet et al. 2014).

During the severe 2003 European heat wave, Saccone et al. (2009) compared the survival of different tree seedlings transplanted into the understory and adjacent gaps of forests from dry ridges and mesic valleys in the Alps. They found that the seedlings of *Abies alba* were facilitated because of increasing shade (and decreasing VPD) below the canopy of adult trees, while those of *Picea abies* and *Acer pseudoplatanus* suffered from competition in the same canopy conditions. Furthermore, they found that facilitation for *Abies* was the highest in the mesic valley forests and collapsed in the drier ridge forests. These results were highly consistent with the ecological requirements of these three species. *Abies alba* can grow in Mediterranean climates with low soil-water availability during the growing season, but not in very continental climates with high VPD (except in shady understories). *Picea abies* and *A. pseudoplatanus* are not found in Mediterranean climates with high soil drought stress during the growing season, but occur in temperate continental climates with high light and high VPD (Michalet et al. 2003).

Most biogeographical studies conducted along drought gradients have not separated the effects of low soil-water availability during the growing season from those of high atmospheric stress related to winter rainfall continentality (Bagnoulous and Gaussen 1953; Stephenson 1990). However, Michalet et al. (2021) analyzed tree-species distribution in western USA in relation to both rainfall continentality (quantified by the amount of winter rainfall, Gams 1932) and summer-water balance. They found that rainfall continentality was the main climatic variable explaining the distribution of tree species in this large geographical area encompassing contrasting climates. In addition, consistent with Saccone et al. (2009) they found that *Abies* species were mostly found in oceanic temperate and Mediterranean climates (characterized by low VPD due to important air humidity) with potentially important soil drought during summer, while *Picea* species were found in continental climates (with potentially high VPD) with higher soil-water availability in summer. Thus, although low rainfall during the growth season has often been proposed as the only important drought factor in bioclimatology and ecophysiology (Bagnoulous and Gaussen 1953; Stephenson 1990), it may be less important for species distribution than rainfall continentality, in particular in mountain areas with strong rainshadow effects.
The apparent consistency existing between tree species distri-
bution at the local and regional scales along climatic gradients
(Michalet et al. 2021) and tree seedling responses to the effects of
neighbors at the neighboring scale (Saccone et al. 2009) suggests
that the direct ecological mechanisms involved in these two
processes are the same. However, to our knowledge, there have
been no studies that have assessed variations in tree-seedling re-
sponses to the effects of canopy trees or understory shrubs along
regional climate gradients. Such studies are worth conducting if
we want to quantify the contribution of plant-plant interactions
to current tree species distributions and how it may influence
species responses to climate change.

The main goal of our study is to quantify the magnitude of changes
in spatial patterns of associations (positive or negative) between
understory shrubs and woody seedlings of several species of con-
trasting geographical distributions along a regional-climatic gra-
dient. We chose the Lebanese Mediterranean coniferous forests
as a model system. In Lebanon, there are a limited number of
relict *Cedrus libani* forests that are distributed along a north-to-
south gradient of drought due to both increasing winter rainfall
continentality and summer drought. In a previous study, Douaihy
et al. (2022) quantified the strength of associations between the
shrub *Juniperus oxycedrus* and woody seedlings below the for-
est canopy and in gaps at the northernmost and wettest Cedar
Reserve of Ehden. They found contrasting responses to the effects
of the shrub depending on the woody species, with (i) important
positive associations for a group of shade-demanding tree species
including *Abies cilicica*, a rare species only located in this wet
northernmost reserve, (ii) weak negative associations for a group
of shade-intolerant species including *C. libani*, found further south
in drier reserves. In the present study, we quantified, with the
same method, associations between the Juniper and tree seedlings
at four sites along the entire drought gradient in order to assess
the magnitude of changes in positive and negative associations
along this complex gradient, with the ultimate goal of quantifying
the relative contribution of facilitation and competition to woody
species distributions. We make the following predictions: (i) con-
sistent with Douaihy et al. (2022), the seedlings of the tree spe-
cies will have contrasting responses to the effects of the Juniper,
(ii) consistent with Saccone et al. (2009) and Michalet et al. (2021)
there should be stronger variation in positive than negative associ-
ations along our treatments, (iii) variation in positive associations
should be highly correlated with variation in VPD, (iv) variation in
positive associations along the drought gradient is likely to explain
the only occurrence of *Abies* in northern Lebanon.

2 | Material and Methods

## 2.1 | Climatic Gradient, Study Sites and Species

Due to the rarity of relict Cedar forests in the country (only eight
Cedar forests), among which four of them are very small (less
than 10 ha and one is highly grazed), we selected in the three
largest Natural Reserves of Lebanon four sites along a drought
gradient driven by several complex factors: latitude, altitude,
climate, topography, and exposure to sunlight (Figures 1-3).
First, the four sites were positioned at three different latitudes.
The two wettest sites were located at the northernmost latitude
(33.2) in the Horsh Ehden Nature Reserve (HENR, hereafter

Ehden), the third site at a slighter southern latitude (33.1) in the
Tannourine Reserve, and the driest site at the southernmost lat-
itude (32.4) in the Shouf Reserve (Figure 1). The three reserves
were protected from grazing and various human activities and
located on a calcareous bedrock.

The three reserves are characterized by a Mediterranean oce-
anic climate with high winter rainfall and long summer drought
(Figure 1). However, Ehden and Tannourine are character-
ized by a slightly shorter summer drought period (4 months,
Figure 1b,c), than Shouf (5 months, Figure 1d), due to the lower
latitudinal position of the latter. Additionally, although Ehden
and Tannourine share very similar macroclimatic conditions due
to their close latitudes (see the position of the sites along the first
axis of the PCA conducted on climate conditions, Figure 3a,b),
the more convex topography and more exposed dominant west
and south exposures of Tannourine (see Figure 2a,b, vs. c and
d, respectively) makes it drier than Ehden. Finally, the two wet-
test sites from Ehden were selected at two different elevations,
1700–2000 m for Ehden High and 1500–1700 for Ehden Low,
with, thus, slightly lower VPD in gaps (Figure 3c) for the for-
mer than the latter and also shorter summer drought season. To
summarize, the four sites were located at four different positions
along a drought gradient that was due to both increasing sum-
mer drought and rainfall continentality.

The four forests were dominated by *C. libani* A. Rich., mixed
with various other woody species, such as junipers (*Juniperus
excelsa* M. Bieb, *J. oxycedrus* L.), oaks (*Quercus calliprinos*
Webb, *Quercus infectoria* Oliv. and other deciduous oaks), wild apple
(*Malus trilobata* Lab. Schneider), maple species (*Acer hermo-
neum* Bornm & Schwer, *Acer tauricolum* Boiss. & Bal.), and sor-
bus species (*Sorbus torminalis*, *Sorbus flabellifolia*). Additionally,
the Cilician fir (*Abies cilicica* Antoine et Kotschy, Carrière
1855) was present (and abundant) only at the two wettest sites,
Ehden High and Low. The Horsh Ehden Nature Reserve is the
southernmost latitudinal limit of the distribution range of the
species, which is more widespread north of Lebanon, in Syria
and Turkey. Additionally, the driest site (Shouf) is also the
southernmost limit of the geographical distribution of *C. libani*
(Figure 2). *Juniperus oxycedrus* was used as a nurse species in
this study because it was present at the four sites, both in forests
and gaps. This is a drought-resistant and light-demanding shrub
species, native to the Mediterranean region. This spiny shrub is
very variable in shape, with either an erect (small tree of maxi-
mum eight meters high) or a prostrate (spreading shrub shape)
phenotype.

## 2.2 | Sampling Design at Each Site

We used an observational approach to quantify associations of nat-
ural recruits of woody species with the Juniper in the four sites
along the regional drought gradient. In early May 2023, we ran-
domly selected at each site five Juniper individuals for each combi-
nation of a Canopy (gap vs. forest canopy) and a Phenotype (erect
vs. prostrate shape) treatment, giving a total of 20 selected shrubs
per site. For each Juniper individual, we selected three plots of 1 m²
at three different Positions: (i) the center, (ii) the periphery, and
(iii) outside the shrub for the open plots (Appendix S1). The open
plots were located in close vicinity (at an approximatively distance
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | Position of the three natural cedar reserves within Lebanon and corresponding ombrothermic diagrams. (a) Map of Lebanon and position of the three natural reserves. (b) Ombrothermic diagram of Horsh Ehden Nature Reserve. (c) Ombrothermic diagram of Tannourine Cedar Forest Nature Reserve. (d) Ombrothermic diagram of Shouf Cedar Nature Reserve. Alt = Altitude, Lat = Latitude, Long = Longitude, P = precipitation (mm), T = mean temperature (°C). Note that ombrothermic diagrams are made at 1600 m above sea level. Climatic data are from ClimateDT tool (https://ibbr.cnr.it/climate-dt/).</figcaption></figure>

of 1.5 m) of the Juniper individuals in similar microhabitats to the Juniper regarding soil depth, stone amount, and exposure.

From May to June 2023, we identified and counted the num-
ber of seedlings of all woody species in each plot. We obtained
a grand total of 2432 counted recruits belonging to 24 different
woody species in the 240 sampled plots (4 Sites × 2 Canopies × 2
Phenotypes and 3 Positions × 5 replicates).

2.3 | Microclimatic Measurements

In each of the four sites, we recorded air temperature (T, °C) and
air relative humidity (RH, %) in the 12 combinations of the three

treatments (2 Canopies × 2 Phenotypes × 3 Positions modalities). We used hygrochron DS1923-F5 I-buttons data loggers (Maxim Integrated) placed directly at ground level. To avoid heating of the loggers due to direct sun radiation, we designed I-button shelters following recommendations by Maclean et al. (2021).

In each site, we could place only 30 I-buttons at 30 plots, with
15 placed in forests and 15 in gaps. We randomly chose five
Junipers in gaps and five Junipers in forests and placed one I-
button at each of the three positions of all selected Junipers (cen-
ter, periphery and open). From early to mid-July to early August,
temperatures and relative humidity were registered at each site
every 2 min and then averaged per hour during two consecutive
dry summer sunny days with no clouds in 2023.
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | Illustrations of the three natural Cedar reserves in Lebanon. (a and b) Horsh Ehden Nature Reserve. (c and d) Tannourine Cedar Forest Nature Reserve. (e and f) Shouf Cedar Nature Reserve.</figcaption></figure>

Since I-buttons were placed at different days between sites, we used the temperature and relative humidity data provided by the weather stations located in each reserve to correct the data from the loggers as follows. For each hour and each reserve, we used the data of the weather stations to calculate the difference of temperature and relative humidity between the hour and date of the hygrochron records (for instance 4th and 5th August in Shouf) and the hour and date of a chosen reference period (7th and 8th of July, dates of the first records in Ehden High). These differences accounted for potential differences of temperature and relative humidity between the different periods of hygro- chron records due to variation of meteorological conditions. Then we applied these differences to correct the data recorded by the loggers, providing temperature and air humidity as if all loggers had been placed in the different sites during the same 2- day period. Then, we calculated average VPD between 11:00 am and 5:00 pm using the following equation:

$$VPD = ((100 - RH)/100) \times 610.7 \times 10^{7}((7.5 \times T)/(237.3 + T)) \quad (1)$$

where, RH is relative humidity and T is temperature.

## 2.4 | Plant Interaction Indices and Response Groups

For each Juniper shrub, we computed the Relative Interaction Intensity index (RII, Armas et al. 2004) to quantify the impact of the Juniper on the woody species recruits at each of the two Positions within the Juniper patch (center and periphery). RIIs were calculated with the following two equations:

$$RII_{\text{center}}: (X_{\text{Juniper center}} - X_{\text{outside Juniper}}) / (X_{\text{Juniper center}} + X_{\text{outside Juniper}}) \quad (2)$$

$$RII_{\text{periphery}}: (X_{\text{Juniper periphery}} - X_{\text{outside Juniper}}) / (X_{\text{Juniper periphery}} + X_{\text{outside Juniper}}) \quad (3)$$

$X_{\text{Juniper center}}$ represents the number of woody seedlings observed at the center of the Juniper, $X_{\text{Juniper periphery}}$ the number of woody

seedlings observed at its periphery, and $X_{\text{outside Juniper}}$ the number of woody seedlings observed at the open position. RII values vary between -1 and +1, with positive values indicating association, and negative values indicating repulsion. The two RII indices (center and periphery) of the five replicates of the 32 treatment combinations (4 Sites × 2 Canopies × 2 Phenotypes × 2 Positions) were calculated at the community level and at the response- group level. For the community level, we calculated RII indi- ces in each plot using the number of all individuals of woody recruits recorded in the plot. For the response-group level, we selected the 13 most common species and calculated mean RII for each species in each treatment combination. Then we per- formed a PCA on this RII of all species (prcomp base R function; fviz_pca_*; factoextra library), followed by a hierarchical cluster analysis (using Ward's method, Thioulouse et al. 1991) based on the Euclidean distance (Appendix S8). This enabled the char- acterization of four different groups of species, having similar responses to the effect of Juniper shrubs in the different treat- ments. Then, we calculated for each response group, position and in each plot, the RII on the pooled number of all recruits belonging to the same species response group.

## 2.5 | Data Analysis and Statistics

We first performed an ANCOVA to analyze the variation of VPD along our climatic gradient (Site) with the other treatments (Canopy, Phenotype and Position) with (Time), and their inter- actions. Mixed models were used (lmer function, lme4 package) because of the dependency of the values coming from the same data loggers located in the same position during the record.

Second, one-sample t-test were used to assess whether RII sig- nificantly differed from zero in the different combinations of treatments. Then, we analyzed the variation of RII according to the different treatments (Site, Canopy, Phenotype and Position) in a four-way ANOVA at the community level (lmer function, lme4 package). We performed a similar ANOVA at the Response group level, including the Response group (and all possible inter- actions) as a fixed factor in the model, and including the Juniper identity as a random factor (lmer function, lme4 package). To
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | Climatic characterization of the four sites. (a) Diagram of a PCA conducted on four climate variables at eight positions in the four study sites (one position at the altitudinal top (triangle symbols) and one position at the altitudinal bottom (round symbols) position of the four sites). Climatic data are from ClimateDT tool (https://ibbr.cnr.it/clima te-dt/). Blue symbols are for Ehden High, green for Ehden Low, yellow for Tannourine and red for Shouf. (b) Variable diagram of the PCA. Continentality: Rainfall continentality calculated on winter rainfall with the Gams index (Gams 1932), Summer Humidity: Ratio of mean precipitation of the three summer months to mean summer max temperature, T. Min: Minimal January temperature, GSL: Growing season length (number of months with mean max temperature above 12.5°C; see Michalet et al. 2021 for precise methods). (c) Mean (±SE, N=15) VPD values in the gaps of the four sites during zenith hours (from 14:00 to 17:00 pm). Colors are the same as in (a) (see Section 2 for VPD calculation). Different letters indicate significant differences between sites (Tukey post hoc multiple comparisons).</figcaption></figure>

reveal the potential variation of RII of the different Response
groups with climate, and more specifically with VPD, we as-
signed the average VPD of each combination of Canopy and
Site (i.e., eight different values of VPD). More precisely, since
Site and Canopy treatments create large variations of VPD, we
replaced these treatments with corresponding eight different
values of VPD for this analysis corresponding to the mean VPD
(from 11:00 am to 5:00 pm) the 2 days of the records. Then, we
conducted an ANCOVA to analyze the variation of RII accord-
ing to VPD, the Response group, the Position in the shrubs, and
their interactions (lmer function, lme4 package). The treatments
Position and Phenotype were kept in the model despite their
influence on VPD because of their potential additional impor-
tance regarding light and water pre-emption and corresponding
effects on plant-plant interactions (Pescador et al. 2014). For
this ANCOVA, a quadratic term of VPD was included in the
model because of the non-linear relationship observed. During

RII analyses, we used mixed modeling (Juniper identity as ran-
dom factor) to control for the dependence between the two RIIs
(RII<sub>center</sub> and RII<sub>periphery</sub>) calculated for the same shrub.

Third, we analyzed the variation of abundance of the different
response groups in our study (lm base R function). However, this
variation may be due to non-controlled effects like higher disper-
sal of seeds in some sites, in some areas within sites and for some
response groups more than others. Therefore, the three values of
abundance of seedlings for the different Response groups in the
three Positions (center, periphery, open) of each Juniper are not
independent. To ease the statistical analysis (and to avoid con-
vergence issues when using mixed modeling) we first built a glm
model (quasi-Poisson distribution) with the abundance (number
or recruits) as response variable and the identity of the Juniper
and the Response group and their interaction as factors. The vari-
ation of the residuals of this first model (hereafter referred to as
"abundance indices") was then analyzed similarly to the variation of RII (see above). A last simple ANCOVA was performed to fit the relationship between the abundance indices and the VPD according to the Response-group level (lm base R function).

All the analyses were conducted using the R software (R Core Team 2024). During all ANOVA analyses, Tukey post hoc multiple comparisons were performed to detect significant differences between different Response groups when the related factors (or interaction between factors) were significant (emmeans function, emmeans package; cld function, multcomp package). When fitting linear models (ANOVA, ANCOVA), the conditions of applications of linear modeling were checked and validated (normal distribution and homoscedasticity of the residuals).

# 3 | Results

## 3.1 | Microclimatic Variation Along the Gradient in Different Treatment Combinations

VPD was higher in Shouf and Tannourine than in Ehden Low and High; higher in gaps than in forests; and higher in the open intermediate in the periphery and lower in the center of plots, respectively (Site, Canopy and Position effects, Appendices S2 and S3). However, there were also highly significant Site by Canopy and Canopy by Position interactions because the Canopy and Position effects were much stronger at Shouf than at the other three sites (Appendices S3 and S4, respectively). There was also a highly significant Site by Canopy by Position interaction because the stronger Position effect at Shouf was also mostly observed in gaps (data not shown). Finally, there were also highly significant interactions between the Time and other treatments because, depending on sites, the buffering effects of the Canopy and Position treatments were stronger either in the morning or the afternoon.

Overall these results show that there were strong buffering effects on VPD by the forest and Juniper canopies inducing a decrease in VPD stress along the drought gradient from the southernmost site Shouf to the northernmost ones Ehden High and from the periphery of the shrub to its center.

## 3.2 | Changes in RII at the Community Level

Only the Position treatment had a significant effect on RII calculated at the community level, with more negative RIIs at the center than at the periphery of the Juniper (Appendices S5 and S6).

## 3.3 | Species Response Groups: PCA and Cluster Analysis

The PCA axes accounted for 23.4% and 20.9% of the explained variance. The 13 most common species had contrasting responses to the Juniper effects depending on treatments (Appendix S7), allowing separation into four response groups with the cluster analysis (Appendix S8). The PCA diagram showed that axis 1 mostly opposed group IV to groups I and II, with group III at an intermediate position. This was likely due to negative RII values for group IV in all treatments except

T.G.P.P. (Tannourine, gap, prostrate, periphery). Axis 2 mostly opposed group I, at its positive end, to group II at its negative end. This was likely due to more positive RIIs for group I in the erect Juniper phenotype and in the prostrate phenotype for group II (Appendix S7).

## 3.4 | Treatments Effects on RII at the Response-Group Level and Its Relation With VPD

There was a significant effect of the Response group by Canopy interaction (Figure 4, Appendix S9), due to more negative RIIs for group IV but more positive RIIs for the other three groups from forests to gaps (Appendix S9, Figure 4). However, there was an additional significant effect of the Response group by Site by Canopy interaction (Appendix S9), because the increasing facilitation from forests to gaps was only significant for group I at Ehden Low (Figure 5a). Finally, there was a significant effect of the Response group by Site × Phenotype interaction, likely due to the occurrence, for group I only, of opposite increasing and decreasing RII variation from Ehden High to Tannourine in the erect vs. prostrate phenotype, respectively (data not shown).

Furthermore, there was a significant VPD² effect on RII showing a unimodal relationship between VPD and RII. However, this relationship was significant only for group I (with positive maximum RII at intermediate VPD and a decrease to null RII with lower and higher VPD), while RII became more negative with increasing VPD for group IV (Figure 5b). This led to a significant VPD by Response group interaction and a marginally significant Response group by VPD² interaction for RII (Appendix S10). There was an additional marginal interaction between Response group and Position (*p* = 0.07). This was due to a significant effect of the Position on the RII when considering group IV only (*p* < 0.001).

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | Means (±SE, N=10) of the relative interaction intensity index (RII) of the four response groups (blue=group I, green=group II, yellow=group III, red=group IV), with filled bars for forest and hatched bars for gaps. Results of one-sample *t*-test are shown when significant (*p*<0.05; ***p*<0.001). Capital letters represent the results of the Tukey test for the Response group effect and small letters the results of the Tukey test for the Canopy by Response group interaction. See Appendix S9 for full statistical results.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

**FIGURE 5** | (a) Means (±SE, N=10) of the relative interaction intensity index (RII) of the four response groups (see Appendix S3) according to the Site, the Response group and the Canopy treatments. Filled and hatched bars are for RII in forest and gaps, respectively. RII values significantly different from zero are shown above or below bars (*p* < 0.05; **p* < 0.01; ***p* < 0.001, one-sample *t*-test). Small letters represent the results of the Tukey test for the Response group by Site by Canopy interaction. See Appendix S9 for full statistical analysis. (b) Relationship between RII and VPD according to Species-response group and Position. A relationship considering the 2nd order polynomial of the VPD was fitted. Since the interaction between the VPD and Response group was significant, the different fitted lines represent the relationships between the RII and VPD for each group. The continuous line indicates a significant relationship (*p* < 0.05 for the considered group only). Dotted lines are for *p* > 0.05, 90% confidence envelopes are shown with the same color code. The left and right vertical dotted lines represent the lower and upper limits of VPD in the forests and in the gaps, respectively. There was also a marginal interaction between Response group and Position (*p* = 0.07). This was due to a significant effect of the Position on the RII for group IV only (*p* < 0.001). Therefore, we chose to represent separately RII response for group IV according to the Position (light color for periphery and dark color for center). See Appendix S10 for full statistical results. The color code used for the different response groups is as in Figure 4.

## 3.5 | Abundance Index and Its Relation With VPD

There was a highly significant effect of the Response group by Position interaction on the abundance index, because the abundance index strongly decreased from the center of the Juniper to its periphery and to open plots for group I, whereas the converse effect was observed for group IV (Appendix S11, Figure 6a). This indicates that for group I, the number of individuals found increased from the open to the center of the Juniper, while the pattern was opposite for group IV.

Additionally, there was a highly significant effect of the Response group by VPD² interaction on the abundance index because, with increasing VPD, the abundance indices of groups I and IV increased and decreased, respectively (Appendix S12, Figure 6b).

## 4 | Discussion

Consistent with our first prediction, we found that the seedlings of the most common woody species had contrasting responses to the effects of the Juniper, with overall negative responses for *Cedrus* and positive responses for a group of species including *Abies*. Consistent with our second prediction, we found stronger variation in positive than negative responses along our treatments. Consistent with our third prediction, we found that variation in positive responses was highly correlated to variation in VPD, responses being the most positive at intermediate positions along the VPD gradient. Finally, consistent with our fourth prediction, variation in positive responses along the drought

gradient may explain the only occurrence of *Abies* in northern Lebanon. Overall, our results highlight the primary importance of VPD for variation in spatial associations likely indicative of plant-plant interactions along this Mediterranean drought gradient and the potential role of facilitation for explaining tree species distribution along spatial climatic gradients.

### 4.1 | The Contrasting Responses of Woody Species to the Effects of Juniper

We found contrasting associations with Juniper at the species level that could be synthesized in four main species groups. There were two extreme responses, with overall positive RIIs for species of group I (including *Abies*) and negative RIIs for *Cedrus*, the only species of group IV. There were also two intermediate groups that had overall weakly positive or null RIIs. This result is consistent with an increasing number of studies showing that communities are formed of species with different functional strategies that have contrasting responses to environmental conditions and the effects of neighbors (Chapin III and Shaver 1985; Choler et al. 2001; Liancourt et al. 2005; Michalet, Chen, et al. 2015; Bai et al. 2021; Bahalkeh et al. 2024).

### 4.2 | Variation in Positive and Negative Responses Along Our Treatments

The species group-specific responses to Juniper were primarily affected by the presence of the forest canopy, with increasing
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 6** | (a) Means (±SE, N=80) of the abundance index of the four response groups are shown at the three positions relative to the Juniper shrubs. Small letters represent the results of the Tukey test for the Response group by Position interaction. See Appendix S11 for full statistical analysis. (b) Relationship between the abundance index and VPD for the different response groups. Relationships were fitted with 2nd order polynomial of the VPD. Continuous lines indicate a significant relationship between RII and VPD (*p* < 0.05); dotted lines are for non-significant relationships (*p* > 0.05). 90% confidence envelopes are shown with the same color code. See Appendix S12 for full statistical analysis. The color code of the different response groups is as in Figure 4.</figcaption></figure>

positive and negative associations with the shrub from forests to gaps for the likely facilitated species (e.g., *Abies*) and the likely outcompeted one (e.g., *Cedrus*), respectively. These results are consistent with those of Douaihy et al. (2022), conducted in the same system in Lebanon but only in the northernmost and wettest Cedar Reserve of the country. In our study including four sites from the wettest northernmost Ehden Cedar Reserve to the driest southernmost Shouf Cedar Reserve, we also found important variation in RII values among sites in interaction with the presence of the forest canopy.

Consistent with our second prediction, we found stronger variation in positive than negative responses along the drought gradient. Specifically, for group I, we found a unimodal variation in positive association along the drought gradient, in particular in the gaps, with the most positive RII at Ehden Low (second wettest site), null RIIs in the most stressed conditions, that is, at Shouf, and intermediate positive RIIs in the wettest site, Ehden High. Therefore, positive responses peaked at intermediate positions along the drought gradient and waned in extreme conditions of stress, consistent with Michalet et al. (2006) and Holmgren and Scheffer (2010).

In contrast, negative responses of *Cedrus* remained relatively stable across sites, only being slightly more intense in gaps than in forests (Figure 5a). A higher sensitivity of facilitation than competition along a drought gradient and the primary importance of canopy effects as compared to the position of sites along a drought gradient have also been observed by Muhamed et al. (2013) and Saccone et al. (2009) in coastal oak forests and mountain forests, respectively.

## 4.3 | Variation in Positive Responses Are Highly Correlated to Variation in VPD

Douaihy et al. (2022) argued that the mechanism of facilitation for *Abies* should be a mitigation of atmospheric stress (high VPD and extreme temperatures) provided by the shade of the shrub, consistent with a number of studies that highlighted that atmospheric stress mitigation by the shade of neighbors is the primary driver of facilitation in dry systems (Gómez-Aparicio et al. 2005; Saccone et al. 2009; Muhamed et al. 2013; Guignabert et al. 2020; Aguirre et al. 2021). Our microclimatic measurements suggest that the driving mechanisms of the positive response of species of group I to the shrub were very likely a mitigation of the atmospheric stress. This result is consistent with the proposition of Michalet et al. (2014) who suggested that the collapse of facilitation is more likely to be observed along environmental gradients when the direct driving mechanism is a non-resource factor such as VPD, temperature, salinity, or disturbance.

The switch back to a competition pattern is more likely to be observed when the driving mechanism is a resource factor such as light, water, or nutrient availability. Species generally compete when the availability of resources decreases (MacArthur and Wilson 1967; Davis et al. 1998; Maestre et al. 2003), but cannot compete for a non-resource factor. Therefore, it is likely that the negative association of Juniper with *Cedrus* was due to competition for either water or light. Although we did not find an increasingly negative RII for *Cedrus* along the drought gradient and did not quantify soil moisture, water is likely more involved than light in this negative association since RII decreased from forests to gaps.
## 4.4 | The Collapse of Facilitation May Explain the Southernmost Limit of Abies in Northern Lebanon

Michalet et al. (2021) showed that winter rainfall continental-
ity was the primary climatic driver of tree species distribution
in the western United States before temperature and summer
water balance. Winter rainfall continentality is a complex cli-
mate factor for plants that is mostly related to cloudiness, VPD,
and minimum temperature. One of the most oceanic species of
the Western United States, *Sequoia sempervirens*, is restricted
to northern coastal California, characterized by a high cloud-
iness, low VPD, and mild winter temperatures (Waring and
Franklin 1979). Simonin et al. (2009) have shown with an ex-
perimental manipulation of air humidity that the seedlings of
the species were able to support the summer drought of the
Californian Mediterranean climate because of the low VPD
values occurring below the forest canopy. In the case of *Abies*,
this is particularly true in Mediterranean climates, where the
low summer water balance limits VPD-sensitive species to the
most oceanic climates (Michalet et al. 2003, 2021; Quézel and
Médail 2003).

Our results showed that the species abundance index of the
group of *Abies* non-linearly decreased when VPD increased,
whereas *Cedrus* abundance showed the converse pattern.
Therefore, it is likely that the restricted occurrence of *Abies* in
the northernmost Cedar Reserve of Ehden could be explained by
the low VPD values occurring only there, due to the combined
effects of latitude, altitude, exposure, and topography (steep and
deep valleys exposed to the North West). Positive association of
*Abies* with the Juniper was correlated to decreased VPD, and
this effect was the strongest in the gaps of Ehden Low, where
stress was intermediate. At Ehden High, and in particular below
the forest canopy, positive associations decreased likely because
the benefit of shade for atmospheric stress was balanced by its
cost for photosynthesis. When drought stress was the highest,
in the gaps of the driest site, Shouf, positive associations sharply
declined likely because atmospheric stress was too high to be
mitigated by the Juniper, which likely stopped transpiring, as
shown by Saccone et al. (2009) for *A. alba* seedlings during the
2003 European heatwave. We conclude that *Abies cilicica* cer-
tainly thrives in the cool, moist valleys of northern Lebanese
forests of Ehden thanks to facilitation, and is absent southward
because of facilitation collapsing. However, manipulative ex-
periments should now be conducted on the same target species,
and in particular *Abies* and *Cedrus*, to provide straightforward
evidence of such interactions. It is worth also assessing indirect
interactions with herbivores, since the positive association of
*Abies* with Juniper might in part also be due to the protection
of seedlings by the spiny leaves of the Juniper (Rousset and
Lepart 2000).

5 | Conclusion

Relict Cedar forests of Lebanon include woody species of dif-
ferent functional strategies that have contrasting spatial asso-
ciations with the dominant shrub *J. oxycedrus*, likely indicative
of positive or negative interactions; however, manipulative ex-
periments are needed to conclude. We found a gradient in the

effect of Juniper on woody seedlings, from negative over *Cedrus*
to facilitative over *Abies*. Positive effects strongly varied along
our treatments, and in particular along the drought gradient and
with the presence of the forest canopy, whereas negative effects
only slightly increased from forests to gaps. Positive effects var-
ied with VPD, being the most positive at intermediate to high
VPD values and decreasing in the wettest and driest conditions.
The southernmost limit of *Abies cilicica* in the Near East is ob-
served at Horsh Ehden Nature Reserve, where positive effects
are the highest. This suggests that the absence of the species
further south might be due to a decline in facilitation, caused
by excessive dry conditions, like in the Tannourine and Shouf
Reserves. *Cedrus libani* appeared more tolerant to high VPD
but less to low light and water availability since it is likely ex-
cluded by Juniper, in particular in the gaps. However, we did
not find evidence that the southernmost limit of *Cedrus* found
in Lebanon can be due to plant-plant interactions, in particu-
lar competition. Our results bring important knowledge for the
management and restoration of the relict coniferous forests of
Lebanon, in particular in the context of climate change. Beyond
this, our study brings crucial insights on the role of plant-plant
interactions for large-scale tree species distribution.

Author Contributions

R.M. and F.D. conceived the research idea and design. A.M. and R.M. conducted the fieldwork and data collection. A.M. and F.D. performed statistical analyses. A.M., R.M., and F.D. wrote the paper.

Acknowledgments

We thank Nayla Darazi Maamary, Georges Maamary, Sophia Maamary, Vanessa Maamary, Jean Beainy, and Ghinwa Kanaan for their help during fieldwork. We also thank David Nemer for his help in statistical analysis. We show our gratitude to Horsh Ehden Nature Reserve, Tannourine Cedar Forest Nature Reserve, and Shouf Cedar Nature Reserve forest managers, as well as their teams, for their welcoming and help in providing us additional data and information; we thank the responsibles of the correspondent weather stations, as well as the LARI (Kfarshakhna—Lebanon) for additional macroclimatic data.

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

The data that support the findings of this study are openly available in
Figshare at https://doi.org/10.6084/m9 figshare.28044908.v1, reference
number 28044908.v1.

References

Aguirre, B. A., B. Hsieh, S. J. Watson, and A. J. Wright. 2021. “The Experimental Manipulation of Atmospheric Drought: Teasing out the Role of Microclimate in Biodiversity Experiments.” *Journal of Ecology* 109: 1986–1999. https://doi.org/10.1111/1365-2745.13595.

Armas, C., R. Ordiales, and F. I. Pugnaire. 2004. “Measuring Plant Interactions: A New Comparative Index.” *Ecology* 85: 2682–2686. https://doi.org/10.1890/03-0650.

Bagnoul, F., and H. Gaussen. 1953. "Saison Sèche et Indice Xérothermique." *Bulletin de La Société D'histoire Naturelle De Toulouse* 88: 193-239.
Bahalkeh, K., M. Abedi, G. A. Dianati Tilaki, and R. Michalet. 2024. "Species-Group Responses Improve Our Understanding of the Effects of Community Dominants on Subordinate Species Along a Grazing Gradient." *Applied Vegetation Science* 27: e12789. https://doi.org/10.1111/avsc.12789.

Bai, Y., R. Michalet, W. She, et al. 2021. "Contrasting Responses of Different Functional Groups Stabilize Community Responses to a Dominant Shrub Under Global Change." *Journal of Ecology* 109: 1676-1689. https://doi.org/10.1111/1365-2745.13588.

Bertness, M. D., and R. Callaway. 1994. "Positive Interactions in Communities." *Trends in Ecology & Evolution* 9: 191-193. https://doi.org/10.1016/0169-5347(94)90088-4.

Bertness, M. D., L. A. Cavieres, C. J. Lortie, and R. M. Callaway. 2024. "Positive Interactions and Interdependence in Communities." *Trends in Ecology & Evolution* 39: 1014-1023. https://doi.org/10.1016/j.tree.2024.09.003.

Chaieb, G., C. Abdelly, and R. Michalet. 2021. "A Regional Assessment of Changes in Plant-Plant Interactions Along Topography Gradients in Tunisian Sebkhas." *Ecosystems* 24, no. 5: 1024-1037. https://doi.org/10.1007/s10021-020-00567-8.

Chapin, F. S., III, and G. R. Shaver. 1985. "Individualistic Growth Response of Tundra Plant Species to Environmental Manipulations in the Field." *Ecology* 66: 564-576. https://doi.org/10.2307/1940405.

Chesson, P. 2000. "Mechanisms of Maintenance of Species Diversity." *Annual Review of Ecology and Systematics* 31: 343-366. https://doi.org/10.1146/annurev.ecolsys.31.1.343.

Choler, P., R. Michalet, and R. M. Callaway. 2001. "Facilitation and Competition on Gradients in Alpine Plant Communities." *Ecology* 82: 3295-3308. https://doi.org/10.1890/0012-9658(2001)082[3295:FACOGI]2.0.CO;2.

Davis, M. A., K. J. Wrage, and P. B. Reich. 1998. "Competition Between Tree Seedlings and Herbaceous Vegetation: Support for a Theory of Resource Supply and Demand." *Journal of Ecology* 86: 652-661. https://doi.org/10.1046/j.1365-2745.1998.00087.x.

De Frenne, P., F. Zellweger, F. Rodriguez-Sánchez, et al. 2019. "Global Buffering of Temperatures Under Forest Canopies." *Nature Ecology & Evolution* 3: 744-749. https://doi.org/10.1038/s41559-019-0842-1.

Douaihy, B., A. Maamary, R. Yammine, et al. 2022. "Forest Canopy and Mistletoe Infestation Alter the Facilitative Effects of *Juniperus oxycedrus* SI on Woody Seedlings on Mount Lebanon (Lebanon)." *Journal of Vegetation Science* 33: e13163. https://doi.org/10.1111/jvs.13163.

Gams, H. 1932. "Die Klimatische Begrebzung Von Pflanzenarealen Und Die Verteilung Der Hygrischen Kontinentalität in Den Alpen." *Zeitschrift Der Geselchafften Für Erkunde* 30: 205-210.

Gómez-Aparicio, L., J. M. Gómez, R. Zamora, and J. L. Boettinger. 2005. "Canopy vs. Soil Effects of Shrubs Facilitating Tree Seedlings in Mediterranean Montane Ecosystems." *Journal of Vegetation Science* 16: 191-198. https://doi.org/10.1111/j.1654-1103.2005.tb02355.x.

Grime, J. P. 1973. "Competitive Exclusion in Herbaceous Vegetation." *Nature* 242: 344-347. https://doi.org/10.1038/242344a0.

Guignabert, A., L. Augusto, M. Gonzalez, C. Chipeaux, and F. Delerue. 2020. "Complex Biotic Interactions Mediated by Shrubs: Revisiting the Stress-Gradient Hypothesis and Consequences for Tree Seedling Survival." *Journal of Applied Ecology* 57: 1341-1350. https://doi.org/10.1111/1365-2664.13641.

Holmgren, M., and M. Scheffer. 2010. "Strong Facilitation in Mild Environments: The Stress Gradient Hypothesis Revisited." *Journal of Ecology* 98: 1269-1275. https://doi.org/10.1111/j.1365-2745.2010.01709.x.

Kitzberger, T., D. F. Steinaker, and T. T. Veblen. 2000. "Effects of Climatic Variability on Facilitation of Tree Establishment in Northern Patagonia." *Ecology* 81: 1914-1924. https://doi.org/10.2307/177281.

Liancourt, P., R. M. Callaway, and R. Michalet. 2005. "Stress Tolerance and Competitive-Response Ability Determine the Outcome of Biotic Interactions." *Ecology* 86: 1611-1618. https://doi.org/10.1890/04-1398.

Liancourt, P., Y. Le Bagousse-Pinguet, C. Rixen, and J. Dolezal. 2017. "SGH: Stress or Strain Gradient Hypothesis? Insights From an Elevation Gradient on the Roof of the World." *Annals of Botany* 120: 29-38. https://doi.org/10.1093/aob/mcx037.

Lortie, C. J., R. W. Brooker, P. Choler, et al. 2004. "Rethinking Plant Community Theory." *Oikos* 107: 433-438. https://doi.org/10.1111/j.0030-1299.2004.13250.x.

MacArthur, R. H., and E. Wilson. 1967. *The Theory of Island Biogeography*. Princeton University Press.

Maclean, I. M., J. P. Duffy, S. Haesen, et al. 2021. "On the Measurement of Microclimate." *Methods in Ecology and Evolution* 12: 1397-1410. https://doi.org/10.1111/2041-210X.13627.

Maestre, F. T., S. Bautista, and J. Cortina. 2003. "Positive, Negative, and Net Effects in Grass-Shrub Interactions in Mediterranean Semiarid Grasslands." *Ecology* 84: 3186-3197. https://doi.org/10.1890/02-0635.

Maestre, F. T., and J. Cortina. 2004. "Do Positive Interactions Increase With Abiotic Stress? A Test From a Semi-Arid Steppe." *Proceedings of the Royal Society of London. Series B: Biological Sciences* 271: S331-S333. https://doi.org/10.1098/rsbl.2004.0181.

Michalet, R. 2006. "Is Facilitation in Arid Environments the Result of Direct or Complex Interactions." *New Phytologist* 169: 3-6. https://doi.org/10.1111/j.1468-8137.2006.01617.x.

Michalet, R., R. W. Brooker, L. A. Cavieres, et al. 2006. "Do Biotic Interactions Shape Both Sides of the Humped-Back Model of Species Richness in Plant Communities?" *Ecology Letters* 9: 767-773. https://doi.org/10.1111/j.1461-0248.2006.00935.x.

Michalet, R., S. Y. Chen, L. Z. An, et al. 2015. "Communities: Are They Groups of Hidden Interactions?" *Journal of Vegetation Science* 26: 207-218. https://doi.org/10.1111/jvs.12226.

Michalet, R., P. Choler, R. M. Callaway, and T. G. Whitham. 2021. "Rainfall Continентality, via the Winter Gams Angle, Provides a New Dimension to Biogeographical Distributions in the Western United States." *Global Ecology and Biogeography* 30: 384-397. https://doi.org/10.1111/geb.13223.

Michalet, R., Y. Le Bagousse-Pinguet, J. P. Maalouf, and C. J. Lortie. 2014. "Two Alternatives to the Stress-Gradient Hypothesis at the Edge of Life: The Collapse of Facilitation and the Switch From Facilitation to Competition." *Journal of Vegetation Science* 25: 609-613. https://doi.org/10.1111/jvs.12123.

Michalet, R., G. Losapio, Z. Kikvidze, et al. 2024. "Assessing the Accuracy of Paired and Random Sampling for Quantifying Plant-Plant Interactions in Natural Communities." *Population Ecology* 66: 39-52. https://doi.org/10.1002/1438-390X.12161.

Michalet, R., J. P. Maalouf, P. Choler, et al. 2015. "Competition, Facilitation and Environmental Severity Shape the Relationship Between Local and Regional Species Richness in Plant Communities." *Ecography* 38: 335-345. https://doi.org/10.1111/ecog.01106.

Michalet, R., C. Rolland, D. Joud, D. Gafta, and R. M. Callaway. 2003. "Associations Between Canopy and Understory Species Increase Along a Rainshadow Gradient in the Alps: Habitat Heterogeneity or Facilitation?" *Plant Ecology* 165: 145-160. https://doi.org/10.1023/A:1022297624381.

Muhamed, H., B. Touzard, Y. Le Bagousse-Pinguet, and R. Michalet. 2013. "The Role of Biotic Interactions for the Early Establishment of Oak Seedlings in Coastal Dune Forest Communities." *Forest Ecology and Management* 297: 67-74. https://doi.org/10.1016/j.foreco.2013.02.023.

Novick, K. A., D. L. Ficklin, P. C. Stoy, et al. 2016. "The Increasing Importance of Atmospheric Demand for Ecosystem Water and Carbon
Fluxes.” *Nature Climate Change* 6: 1023–1027. https://doi.org/10.1038/
nclimate3114.

Pescador, D. S., J. Chacon-Labella, M. de la Cruz, and A. Escudero. 2014.
“Maintaining Distances With the Engineer: Patterns of Coexistence in
Plant Communities Beyond the Patch-Bare Dichotomy.” *New Phytologist*
204: 140–148. https://doi.org/10.1111/nph.12899.

Quézel, P., and F. Médail. 2003. *Ecologie Et Biogéographie Des Forêts Du Bassin Méditerranéen.* R Core Team.

R Core Team. 2024. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing. https://www.R-project.org/.

Rousset, O., and J. Lepart. 2000. “Positive and Negative Interactions at Different Life Stages of a Colonizing Species.” *Journal of Ecology* 88: 401–412. https://doi.org/10.1046/j.1365-2745.2000.00457.x.

Saccone, P., S. Delzon, J. P. Pagès, J. J. Brun, and R. Michalet. 2009.
“The Role of Biotic Interactions in Altering Tree Seedling Responses to
an Extreme Climatic Event.” *Journal of Vegetation Science* 20: 403–414.
https://doi.org/10.1111/j.1654-1103.2009.01012.x.

Simonin, K. A., L. S. Santiago, and T. E. Dawson. 2009. “Fog Interception by *Sequoia sempervirens* (D. Don) Crowns Decouples Physiology From Soil Water Deficit.” *Plant, Cell & Environment* 32: 882–892. https://doi.org/10.1111/j.1365-3040.2009.01967.x.

Soliveres, S., C. Smit, and F. T. Maestre. 2015. “Moving Forward on Facilitation Research: Response to Changing Environments and Effects on the Diversity, Functioning and Evolution of Plant Communities.” *Biological Reviews* 90: 297–313. https://doi.org/10.1111/brv.12110.

Steinbauer, M. J., C. Beierkuhnlein, M. A. S. A. Kahn, et al. 2016.
“How to Differentiate Facilitation and Environmentally Driven Co-
Existence.” *Journal of Vegetation Science* 27: 1071–1079. https://doi.org/
10.1111/jvs.12441.

Stephenson, N. L. 1990. "Climatic Control of Vegetation Distribution: The Role of the Water Balance." *American Naturalist* 135: 649-670.
https://doi.org/10.1086/285067.

Thioulouse, J., J. Devilliers, D. Chessel, and Y. Auda. 1991. "Graphical Techniques for Multidimensional Data Analysis." In *Applied Multivariate Analysis in SAR and Environmental Studies*, edited by J. Devillers and W. Karcher, 153-205. Kluwer Academic.

Tielbörger, K., and R. Kadmon. 2000. “Temporal Environmental Variation Tips the Balance Between Facilitation and Interference in Desert Plants.” *Ecology* 81: 1544–1553. https://doi.org/10.2307/177305.

Waring, R. H., and J. F. Franklin. 1979. “Evergreen Coniferous Forests of the Pacific Northwest: Massive Long-Lived Conifers Dominating These Forests Are Adapted to a Winter-Wet, Summer-Dry Environment.” *Science* 204: 1380–1386. https://doi.org/10.1126/science.204.4400.1380.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section. **Appendix S1:** Detailed design of a Juniper patch. **Appendix S2:** Results of the four-way ANOVA on the Vapor Pressure Deficit (VPD). **Appendix S3:** VPD variations according to the Site, Canopy and Time treatments. **Appendix S4:** VPD variations according to the Canopy, Position and Time treatments. **Appendix S5:** Results of the four-way ANOVA conducted on RII at the community level. **Appendix S6:** Effect of the Position on RII at the community level. **Appendix S7:** Results of the PCA conducted on the responses of the 13 most common species. **Appendix S8:** Dendrogram of the cluster analysis conducted on 13 most common species. **Appendix S9:** Results of the five-way ANOVA conducted on RII at the response-group level. **Appendix S10:** Results of ANCOVA test for the effects of the VPD (and its 2nd order polynomial), Response group, Position and their interactions on RII at the response-group level. **Appendix S11:** Results of the ANCOVA test for the effects of the VPD (and its 2nd order polynomial), Response group, Position and their interactions on the abundance index.

**Appendix S12:** Results of the ANCOVA test for the effects of the VPD (and its 2nd order polynomial) and its interaction with the Response group on the abundance index.