Introduction

An understanding of the spatial distribution of diversity in plant communities constitutes essential information for conservation planning (Socolar et al. 2016). At local scales, community assembly and species coexistence are determined primarily by biotic interactions among neighboring individuals (Pearson and Dawson 2003; Götzenberger et al. 2012; Mod et al. 2020; King et al. 2021). Traditionally, competition for resources and space has been considered the most important interaction involved in community organization (Grime 1973; Tilman 1982; Chesson 2000). In this respect, only the best-adapted and competitive species are able to persist in the community, leading to the exclusion of the less competitive ones. Yet competition can also favor species coexistence, as it is a strong force in species differentiation (Chesson 2000; Bruno et al. 2003). However, over the last two decades, a growing body of studies has acknowledged the importance of facilitation in shaping the structure and composition of plant communities, especially those living under more severe environmental conditions (Bruno et al. 2003; Brooker et al. 2008; McIntire and Fajardo 2014; Michalet and Pugnaire 2016). The presence of species with a predominance of facilitative interactions (i.e., nurse species) in the community enables the persistence of species that usually are poor competitors, less-well adapted to harsh abiotic conditions or have a low reproductive performance (Soliveres et al. 2015a; O’Brien et al. 2019; Calatayud et al. 2020); thus, promoting and sustaining local diversity and species coexistence (Le Bagousse-Pinguet et al. 2014; Cavieres et al. 2014). More specifically, nurse species facilitate other species by expanding their niche through microhabitat amelioration (Bruno et al. 2003; Callaway 2007; O’Brien et al. 2019) and by enhancing their populations (Soliveres et al. 2015a). Therefore, plant–plant interactions constitute the so called “biotic filter” (Götzenberger et al. 2012) that selects which species can appear in the vicinity of each species in the community, shaping the local pattern of diversity and the subsequent spatial heterogeneity in species composition.

Accordingly, one of the most common methods for identifying interspecific plant–plant interactions in natural communities is to test whether plant species are spatially independent from each other (Soliveres and Maestre 2014; López et al. 2016; Alados et al. 2017; Saiz et al. 2017; O’Brien et al. 2019). Specifically, a positive net interaction (i.e. facilitation) between two or more species can be inferred when they aggregate in space more than could be expected by chance. Conversely, a negative net interaction (i.e. competition or allelopathy) between two or more species can be inferred when they segregate in space more than expected. However, spatial patterns in plant species are also determined by other ecological factors apart from plant–plant interactions, including fine-scale environmental and resources heterogeneity (Pescador et al. 2020), herbivory (Saiz and Alados 2012), seed dispersal limitation and other stochastic processes (Alados et al. 2010; Chacón‐Labella et al. 2017; D’Amen et al. 2018). Thus, inferring the magnitude and direction of biotic interactions from the analysis of species patterning requires the use of suitable null models that can separate processes operating at different scales (Wiegand and Moloney 2014). This is particularly relevant when multiple plant communities are compared. However, the few studies evaluating the effect exerted by each species in the community in shaping the local pattern of diversity (i.e., diversity aggregation or segregation in their vicinity) as a result of plant–plant interactions are often performed on one or two plant communities (Wiegand et al. 2007; Rayburn and Wiegand 2012; e.g., Arroyo et al. 2015; Chacón‐Labella et al. 2016), rather than including multiple plant communities along a regional environmental gradient.

