Aim

Progress has been made in understanding the relationship between biodiversity and ecosystem functioning (BEF) in both experimental and real-world ecosystems. Yet, we have a limited understanding of the extent to which biodiversity affects ecosystem functioning in heterogeneous environments and whether variation in ecosystem functioning between communities is related to variation in species richness or turnover. Here, we quantify the relative contribution of variation in species richness and species turnover to variation in ecosystem functioning between communities (i.e., the diversity effect) along two tropical elevational gradients.

Location

Andes (Ecuador) and Mt. Kilimanjaro (Tanzania).

Taxa studied

Woody plants, springtails, soil arthropods, ants, and frugivorous birds.

Methods

We collected data on seven ecosystem functions, including biomass and process rates, across six ecosystem types along the two elevational gradients. We then combine the ecological Price equation with the concept of β-diversity to quantify how the diversity effect is shaped by environmental heterogeneity within and across ecosystem types, and whether the effect of environmental heterogeneity is primarily mediated by variation in species richness or species turnover.

Results

The diversity effect on ecosystem functioning increased consistently with environmental heterogeneity on both mountains. Species richness and turnover, on average, contributed similarly to the diversity effect on ecosystem functioning in both mountain regions, but effect sizes varied across functions. The increase in the diversity effect with environmental heterogeneity was primarily mediated by species richness, while species turnover played a secondary role in mediating the effects of environmental heterogeneity.

Main Conclusions

Our study reveals that the diversity effect on ecosystem functioning increases with environmental heterogeneity and that species richness, rather than species turnover, primarily drives this relationship. The dominant role of species richness in mediating the effect of environmental heterogeneity indicates that BEF relationships along environmental gradients are strongly influenced by environmental filters that limit local species coexistence.

Methods: Study area

The Ecuadorian Andes have been formed by tectonic activity whereby the subduction of the Pacific plate under the continental plate uplifted the mountain cahin with its origin in the early Miocene or earlier (Evenstar et al. 2015).

In the Ecuadorian Andes, the data was collected at three sites in southern Ecuador located in and next to the Podocarpus National Park (PNP) (Beck et al., 2019, Bendix et al. 2021). Here, the annual mean temperature ranges between 19°C and 13°C (Richter 2003; Bendix et al. 2006). The study sites experience a tropical humid climate with annual precipitation ranging from approximately 2,300 mm at 1,850 m a.s.l. to over 6,000 mm at 3,100 m a.s.l (Bendix et al. 2006). The wet season, which is exceptionally moist, occurs from April to July, while a period with less precipitation extends from September to December. The central research station is located in the valley of Rio San Francisco (Estación Científica San Francisco; Lat. 3°58’18’’ S (-3.971667), Long. 79°4’45’’ W (-79.079167); Beck et al., 2019). Study sites with premontane rainforest (at elevations of ~1,000 m a.s.l.) were located at the Bombuscaro area within PNP. Study sites with montane rainforest (at elevations of ~2,000 m a.s.l.) were located at the humid eastern slopes of the main Cordillera, in the valley of the Rio San Francisco at the border of the PNP (Reserva Biológica San Francisco). Study sites at the upper montane sites (elevations of ~3,000 m a.s.l.) were located at the Cajanuma area of PNP. Study sites in disturbed forests were located on private land next to the PNP, mostly embedded in pastures, pasture successions with southern bracken fern or exotic tree plantations (pines) (Knoke et al. 2014, 2016; Curatola Fernández et al. 2015).

