Frontiers of Biogeography Co-occurrence frequency in vegetation patches decreases towards the harsh edge along an arid volcanic elevational gradient

co-occurrences helps to improve general hypotheses about biotic interactions. Abstract Positive plant–plant interactions are thought to drive vegetation patterns in harsh environments, such as semi-arid areas. According to the stress-gradient hypothesis (SGH), the role of positive interactions between species ( facilitation ) is expected to increase with harshness, predicting associated variation in species composition along environmental gradients. However, the relation between stress and facilitation along environmental gradients is debated. Furthermore, differentiating facilitative interactions from other underlying mechanisms, such as microtopographic heterogeneity, is not trivial. We analysed the spatial co-occurrence relationships of vascular plant species that form patchy vegetation in arid lapilli fields (tephra) from recent volcanic eruptions on La Palma, Canary Islands. We assume a harshness gradient negatively correlated with elevation because of more arid conditions at lower elevations where water availability is considered the most limiting resource. Based on the SGH we expect a greater degree of co-occurrence at lower elevations, as an outcome of facilitation is plants co-occurring in the same patch. We tested this at both the species and the individual plant level. We analysed the species composition of 1277 shrubby vegetation patches at 64 different sampling points, ranging from the coast to around 700 m a.s.l. Patch morphology and microtopographic heterogeneity variables were also measured, to account for their potential effects on the species composition of patches. We used generalized linear models and generalized mixed-effects models to analyse species richness, number of individuals in patches and percentage of patches with positive co-occurrences, and a pairwise co-occurrence analysis combined with a graphical network analysis to reveal positive links between 13 of the species. We found that the percentage of patches with positive co-occurrences increased at higher elevations, in contrast to the predictions of the SGH, but in accordance with a refined stress-gradient hypothesis for arid sites, generalized explaining We according

• We investigate an elevational gradient on La Palma, Canary Islands, where harshness inversely correlates with elevation because of aridity, using a novel approach combining a co-occurrence and network analysis.
• We find a decrease of positive co-occurrences with increasing harshness, and different patterns for species richness and abundance.
• Co-occurrence patterns and their drivers appear to differ depending on whether temperature or aridity determine the underlying harshness gradient.
• Addressing the effects of harshness on co-occurrences helps to improve general hypotheses about biotic interactions.

Abstract
Positive plant-plant interactions are thought to drive vegetation patterns in harsh environments, such as semi-arid areas. According to the stress-gradient hypothesis (SGH), the role of positive interactions between species (facilitation) is expected to increase with harshness, predicting associated variation in species composition along environmental gradients. However, the relation between stress and facilitation along environmental gradients is debated. Furthermore, differentiating facilitative interactions from other underlying mechanisms, such as microtopographic heterogeneity, is not trivial. We analysed the spatial cooccurrence relationships of vascular plant species that form patchy vegetation in arid lapilli fields (tephra) from recent volcanic eruptions on La Palma, Canary Islands. We assume a harshness gradient negatively correlated with elevation because of more arid conditions at lower elevations where water availability is considered the most limiting resource. Based on the SGH we expect a greater degree of co-occurrence at lower elevations, as an outcome of facilitation is plants co-occurring in the same patch. We tested this at both the species and the individual plant level. We analysed the species composition of 1277 shrubby vegetation patches at 64 different sampling points, ranging from the coast to around 700 m a.s.l. Patch morphology and microtopographic heterogeneity variables were also measured, to account for their potential effects on the species composition of patches. We used generalized linear models and generalized mixed-effects models to analyse species richness, number of individuals in patches and percentage of patches with positive co-occurrences, and a pairwise co-occurrence analysis combined with a graphical network analysis to reveal positive links between 13 of the species. We found that the percentage of patches with positive co-occurrences increased at higher elevations, in contrast to the predictions of the SGH, but in accordance with a refined stress-gradient hypothesis for arid sites, in which characteristics of the interacting species are incorporated.

