Fine- and broad-scale distribution of a cushion plant species: Patterns and predictors for Euphorbia clavarioides

ABSTRACT

Cushion plants are a key Arctic, Antarctic, and alpine growth form, with many cushion-forming species strongly affecting community structure in abiotically stressful environments. Despite their ecological importance, there is little information about what drives the distribution of species exhibiting this growth form. This study investigates the determinants of the distribution of a cushion plant species, Euphorbia clavarioides at (1) a fine scale, using field-collected predictors from an alpine landscape; and (2) a broad scale, using distribution records and climate data across the species’ distributional range. At the fine scale, the species was locally rare (occurring in about 4 percent of samples) and may be limited to specific microsites by interspecific competition with taller-growing species. Broad-scale species distribution modeling showed that both temperature and rainfall are important in predicting the distribution of E. clavarioides with a higher probability of occurrence in areas with higher annual precipitation and mean annual temperatures < 15°C. Given the species’ sensitivity to competition and abiotic conditions (high temperatures and low precipitation), E. clavarioides may be vulnerable to environmental changes. Therefore, assuming that other cushion plant species exhibit similar patterns, species of this growth form may be particularly useful indicators of change in alpine, Arctic, and Antarctic areas.

KEYWORDS:

• Cushion growth form
• geographic range
• spatial scales
• species distribution modeling

Introduction

The distribution of species, growth forms, biomes, and other biological units reflect their ecology, physiology, and evolutionary history (Rapoport 1982). For many plant growth forms, the factors determining their local and global distributions are well known. For example, the distribution of woody species is limited at both fine and broad scales by low temperatures, precipitation extremes, and the occurrence of disturbances (Körner 2007; Woodward, Lomas, and Kelly 2004). In contrast, annual plants are limited by dense vegetation and require sites where conditions are favorable for germination (Freas and Kemp 1983; Maron and Crone 2006; Turnbull et al. 2004). However, the factors driving the distribution of several other growth forms are still largely undocumented.

Cushion-forming plants are a key alpine, Arctic, and sub-Antarctic growth form (Aubert et al. 2014; Körner 2003), but the drivers of cushion plant species’ distributions are poorly understood. The cushion growth form is defined by compact, prostrate, herbaceous, perennial plants that form a smooth hemispherical shape (Billings and Mooney 1968; Cabrera, Rada, and Cavieres 1998; Ramsay 1992). Cushion plants strongly modify the local abiotic environment as their dome-shaped growth form and compact canopy slows heat loss and traps organic matter, which results in increased water-holding capacity (Aubert et al. 2014; Kleier and Rundel 2004; Ralph 1978; Ramsay 1992; Sklenář 2007), further enhancing cushion plants’ buffering of temperatures due to water’s high thermal inertia. Due to the impacts that cushion plants have on the abiotic environment, they often increase the local abundance, richness, and diversity of vascular plant and arthropod species (Badano et al. 2006; Barendse and Chown 2001; Cavieres et al. 2002; Molenda, Reid, and Lortie 2012; Sklenář 2009); they may enhance phylogenetic uniqueness of communities (Butterfield et al. 2013), and alter community structure by alleviating stress and providing resources for beneficiary species (Antonsson, Björk, and Molau 2009; Arredondo-Nuñez, Badano, and Bustamante 2009; Arroyo et al. 2003; Brancaleoni, Strelin, and Gerdol 2003). Indeed, in some systems they are considered keystone species (e.g. Lortie and Reid 2012; Phiri, McGeoch, and Chown 2009).

However, despite the importance of cushion plants for community structure in high-elevation and high-latitude systems, studies have not explicitly documented the environmental factors determining the distribution of plant species with cushion growth forms; reports on the drivers of cushion plant distributions have mainly been anecdotal (Table 1). The global distribution of cushion-forming species was recently reviewed (Aubert et al. 2014; see also Hauri and Schröter 1914), identifying regions with the highest richness of cushion plant species (e.g. South America and Europe; see also Boucher et al. 2016). However, this study did not identify the environmental drivers of these patterns. In a fairly coarse analysis, Boucher et al. (2016) showed that cushion richness within large political land units was correlated positively with low temperature and aridity (see also Boucher et al. 2012), while at finer scales, occurrence patterns for some species have been reported and some generalities emerge (Table 1); but the influence of almost none of the purported drivers have been statistically verified (although see e.g. Kleier and Rundel 2004). These anecdotal observations suggest that cushion plants occur in harsh abiotic environments (both locally and globally), but quantitative analyses of the distribution of cushion plants across any environmental gradients are still lacking.