Mt. Kilimanjaro is a dormant volcano from the Pleistocene or earlier covering multiple climatic and cultivation zones (Wilcockson 1956; van der Plas et al. 2021). Here, data was collected at the southern and southeastern slopes (Tanzania, East Africa; 2° 45′–3° 25′ S, 37° 00′–37° 43′ E) (Peters et al. 2019). Mt. Kilimanjaro covers lowlands with elevations of 700 m a.s.l. to a snow-capped summit with an elevation of 5,895 m a.s.l. However, for this study, we focused on elevations up to ~3,000 m a.s.l. The mean annual temperature ranges within the studied elevational range between 4-18°C at the lowlands to 8.8-10°C at 2,770-3,060 m a.s.l (Peters et al. 2019). The mean annual precipitation peaks at ~2,200 m a.s.l. with 2700 mm (Hemp 2006a). Mt. Kilimanjaro is characterized by a rainy season that lasts from March to May and short rains in November (Hemp 2006a; Peters et al. 2019). Mt. Kilimanjaro has been characterized for a long time by land use and thus natural ecosystem types (in particular in the lowlands) have been degraded by fire, wood extraction, and agroforestry practices (Hemp 2006a, b). Only the area above 1,800 m a.s.l. is protected and has been designated as a national park (Mt. Kilimanjaro National Park; (Hemp 2006b; Peters et al. 2019).

Methods: Study design

Data was collected on 15 to 67 study sites (median = 18 study sites) in the Ecuadorian Andes, south-eastern Ecuador, and 12 to 30 study sites (median = 29 study sites) on Mt. Kilimanjaro, Tanzania. In both mountain regions, we covered an elevational gradient of approximately 1,000-3,000 m.a.s.l. In the Ecuadorian Andes, data was collected in undisturbed premontane forest and disturbed premontane forest at elevations of 960–1,268 m a.s.l., in undisturbed lower montane forest and disturbed lower montane forest at elevations of 1,850–2,450 m a.s.l. as well as in undisturbed upper montane forest and disturbed upper montane forest at elevations of 2,679–2,931 m a.s.l. (Supplementary Methods Table 1). At Mt. Kilimanjaro, Tanzania, sampling was performed in undisturbed and disturbed lower montane forest (i.e., chagga home gardens) at elevations of 1,171–2,097 m a.s.l., undisturbed and disturbed Ocotea forest at elevations of 2,150–2,741 m a.s.l., Podocarpus forest at elevations of 2,720-2,970 m a.s.l. as well as in undisturbed and disturbed podocarpus forest at elevations of 2,753–3,009 m a.s.l. (Supplementary Methods Table 1).

Methods: Environmental variables

The methods for measuring the environmental variables in this study were previously explained in detail. We offer a concise summary with a reference to the respective publications (Appelhans et al. 2016; Peters et al. 2019).

Methods: Soil variables

In the Ecuadorian Andes, soil properties were determined to represent the upper 0.5 m of the mineral soil. Two sampling designs were used depending on the location: Soils were sampled from three profiles on 18 study sites along the elevation gradient from 1,000-3,000 m a.s.l. and at 8-10 profiles along transects in four micro-catchments at 1,900-2,200 m a.s.l. Soil profiles were excavated down to 1 m depth and representative samples were taken from each horizon of the mineral soil. The samples were dried at 40 °C, sieved to <2 mm and an aliquot milled in a planetary ball mill (PM400, Retsch, Haan, Germany). Total C and N content were measured by elemental analysis (Flash 2000 HT Plus, Thermo Scientific, Bremen, Germany) via thermal combustion at 1020 °C. The total P content was determined after microwave digestion (Mars 6, CEM, Kamp-Lintfort, Germany) with HNO3, H2O2, and HF by ICP-OES (5100 OES VDV, Agilent Technologies, Waldbronn, Germany). We calculated a weighted average of the measured element contents to 0.5 m mineral soil depth by considering the thickness and density of each soil layer. Because the soils were free of carbonates, total C corresponds to organic carbon (Corg). To determine the content of soil organic carbon at Mt. Kilimanjaro, samples of the organic horizon and the mineral soil (0-5cm) were collected at five locations per study site. Mineral soil and organic horizon materials were air-dried until constant weight. Soil was sieved to 2 mm with visible root fragments being further removed, while the organic horizon material was shredded prior to grinding with a mixer mill (MM200, Retsch, Haan, Germany) (Peters et al. 2019). All samples were analyzed using dry combustion elemental analyzer (Flash EA1112, Thermo Scientific, Bremen, Germany) to determine the C content (i.e., Corg, because all soils were carbonate-free) at 950 °C. Individual element contents of soil and the organic layer were averaged across the five locations for each study site. To quantify the C/N and N/P ratios, mineral soil samples were collected from soil pits at 10 cm intervals using a standard soil auger. Mineral soil samples were dried at 60 °C for 24 h and root and plant materials removed by sieving to < 2 mm before grinding them for further analysis. The total C and N concentrations of the soils were measured using an elemental analyzer (Vario EL, Elementar) at 950 °C. The total P content was determined using inductively coupled plasma optical emission spectrometry (Spectro Analytical Instruments) after pressure digestion with concentrated HNO3.

Methods: Climate variables

In the Ecuadorian Andes, temperature and precipitation information was gathered for each 1-ha site. The average monthly temperature (i.e., the monthly mean of daily mean temperatures) was acquired using an air temperature regionalization tool developed specifically for the study region (Fries et al. 2012). Monthly mean precipitation (i.e., the average of the monthly precipitation sum) was gathered with a hybrid approach blending ground-based and space-born remote sensing data (using local area weather radar and satellite imagery) with observation data of the meteorological gauge network (Rollenbeck & Bendix, 2011, Bendix et al. 2017).

To assess the climatic conditions at the study sites at Mt. Kilimanjaro, mean annual temperature and mean annual precipitation was assessed. To measure annual temperature, we installed temperature sensors approximately 2 meters above the ground at all study sites (Peters et al. 2019). The sensors recorded temperature at 5-minute intervals for about 2 years. The mean annual temperature (°C) was calculated by averaging all the measurements per study site. Mean annual precipitation (mm yr⁻¹) was estimated by interpolating a 15-year dataset from a network of about 70 rain gauges on Mt. Kilimanjaro using a co-kriging approach (Appelhans et al. 2016; Peters et al. 2016).

Methods: Ecosystem functions

In each mountain region, data on seven ecosystem functions belonging to either biomass stock or process rates were collected (Supplementary Methods Table 2). We were able to directly measure the species-specific functional contributions for 5 out of 7 functions in each community. However, we did not have direct estimates for the functional contributions of litter decomposition by oribatid mites and springtails. Instead, we assumed that the functional contribution of a species was proportional to its relative abundance or biomass at each site.

Methods: Biomass stocks

To determine the biomass stocks of springtails, oribatid mites, ants, and birds in the **Ecuadorian Andes and at Mt. Kilimanjaro at each study site, we combined data on species abundances with data on species-specific per capita mass. Species abundances at each study site were assessed using standardized methods with taxon-specific sampling techniques, which have been detailed in previous studies (for Ecuador: Marian et al., 2018; Santillán et al., 2018; Tiede et al., 2017; for Kilimanjaro: Peters et al., 2016). To estimate the per capita mass for birds, we used existing literature (Dunning 2008; Wilman et al. 2014). For ants and springtails, we applied allometric equations to morphometric measures (i.e., head length for ants, and total body length for springtails) to derive species-specific estimates of per capita mass. Morphometric measurements were taken from up to ten randomly selected individuals per species using a binocular microscope with a calibrated ocular micrometer. The individual biomass estimated using allometric equations closely matched the true biomass determined with a precision scale.

In the Ecuadorian Andes, we included only trees (no shrubs) with a diameter at breast height greater than or equal to 10 cm (Homeier & Leuschner, 2021). We used basal area as a proxy for aboveground biomass (AGB). We calculated the species-specific basal area by summing the individual basal area measurements for each species at each study site. To determine the aboveground biomass stocks of woody plants at Mt. Kilimanjaro, we included all woody plant individuals that were taller than 1.3 m and had a diameter at breast height greater than or equal to 10 cm (trees) or less than 10 cm (shrubs) (Ensslin et al. 2015). We applied pantropical allometric equations to measures of plant height, diameter at breast height, and wood density (Ensslin et al. 2015). We calculated species-specific aboveground biomass stocks by summing individual biomass estimates for each species at each study site.

Methods: Process rates

To measure seed dispersal by birds, resource use by ants, and litter decomposition by microorganisms, process-specific protocols were used. The approach for seed dispersal measurement by birds followed standardized methods, which have been previously described in the literature (Ecuador: Quitián et al., 2018; Mt. Kilimanjaro: Albrecht et al., 2018). In the Ecuadorian Andes as well as at Mt. Kilimanjaro, bird-fruit interactions were monitored on a site measuring 30 m x 100 m at each study site. Birds were observed over four consecutive days, for a total of 25 hours, using binoculars to record interactions with fruiting plants. For each bird species, the number of visits to each fruiting plant species was recorded, as well as their behavior. We calculated species-specific contributions to seed dispersal as the number of visits to all fruiting plants by each bird species, considering only those visits that involved legitimate seed removal events, such as swallowing or carrying away fruits from the mother plants.

To assess resource use by ants, bait experiments were conducted at each study site. In the Ecuadorian Andes, six 50 ml Falcon tubes were used at five subplots, respectively (Tiede et al. 2017). At Mt. Kilimanjaro, thirty 50 ml Falcon tubes were placed along three 50-m transects at ground level (Peters et al. 2014). In the Ecuadorian Andes, each tube contained 15 mL and at Mt. Kilimanjaro 10 mL of nutrient solutions, such as sugar, sugar-protein, protein, water, salt, or oil, with five replicates per nutrient. In the Ecuadorian Andes, ant data was collected at 1,000, 1,500, 2,000, 2,500, and 3,000 m a.s.l. after 2, 3, 4, 4.5, and 5 hours, respectively. Based on the occurrence of ants at the baits, the species-specific contributions to resource use were calculated as the proportion of baits detected by each ant species. In cases where more than one ant species was recorded at a bait, we assigned species-specific contributions to resource use in proportion to the relative abundance of each ant species at that bait. At Mt. Kilimanjaro, the baits were collected after two hours and the number of individuals of each ant species that were present at each bait were recorded.

To study net litter decomposition rates, standardized litter bags with leaves or roots (with a 4 mm mesh and containing 10 g of leaves or roots) were utilized in the** Ecuadorian Andes**. 120 litterbags were placed at each study site, with one at each elevational level, and collected after 6, 12, 24, 36, and 48 months. We followed an established protocol to process the leaves and roots (Marian et al. 2017). Decomposition rates were calculated based on the remaining carbon (CR) in the litterbags at the sampling dates (n) expressed as a percentage of the initial amount of carbon placed in the litterbags (C0). Changes in the remaining nitrogen (NR) were similarly expressed as a percentage of the initial amount of nitrogen placed in the litterbags (N0). We used the following formulas: CR [%] = (Cn/C0) × 100 and NR [%] = (Nn/N0) × 100, with Cn and Nn being the amount of carbon and nitrogen remaining at each sampling date n. We calculated the contents of carbon (CC) and nitrogen (NC) in the litter using the following equations: CC [%] = (Cn/DWn) × 100 and NC [%] = (Nn/DWn) × 100, with DWn being the dry weight of litter remaining at sampling date n.

At **Mt. Kilimanjaro **litter decomposition rates were assessed using litterbags filled with dried maize straw (10 cm x 15 cm, 20 μm x 20 μm mesh size, and containing 5 ± 0.05 g of maize husks; Peters et al., 2019). Three bags were placed at each study site and collected after 69 to 86 days. To process the leaves, established protocols were applied. Due to logistical reasons, bags at lower elevations were exposed for a longer time than bags at higher elevations. To adjust for these differences, decomposition rates per day were calculated using the equation k = −ln(m_LOI / m_OAF) / t, where m_LOI is the weight after loss-on-ignition, m_OAF is the original ash-free weight, and t is the number of days the bags have been exposed. The decomposition rates per study site were averaged.

In both the Ecuadorian Andes and at Mt. Kilimanjaro, we did not have direct measures of species-specific contributions to litter decomposition. Thus, we estimated the specific contribution of each species to decomposition at each site based on the relative abundance of each species at that site. In the Ecuadorian Andes, decomposition rates were related to the abundance of oribatid mites (Marian et al. 2018) and at Mt. Kilimanjaro to springtails (Peters et al. 2016).

Methods: Analytical approach

Quantifying the effect of diversity on ecosystem functioning. The framework is generally applicable to any ecosystem function that is determined by the summed functional contributions of individual species. It is based on a community matrix F (n × s), which summarizes the contribution of s species to a given ecosystem function in n study sites (which we refer to as communities; Albrecht et al., 2021). Each element *f_ij *of matrix F represents the contribution of species j to a given ecosystem function at study site i. Using matrix F, we constructed three binary matrices P, Q, and O. Matrix P is a binary species incidence matrix (with the same dimensions as F, in which *p_ij *= 1 if f_{ij *> 0 and *pij} = 0 otherwise). Further, matrix Q is the complement of matrix P (Q = 1 – P), so that q_ij = 1 if f_{ij *= 0 and *qij} = 0 otherwise. And finally, matrix O is given by the sum of matrices P and Q, so that all elements o_ij = 1. The total difference in the magnitude of a given ecosystem function between two communities i and j is given by the element ∆f_ij of the n × n square matrix ∆F = FO^T – OF^T. The difference in the magnitude of an ecosystem function between two communities i and j caused by changes in the functional contributions of shared species is given by the element ∆s_{ij of the n × n square matrix ∆S = FP^T – PF^T. Lastly, the difference in the magnitude of a given ecosystem function between two communities i and j caused by differences in species richness and species turnover between communities is represented by the element ∆*dij} of the n × n square matrix ∆D = FQ^T – QF^T.

Equation (1) comprises three components that possess the dimension of the corresponding ecosystem function and an expected value of zero if there is no overall difference in this function between communities i and j (∆f_ij = 0), no difference in function because of differences in species richness and species turnover (∆d_ij = 0), or no difference in function because of differences in the functional contributions of shared species (∆s_ij = 0). The components can be positive or negative, and the diversity and shared species components may cancel each other out, resulting in a net difference in the function of zero. Note that the formula presented in equation (1) represents a simplified version of the of the 5-part and 3-part Price equation presented in Fox and Kerr (2012) and Bannar-Martin et al. (2018). However, the formulation described here overcomes an important limitation of previous frameworks (Fox and Kerr 2012; Bannar‐Martin et al. 2018), as it can be applied to any pair of communities, regardless of whether they have species in common.

Using equation (1), we calculated the relative contribution of the diversity of the focal taxa to variation in their ecosystem functions between communities, which we refer to as the diversity effect (Y).

Where the absolute values of ∆d_ij and ∆s_ij are denoted as |∆d_ij| and |∆s_ij|, respectively. Absolute values are used because ∆d_ij and ∆s_ij can either be positive or negative, thus computing the mean of these values would mask the total contribution of both factors. The diversity effect is a dimensionless measure that ranges from 0 to 1. The diversity effect equals 0 if all differences in ecosystem functioning between communities result from differences in the functional contributions of shared species, which can be attributed to differences in abundance or individual performance. Conversely, the diversity effect equals 1 if all differences in ecosystem functioning result from the combined effects of differences in species richness and turnover between communities. For worked example scenarios of how changes in species richness, species turnover, or the functional contribution of shared species are related to the diversity effect see Fig. S1.1 in Supporting Information.

Because our metric for the diversity effect is an absolute measure, we also quantified the correlations of the raw values of ∆f_i with ∆d_i and ∆_si to assess how strongly total variation in ecosystem functioning between communities (∆f_i) is related to variation in ecosystem functioning, due to differences in species richness and turnover (∆d_i), and to variation in ecosystem functioning, due to changes in the contributions of shared species (∆s_i; Fig. S1.3). Across all 14 functions as well as for individual functions, we observed a strong positive correlation between ∆f_i and ∆d_i~ (Fig. S1.3 and Table S1.11). Moreover, across the 14 ecosystem functions, the diversity effect was strongly positively related to the correlations between ∆fi and ∆di of individual functions (r = 0.78; Fig. S1.3) and strongly negatively related to the correlation between ∆*fi~ *(and ∆s_i~ ~(r = -0.89; Fig. S1.3). This indicates that our metric of the diversity effect is able to quantify the contribution of differences in species richness and turnover to the total variation in ecosystem functioning between communities.

Quantifying variation in species richness and turnover. We used the Jaccard index to partition the variation in species composition (β) into variation resulting from differences in species richness (R) and species turnover (T) between n communities (Albrecht et al. 2021). The species richness and turnover components form an additive partition of the total variation in species composition, meaning that β = R + T (Legendre 2014). To estimate the total variation in species composition and its two components, we calculated the number of shared species, unique species in community* i*, and unique species in community* j*, using the species incidence matrix P and its complement Q. The number of shared species between communities* i* and j is represented by element a_ij of the n × n matrix A = PP^T, the number of species unique to community i by element b_ij of the n × n matrix B = PQ^T, and the number of species unique to community j by element c_ij of the n × n matrix C = QP^T. Using these matrices, we can express the total variation in species composition (β), the richness component (R), and the turnover component (T).

Quantifying environmental heterogeneity. We used the Gower distance to determine the environmental distance h between sites based on a set of environmental variables. The Gower distance between two sites* i* and j equals the mean difference in environmental variables across all environmental variables after standardizing the environmental variables by their ranges.

In equation (6), n is the number of environmental variables, x_ik and x_jk are the values of variable k on study sites i and* j.* The Gower distance is preferred over the Euclidean distance as it is less sensitive to extreme values and facilitates the inclusion of categorical measures. Moreover, the range standardization ensures that each environmental variable contributes equally to the distance metric, and the maximum value of the distance function is 1. As some study sites had missing data for some environmental variables, we calculated the pairwise distances by using a pairwise deletion of missing observations. We defined environmental heterogeneity (H) as the mean environmental distance between n combinations of study sites.

Quantifying how environmental heterogeneity modulates the effect of diversity on ecosystem functioning. To quantify the direct effect of environmental heterogeneity on the diversity effect as well as the indirect effects that are mediated by variation in species richness and species turnover within and across ecosystem types, we segmented the distance matrices for Y, R, T and H into submatrices that contained comparisons within single and across multiple ecosystem types (Fig. 1d). Then, we averaged these submatrices to obtain estimates for Y, R, T and H for comparisons between sites based on different numbers of ecosystem types (Fig. 1e).