Our study was carried out on Tenerife, the largest of the Canary Islands, situated in the mid-Atlantic between 27° 38′ and 29° 25′ off the West-African coast. The highest point on the island is the summit of the volcano Mount Teide at 3,718 m a.s.l. The crossover of different weather system causes the humid air masses of the north-eastern trade winds to remain on the northern slopes. Consequently, the southern slopes of Tenerife, where this study’s elevational gradients are located, are characterized by an arid and warm climate (Fernández-Palacios, 1992). Temperature decreases along the elevational gradient with a lapse rate of approximately 0.5°C per 100 m (Morales Matos and Pérez González, 2000). Precipitation increases from approximately 100 mm per year at the coast to a maximum of approximately 500 mm at mid-elevations and moderately decreases again toward the summit, but with high variability due to the island’s complex topography (Morales Matos and Pérez González, 2000; AEMET, 2012). This characteristic distribution of precipitation and temperature, the history of volcanic activity, and the specific topography of the island, lead to the formation of many different types of soil (Morales Matos and Pérez González, 2000).

Data collection took place in April and May 2018. All data and samples were retrieved from the permanent plots which were established in 2008 following the design of the standardized mountain road survey of the Mountain Invasion Research Network (MIREN) (Bacaro et al., 2015; Haider et al., 2021). Thereby, three roads on the island’s southern slopes were selected, spanning from the coast to the crater of Mt. Teide. By dividing the elevational range of each road into 19 equal portions, 20 sampling locations were chosen evenly over an ascent of approximately 2,350 m. At each sampling location, two plots with 2 m × 50 m each were established, one directly adjacent and with the long side parallel to the road (hereinafter roadside plot), and a second one perpendicular to the roadside plot, at 50–100 m away from the road (hereinafter distant plot). Steep topography and private property prevented sampling in some locations, resulting in a total of 111 plots.

In each plot, all vascular plant species were identified, and their cover was visually estimated and recorded via classes from 1 to 7 (where class 1 corresponds to a cover of up to 1%; 2 = 1–5%; 3 = 6–25%; 4 = 26–50%; 5 = 51–75%; 6 = 76–95% and 7 = 96–100%). For the leaf trait sampling, we primarily aimed to collect the most abundant and most frequent species, both native and non-native species, and all life forms (grasses, forbs, shrubs, and trees). For each species, we collected 10–50 leaves from one to three typical populations, so that the sampling covered the species’ elevational distribution. Later, trait values of all populations of a species were averaged to get a single value for each trait of each species. With this approach, we obtained for each plot trait information of species making up at least 80% of the relative plot cover (except for five plots); an established threshold for a reliable characterization of the functional properties of a community (Pakeman and Quested, 2007). In total, 332 leaf samples for 189 species (out of 270 species recorded in the survey) were collected.

For a subset of the leaf samples, we conducted laboratory analyses of five traits associated with an acquisitive (SLA, leaf N, and leaf N:P ratio) and a conservative growth strategy (LDMC and leaf C:P ratio). In a second step, we used this subset of samples as “calibration samples” to predict the trait values for the remaining samples via the high-throughput method of near-infrared reflectance spectroscopy (NIRS) as described for example by Foley et al. (1998) and applied in previous projects (Kühn et al., 2021; see Supplementary Material for more details).

In the field, the calibration samples were kept in bags with a wet tissue to prevent water loss during the sampling day, then weighed and scanned with a flatbed scanner on the same day. Leaf area was calculated using WinFOLIA (Regent Instruments, Quebec, CA). The samples were subsequently dried for 72 h at 80°C in a drying oven and weighed again. LDMC was calculated as the ratio of leaf dry weight to leaf fresh weight, and SLA as the ratio of fresh leaf area to leaf dry weight (Kleyer et al., 2008; Pérez-Harguindeguy et al., 2013). The dry samples were milled into fine powder using a ball mill (MM 400, Retsch, Haan, DE). Following a nitric acid digestion and adding ammonium heptamolybdate and ascorbic acid, the concentration of leaf phosphorus was measured by visual assay with a Photometer SmartSpec^TM 3000 (BIO-RAD). The total carbon and nitrogen concentrations were measured via gas-chromatography with the Dumas method (Vario EL Cube, Elementar Analysensysteme, Langenselbold, DE). From these single element values, we calculated leaf N:P and leaf C:P ratio.