Table 1. Summary of published information on the determinants of the fine-scale distribution of cushion plant species. Most of these data are observational, and only studies indicated with an asterisk statistically quantified the purported patterns

Species	Description	References
Androsace alpine (L.) Lam.	Confined to sparse vegetation	Schönswetter, Tribsch, and Niklfeld 2003
Androsace tapete Maxim.	Reaches maximum coverage at a specific suite of optimal climatic conditions (growing season temperatures 2.6–11.1°C, minimum soil temperature 5.6°C, moisture 0.7 mm/°C)	Li et al. 2013*
Arenaria tetraquetra L.	No obvious preference for specific microsites	Schöb et al. 2013
Armeria caespitosa (Gómez Ortega) Boiss.	Found on granite and gneiss outcrops and in fellfields and pastures with low temperatures	García-Camacho and Escudero 2008
Azorella compacta Phil.	Excluded from sandy soil, found on equator-facing slopes, grow next to large boulders	Kleier and Rundel 2004*
Azorella macquariensis Orchard	Found more frequently at higher elevations, on ridges, in areas that receive high amounts of solar radiation, and further from the coast	Bricher et al. 2013
Azorella madreporica Clos	Individuals growing in open or rocky areas	Badano and Marquet 2009, Badano, Marquet, and Cavieres 2010
Azorella monantha Clos	Grows on rocky, bare ground	Badano et al. 2006, Badano and Marquet 2008
Azorella selago Hook.f.	Windy areas, not found in areas of reduced light availability	Hughes 1987, Phiri et al. 2009
Dianthus lumnitzeri Wiesb.	Only grows in rocky areas of dry calcareous grassland	Kampichler 1990
Diapensia lapponica L.	Higher densities in exposed areas	Molau 1996
Distichia sp.	Abundant on damp soil at valley bottoms	Cleef 1978
Donatia novae-zelandiae Hook.f.	Found in bogs, often in depressions. Grow where the water table is 10–20 cm below the surface	Lough et al. 1987, Dobson 1975
Gaimardia setacea Hook.f.	Occurs in wet sites where the water table is close to the surface	Dobson 1975
Laretia acaulis (Cav.) Gillies & Hook	Individuals surrounded by open or rocky areas	Badano, Marquet, and Cavieres 2010
Lyallia kerguelensis Hook.f.	Grows in exposed, humid areas on gravel, moraine deposits or syenite-derived soils	Wagstaff and Hennion 2007
Oreobolus cleefii L.E.Mora	Occurs on thin soil peat layer, a high concentration of NO2 and a low soil pH	Bosman et al. 1993*
Oreobolus pectinatus Hook.f.	Found on less fertile soils	Dobson 1975
Plantago rigida Kunth	Low soil concentrations of NO3 and PO4, and higher concentrations of Fe. Abundant on damp soil at valley bottoms	Bosman et al. 1993*, Cleef 1978
Pycnophyllum aristatum Mattf.	Found in rocky areas and on cliffs.	Weberbauer 1931
Saxifraga oppositifolia L.	Common on exposed ridges where the soil is well drained	Grulke and Bliss 1985
Schoenus calcatus K.L.Wilson	Occurs on sand	Parsons and Gibson 2009
Sibbaldia tetrandra Bunge	Found on slopes covered in rock and scree	Pyšek and Liška 1991
Silene acaulis (L.) Jacq.	Grows on ridges with a high wind exposure, slopes with high rock cover and in open grassland	Molenda, Reid, and Lortie 2012
Thylacospermum spp.	Exclusively found in areas of high soil moisture	Qin and Xie 1980
Xenophyllum humile (Kunth) V.A.Funk	Found in sandy soil	Sklenář 2009
Cushion plants in general	Grow in windy areas	Körner 1993
	Individual cushion plants are surrounded by open rocky areas	Badano and Cavieres 2006
	Grow in higher densities on wind-exposed ridges	Bliss, Svoboda, and Bliss 1984
	Occur in-between rocks and pebbles	Hauri 1912
	Damp, salty, coastal cliffs, moors, swamps, mires, snowbeds, springs, and wet alpine areas. Restricted to areas with high light availability	Hauri and Schröter 1914
	Exposed elevations, often also on rocky ground	Hodge 1946
	Limited by high summer temperatures and drought	Holten 2003
	Wind-exposed areas that are free from snow cover	Johnson and Billings 1962
	High-altitude areas with a relatively high rock cover, high exposure and low disturbance	Ramsay and Oxley 1997*