Importantly, the net balance between positive and negative interactions depends on environmental conditions. In particular, the frequency of positive interactions is expected to rise as the level of environmental severity increases, while under milder conditions communities are expected to be governed by negative interactions (Bertness and Callaway 1994; He et al. 2013). Nevertheless, some studies have found a collapse of facilitation under extreme environmental conditions (Maestre et al. 2009; Holmgren and Scheffer 2010; Michalet et al. 2014; Alados et al. 2017). Future climate change scenarios foresee an increase in aridity, especially in arid and semiarid ecosystems (IPCC 2014), where positive plant interactions are closely related to the structure of plant diversity (Saiz et al. 2018) and their decline can precede a diversity loss (Valiente‐Banuet et al. 2015; Berdugo et al. 2020). In addition, biodiversity is an important biotic attribute that enhances key ecosystem functions, such as productivity and nutrient cycling (Hooper et al. 2005; Maestre et al. 2012, 2016). Therefore, ecosystem management strategies should focus primarily on identifying species involved in positive interactions in order to prevent biodiversity loss. Usually, shrubs act as nurses for other species in the community by creating and maintaining diversity islands, and have been a successful tool in restoration practices (Gómez-Aparicio et al. 2004; Padilla and Pugnaire 2006; Gómez-Aparicio 2009; Liczner et al. 2019). However, most research has focused on a few shrub species of interest and their species-specific effects, while there are much fewer studies taking a multi-species approach to assess their general role in structuring diversity at community level (Rayburn and Wiegand 2012; Soliveres and Maestre 2014; Chacón‐Labella et al. 2016).

Beyond traditional local diversity metrics, beta-diversity measures the spatial variation in the composition of species among assemblages (Whittaker 1960). Recent developmental framework proposes that beta-diversity can be decomposed into two additive components, namely species turnover and nestedness, which describe two different ecological processes causing spatial heterogeneity in species composition (Baselga 2010; Legendre 2014; Baselga and Leprieur 2015). Specifically, species turnover measures the replacement of some species by others, while nestedness measures the loss (or gain) of species as a consequence of some assemblages being subsets of richer assemblages (Baselga 2010; Legendre 2014; Socolar et al. 2016). Beta-diversity metrics are normally used to assess the spatial heterogeneity in species composition at regional scale (Umaña et al. 2020; Fontana et al. 2020). Nevertheless, quantification of how the local patterns of beta-diversity and its components change within natural plant communities may provide further insight into the effect of plant–plant interactions in local diversity organization. For instance, although species turnover is usually the largest component of beta-diversity (Soininen et al. 2018), in plant communities where nurse species are present an increase in the contribution of the nestedness component to the overall local beta-diversity can be expected, due to beneficiary species appearing mostly in microsites created by nurse plants. However, to our knowledge, the local pattern of beta-diversity has seldom been explicitly addressed within plant communities.

In this study we aimed to investigate, using a multi-species approach, the spatial structure of diversity in a set of Mediterranean plant communities, ranging from semi-arid to subalpine habitats in Spain, in order to identify species playing a key role in shaping local diversity through plant–plant interactions. Specifically, we assessed (1) the spatial pattern of diversity (i.e., diversity aggregation or segregation) in the local neighborhood of perennial plant species and (2) the local pattern of beta-diversity (i.e., spatial heterogeneity in species composition among local assemblages), including the contribution of species turnover and nestedness components to overall local beta-diversity. We expected that (i) the proportion of species aggregating diversity in their local neighborhood would be higher than the proportion of species segregating diversity in plant communities subjected to more severe environmental conditions, as a consequence of the prevalence of facilitative interactions in these communities. Furthermore, we expected that species with different life forms, and thus with different biological traits, would exhibit different effects in shaping the local spatial pattern of diversity. We particularly expected that (ii) shrubs and dwarf shrubs would aggregate diversity in their local neighborhood due to their prominent role as nurse species, while other plant groups would mainly act as diversity segregators. Finally, we expected that (iii) in those plant communities under more severe environmental conditions, where positive plant–plant interactions prevail, local species assemblages would be characterized by high nestedness due to the presence of key nurse species that promote diversity islands.

Materials and methods

Study area and vegetation surveys