Species richness was defined as the total number of species recorded per plot. Community-weighted trait means (CWMs) were calculated for each of the five leaf traits (SLA, LDMC, leaf N, leaf N:P ratio, leaf C:P ratio) and for each plot using the SYNCSA package (Debastiani and Pillar, 2012) in the R statistical environment (R Core Team, 2020):

where p_i is the relative cover of species i in the community, and trait_i is the trait value of species i.

To quantify functional diversity (FD), we calculated Rao’s quadratic entropy (Rao’s Q; Rao, 1982; Botta-Dukát, 2005) for each plot:

where the species’ trait distances are weighted by their relative abundance in the community. N is the number of species in the plot, p_i and p_j are the relative cover of species i and j, and d_ij is the trait distance between species i and j. Thus, Rao’s Q gives the mean functional distance between any two species in a plot. We included all five leaf traits in the calculation of Rao’s Q (multi-trait functional diversity). The calculation of functional diversity was done with species trait values scaled to unit variance, in the FD package in R (Laliberté and Legendre, 2010). Elevational variation of CWMs, functional diversity (Rao’s Q) and species richness and the effect of road proximity are shown in Supplementary Figure 1 and Supplementary Table 1. See Supplementary Table 2 for an overview of all the leaf trait, soil and temperature variables used in this study.

From three locations in each plot, we collected a bulked soil sample from the upper 10 cm of the mineral soil layer. The samples were sieved (<2 mm) and remaining, visible plant particles were removed. Soil pH was measured in 1 molar KCl. For gas-chromatographically measuring total soil carbon and nitrogen concentrations, fresh soil was dried for 72 h at 105°C and then milled to fine powder. To determine the concentrations of the base cations calcium, magnesium and potassium, fresh soil was percolated in a solution of 0.1 molar barium chloride. Afterward, the cation concentrations were measured by atomic absorption spectrometry (ContrAA 300 AAS, Analytik Jena, Jena, Germany). The concentration of hydrogen ions was measured by determining the pH of the percolates, and the cation exchange capacity (CEC) was calculated as the sum of ion equivalents of all measured exchangeable cations (Ca, K, Mg, and H). For the determination of bioavailable phosphorus, the soil samples were analyzed using the method proposed by Olsen (1954). Like for the leaf samples, the measurement of soil P was done photometrically after adding a molybdate reagent solution and an ascorbic acid solution. The ratios of N:P and of C:P were calculated.

Near-surface soil temperature was measured using miniature loggers (HOBO 64K Pendant, with 0.53°C accuracy, Bourne, MA, United States). They were installed in the center of the roadside plot and 50 m away from it, at the start of the distant plot, at a depth of 3 cm below the soil surface. They were programmed to log at 1.5 h intervals from May 2018 until April 2019. Data quality control was done in R (R Core Team, 2020), and days with less than 13 measurements per day (out of 16) were excluded when calculating the daily averages, as well as months with less than 27 days with enough measurements. After data cleaning, daily and monthly averages and the bioclimatic variables Bio10 (mean temperature of the warmest quarter) and Bio11 (mean temperature of the coldest quarter) were calculated using the climvars package (Maclean and Gardner, 2021), which follows the definitions of WorldClim (Hijmans et al., 2005). By taking into account climate of different seasons and including the periods where temperature stress might be strongest for plant growth, survival and reproduction, these two bioclimatic variables are ecologically and physiologically relevant and are widely used in species distribution studies (Lembrechts et al., 2019). We did not choose to use annual means, for example, since they cover periods that have less influence on foliage. Temperature variation with elevation and proximity to roads is shown in Supplementary Figure 2.

All statistical analyses were carried out in R version 4.1.0 (R Core Team, 2020). To get an overview about how the measured chemical soil properties are correlated, a principal component analysis (PCA) was performed (function “prcomp”) and visualized with the R packages “FactoMineR” (Lê et al., 2008) and “factoextra” (Kassambara and Mundt, 2020).