Introduction
Plant-plant interactions, whether positive or negative or reflecting complex dynamics, directly affect the distribution and spatial organization of species at local and regional scales (McIntire and Fajardo 2014). Negative plant-plant interactions (competition) have for a long time been thought to act as a biotic filter (White and Jentsch 2004) and thereby an important factor shaping plant communities (Brooker et al. 2008). However, numerous studies in the last few decades have provided evidence that positive interactions (facilitation) play a crucial role as well (Bertness and Callaway 1994, Bruno et al. 2003, Michalet 2006. Facilitation describes positive interactions between species within a single trophic level, for which proposed underlying mechanisms are diverse. Facilitative interactions between plants can be direct (e.g. nurse plants enabling the survival of a seedling that would otherwise die) or indirect (e.g. enhanced reproductive fitness through pollinator attraction by neighbouring plants) (Brooker et al. 2008, Tielbörger andKadmon 2000), and change over time. It is now widely considered that positive interactions often contribute substantially to species distributions, and enhance local biodiversity, especially on sites where species have to cope with unfavourable conditions for plant growth and meet their environmental limits Badano 2009, McIntire andFajardo 2014).
Stressful conditions usually prevail at the extremes of environmental gradients due to resource restrictions. Attempts to understand how the balance of facilitative and competitive plant-plant interactions change in increasingly harsh environments have led to the development of different conceptual models. Bertness and Callaway (1994) proposed the stress-gradient hypothesis (SGH), according to which facilitative interactions among plant species increase when environmental conditions get harsher or disturbance increases, while the role of competition between plants plays a more important role at benign sites. The SGH has received empirical support from a number of studies along different stress gradients, such as salinity (He et al. 2011, Zhang andTielbörger 2019) and especially along elevational gradients in cold alpine environments (Callaway et al. 2002. However, the SGH is still debated, with several recent studies reporting contradictory findings, such as a decreasing (Bonanomi et al. 2016) or hump-shaped relationship (Maestre and Cortina 2004) between stress and facilitation, and even a lack of any facilitative interaction  or a shift to enhanced competition Kadmon 2000, Armas andPugnaire 2005) at the most extreme sites. Because most studies did not find enhanced facilitation in extremely dry areas (e.g., Tielbörger and Kadmon 2000, Berdugo et al. 2018, Filazolla et al. 2020), a refined SGH was proposed for aridity gradients , Maestre et al. 2009). This revised hypothesis incorporates the life histories of the interacting species, as well as the respective stress factor in terms of environmental conditions and resource availability. Indeed, a collapse of positive plant-plant interactions towards extremely arid sites may be common in ecosystems where water availability is the limiting factor.
Plant patchiness is a characteristic pattern for waterlimited sites (Lejeune andTlidi 2002, Rietkerk et al. 2004). The drivers of species' aggregation within patches can include environmental heterogeneity, seed trapping, vegetative propagation, dispersal limitation and true biotic interactions (López et al. 2016). Disentangling these reasons for species co-occurrences is not trivial because microsite particularities that favour the establishment of one or more individual plants (or at least increase survival) will gain enormously in importance at harsh sites compared with more benign sites (Steinbauer et al. 2016). As a result, the effects of both microtopographic heterogeneity and positive biotic interactions may increase in a similar way. In other words, patchy vegetation cover might reflect patchy environmental suitability, facilitative biotic interactions or a combination of both. Additionally, species tend to interact in complex networks rather than pairwise combinations, which results in indirect effects among plants through diffuse or intransitive interactions (Losapio et al. 2019). Furthermore, patterns of biotic interactions may change over time, with a changing hierarchy of resource ratios, changing stoichiometric requirements of individual species and thus of limiting factors, as suggested by the theory of pulse dynamics and disturbance in ecology (Jentsch and White 2019). However, there is increasing interest in inferring plant-plant interactions from co-occurrence field data and promising statistical methods to analyse complex co-occurrence data have been proposed. These include pairwise interaction matrices (Gotelli 2000, Veech 2013), uni-and bipartite network approaches (Morueta-Holme et al. 2016, Montesinos-Navarro et al. 2018, Losapio et al. 2019, spatial pattern analysis (Wiegand et al. 2013, López et al. 2016) and combinations of these (e.g., Fournier et al. 2016), which can be applied to co-occurrence data derived from field observations. Plants of the lapilli fields in the south of the volcanic Canary Island La Palma spatially organise in a characteristic mosaic of vegetation patches and bare (tephra) soil. Here, we study an elevational gradient on these lapilli fields ranging from the coast towards the higher-elevation San Antonio volcano. This elevational study gradient represents an inverse harshness gradient, from most arid (harshest) near the coast to more mesic with increasing elevation Beierkuhnlein 2011, Irl et al. 2020). Harshness towards the coastal edge of the gradient may also be increased through increased salinity from salt spray exposure. We therefore expect patches at lower elevations to reflect the greater harshness by smaller patch size and sparser patch organisation. According to the SGH, positive plant-plant interactions should increase in preponderance at harsher sites. Here, several species might act as facilitators (nurse plants), which should be reflected in high numbers of positive links with other plants of the same or different species. From this Frontiers of Biogeography 2021, 13.3, e49743 © the authors, CC-BY 4.0 license 3 hypothesis we therefore expect a greater degree of species co-occurrence in patches at lower elevations, which should result in increased a) species richness, b) number of individuals and c) percentage of patches with positive co-occurrences.