This study was conducted in 14 well-preserved Mediterranean plant communities distributed across four regions in Spain (Cabo de Gata-Níjar, middle Ebro Valley, Sierra Nevada and Sierra de Guara; Fig. 1). Study sites had varying environmental conditions (Table1). For instance, mean annual rainfall ranged from 250 to 1150 mm, mean annual minimum temperature ranged from 1.6 to 13.8 °C and elevation ranged from 125 to 2270 m.a.s.l. (Table 1). Productivity, more precisely aboveground biomass, is considered to be closely related to the overall environmental severity (Grime 1973; Michalet et al. 2006, 2021). In fact, most studies use biomass as surrogate for productivity (Fraser et al. 2015; Guo et al. 2019; Li et al. 2020; Michalet et al. 2021). Accordingly, we defined a productivity gradient based on perennial plant species biomass (but see Gillman and Wright 2006; Jenkins 2015). Specifically, for each study site, aboveground biomass was computed as the sum of every perennial plant species height × cover (i.e., species volume). Thus, the study site with the lowest plant biomass represented the extreme of the gradient with the highest level of environmental severity. Conversely, sites with higher biomass represented the milder extreme of the gradient (Michalet et al. 2006, 2021).

Fig. 1
figure 1

Locations of study sites in Spain

Full size image
Table 1 Location, year of field survey, topographic (elevation, slope and soil type), climatic (mean annual rainfall, annual mean of weekly maximum and minimum temperatures and aridity index) and biotic (species richness, perennial plant species cover and Pielou’s evenness index) parameters for each study site
Full size table

Vegetation in Cabo de Gata-Níjar and the middle Ebro Valley sites consisted of relatively open steppes dominated by dwarf shrubs and perennial grasses. In the Sierra Nevada sites, vegetation consisted of subalpine shrublands dominated by cushion plant species, and in Sierra de Guara vegetation consisted of a Mediterranean deciduous forest dominated by Quercus faginea. At each study site, vegetation was surveyed by performing three 500-m linear transects following the point-intercept method (Goodall 1952). Specifically, we recorded the presence of all species found at evenly located points along the transects, without distinguishing plant ontogenetic stages. Points were located every 20 cm in steppes and shrublands (11 sites, 2500 points per transect), and every 40 cm in forests (three sites, 1250 points per transect).

Local pattern of diversity

We analyzed the spatial pattern of diversity in the local neighborhood of all perennial plant species occurring more than 75 times within the set of transects carried out at each study site (i.e., perennial species with a cover greater than 1% in steppes and shrublands or 2% in forests; referred to as target species). Specifically, the ISAR method (individual species–area relationship; Wiegand et al. 2007) was used to evaluate the expected number of species j found within a distance d from the target species t,

$${\mathrm{ISAR}}_{(d)} = \sum_{j=1}^{\mathrm{S}}\left[1- {P}_{tj}\left(0, d\right)\right],$$

where 1 − Ptj(0, d) is the probability that species j was present in the vicinity (d) of individuals of target species t. Thus, the ISAR(d) value for such target species would be the sum of the probabilities for all species j in the plant community. The ISAR(d) was computed along the 1500-m transect resulting from combining the three 500-m transects performed at each study site to ensure enough statistical power, and to a maximum distance d of 4 m (i.e., distance within plant–plant interactions are expected to occur; Arroyo et al. 2015).

To determine whether the ISAR(d) value for a given target species was different from that expected by chance, we computed confidence envelopes by performing a Monte Carlo test with 199 heterogeneous Poisson null model simulations (Wiegand et al. 2007). The heterogeneous Poisson null model accounts for the effect of large-scale habitat association (Wiegand et al. 2007; Wiegand and Moloney 2014), thus controlling for potential effects of soil heterogeneity driving plant spatial patterns. Each simulation randomized the position of target species in the transect within the maximum distance d (i.e. 4 m), while maintaining their spatial distribution at longer distances. When the ISAR(d) value was greater than the fifth greatest simulated value, the local neighborhood of target species was more diverse than expected by chance (p ≈ 0.05). Conversely, when the ISAR(d) value was lower than the fifth lowest simulated value, the local neighborhood was less diverse than expected (p ≈ 0.05). If the ISAR(d) values fell within simulated envelopes, then the local pattern of diversity was as diverse as expected by the null model (Wiegand et al. 2007). Thus, diversity aggregation or segregation in the local neighborhood of target species can be regarded as a consequence of the predominance of positive or negative plant–plant interactions respectively, whereas a predominance of weak, or a neutral net balance of positive and negative interactions, leads to a local neighborhood as diverse as expected by chance (Wiegand et al. 2007; Arroyo et al. 2015; Chacón‐Labella et al. 2016).