To analyze elevational changes of chemical soil properties and the influence of road proximity, we used linear mixed-effects models (R package “lmerTest”; Kuznetsova et al., 2017). Models were fitted with each of the eight investigated soil variables as response, and elevation, proximity to roads (roadside or -distant plot) and their interaction as fixed effects. Sampling location nested in road identity was added as random factor. Elevation was scaled prior to analyses. Because the interaction term was never significant, a simplified model without the interaction was fitted. P-values were calculated from F-statistics of type III sum of squares with Satterthwaite approximation to estimate the denominator degrees of freedom. We visually checked, if model assumptions were met, and transformed the response variables if necessary. This was the case for soil N and soil Mg concentrations, soil C:P and N:P ratios (log-transformed) and for soil Ca and soil K concentration and CEC (square-root-transformed). Additionally, we fitted models with the same structure, but including the quadratic term of elevation, using its second-order polynomial. We compared the AIC (Akaike’s Information Criterion) of linear and quadratic models and considered those with lower AIC as the better models. For soil pH, C:P ratio and K, the quadratic model had a lower AIC than the linear one.

To test and compare the effects of soil properties, temperature and road proximity, we used a multivariate generalized linear mixed model based on a Markov Chain Monte Carlo sampler (“MCMCglmm” function in the R package “MCMCglmm”; Hadfield, 2010). Response variables were the CWMs of SLA, LDMC, leaf N, leaf N:P, and leaf C:P, Rao’s Q of these traits (multi-trait functional diversity) and species richness. All response variables were scaled and, if necessary, transformed prior to scaling. This was the case for SLA, leaf C:P ratio, Rao’s Q, and species richness (log-transformed). As predictors we included temperature of the warmest and of the coldest quarter, proximity to roads (roadside or distant plot), and chemical soil properties, summarized by using the values of the first three axes of the PCA described above. We also permitted the interactions between road proximity and soil variables. All predictors were scaled, which allows to directly compare their effect sizes. The random effect of sampling location nested in road was added to control for the fact that all response variables were measured in the same plots. The model was fitted with a Gaussian distribution for all response variables. We ran the model with 100,000 iteration steps, discarded the first 10,000 samples as burn-in and thinned the remaining samples by storing every 25th iteration step.

We found three major axes of variation in chemical soil properties. The first PCA axis summarized cation concentrations in the soil (hereafter referred to as “cation axis”), and explained 41.1% of the total variation (Figures 2A,C). The soil cations Ca, Mg and the CEC were strongly negatively correlated with the cation axis (i.e., high positive values along this axis mean a low concentration of soil Ca and soil Mg, and low CEC). The second PCA axis summarized the concentration of soil N as well as soil C:P and N:P ratios (hereafter referred to as “nutrient axis”). This axis explained 28.6% of the total variation (Figures 2A,B). Soil N and both soil C:P and N:P ratios were negatively correlated with the nutrient axis. The third PCA axis summarized soil acidity and the concentration of soil K (hereafter referred to as “soil acidity axis”), explaining 14.2% of the total variation (Figures 2B,C). Soil pH and K were negatively correlated with the soil acidity axis.

Most of the investigated chemical soil properties (pH, Ca, K, Mg, and CEC) decreased significantly with increasing elevation (Figures 3A,E–H and Table 1). While soil Ca, Mg and CEC decreased continuously, pH and soil K decreased more steeply in the lower part of the elevational gradient. Soil C:P ratio peaked at mid elevations, while while soil N:P ratio did not respond to elevation (Figures 3C,D and Table 1). Soil pH values were significantly higher in roadside plots compared to distant plots (Figure 3A and Table 1). For concentrations of N and Mg, the opposite was observed (Figures 3B,G and Table 1).

Mean temperature of the warmest quarter had a positive effect on leaf N and a negative effect on leaf C:P ratio (Figure 4 and Supplementary Table 3). Mean temperature of the coldest quarter had a positive effect on SLA and species richness (Figure 4 and Supplementary Table 3). SLA increased with increasing CEC and soil cation concentrations (negative effect of the “cation axis,” Figure 4 and Supplementary Table 3). Species richness was higher in roadside plots compared to distant plots. Among the environmental variables tested, mean temperature of the coldest quarter had the strongest effects on vegetation characteristics (strongest on species richness, followed by the effect on SLA). The second strongest influence was exerted by mean temperature of the warmest quarter (on leaf N and second on leaf C:P ratio), followed by road proximity (on species richness) and soil cations (on SLA).