Other aspects of the biology and ecology of cushion plants, in contrast to the drivers of cushion plant species distributions, have been studied in more detail. For example, recent research on cushion plants has focused on cushion plants as ecosystem engineers, facilitators, and foundation species (e.g. reviewed by Kikvidze et al. 2015; Molenda, Reid, and Lortie 2012; Reid, Lamarque, and Lortie 2010), while earlier reviews focused on convergent morphology and anatomy across cushion plant species (e.g. Hauri 1916; Rauh 1939; Weberbauer 1931). Despite this extensive research on cushion plants, succulent cushion-forming species have been relatively poorly studied and less is known about cushion plant species from African ecosystems than the rest of the world.

While cushion plants appear to be most common in abiotically extreme environments, they are often absent from lower elevations and lower latitudes (see e.g. Aubert et al. 2014). This likely reflects that cushion plants cannot tolerate high temperatures (see Buchner and Neuner 2003; Körner 2003). Additionally, cushion plants may be excluded from lower elevations and latitudes by interspecific competition (particularly for light), to which cushion plants are thought to be poorly adapted (Ramsay 1992; see also Armesto, Arroyo, and Villagran 1980). Therefore, cushion plants may be expected to occur in cool and wind-exposed microsites, with low and/or sparse vegetation (i.e. where conditions are abiotically unfavorable and productivity is low), and in regions of stressful climatic conditions (including high wind and low temperature).

The relative influence of biotic interactions and abiotic factors on species’ distributions are thought to vary with the spatial scale of investigation (McGill 2010). While it has been suggested that biotic factors will be more important in determining species occurrence patterns at fine scales, with abiotic factors being the dominant drivers of broad-scale species distributions (Brooker and Callaghan 1998; McGill 2010; Pearson and Dawson 2003), there is no consensus on the relative importance of abiotic and biotic factors for the distribution of species across different scales (Svenning et al. 2014). Therefore, both abiotic and biotic factors should be investigated when aiming to quantify the distribution of species and, where possible, studies should additionally examine species occurrence patterns at different scales to identify all of the key variables driving species distributions.

This study investigates the drivers of a cushion plant species’ distribution at broad and fine scales, using the succulent cushion plant, Euphorbia clavarioides Boiss. (Euphorbiaceae), as an exemplar species. Euphorbia clavarioides occurs in southern Africa (Aubert et al. 2014) and is widespread in the eastern parts of South Africa (Dyer 1947), but little else is known about its distribution. Therefore, the aims of this study were to document the distribution of E. clavarioides, at (1) a fine spatial scale, correlating occurrence patterns with field-quantified biotic and abiotic predictors; and (2) a broad scale, by modeling the species’ broad-scale distribution with climate and topography-related variables. This is the first study quantitatively assessing the drivers of distribution of a cushion plant species.

Methods

Local distribution

Study species and study area

Euphorbia clavarioides (Figure A1) is a poorly studied succulent cushion plant species (see e.g. Archer, Victor, and von Staden 2014). Several other cushion-forming species occur in the genus, including E. acanthothamnos, E. aggregata and E. pulvinata (Aubert et al. 2014), but their ecology is also largely unstudied.

Figure A1. (a) A Euphorbia clavarioides cushion; (b) E. clavarioides growing on an eroded bank with parts of its root system exposed; (c) a flower of E. clavarioides; (d) milky latex exuding from damaged stems; (e) Senecio rhomboideus rooted within an E. clavarioides cushion; (f) E. clavarioides growing in a rocky crevice

The study was conducted during January 2016 and January 2017 in the Golden Gate Highlands National Park (Figure A2), which lies at the foothills of the Maluti Mountains, South Africa (28°31ʹS, 28°37ʹE). The mean daily minimum and maximum temperatures from the closest weather station (Bethlehem, elevation c. 1700 m, 45 km away) are 6.4°C (1980–2009; range: −6.1 to 14.3°) and 22.2°C (12.3 to 29.7°) respectively, with the hottest months generally being December and January, and the coldest June and July (1980–2009; unpublished data, South African Weather Service). Most of the rainfall occurs in spring and summer (between September and March), although precipitation does occur throughout the year, with an average annual precipitation of 700 mm per year (1980–2009; unpublished data, South African Weather Service). The Golden Gate Highlands National Park ranges in elevation from c. 1900 to 2800 m a.s.l. and is probably somewhat colder and wetter than Bethlehem. The national park is chiefly comprised of grassland, with isolated areas of shrubland and Afromontane forest (Daemane, van Wyk, and Moteetee 2010).