Later, we computed the proportion of target species with a local neighborhood more diverse, less diverse and as diverse as expected by the null model, per study site and distance d. In addition, each target species was arranged as a diversity aggregator, segregator or neutral based on the observed pattern of diversity below d = 1 m, since most departures of the ISAR curve from the confidence envelopes were found below this distance. Also, each target species was classified into one of six groups based on their biological traits and life form: trees, climbing phanerophytes, shrubs, dwarf shrubs (including cushion species), graminoids (i.e., Poaceae, Juncaceae and Cyperaceae) and forbs (i.e., the rest of perennial herbaceous species). Thus, the local pattern of diversity for each group was summarized as:

$$ {\text{Local pattern of diversity }}(d < 1\;{\text{m}}) = { }\frac{{{\text{n}^\circ}_{{{ }diversity \,\, aggregators}} { } - {\text{ n}^\circ}_{{{ }diversity \,\, segregators}} }}{{{\text{n}^\circ}_{{{ }target \,\, species}} }}, $$

This index ranges between 1 and − 1. Positive values indicate that species within the group predominantly aggregated diversity, while negative values indicate that they predominantly segregated diversity. Values around zero indicate that species within the group had a neutral net effect on diversity, due to a predominance of weak interactions or a similar proportion of diversity aggregators and segregators.

ISAR analysis and Monte Carlo simulations were performed with MATLAB R2013a.

Local pattern of beta-diversity

At each study site, we measured local beta-diversity as the spatial variation in species composition among local assemblages (i.e., species assemblages co-occurring with target species). More specifically, we computed incidence based beta-diversity, species turnover (i.e., compositional heterogeneity due to species replacement among local assemblages) and nestedness (i.e., compositional heterogeneity due to differences in species richness among local assemblages) metrics following the framework in Baselga (2010). To control for the effect of target species abundance on the composition of local assemblages, we built 999 co-occurrence matrices (target species in rows and presence-absence of co-occurring species in columns) per study site using a resampling procedure. Each matrix was built from 75 randomly selected records per target species, according to the minimum cover used to define a target species. A final matrix of co-occurrences for each study site was obtained by including only those species co-occurring with target species in more than 95% of cases.

When the attribute of interest is the overall heterogeneity in species composition among a pool of assemblages (e.g., target species), multiple site dissimilarity measures are recommended (Baselga 2013). Therefore, beta-diversity and its components were computed using the beta.multi function in the betapart package (Baselga et al. 2018) in R (R Core Team 2020). This function quantifies overall beta-diversity as Sørensen-based multiple-site dissimilarity (βSOR), species turnover as Simpson-based multiple-site dissimilarity (βSIM) and nestedness as βSOR − βSIM (βNES). These indices range between 0, meaning a lack of heterogeneity in species composition, and 1, indicating maximum heterogeneity in species composition among local assemblages. Values computed for sites with a different number of target species were made comparable by using a resampling procedure, taking 8 (according to the site with the fewest target species) random target species per site 999 times and calculating the averaged value (beta.sample function; Baselga 2010).

Finally, the relative contribution of species turnover and nestedness to the overall local beta-diversity at each study site was calculated as:

$$ {\text{Relative Contribution}} = { }\frac{Species\;\; Turnover - Nestedness}{{Beta \;\;Diversity}} $$

This index ranges between 1 and − 1. Positive values indicate a higher contribution of the species turnover component to local beta-diversity, while negative values show a predominance of the nestedness component.

Statistical analyses

The effect of environmental conditions and target species life form, as well as their interaction, in the local pattern of diversity was evaluated by fitting a linear mixed model (LMM). Study site was included as a random factor in the model. Selection of the best random effect structure (i.e., random intercept effect vs random slope effect) was based on Akaike´s Information Criterion (> 2 AICc points difference; Zuur et al. 2009). Aboveground biomass was used as the indicator of environmental conditions. Post hoc Tukey´s honest significant difference (HSD) test was used to detect significant differences between pairs of plant life form groups.