Soil properties, such as soil acidity and the amount of nutrients, vary with environmental conditions and can reflect e.g., climatic differences, such as lower temperatures or aridity. This effect can be seen in particular along elevational gradients due to the strong environmental variation. In our study, the CEC, base cation concentrations (Ca, Mg, K) and pH decreased significantly along the elevational gradient on the south of Mt. Teide, on Tenerife. This elevational response of exchangeable cations is overall consistent with the results of another study conducted from mid to high elevations (1,400–3,100 m a.s.l.) on the northeastern slopes of Mt. Teide (Köhler et al., 2006), and was also found, for example, on the Galápagos Islands (Dinter et al., 2020). In contrast, a study conducted in the equatorial Andes did not observe clear elevational changes of soil base cations (Unger et al., 2012), and in Costa Rica, CEC, Ca and Mg concentrations were found to be highest at the highest elevation sampled (2,600 m a.s.l.; Grieve et al., 1990). Soil nutrient concentrations depend on multiple factors which do not necessarily change linearly with elevation (Körner, 2007). Differences in bedrock, soil age, climatic conditions or vegetation between regions result in region-specific soil property changes along elevational gradients (Köhler et al., 2006; Dinter et al., 2020). The distinct vegetational cover along the elevational gradient and thus input of litter between the study site in the equatorial Andes (Unger et al., 2012) and our study makes it difficult to sufficiently compare the results. Although soil age and bedrock are similar on the study site in Costa Rica, the investigated soil shows much greater depth and different structure. In contrast to our results, soil pH was found to increase from mid-elevations toward the tree line (around 2,100 m a.s.l.) on the northeastern slopes Mt. Teide (Köhler et al., 2006). These opposing responses of soil pH in the same mountain system are likely to be caused by different elevational patterns of precipitation: annual precipitation mainly increased within our studied elevational gradient from the coast (∽100 mm/y) to the Canary pine forest belt (∽500 mm/y), plateauing in the highest study plots. Vice versa, the lowest part of the elevational gradient studied by Köhler et al. (2006) was located in the most humid belt at mid-elevations (∽1,050 mm/y), and precipitation decreased from there toward higher elevations (480 mm/y measured at 2,160 m a.s.l., in short distance to our highest study plots; Köhler et al., 2006). Hence, pH seems to be negatively correlated to precipitation and thus soil moisture—a pattern also reported from elevational gradients, for example, on Hawaii (Vitousek and Chadwick, 2013), on the Galápagos Islands (Dinter et al., 2020), and in Italy (Dixon et al., 2016).

In our study, soil C:P ratio peaked at mid-elevations, within the belt of Canary pine forest. A high ratio of soil C:P might indicate an efficient use of phosphorus by plants and microorganisms (Zhang et al., 2016), which in our case is likely to be supported by moister conditions in this elevational belt. On the other hand, the decomposition of pine needles is a slow process, potentially slowing down the relocation of phosphorus from dead plant material into the soil. At higher elevations, colder temperatures might further reduce decomposition rates and thus also increase the soil C:P ratio. The latter was found for example in a study in Northern China, where the soil C:P ratio significantly increased with elevation (Xu et al., 2019).

Soil N and N:P ratio showed no significant variation along our studied gradient. In our study, soil analyses were only conducted on mineral soil which has lower N concentrations than the organic layer. It was found that N concentrations were highest in the organic layers due to decomposing litter (Köhler et al., 2006). Thus, the increased litter in the pine forest resulted in higher N concentrations in mid-elevation soil samples which were, however, not significantly different to other samples. The lack of significant elevational changes of N and P coincides with the results of the studies conducted on the northeastern slope of Mt. Teide on Tenerife (Köhler et al., 2006), in the equatorial Andes (Unger et al., 2012), and on Mt. Gongga in China (Tan and Wang, 2016). Overall, young soils of volcanic origin—as found on Tenerife—are usually poorly developed and have lower nutrient concentrations and lower CEC compared to older soils (Köhler et al., 2006; Tan and Wang, 2016).