Figure A2. Examples of sites where Euphorbia clavarioides was present (a, b) and absent (c, d)

Field methods

To quantify the occurrence of Euphorbia clavarioides, 4 × 4 m quadrats were surveyed at 100-m intervals (or 250-m intervals in relatively flat or homogenous areas) along hiking trails, management tracks, and roads, spanning an elevational range and topographic variation representative of the Gladstone section of the reserve. Quadrats were located 10 m away from hiking paths, roads, access tracks and preplanned routes in inaccessible areas, positioned on alternating sides of the trail, road, track or route. Because of the localized distribution of E. clavarioides, this systematic sampling yielded few presence records and ad hoc presences for the species were thus additionally recorded in the areas where systematic sampling was conducted (similar to le Roux et al. 2013), with a 4 × 4 m quadrat centered over the first E. clavarioides individual observed during sampling. A minimum distance of 100 m was maintained between ad hoc quadrats.

The presence or absence of live E. clavarioides cushion plants within each quadrat was recorded. At each quadrat geographical coordinates (from a handheld GPS unit; Garmin Dakota 20, 3 to 5 m accuracy) and aspect were recorded, along with categorical estimates of quadrat slope (shallow, intermediate, steep), topography (convex, flat, concave) and wind exposure (low, medium, high; estimated visually based on the surrounding topography). Soil depth was measured at each corner of the quadrat using a soil depth probe with a maximum depth of 60 cm (deeper soils were assigned a value of 70 cm). Based on texture, soil type at each quadrat was classified as clay, loam, or sand. The vegetation type within the quadrat was classified as either grassland (less than 10 percent woody plant cover) or shrubland (which included forest). The percentages of rock cover, woody plant cover, and non-woody plant cover were visually estimated. The tallest woody and non-woody plants in each quadrat were measured (excluding reproductive structures). Any visible erosion within the quadrat was noted and the occurrence of mammalian disturbance was recorded based on evidence of trampling, burrowing, grazing, or the presence of dung.