Linear models were fitted to assess the significance of the effect of environmental conditions on beta-diversity metrics (i.e., local beta-diversity, species turnover, nestedness and the relative contribution of species turnover and nestedness to local beta-diversity). Aboveground biomass was used in the models to define environmental conditions. On the other hand, the significance of the relationship between the local pattern of diversity and beta-diversity metrics was tested, for each target species life form group, by fitting separate linear models. The region where study sites are located was not included as a random factor in any of these models, since no significant improvement was found (< 2 AICc points difference; Zuur et al. 2009).

Trees and climbing phanerophytes were excluded from data analysis because they were only present in few study sites. All statistical analyses were performed in R 4.0.3 (R Core Team 2020).

Results

Spatial pattern of diversity in the local neighborhood of target species

The local pattern of diversity was analyzed for 190 target species (Supplementary Figure S1), ranging between 8 (Planerón, middle Ebro Valley) and 19 (Rapun, Sierra de Guara) per study site. We found target species with a diversity pattern in its local neighborhood significantly different than expected by the null model in all study sites (Supplementary Figure S2; Table 2). The distance d at which such departures were observed (i.e., scale at which plant interactions occur) differed among target species, but, overall, they were found at d < 1 m, while significant departures of the ISAR curve at d > 1 m were rare (Supplementary Figure S2). Accordingly, the prevalence of diversity segregation was particularly relevant at d = 0 (i.e., at those points at which different species co-occur; Supplementary Figure S2). Conversely, the proportion of target species with a local pattern of diversity as diverse as expected by the null model prevailed at d > 0.5 m for all study sites (Supplementary Figure S2).

Table 2 Local pattern of diversity (d < 1 m) and relative abundance (%) of target species at each study site
Full size table

We found a significant relationship between the spatial pattern of diversity in the local neighborhood of target species and the productivity gradient (F1,12 = 12.81, p < 0.01). Specifically, a predominance of diversity segregation was observed in more productive environments, whereas the local pattern of diversity tended to be neutral or slightly aggregated under less productive conditions (Fig. 2).

Fig. 2
figure 2

Relationship between the local pattern of diversity and the productivity gradient defined as aboveground biomass (log-transformed for better visualization). Positive values indicate a predominance of diversity aggregation. Negative values indicate a predominance of diversity segregation. Different letters indicate significant differences among plant life forms (Tukey´s HSD test). Trees and climbing phanerophytes were not included in data analysis

Full size image

Significant differences in the local pattern of diversity were found among plant life forms (F3,31 = 4.14, p < 0.05). Post hoc Tukey’s HSD test showed that diversity was more segregated in the local neighborhood of graminoids than in the surroundings of forbs species, while shrubs and dwarf shrubs had an intermediate behavior with respect to diversity (Fig. 2). On the other hand, we did not find a significant interaction between plant life form and environmental conditions (F3,31 = 1.40, p = 0.26), indicating that all plant groups showed similar responses along the productivity gradient (Fig. 2).

Compositional heterogeneity of species assemblages associated with target species

We found that neither local beta-diversity, nor species turnover, nor nestedness varied significantly along the productivity gradient (Fig. 3A–C). Species turnover was the beta-diversity component with the largest contribution to the overall spatial heterogeneity in species composition among local assemblages, and this effect did not vary significantly with environmental conditions (Fig. 3D).

Fig. 3
figure 3

Relationship between A local beta-diversity, B species turnover, C nestedness and D relative contribution of species turnover and nestedness to the overall local beta-diversity, and the productivity gradient defined as aboveground biomass (log-transformed for better visualization). In D, positive values indicate a predominance of species turnover over the nestedness component

Full size image

On the other hand, we found a significant relationship between the nestedness component of beta-diversity and the pattern of diversity in the local neighborhood of shrub species (F1,9 = 6.13, p < 0.05). Specifically, a higher nestedness among local assemblages (i.e., higher richness differences among species assemblages co-occurring with target species) was associated with more aggregated (or less segregated) patterns of diversity in the local neighborhood of shrub species (Fig. 4C). Non-significant relationships between beta-diversity metrics and the local pattern of diversity were found for other beta-diversity components and plant life forms (p > 0.05).