Study area
La Palma is one of the seven main islands of the Canary archipelago in the Atlantic Ocean. The island originates from eruptions of several volcanoes. Our study was conducted in the southern part of the island in the Fuencaliente municipality ( Fig. 1), where the two most recent eruptions of the San Antonio volcano (1677) and the Teneguía volcano (1971) formed an arid, black landscape consisting of lava and lapilli fields (Klügel et al. 1999). Lapilli describes pyroclastic fragments (tephra) with grain sizes ranging from 2 to 64 mm, which are emitted during mafic and silicic eruptions (Troll et al. 2017). The usually high porosity of interconnected voids and a lack of reactive silica characterize lapilli fields as nutrient poor and dry environments with low water storage capacity (Lomoschitz et al. 2006). Sampling only on lapilli ensures consistency of substrate type, controlling substrate variation. La Palma is characterized in general by a subtropical climate with humid winters and dry summers. However, precipitation varies locally with high amounts on the north-eastern side and lower amounts on the western and especially on the southern part of the island, and semi-arid coastal areas. Mean annual precipitation ranges from 1,500 mm at the north-eastern part to 200 mm at the south-western part of the island (Irl et al. 2020). Fog drip can supplement lower amounts of precipitation during the drier summer months on a local scale, depending on the topography and vegetation type (Irl et al. 2015). Mean annual temperature decreases from around 22 °C at the coast to around 9 °C at highelevation summits (Irl et al. 2015). On the lapilli fields of the Fuencaliente municipality, the arid conditions of both climate and soil are associated with very patchy vegetation ( Fig. 1).

Sampling Design and Methods
In March 2018, we sampled 1277 vegetation patches at 64 different sampling sites along an elevational gradient of approximately 650 metres, ranging from the coast (17 m a.s.l.) to the top of the San Antonio volcano (675 m a.s.l.) thereby collecting information on a total of 31 plant species. The 64 sampling sites were chosen randomly along the elevational gradient at intervals of roughly 15 m elevation. We avoided sites in close proximity to roads and buildings as well as sites that showed any signs of human disturbance or intensive land use. We recorded elevation (m a.s.l.), slope (%) and aspect (°) at the centre of each site. Starting from the central point, we marked the 20 closest vegetation patches and measured the distance to the closest as well as the most remote patch. Patches with vegetation entirely < 0.10 m in height were excluded from sampling, to avoid recording seedlings and annual species, as we expect high seedling mortality rates and annual species to be more prone to react to short-term weather fluctuations than to longer-term biotic interactions (Morris et al. 2008). The sampled patches met our criteria for a patch according to the definition by Tongway and Hindley (2004), which defines a patch as a single plant or an association of plants, connected by living plant matter overlapping each other physically above ground or at least via a litter bridge. For each vegetation patch, we recorded species richness, identity and abundance (number of individuals) of perennial species. We used the flora of the Canary Islands (Muer et al. 2016) for species identification and taxonomic nomenclature. We also measured the patch height, width and length (all in m), estimated the percentage of dead wood and litter, the orientation of the longest patch side and the distance (m) to the nearest neighbouring patch. To account for potential micro-topographic specificities and micro-availability of nutrients, we recorded the presence or absence of rocks, hills, depressions and rabbit excrement for each patch.
For each patch, we calculated the patch volume (m 3 ) by its measured height, length and width, respectively. For this, we used the formula 2 * pi *a * b *c 3 , where a, b and c are patch height, length and width, respectively. The formula describes a dome-shaped semi-ellipsoid that we considered the most proximate geometric form for the patch biomass above ground. We calculated a density proxy for each of the 64 sites by dividing the distance from the site centre to the closest and the most remote patch at each site. Density values closer to 1 therefore indicate patches in close proximity to each other (while still being small patches rather than continuous vegetation), while values close to 0 represent sites with very sparse patch organization. The amount of mean annual precipitation (mm) was extracted for each elevational site (n=64) from the interpolated precipitation dataset published by Irl et al. (2020).
To approximate normality of residuals, the variables patch volume, site density, amount of dead wood and litter were log transformed (natural logarithm), and the variables slope and distance to nearest neighbour were square-root transformed.