We found significantly higher concentrations of soil N and Mg in plots away from the road compared to plots adjacent to the road. This contradicts our expectation that road verge soil might contain higher nutrient concentrations as a consequence, for example, of road dust, and moister conditions from road drainage (Forman and Alexander, 1998; Müllerová et al., 2011). Our results might be explained by higher bulk density caused by vehicles accessing the road verge which might lead to lower soil permeability and reduced water infiltration (i.e., greater surface runoff) and nutrient flows (Johnston and Johnston, 2004; Makineci et al., 2007; Deljouei et al., 2018; Yang et al., 2018). Other studies which investigated the effect of road proximity on soil nutrients also found higher soil N and soil Mg concentrations away from the road. They explained this increase of nutrients by higher vegetation density and higher N mineralization rates (Johnston and Johnston, 2004; Deljouei et al., 2018). Moreover, along roads, due to mechanical disturbance mainly through vehicles, finer soil material is removed and the more coarse-grained mineral soils have lower concentrations of the major nutrients P and N as well as of exchangeable cations (Johnston and Johnston, 2004).

Besides the resulting disturbance along roadsides, the construction of roads might have severe impacts on soil properties. These are especially strong, if alkaline construction material is used in naturally acidic environments (Müllerová et al., 2011). As a consequence, soil pH might increase next to roads—as found in our study as well as in other regions (Müllerová et al., 2011; Deljouei et al., 2018).

The main objective of our study was to evaluate how two important components of the environment—chemical soil properties and temperature—shape vegetation characteristics along the elevational gradient on Tenerife. Moreover, we analyzed the influence of roads on vegetation, and how the effects of soil are modified through road proximity. The strongest effects detected on different vegetation aspects analyzed in this study was that of temperature of the coldest and the warmest quarter, both variables decreasing with elevation. Species richness and SLA (acquisitive trait) significantly increased with mean winter temperature, and leaf N (acquisitive trait) increased with mean summer temperature, while leaf C:P ratio (conservative trait) decreased with mean summer temperature (Figure 4). The effect of temperature on vegetation has been widely studied (Kreft and Jetz, 2007; Šímová et al., 2011) and shifts of functional traits toward more conservative trait values with decreasing temperature and increasing elevation have been observed (Read et al., 2014; Schellenberger Costa et al., 2017). Cold temperatures and frost affect plants by causing physiological stress via cellular membrane damage on leaves and lower metabolic activities, which ultimately affect growth and reproduction (Yadav, 2010). At the high elevations of Mt. Teide, where we measured winter temperatures in the soil as cold as 1.7°C (Supplementary Figure 1), the specialist sclerophyllous shrubs dominate the communities. They seem to cope with the extreme conditions by having smaller, thicker leaves (low SLA) with high carbon content and low nutrient concentrations (high leaf C:P ratio). Because only specialist species are able to survive in such harsh conditions at high elevations, temperature also has a filtering effect on species richness toward colder sites, as observed as well on the neighboring island La Palma (Irl et al., 2015), but also at the global scale (Wright et al., 2017; Bruelheide et al., 2018) and e.g., for Kilimanjaro (Peters et al., 2019).

In our study, soil cation concentrations were the only chemical soil properties exerting a significant effect on the vegetation. SLA increased with increasing soil Ca and Mg concentrations and higher CEC. More favorable soil conditions at low elevations (with higher Ca, Mg and CEC) seem to be an important factor—in addition to temperature—which favors species with acquisitive strategies in these communities. As soil conditions get less favorable (i.e., with increasing elevation), the communities are characterized by a more conservative strategy with smaller SLA. Our results deviate from those found at Mt. Kilimanjaro, where soil nutrients (N, P, Ca, Na, K) had a negative effect on LDMC, and a positive effect on leaf N and leaf P, but not on SLA (Schellenberger Costa et al., 2017). Despite the different single results in terms of which traits were affected by soil nutrients, the overall result from Schellenberger Costa et al. (2017) also showed soil nutrients to favor acquisitive strategies. The lack of an effect from our soil nutrient axis (determined by soil N concentration, soil C:P, and N:P ratio) might be due to the very low concentrations of N and P in the poorly developed soils at Mt. Teide, as also found in another study on Tenerife (Köhler et al., 2006) and in general for Mediterranean soils (Torrent, 2005). Even though the N concentrations we found in the soil along the southern slopes of Tenerife spanned over a wide range (0.02–1.57%), comparable to values found in alpine grasslands (0.2–1.2%, Körner, 2003), only a few sites contained more soil N than found by Köhler et al. (2006) on the northern slopes (0.02–0.26%). The overall poorly developed soils in our study area might also explain why we could not detect a positive effect of soil quality on species richness, although soil quality, and species richness almost paralleled each other along the elevational gradient.