Fig. 4
figure 4

Relationship between A local beta-diversity, B species turnover, C nestedness and D relative contribution of species turnover and nestedness to the overall local beta-diversity, and the pattern of diversity in the local neighborhood of shrubs species. In D, positive values indicate a predominance of species turnover over the nestedness component

Full size image

Discussion

The purpose of this study was to investigate, from a multi-species perspective, how biotic interactions of different plant species and plant life forms determine the spatial structure and composition of local diversity in a range of Mediterranean plant communities. We found that diversity segregation in the local neighborhood of target species decreased in the less productive environments. This finding is in line with a reduction in competition in plant communities under more severe environmental conditions (Bertness and Callaway 1994; He et al. 2013). In fact, we found that in some plant communities at the less productive extreme of the gradient (i.e., Planerón and Leciñena) the proportion of target species that aggregated diversity in their vicinity was higher than those that segregated it, in agreement with results from previous research (Chacón‐Labella et al. 2016; Foronda et al. 2019). However, we did not observe a general shift towards local diversity aggregation (i.e., towards a net positive interaction outcome) with increasing environmental severity. This kind of shift is frequently found in studies testing how the magnitude of interactions between a few pairs of species changes with the environmental conditions (Soliveres et al. 2015b; Liancourt et al. 2017). Unlike those studies, our diversity spatial pattern analysis allowed us to evaluate biotic interactions at community level, including interactions among all plant species in the community and giving a more realistic perspective on the prevalence of positive interactions under different conditions (He and Bertness 2014; Cavieres et al. 2014; Soliveres et al. 2015b; Liancourt et al. 2017; Liancourt and Dolezal 2020). At this level, a collapse in facilitation can be found at the least productive extreme of the gradient (Soliveres et al. 2015b; Liancourt et al. 2017), especially when communities along the gradient are composed by different species (Liancourt et al. 2017). However, we did not find this pattern either. One possible explanation is that our productivity gradient was not wide enough to include the most severe extreme, which would have prevented us from finding the expected shift (He et al. 2013; Liancourt et al. 2017). Nonetheless, our results emphasize that, although the relative importance of positive interactions increases with environmental severity, it is not always enough to overcome the negative ones at community level, while further studies are necessary to explain the exact causes.

Interestingly, in our study, the amount of target species with a neutral effect on diversity was notably high, especially in the less productive study sites (approx. 50% of target species; Supplementary Figure S3A). This finding contrasts with those reported by Chacón-Labella et al. (2016), who found that less than 20% of perennial plant species present in a semiarid Mediterranean plant community had a local neighborhood as diverse as expected. It could be argued that a lack of statistical power might cause scarce or rare species to appear as neutral despite having strong interspecific interactions (Wiegand et al. 2007; Rajala et al. 2019). However, in our study sites, the proportion of neutral species was not associated with species richness (Supplementary Figure S3B), assuming that richer plant communities have more rare species (Rajala et al. 2019). Furthermore, we found target species with very low relative abundance (around 1%) and a diversity pattern in their local neighborhood different than expected by the null model. On the other hand, some authors have argued that neutral or random biological processes govern the organization of plant communities, or that plant–plant interspecific interactions are weak (Hubbell 2001; Chacón‐Labella et al. 2017; Perry et al. 2017). In this line, other studies have also found that only a few species aggregate or segregate diversity in comparison with neutral species (Perry et al. 2017; Foronda et al. 2019). However, interspecific positive and negative interactions occur simultaneously (Holzapfel and Mahall 1999; Schöb et al. 2013; Poulos et al. 2014; Saiz et al. 2017). In this sense, it is likely that target species had a net positive effect on a set of beneficiary species, while being segregated from others, eventually resulting in a neutral effect on diversity (Wiegand et al. 2007; Arroyo et al. 2015; Foronda et al. 2019). Thus, taken as a whole, our results suggest that the observed shift from a segregated pattern of diversity to a local pattern as diverse as expected would respond to an increase in positive interactions in less productive environmental conditions, counteracting negative ones, and highlighting the importance of facilitation for diversity maintenance in these plant communities.