Statistical analysis
All statistical analyses were performed using R version 4.0.2 (R Core Team, 2020). The threshold for determining significance was alpha = 0.05 for all analyses.
We used single-predictor generalized linear models with a Gaussian error distribution to test how mean annual precipitation, patch volume, the distance to the nearest neighbour as well as the site density proxy changed with elevation. We used single-predictor and multi-predictor generalized linear models with a Poisson error distribution (natural logarithm link function) to analyse changes in species richness and number of individuals within the patches (n patches = 1277, n species = 31) along the elevational gradient. Single-predictor models were calculated for each explanatory variable, respectively. For the multi-predictor models, environmental parameters were tested for collinearity by using Spearman's rank correlation. For all pairs of environmental parameters with correlation coefficient > 0.6, one was excluded from further analysis based on strength of correlation with the response variable (species richness or number of individuals in the patches) in the single-predictor models. The remaining variables were scaled to yield estimates in standard deviation without centering by the scale() function and then included in a first multi-predictor model. We excluded non-significant parameters by manual stepwise model simplification, and compared models using the Akaike information criterion (AIC) and pseudo-R 2 .
Single-predictor and multi-predictor generalized linear models were performed using the glm() function in R. Pseudo-R 2 (pseudo-r 2 in the case of single-predictor models) was calculated using the Nagelkerke method, with the PseudoR2() function in the "DescTools" package (Signorell et al. 2020). To quantify the relative contribution of the final explanatory variables, we performed hierarchical partitioning with the hier. part() function in the "hier.part" package (Mac Nally and Walsh 2004), with root-mean-square prediction error as the goodness-of-fit measure for Poisson-error glms. Categorical microsite specificities (Rocks, Hills, Depressions, Rabbit excrements, None, Several) were included as random effects in a subsequent generalized mixed-effects model with Poisson error distribution, to account for microsite heterogeneity along the elevational gradient. We performed the generalized mixed-effects models with the glmer() function in the "lme4" package (Bates et al. 2015) and used AIC to compare this model with the final multi-predictor generalized linear model, as Nagelkerke's pseudo-R 2 is not available for generalized mixed-effects models.
We analysed the species associations of the multispecies patches by calculating a pairwise probabilistic co-occurrence model, using the approach of Veech et al. (2013). We implemented this in the "cooccur" package in R (Griffith et al.2016) in combination with the R packages "vegan" (Oksanen et al. 2019) and "MASS" (Venables and Ripley 2002). To reduce the impact of rare species, we excluded all species occurring in fewer than 5 of the 1277 patches from the co-occurrence and network analysis. Number of individuals within patches along the study gradient is best explained by a weak U-shaped relationship. C) Species richness increases with increasing number of individuals in patches (please note that the number of individuals in patch is log-transformed by natural logarithm transformation). D) Mean annual precipitation significantly increases with elevation, representing an inverse harshness gradient. E) The number of individuals in a patch significantly increases with increasing patch volume (please note that the x-axis in E is log-transformed by natural logarithm transformation). Pseudo-r 2 by Nagelkerke. In plots A, B, C and E, n=1277 patches. Please note that explained variation is low for models shown in A-C.
By using the cooccur() function ("cooccur" package), we calculated a pairwise species co-occurrence matrix. This matrix was then compared to a probability matrix, containing the expected co-occurrences of the species, based on the probability mass function of a hypergeometric distribution (Griffith et al. 2016). A final output matrix was produced with positive (1), neutral (0) or negative (-1) values for each species pair, depending on whether the observed co-occurrences are higher, equal or lower than the frequencies expected from the probability matrix. In the cooccur() function, we chose to activate a threshold that removed species with less than one expected co-occurrence, and thus reduced the amount of potential species pairs. A truerandom classifier in the function has to be defined to statistically distinguish truly random from unclassifiable associations, which we set to a stricter value (0.05) than the default (0.1). This means that the deviance from the expected co-occurrences should not be more than 0.05 times the total number of included sites to classify a species association as truly random. As we were mainly interested in positive plant-plant associations, we subsequently excluded all random associations to build an undirected species co-occurrence network, where random or negative links between species are not shown. We assessed the modularity of our final network graph by calculating the leading non-negative eigenvector of the community matrix according to the method developed by Newman (2006) for undirected unipartite networks. The network was built and its modularity analyzed using the "igraph" package (Csardi & Nepusz 2006); illustration was done in combination with the free software "Inkscape". Based on the results of the co-occurrence and network analyses, we identified all patches that contained at least one positive co-occurrence pair. We calculated the percentage of patches with positive co-occurrences (number of patches with positive cooccurrences/ number of patches at the site), as well as the percentage of multi-species and multi-individuals patches at each site. We used generalized linear models with binomial errors to test these variables against elevation.