After field sampling, a digital elevation model (DEM, at 30 m resolution, with vertical height accuracy < 16 m) from NASA’s Space Shuttle Radar Topography Mission was downloaded from the United States Geological Survey Earth Explorer (http://earthexplorer.usgs.gov/). A topographic wetness index (TWI) was calculated from the DEM as a proxy for soil moisture in SAGA-GIS version 2.2.7 (Conrad et al. 2015). The TWI combines the local slope angle and the upslope area to estimate the influence of topography on the hydrology of an area (Beven and Kirkby 1979; Sörensen, Zinko, and Seibert 2006).

Analyses

Local distribution

For fine-scale analyses, collinearity between predictor variables was assessed using the generalized variance inflation factor (GVIF), where GVIF values > 3 were considered to be strongly collinear (see e.g. Mizerek et al. 2016). Non-woody cover and rock cover were the only variables, with a GVIF > 3 (Pearson’s correlation coefficient = −0.730, P < .001), and rock cover was excluded from further analyses since non-woody cover better explained the occurrence of E. clavarioides in univariate analyses.

The relationship between the occurrence of E. clavarioides and each predictor variable was initially assessed using univariate generalized linear models (GLMs), assuming a binomial distribution. For continuous predictor variables, the GLMs were run as both linear and quadratic models to determine which model provided a better fit to the data (using a log-likelihood ratio test to test for significant differences between the linear and quadratic models). Model significance was determined using a log-likelihood ratio test to compare the performance of the best fit model (quadratic or linear) for each predictor variable against a null model. Quadrat aspect was converted to northness prior to analyses, where the cosine of aspect (0–360°) is calculated to produce values between 1 (north-facing slope) and −1 (south-facing slope; following e.g. Bader and Ruijten 2008).

The data were then analyzed following a best-subsets approach, modeling the occurrence of E. clavarioides as a function of all combinations of predictor variables, including polynomial terms for those variables that performed significantly better as a quadratic model during the univariate analyses. The best model was identified based on the lowest Akaike Information Criteria (AIC) value (Burnham and Anderson 2004). These statistical analyses were conducted in R statistical software (R Core Team 2016).

Regional distribution

For broad-scale distribution modeling, presence records for E. clavarioides (including its synonyms, see The Plant List 2013) were collected from the Global Biodiversity Information Facility, the South African National Biodiversity Institute, local and online herbaria (J, LYD, PRU), and iSpot nature (http://ispotnature.org, with identifications verified from photographs).

Climatic variables that were thought to be biologically important (annual mean temperature, maximum temperature of the warmest month, minimum temperature of the coldest month, annual precipitation, precipitation seasonality (coefficient of variation of weekly precipitation), and precipitation of the coldest quarter) were downloaded from WorldClim, version 1.4 at a resolution of 5 minutes (current data interpolated for 1960–1990; Hijmans et al. 2005). Some of these variables represent extremes (e.g. minimum temperature of the coldest month, maximum temperature of the warmest month), while others represent total resource availability (e.g. annual mean temperature, annual precipitation). Two additional predictor variables were calculated from the SRTM DEM, version 4.1 (Jarvis et al. 2008), which has a resolution of 90 m. The range and standard deviation of elevation were calculated for each 1 km² grid cell, as measures of elevational range and ruggedness, respectively. All predictor variables were clipped to southern Africa south of 20° S, and all data were reprojected to the Albers equal area conic projection with a resolution of 1 km².

Although repeating analyses of the species’ broad-scale distribution with the same predictor variables that were recorded in the fine-scale analyses would have been valuable, these data are not available at broader scales and the variables considered in the broad-scale analyses show very little variation at small scales. Additionally, since the dominant drivers of species distributions are expected to shift from biotic factors at fine scales to abiotic factors (predominantly climatic factors) at broad scales, we consider our choice of predictor variables appropriate (Brooker and Callaghan 1998; McGill 2010; Pearson and Dawson 2003). All map modifications were conducted in ArcMap v. 10.3.

Occurrence records were modeled using the climatic predictors and elevational range and standard deviation, following a maximum entropy approach in MaxEnt v. 3.3 (Phillips, Dudík, and Schapire 2004), using 5-fold cross-validation (an average of 10049.6 background points from across southern Africa was used for each of the five replicates, with an average of 88.8 training samples; see e.g. Speed and Austrheim 2017). The percentage contribution for each variable was calculated based on the order in which variables were added to the model, while permutation importance was also calculated to reflect the contribution of variables based only on the final model (irrespective of the order in which variables were added; Phillips 2005). A presence/absence map was created from the probability of occurrence map using the equal training sensitivity and specificity threshold (see Liu et al. 2005).

Results

Local distribution

A wide variety of habitat conditions were surveyed, with 367 quadrats surveyed across a broad range of elevations (1836 to 2637 m). Euphorbia clavarioides was recorded in about 4 percent (15 quadrats) of the systematically surveyed quadrats, but was also recorded in an additional 28 ad hoc quadrats.

Of the categorical predictor variables, wind exposure, soil type, and vegetation type significantly predicted the occurrence of E. clavarioides in univariate models, while slope was a marginally significant predictor (Figure 1, Table A1). The species only occurred at sites with medium or high wind exposure, and only in grassland sites. The probability of occurrence of E. clavarioides was lowest in loam soils, relative to clay and sandy soils. Other categorical predictors were not significant in univariate analyses.

Table 2. Best fit generalized linear model for the occurrence of Euphorbia clavarioides. Deviance explained = 60.53 percent, Likelihood ratio test = 160.5, P < .001. See Table A3 for all models differing in AIC < 2 from the best model

Variable	Effect^a	P-value
Slope	shallow < intermediate < steep	< .001
Wind exposure	low < medium < high	.003
Woody height	−0.128 ± 0.054	.001
Woody cover	3.178 ± 1.619	.003
Woody cover^2	−0.375 ± 0.238	.002
Non-woody height	−0.039 ± 0.018	.016
Rock cover	0.138 ± 0.035	< .001
Rock cover^2	−0.001 ± 0.001	.009
Elevation	−0.006 ± 0.002	.003
Northness	3.207 ± 0.792	< .001

^aCoefficient ± SE for continuous predictors; ranking of levels for categorical predictors.

Table A1. Results from the univariate analyses of the occurrence of E. clavarioides in response to each categorical predictor variable (n = 367)

			Null model vs. univariate model
Variable	Order of factor levels	Percentage deviance explained	∆ deviance	P
Slope steepness	shallow < intermediate < steep	2.23	5.913	.052
Topography	flat < convex < concave	0.36	0.962	.618
Wind exposure	low < medium < high	7.56	20.045	<.001
Soil texture	loam < clay < sandy	8.39	22.241	<.001
Vegetation type	other < grassland	3.16	8.371	.004
Presence of erosion	no < yes	1.02	2.711	.100
Presence of mammal disturbance	yes < no	0.05	0.141	.707

Figure 1. Spineplots indicating the proportion of quadrats where Euphorbia clavarioides was present and absent, split by: (a) low (L), medium (M), and high (H) wind exposure; (b) soil type; (c) vegetation type; and (d) slope (Int = intermediate). The width of the bars corresponds to the number of samples in each category

All of the continuous predictor variables were significantly related to the occurrence of E. clavarioides in univariate models (Figure 2, Table A2). Soil depth, non-woody plant height, and TWI were significantly negatively related to the probability of occurrence for E. clavarioides, while northness was significantly positively related (Figure 2). Elevation, non-woody plant cover, rock cover, woody vegetation height, and woody vegetation cover all showed a significant hump-shaped relationship with the occurrence of E. clavarioides (Figure 2).

Figure 2. The relationship between continuous predictor variables and the probability of occurrence of Euphorbia clavarioides, with the best fit regression line indicated. Quadratic models performed significantly better than linear models when explaining variation in the occurrence of E. clavarioides in response to elevation, woody plant height, non-woody plant cover, and rock cover. Ticks on the upper x-axis indicate the values for occurrence records, and on the lower x-axis the values for absence records. Percentage deviance explained (%DE; estimating the proportion of the variability in occurrence explained by the predictor variable) and the significance of each relationship is indicated

The best model of slope, wind exposure, woody plant height, woody plant cover, woody plant cover², non-woody plant height, rock cover, rock cover², elevation, and northness, identified by the best subsets modeling approach, performed well, explaining 61 percent of the deviance in the data (Likelihood ratio test = 160.5, P < .001, Table 2). The direction of the effect for each predictor variable was consistent with their effects in the univariate models (Table A1 and A2) and in all the other models with AIC differences < 2 (Table A3).

Table A3. Results for the top-ranked models (i.e. all models with AIC values that differ < 2 from the best model) predicting the fine-scale occurrence patterns of Euphorbia clavarioides. Log-likelihood ratio tests indicate that all 12 models perform significantly better than expected by chance (all P < .001). Slope: sh (shallow), i (intermediate), st (steep). Wind exposure: L (low), M (medium), H (high). Soil type: L (loam), C (clay), S (sandy). Vegetation type: O (other), G (grassland). Erosion: N (no), Y (yes)

Model rank	AIC	DeltaAIC	Number of parameters	Slope	Wind exposure	Woody height	Woody cover	Non-woody height	Rock cover	Elevation	Northness	Soil type	Vegetation type	Soil depth	Mammal disturbance	Erosion	TWI	Woody cover^2	Rock cover^2	Elevation^2	Woody height^2
1	130.65	0	10	sh<i<st	L<M<H	−0.128	3.178	−0.039	0.138	−0.006	3.207							−0.375	−0.001
2	131.59	0.94	11	sh<i<st	L<M<H	−0.103	2.317	−0.038	0.143	−0.005	3.236	L<C<S						−0.254	−0.001
3	131.64	0.99	11	sh<i<st	L<M<H	−0.115	2.880	−0.037	0.137	0.049	3.248							−0.338	−0.001	0.001
4	131.89	1.24	11	sh<i<st	L<M<H	−0.128	3.134	−0.038	0.140	−0.006	3.079		O<G					−0.369	−0.001
5	132.21	1.56	11	sh<i<st	L<M<H	−0.124	3.191	−0.042	0.150	−0.005	3.228			0.023				−0.380	−0.001
6	132.44	1.79	11	sh<i<st	L<M<H	−0.196	4.046	−0.039	0.140	−0.006	3.152							−0.494	−0.001		0.001
7	132.55	1.90	11	sh<i<st	L<M<H	−0.126	3.161	−0.04	0.139	−0.005	3.254				−0.177			−0.376	−0.001
8	132.57	1.92	11	sh<i<st	L<M<H	−0.129	3.22	−0.039	0.136	−0.006	3.195					N<Y		−0.379	−0.001
9	132.61	1.96	11	sh<i<st	L<M<H	−0.129	3.197	−0.04	0.139	−0.006	3.226						0.042	−0.377	−0.001
10	132.62	1.97	12	sh<i<st	L<M<H	−0.100	2.264	−0.038	0.138	−0.005	3.216	L<C<S				N<Y		−0.240	−0.001
11	132.88	2.23	12	sh<i<st	L<M<H	−0.109	2.871	−0.04	0.152	0.057	3.296			0.03				−0.342	−0.001	0.001
12	132.93	2.28	12	sh<i<st	L<M<H	−0.096	2.258	−0.042	0.157	−0.004	3.218	L<C<S		0.029				−0.25	−0.001

Regional distribution

A total of 226 occurrence records for E. clavarioides were collected (with 43 records coming from this study’s fine-scale surveys), representing 111 grid cells at a 1 km² resolution. The MaxEnt model for E. clavarioides performed well (mean AUC ± SD across fivefold cross-validation = 0.92 ± 0.01). The predicted species distribution is constrained to higher elevations inland of the Great Escarpment (Figure 3(a)), and falls within the boundaries of the grassland biome (Figure 3(b)).

Figure 3. The predicted distribution of Euphorbia clavarioides in southern Africa showing (a) probability of occurrence values, and (b) modeled distribution, with green indicating predicted presence (as calculated using the equal training sensitivity and specificity threshold). Provincial and national boundaries are indicated in (a), while the extent of the grassland biome is shown in (b). Black dots indicate occurrence records

Annual mean temperature had the highest percentage contribution to the model, while minimum temperature of the coldest month had the highest permutation importance (Table A4, Figure A3).

Figure A3. Response curves for the individual predictor variables included in the MaxEnt species distribution model for Euphorbia clavarioides. Blue lines indicate the effect of each predictor variable in isolation, while red lines indicate the effect of each predictor variable after accounting for all the other predictor variables included in the model. Max. Temp. and Min. Temp. refer to maximum and minimum temperature, respectively

Discussion

The distribution of Euphorbia clavarioides was explained well at the fine scale, along with a good prediction of the distribution at the broad scale, with multiple important drivers identified at both scales. E. clavarioides was locally scarce within the Golden Gate Highlands National Park (i.e. a small area of occupancy), but widespread across the grassland biome of southern Africa (i.e. it has a large extent of occurrence). Therefore, this species appears to be limited to specific microsites at fine scales (e.g. habitats with low vegetation cover) within regions with a moderately cool and wet macroclimate.

Our results suggest that, at fine scales, E. clavarioides appears to be intolerant of competition, but able to tolerate abiotic environmental extremes. The species was found in areas of low shading (i.e. areas with short woody plants and low vegetation cover, as seen for another cushion plant species, Androsace alpina, Schönswetter, Tribsch, and Niklfeld 2003), and thus one possibility is that its occurrence could be limited by taller vegetation (following e.g. Ramsay 1992). This is in agreement with our hypothesis that the distribution of E. clavarioides will be limited by interspecific competition due to its prostrate growth form (see also Armesto, Arroyo, and Villagran 1980, although noting that our data are not direct measures of competition). These results also highlight that E. clavarioides is able to tolerate a variety of abiotically stressful conditions. For example, E. clavarioides was more likely to occur in areas with lower soil moisture (in contrast to Cleef 1978; Dobson 1975), where competition with other plants that require higher soil moisture is likely to be lower. E. clavarioides has two attributes allowing it to survive drier conditions: an extensive root system (personal observation; Figure A1), which is a common attribute of cushion plants (e.g. Azorella selago, Frenot et al. 1998; Ternetz 1902), and succulent stems. Furthermore, E. clavarioides was also more common in shallower soils (although not included in the final model), probably reflecting that more competitive taller and/or faster-growing species are excluded from these less favorable areas (Bernard-Verdier et al. 2012). E. clavarioides’ disproportionate occurrence in areas with high wind exposure (see e.g. Grulke and Bliss 1985; Körner 1993; Molenda, Reid, and Lortie 2012), and its complete absence from wind-sheltered sites in the fine-scale study, are also consistent with Körner’s (1993) observation that true cushion plants are absent from areas with low wind exposure. Therefore, low and sparse vegetation, driven potentially by various abiotic mechanisms (e.g. thin, dry soils), appears to be the key driver of fine-scale occurrence patterns for this cushion plant species.

At a broad scale, the species’ distribution was more strongly related to climate than to topographic variation, with E. clavarioides occurring disproportionately frequently under lower mean, maximum, and minimum temperatures. In contrast, annual precipitation was positively related to the occurrence of E. clavarioides. This result likely reflects that E. clavarioides is confined to mountainous regions, which tend to receive higher rainfall. Furthermore, the species’ distribution largely overlaps with that of the grassland biome in South Africa, which typically has relatively short vegetation with a low cover and richness of woody species (Mucina et al. 2006), corresponding to the sensitivity of E. clavarioides to vegetation height and cover at fine scales.

While there is little quantitative information about what drives the distribution of other cushion species, several variables show consistent impacts on the distribution of cushion plant species (e.g. comparing results within Table 1 and between Tables 1 and 2), suggesting that what has been observed in this and previous studies may be representative of other cushion plant species. Therefore, the variables most strongly correlated with the distribution of E. clavarioides are likely to be similar for other cushion-forming species due to the convergent evolution of the growth form (Boucher et al. 2016). Therefore, globally, cushion plants may be expected to be excluded from areas of tall and/or dense plant cover due to competition, and confined to less favorable (i.e. abiotically stressful) habitats, such as areas with high rock cover, shallow soils, low soil moisture, and high wind exposure. Because E. clavarioides occurs in relatively mild environments, other abiotic stresses and disturbances may additionally be important in colder areas. For example, periglacial processes (such as frost heave and freeze–thaw cycles) might be important processes influencing the distribution of cushion plants in colder areas (e.g. via frost heaving and increased substrate instability).

At broad spatial scales, the distribution of E. clavarioides is strongly linked with temperature and precipitation parameters, suggesting that habitat suitability will change directly for this species in response to shifting climatic conditions. At finer scales, while the direct effects of warming temperatures and lower rainfall for cushion plants have been demonstrated experimentally (Cranston et al. 2015; le Roux et al. 2005) and observed through long-term monitoring (Bergstrom et al. 2015; Doak and Morris 2010; see also Frenot, Gloaguen, and Tréhen 1997; Fowbert and Smith 1994), and the sensitivity of cushion plants to overheating has been suggested (due to slow rates of heat loss) (Buchner and Neuner 2003), the indirect effects of warmer conditions may also be very important. Indeed, cushion plants may be indirectly affected through the upslope range expansion of taller and faster-growing species (assuming these species are able to establish in sites with, e.g., thinner soils once thermal constraints are relaxed), leading to an increase in the intensity of interspecific competition (Hughes 2000; Thuiller 2007), which appears to be a strong determinant of cushion species occurrence patterns. Therefore, climate change may have important direct and indirect implications for cushion plants at both coarse and fine spatial scales, potentially having knock-on effects on the community of species that benefit from positive interactions with cushion plants (Molenda, Reid, and Lortie 2012). If cushion plant species generally follow these same trends, then species of this growth form may be particularly useful to monitor as indicators of environmental change in alpine areas.

Acknowledgments

We thank South African National Parks, and especially the management and rangers of the Golden Gate Highlands National Park, for allowing and supporting this research (Permit LEROP 1314 and LEROA 1319). We thank I. Buyens, D. Mogashoa, T. Strydom and J. Potgieter for field assistance. We thank Catherine Kleier and one anonymous reviewer for their valuable comments on the manuscript.

Additional information

Funding

We acknowledge funding from the University of Pretoria’s Research Development Programme and the National Research Foundation (Incentive Funding for Rated Researchers, grant number 95768). MM received additional funding from the South African Association of Botanists and the National Research Foundation’s Scarce Skills Programme.

Appendix

Table A2. Results from the univariate analyses of the occurrence of E. clavarioides in response to each continuous predictor variable (n = 367). Results for quadratic models were only included for those predictor variables where the quadratic model performed significantly better than the linear model

						Null model vs. Linear model		Null model vs. Quadratic model
Variable	Linear coefficient	P	Quadratic coefficient	P	% deviance explained	∆ deviance	P	∆ deviance	P
Soil depth	−0.115	<.001	–	–	14.19	37.630	<.001	–	–
Non-woody plant height	−0.038	<.001	–	–	9.79	25.947	<.001	–	–
Topographic Wetness Index	−0.780	<.001	–	–	11.89	31.534	<.001	–	–
Northness	2.134	<.001	–	–	14.70	38.967	<.001	–	–
Woody plant cover	−0.145	.0227	–	–	7.82	20.724	<.001	–	–
Elevation	0.103	.003	−0.000	.004	8.30	9.801	.002	12.194	<.001
Woody plant height	0.068	.040	−0.001	.033	11.24	20.246	<.001	9.546	<.001
Non-woody plant cover	.035	.276	−0.001	.010	28.10	66.258	<.001	8.237	<.001
Rock cover	.120	<.001	−0.001	<.001	32.37	73.475	<.001	12.347	<.001

Table A4. Variables’ contribution to the species distribution model. Response curves for individual variables are presented in Figure A3

Variable	Percentage contribution	Permutation importance
Annual mean temperature	39.8	10.2
Minimum temperature of coldest month	21.2	41.4
Maximum temperature of warmest month	16.5	9.6
Annual precipitation	14.1	17.7
Precipitation of coldest quarter	5.0	9.4
Altitudinal range	1.5	4.2
Altitude standard deviation	1.3	3.6
Precipitation seasonality	0.7	3.9