As expected, we found differences among plant life form groups in shaping the local pattern of diversity. Overall, graminoids acted as diversity segregators, creating less diverse neighborhoods than predicted by the null model (e.g., Brachypodium retusum, Lygeum spartum or Stipa lagascae). Graminoids are known to exert a strong belowground competition for resources, especially in water uptake, due to their dense root system in the upper soil layers (Jackson et al. 1996; Armas and Pugnaire 2011; Pueyo et al. 2016). In addition, species with a clonal growth form, as is commonly the case for grasses, are strong competitors for space (Pottier and Evette 2010; Benot et al. 2013). Furthermore, it has been documented how some grass species are facilitated in their early stages of development, to eventually replace the nurse plant when becoming adults (Armas and Pugnaire 2005). Thus, the high competitive capacity of graminoids prevents other species from being stablished in their surroundings, leading to the observed pattern of diversity segregation. On the contrary, we found that forbs generally had a more diverse local neighborhood than expected in less productive study sites (e.g., Asperula aristata, Polygala rupestris or Psoralea bituminosa). It is probable that plant species falling within this group were mostly facilitated species that appeared, together with annuals, in microhabitats with less severe environmental conditions (Holzapfel et al. 2006; Arroyo et al. 2015; Gonzalez and Ghermandi 2019). As these microhabitats in Mediterranean plant communities are usually highly diverse, those facilitated forb species were ultimately observed as diversity aggregators. Hence, they could be viewed as indicators of diversity islands in these plant communities.

Interestingly, the effect of shrub and dwarf shrub species on diversity patterning was intermediate. Shrubs are usually presented as good candidates to act as nurse species. In Mediterranean ecosystems, shrubs are relatively long-lived species adapted to local conditions with a dense canopy and deep roots that enhance microclimatic conditions for other species in the community (Callaway 2007; Pueyo et al. 2016). Accordingly, we found a number of shrub and dwarf shrub species that aggregated diversity in their local neighborhood, especially in plant communities at the less productive extreme of our gradient (Table 2). However, we also found many other shrubs and dwarf shrubs with neutral and negative effects on the local diversity pattern. On the one hand, the net interaction outcome between shrubs and other species in the community shifts to neutral (e.g., Teucrium charidemi) or negative (e.g., Plantago albicans, Ulex parviflorus or Artemisia herba-alba) in milder environmental conditions (Holzapfel et al. 2006), thus explaining their behavior as diversity segregators in more productive plant communities. On the second hand, some shrub and dwarf shrub species present key traits that can also explain discrepancies in the net interaction outcome within this group. For instance, previous studies have reported the presence of shrub and dwarf shrub species with a negative net effect on diversity in low-productive plant communities, due to their allelopathic activity (e.g., Artemisia herba-alba; Arroyo et al. 2015) or the ability to change soil salt concentration (Zhang et al. 2016). Moreover, legumes are often assumed to have a positive effect on neighboring plants and diversity due to their ability to increase nitrogen availability (Gómez-Aparicio et al. 2004; Wright et al. 2017). However, our results revealed that most woody legumes did not act as diversity aggregators (e.g., Anthyllis cytisoides, Astragalus nevadensis, Dorycnium pentaphyllum, Erinacea anthyllis, Genista sps., Hippocrepis emerus, Ononis sps., Ulex parviflorus, etc.), suggesting that only a small fraction of the fixed nitrogen ends enriching the surrounding soil (Ulm et al. 2017). Alternatively, it is likely that nutrient enrichment in the surroundings of legume species promoted competitive processes between facilitated species, reducing their local diversity (Schöb et al. 2013). Altogether, our results show the existence of key nurse species in shaping the local pattern of diversity in Mediterranean plant communities, and highlight the importance of a multi-species approach in order to clarify their general role at community level.