Co-occurrence pattern along the elevational gradient
Within the 1277 patches among 64 different sampling sites, we recorded 31 plant species in total. Overall, 878 patches consisted of only one species (single-species patches) and 399 patches consisted of more than one species (multi-species patches). Regarding the number of individuals (irrespective of species), 613 patches consisted of only one individual, while 664 patches consisted of more than one individual.
Species richness per patch slightly increased with increasing elevation, our inverse proxy for harshness (pseudo-r 2 = 0.034; Fig. 2A; Table 1), while the number of individuals per patch showed a weak U-shaped relationship with elevation (pseudo-r 2 = 0.05) (Fig. 2B, Table 1). Species richness significantly increased with increasing log-transformed number of individuals per patch (Fig. 2C, Table 1). In the single-predictor Frontiers of Biogeography 2021, 13.3, e49743 © the authors, CC-BY 4.0 license 6 generalized linear models, the number of individuals per patch was best explained by the patch volume (pseudo-r 2 = 0.84) with more individuals in bigger patches (Fig. 2E, Table 1). Mean annual precipitation significantly increased with elevation (r 2 = 0.967) (Fig. 2D, Table 2). Mean annual temperature decreased significantly with elevation (see Fig. S1, Supplementary material). The final multi-predictor model for species richness per patch (pseudo-R 2 = 0.187, AIC = 3220.4) included five explanatory variables (Table 3). Adding microsite specificities as a random factor slightly improved the model (AIC = 3161.1). Hierarchical partitioning revealed that the relative contribution of patch volume and elevation was highest (Fig. S2, Supplementary  material). The final multi-predictor model for the number of individuals per patch (pseudo-R 2 = 0.914, AIC= 5116.6) included seven explanatory variables and was improved by the addition of microsite specificities as random effects (AIC= 4992.5). Patch volume had the highest contribution to the total explained deviance (Fig. S2, Supplementary material).