We expected species richness to be higher next to roads compared to non-road habitats, because of the disturbance in the roadside habitats, which should weaken competition among plants through an increase in bare soil, light and nutrients and thus allow more species to co-exist (Catford et al., 2012). Indeed, communities at the road verge harbored more species than the communities away from the road in our study region. However, this was likely primarily because of the heterogeneity generated by the human activities associated to roads and more available space and light than increased nutrient availability (lower soil N and Mg concentrations next to the road). The removal of biomass during road maintenance in Tenerife, for example, prevents continuous cover of highly competitive species by creating empty spots and opportunities for other species to grow. We also expected disturbance to favor species with acquisitive traits, such as high SLA (McIntyre et al., 1999; Bernhardt-Römermann et al., 2011), which allow for a quick uptake and use of the available resources. However, in our study, CWMs of the functional traits investigated did not differ between communities next to and away from the road. This finding, together with the absence of an increase in functional diversity next to roads, is surprising and difficult to explain. A possible explanation might be that during the warmest quarter, temperature in roadside plots was significantly higher than in distant plots (Supplementary Material). Therefore, the effect of road disturbance might have been captured by the overall stronger effect of temperature. Alternatively, the new species favored by disturbance at roadsides might be functionally similar and thus do not significantly contribute to functional diversity or changes in CWMs.

We found no support for the expected interaction between soil properties and road disturbance. We expected that the same amount of soil nutrients or cations, for example, would have a different impact on roadside vegetation compared to the distant vegetation, as these nutrients might not be available for plant uptake due to higher soil compaction from vehicle traffic. However, these interacting effects were not detected in any of the vegetation characteristics analyzed here. This might be explained by a less pronounced difference in soil bulk density or grain size distribution in roadside and distant plots than we presumed, and therefore a similar nutrient availability for plants.

Surprisingly, functional diversity was not affected by any of the environmental factors considered. While other studies found functional diversity to decrease through environmental filtering (but see Schellenberger Costa et al., 2017; López-Angulo et al., 2020; Midolo et al., 2021) and to increase through competition (Mayfield et al., 2010), we could not detect any differences. While our selection of traits was done carefully in order to encompass the aspects of vegetation relevant to our questions, we also acknowledge that the leaf traits considered in this study do not capture all dimensions of trait variation. For example, seed and dispersal traits may also be important for community assembly (although not independent of life form, which make them less suitable for this study). Since the studied environmental variables are all related to competition intensity, we consider the choice of leaf traits related to the species’ competitive and growth ability as appropriate. Thus, regarding the lack of a response of functional diversity in the study presented here, we emphasize that although elevation is often criticized of being merely an umbrella term, it also combines the effects of several environmental factors, thus revealing ecological patterns which cannot be detected with individual analyses of single variables (Körner, 2007).

In our study, temperature had the strongest effects on vegetation characteristics along the elevational gradient of Tenerife—stronger than chemical soil properties and disturbance. As decreasing temperature gradients are general patterns in mountains, our findings align with observations from elevational gradients in other regions. In addition, our results are in accordance with studies explaining broad-scale patterns of species richness and functional trait variation, particularly those looking at continental and global geographic patterns (Hawkins et al., 2003; Kreft and Jetz, 2007; Šímová et al., 2011; Bruelheide et al., 2018). While elevational temperature patterns are quite linear, soil patterns are more idiosyncratic, and thus not simply transferrable across regions. At the species level, a global compilation study found soil measured at the same site as the traits to be a more important driver of trait variation than temperature data extracted from smaller resolution sources (Ordoñez et al., 2009; but see Liu et al., 2017). Together with the results from our regional study, this emphasizes the importance of including data from local environmental conditions for climate, disturbance and fine-scale soil conditions for understanding biodiversity patterns.