Contrary to our hypothesis, we found that the spatial heterogeneity in species composition among local assemblages (i.e., local beta-diversity) did not change significantly along the productivity gradient. Nevertheless, we still found that an increase in the nestedness component of beta-diversity was associated with less segregated diversity patterns in the local neighborhood of shrub species. This finding suggests that an increase in positive plant–plant interactions in less productive environmental conditions promotes spatial heterogeneity among local species assemblages due to differences in plant species richness, and points to key nurse species maintaining overall community diversity by creating diversity islands. This result agrees with previous studies on the role played by nurse species and facilitation as biodiversity drivers in severe environments (McIntire and Fajardo 2014; Le Bagousse-Pinguet et al. 2014; Cavieres et al. 2014; Gonzalez and Ghermandi 2019), supporting the fact that their protection should be a primary focus for management and conservation strategies in order to prevent the loss of diversity and functioning of these ecosystems (Valiente‐Banuet et al. 2015; Maestre et al. 2016). A high nestedness of diversity could imply to focus on the conservation of only the more diverse sites or patches (Socolar et al. 2016; Oksuz et al. 2020). However, our results revealed that species turnover was still higher than the nestedness of assemblages in all study sites. This means that nurse species were not such generalist facilitators, in line with the finding of few species with a net positive effect on diversity, and different target species harbored different sub-communities regarding species composition. In fact, patch-forming species with different functional traits modify their surroundings in very distinct ways, promoting dissimilarities in species assemblages (Macek et al. 2016; Saiz et al. 2018; Tormo et al. 2020). Hence, protecting multiple plant species in Mediterranean plant communities, and not only focusing on those more diverse patches, would be a more suitable strategy to preserve diversity in face of the expected consequences of climate change.

Our results can also be directly applied to restoration practices, as we offer a list of suitable species that can be used to enhance the success of diversity restoration projects performed in Mediterranean plant communities (Padilla and Pugnaire 2006; Gómez-Aparicio 2009; Hulvey et al. 2017). In this respect, the best candidates would be the abundant shrub species acting as nurses in aggregating diverse plant species assemblages (Table 2). For instance, Ononis tridentata, Salsola vermiculata and Gypsophila struthium subsp. hispanica would be good candidates in plant communities of the Middle Ebro Valley, Lavandula multifida, and Teucrium charidemi would be suitable in Cabo de Gata, while Arenaria armerina and Rosa sp. should be chosen for Sierra Nevada communities. Moreover, neutral shrub species might be also suitable candidates provided that they have a positive net effect on certain species of interest. This would be the case of Anthyllis cytisoides and Lycium intrincatum in Cabo de Gata or Berberis vulgaris subsp. hispanica and the cushion Hormathophylla spinosa in Sierra Nevada. On the contrary, we did not find any species that aggregated diversity in its local neighborhood in Sierra de Guara sites, with the exception of Quercus faginea in Ibort (Table 2). This finding might be explained by the fact that in Ipies and Rapun plant species establish predominantly in gaps where Q. faginea is absent, supporting its negative role on local diversity, while in Ibort, where Q. faginea is quite abundant, species in the community can appear within its local neighborhood too, explaining its positive role on diversity.

In conclusion, in our study we found that diversity segregation in the local neighborhood of plant species decreases in the less productive plant communities, according to an increase in the relative importance of positive interactions. However, we could not find a clear shift towards net diversity aggregation at the most severe environmental conditions. On the other hand, despite the general view that shrubs are good candidates to act as nurse species, our results emphasize that more studies taking a multi-species approach should be performed to clarify their general role at community scale. In fact, the number of species with a net positive effect on diversity was relatively low. Nonetheless, these nurse shrub species seem to play an essential role in low productive plant communities by creating and maintaining diversity islands. In addition, the dominance of species turnover over the nestedness component of beta-diversity suggests that each nurse could interact with different species, highlighting the presence of contrasted facilitated species assemblages. Overall, our results provide further insights into the effect of plant–plant interactions in shaping the local structure and composition of diversity in Mediterranean plant communities, and emphasize plant species and plant life form groups that should be the focus of management and conservation strategies in order to prevent diversity loss.