Co-occurrence network analysis of main species
Our co-occurrence matrix analysis included 1277 patches and 15 species, and led to 105 potential species pairs, of which 28 pairs were removed as the expected co-occurrence was less than one. The analysis found 45 of the remaining co-occurrences to be non-random, of which 22 were classed as positive co-occurrences and 23 as negative. All 15 species were included in the network produced from these data (Fig. 3A). Due to our focus on positive species associations, we chose an undirected network, where links between species represent positive co-occurrences, while random or negative links are not shown (Fig. 3A). For 13 out of 15 species we detected at least one positive co-occurrence with other plant species. Kleinia neriifolia (Asteraceae) showed the most co-occurrences with other species (n partners = 9), followed by Micromeria herpyllomorpha (Lamiaceae, n partners = 7) and Todaroa aurea (Apiaceae, n partners = 4) (Fig. 3, Table S1 Supplementary material). The final network had a relatively low modularity (0.24), with five modules. The two species without any positive co-occurrence with other plant species (Astydamia latifolia and Bystropogon origanifolius) constituted single-species modules.
The percentage of multi-species patches and multi-individual patches did not change significantly with elevation (Table 4). The percentage of patches containing positive co-occurrences significantly increased with increasing elevation (pseudo-r 2 = 0.183) Table 1. Model statistics of best fit single-predictor generalized linear models explaining species richness (patch richness) and number of individuals (patch individuals) within patches (log = natural log-transformation, sqrt = square root transformation). Pseudo-r 2 according to Nagelkerke. For all models, sample size was 1277 and the Poisson error distribution was used. For most models, high sample size led to statistical significance, while explained variance (pseudo-r 2 ) was very low.  (Fig. 4, Table 4). The percentage of patches containing negative co-occurrences did not change significantly along the elevational gradient (see Fig. S3, Supplementary material).

Discussion
In our study site along an elevational harshness gradient in an arid volcanic ecosystem, we found plant species growing in a patchy vegetation, as is characteristic for dry sites. Such arid and semi-arid sites comprise characteristic vegetation patterns with scarce biotic islands within a matrix of bare substrate; in the islands, perennial plants ameliorate the microconditions and enable plant growth (Pugnaire et al. 2011, Verwijmeren et al. 2013. Despite the significant increase in precipitation with elevation, we did not find the expected pattern of smaller and more sparsely growing patches toward the harsh edge at the coast. Against our expectation from the stress-gradient hypothesis (SGH), we observed no increase in species richness within those patches or percentage of positive Table 2. Model statistics of best fit single-predictor generalized linear models explaining site parameters (distance to nearest neighbor, patch volume and site density) as well as mean annual precipitation along the study gradient. (log = natural log-transformation, sqrt = square root transformation). All models used Gaussian error distribution. For the first three models, high sample size led to statistical significance, while explained variance (r 2 ) was very low.  Table 3. Model statistics of best fit multi-predictor generalized linear models explaining species richness within patches (patch richness) and number of individuals (patch individuals) in patches (log = natural log-transformation, sqrt = square root transformation). Pseudo-R 2 according to Nagelkerke. All models used Poisson error distribution and had a sample size of n = 1277. The undermost row (1|microsite specificity) shows the model statistic of a generalized mixed-effects model, where microsite specificities (Rocks, Hills, Depressions, Rabbit excrements, None, Several) were added to the final generalized linear model. We compared the models according to their AIC values. dependent variable model pseudo-R 2 AIC p-value patch richness~elevation 0.187 4220.4 0.000 *** +slope (sqrt) 0.000 *** +patch volume (log) 0.000 *** +distance to nearest neighbour (sqrt) 0.000 *** +litter (log) 0.000 *** +(1|microsite specificity) -3161.1 patch individuals~elevation 0.914 5116.6 0.000 *** +slope (sqrt) 0.000 *** +patch volume (log) 0.000 *** +orientation of longest patch side 0.000 *** +distance to nearest neighbour (sqrt) 0.000 *** +dead wood (log) 0.000 *** +litter (log) 0.000 *** +(1|microsite specificity) -4992.5 Table 4. Model statistics of best fit single-predictor generalized linear models explaining the percentage of multi-species patches, multi-individual patches and patches containing positive co-occurrences against elevation. Pseudo-r 2 according to Nagelkerke. For all models, sample size was 64 and a binomial error distribution was used. co-occurrences towards the harsh edge of the gradient at low elevation. In contrast, there was a significant, but weak increase of species richness per patch and percentage of patches with positive co-occurrences with increasing elevation (i.e. towards the more benign gradient section). This contradicts the findings of several studies along elevational gradients, which reported a notable shift from negative to positive plant-plant interactions with increasing harshness, as predicted by the original stress-gradient hypothesis of Bertness and Callaway (1994). However, most of these studies were conducted in high-elevation (Choler et al. 2001, Callaway et al. 2002, Cavieres et al. 2002 or high-latitude ecosystems (Dorrmann and , Gavini et al. 2019, where changes in elevation primary reflect a temperature gradient (Körner 2003), whereas the elevational gradient in this study is assumed to reflect a water limitation gradient. Facilitative plant interactions seem to play an important role when environmental conditions get colder, enhance local biodiversity and might even represent a future protection mechanism for alpine ecosystems (Cavieres et al. 2016). As we do not reach comparable elevations, and precipitation amount increases with elevation while evaporative demand decreases (Irl and Beierkuhnlein 2011), our sampling gradient reflects an aridity-stress rather than a temperature-stress gradient, being most stressful at the lowest elevations. Our findings indicate that different abiotic stress factors act when not sampling along a classic high-elevational gradient, resulting in a very different pattern of plant-plant interactions. The fact that the nature of the stress factor might affect the outcome of studies analysing biotic interaction was recognized by previous studies. Maestre et al. (2009) proposed a refinement of the original SGH, recommending the incorporation of species life history and stress factor characteristics. Our results agree with those of other studies at arid sites, where a peak of facilitative interactions between plants at most extreme aridity seems to be an exception and a reverse shift from positive to negative or neutral interactions was frequently observed (Berdugo et al. 2019, Lucero et al. 2020. Species-specificity is probably one explanation for this collapse of facilitative interactions at the most extreme sites (Soliveres et al. 2012, Filazolla et al. 2020) and also likely one of the reasons at our coastal gradient edge (see Fig. 3B). The interpretation of co-occurrence patterns is further complicated by potential additional or overlapping stress gradients and species turnover. Both factors might play a role when sampling across environmental gradients or vegetation communities (Gallien et al. 2018). An overlapping stress gradient imposed by enhanced salinity probably applies to the low elevations of our study gradient at the coast, where we did not detect any patches with positive cooccurrences (Fig. 4). Here, patches were mostly singleindividual patches of Astydamia latifolia, a compact, salt-tolerating, succulent shrub that is typically found in rocky habitats along the Mediterranean coast (Pérez-Alonso et al. 1999). This species is evidently well adapted to the arid and saline conditions in coastal areas. It showed no significant positive co-occurrence link to other species of the study area, which indicates that the effect of adaptive traits might play a more important role compared to plant-plant interactions to establish and survive in extremely dry habitats such as these (Soliveres et al. 2012). Hence, species identity and specificity should be considered when interpreting co-occurrence patterns.
Patch morphology had an additional impact on species co-occurrences along the gradient. Small and compact patches at the coast probably do not offer similarly favourable conditions to establish as do the larger and less dense canopies of patches at higher elevations. This could have notable impact on the potential to act as a facilitator or nurse plant, as was shown for cushion plants with varying morphology along an elevational gradient (Bonanomi et al. 2016). The number of individuals per patch did not increase towards the harsh gradient edge, but significantly increased with elevation and patch volume. As the species richness per patch did not increase in a similar proportion, we assume different drivers for these cooccurrence patterns. Plant aggregation in a patchy mosaic can be driven by biotic interactions, but as patch volume explained a large part of the number of individuals in our sampled patches, seed trapping and seedling survival in the shadow of adults might be the major driver. In sparsely vegetated environments, patches probably act as key accumulation objects for wind-dispersed seeds (López et al. 2016). This explanation is supported by the fact that including the percentage of litter and dead wood as well as microsite heterogeneity improved the explanatory power of the final model predicting the number of individuals per patch. Hence, patch morphology, and microsite heterogeneity and sun shelter by tall adults probably facilitate the establishment and survival of plants within the patch, rather than speciesspecific facilitation mechanisms, as proposed by Steinbauer et al. (2016). The limited dispersal ranges of island species (Carlquist 1974) might furthermore contribute to seed dispersal within or close to the parental patch. Additionally, woody shrubs in arid sites often disperse by vegetative propagation (Bond and Midgley 2001). We determined the number of individuals within a patch by counting the above ground distinguishable present plants. Regarding the strong slope of most sites and the easily erodible lapilli, true differentiation of individuals is therefore only possible by genetic analyses or (destructive) sampling of the rooting network. Our results emphasize the importance of including information on morphology, microsite topography and microsite soil water availability as well as testing co-occurrence patterns on different levels (species richness, number of individuals and relative contribution of positive interactions) to disentangle potential drivers of species aggregation in patches.
Interactions between species are complex, and differentiation becomes more difficult the more species a study area contains (Wisz et al. 2013). As a consequence, a lot of surveys focus on predefined target species or few species pairs and neglect further explanatory variables (Steinbauer et al. 2016). However, facilitative interactions are frequently concluded based on species co-occurrences. We wanted to test whether positive species interactions could be identified by co-occurrence data without previous knowledge of potential facilitator species. Across all 31 species along the studied gradient, we found a network of 13 species with at least one positive co-occurrence link to another species. Most of the species co-occurring with others are common shrubs that are endemic to La Palma or the Canary Islands. The species most often occurring with other species was Kleinia neriifolia, a shrub from the Asteraceae family with archipelago-wide distribution in the succulent scrub vegetation of the Canary Islands (del Arco-Aguilar et al. 2010); it rarely grew alone. This indicates that the patchy vegetation pattern here could be at least partly explained by biotic interactions and that several species might either act as facilitators or benefactors. As it is difficult to infer true facilitation exclusively by observed co-occurrence data, we here suggest a potential interaction network without defining the mechanism, the direction or the intensity of the respective plant-plant interactions. Moreover, quantifying facilitative interactions is complicated by the complexity of potential interactions in vegetation communities. It is difficult enough to measure pairwise interactions, but understanding their direction and effects becomes even more challenging with increasing species richness. More species result in more complex networks, where indirect or intransitive positive or negative effects of different species on each other might mask the underlying biotic interactions (Losapio et al. 2019). Analyzing co-occurrences at a community level, as we have done, can serve as the first step towards identifying potential biotic drivers and relevant plant traits behind vegetation patterns of co-occurrence, especially in a study system where knowledge of underlying interactions is lacking, as are a priori hypotheses on which species may be facilitator or benefactor species.

Conclusions
We found that the percentage of positive cooccurrences slightly increased along an elevational gradient in vegetation patches in an arid volcanic environment and thus, did not increase towards the harsh edge of the sampled elevational gradient. This is consistent with the refined stress gradient hypothesis, where positive plant-plant interactions are predicted to decrease at most extreme sites, which may be most applicable to aridity as a stressor. The number of individuals per patch revealed a different pattern and increased with increasing elevation towards the more benign end of the gradient. Number of individuals was mainly correlated with patch morphology and microsite heterogeneity. Studying co-occurrences at the community scale can help to reveal potential interaction networks, but finds its limit when it comes to quantifying the intensity, the direction and the mechanism of facilitation.

Acknowledgements
CB was supported by the European Horizon 2020 project 64176 ECOPOTENTIAL: Improving future ecosystem benefits through Earth Observations. We kindly thank the Caldera de Taburiente National Park Directorate and Felix Medina from the Consejería de Medio Ambiente, Cabildo de La Palma, for field work permissions and their expertise on the study sites and flora. We thank the technical staff of the Biogeography department for their help during field work preparation. We thank all participants of the La Palma Science School 2018 for fruitful discussions, Dagmar Hanz for her valuable review of the first manuscript version and support during the statistical analysis and David Kienle for his support regarding the statistical analysis and the figure illustration. We thank Robert Whittaker, Lea de Nascimento and one anonymous reviewer whose insightful comments and suggestions have substantially improved the manuscript.

Author Contributions
PE, JE, TK, CB, AJ and SI designed the concept of this study. PE, JE, TK conducted the field work and sampled the data with help and input from SI, CB, AC, OV and RF. PE analysed the data with contributions from JE. PE wrote the manuscript with contributions from all authors.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Table S1. Species of the 15 main species of the cooccurrence analysis. Figure S1. Mean annual temperature along the elevational study gradient. Figure S2. Hierarchical partitioning of the final multipredictor glms. Figure S3. Patches containing negative co-occurrences along the study gradient.