Area modulates the effect of elevation but not of land use or canopy on tropical plant species richness

One of the few general patterns in ecology is the increase of species richness with area. However, factors driving species-area relationship (SAR) are under debate, and the role of human-induced changes has been overlooked so far. Furthermore, SAR studies in tropical regions, in particular in multilayered rain forests are scarce. On the other side, studies of global change-induced impacts on biodiversity have become increasingly important, particular in the tropics, where these impacts are especially pronounced. Here, we investigated if area modulates the effect of land use, elevation and canopy on plant species richness. For the first time we studied SAR in multilayered tropical forests considering all functional groups. We selected 13 natural and disturbed habitats on Kilimanjaro in Tanzania, distributed over an elevational range of 3700 m. In each habitat type, we set up three to six modified Whittaker plots. We recorded all plant species in 64 plots and 640 subplots and described SAR using the power function. Area consistently modulated effects of elevation on plant species richness, partly effects of land use but not effects of plant canopy. Thus, area needs to be taken into account when studying elevational plant species richness patterns. In contrast to temperate regions open and forest habitats did not differ in SAR, probably due to a distinct vertical vegetation zonation in tropical forests. Therefore, it is important to consider all vegetation layers including epiphytes when studying SAR in highly structured tropical regions.


Introduction
A fundamental goal in ecology is describing and explaining the distribution and abundance of species on earth (Gaston 2000). One of the few very well documented species richness patterns is the species-area relationship (SAR) (Rosenzweig 1995). It describes the overall diversity and turnover between small and larger areas and thus allows quantifying and comparing the spatial distribution of species (Connor and McCoy 1979;Drakare et al. 2006;Cencini et al. 2012). It can also be useful to predict species and habitat loss and a reduction of speciation in response to global environmental change (Rosenzweig 2001;Thomas et al. 2004).
When studying species-area relationships the perpetual observation is that species richness increases with area. This holds true for all continents (Storch et al. 2012), for all types of taxa (Lomolino 2000;Storch et al. 2012;Patino et al. 2014) and for all spatial scales (Crawley and Harral 2001).
There are neutral, ecological and sampling-related explanations for the constantly observed increase of species richness with area. The underlying idea of the neutral model is that greater areas have lower extinction and higher immigration and speciation rates and thus harbor a greater number of species than smaller areas (MacArthur and Wilson 1963;Hubbell 2001). The ecological explanations involve environmental heterogeneity, i.e. the fact that larger areas include more habitats, as well as spatial aggregation and segregation of species (Rosenzweig 1995;Crawley and Harral 2001). Finally, sampling-related explanations argue that in larger areas rare species are more likely recorded than in smaller areas due to the higher number of individual plants (Connor and McCoy 1979;Rosenzweig 1995).
The first mathematical description of the SAR was proposed nearly a hundred years ago: the power function model (Arrhenius 1921). Up to date, this model is the most frequently used and considered best to describe species-area curves (Drakare et al. 2006;Dengler 2009). It states that the number of species (S) equals a constant c multiplied by area (A) to the power of another constant z. The two parameters in this model, c and z, correspond to the intercept and slope of the linearized version of this function, respectively. While the biological interpretation of the intercept is the average species richness (alpha diversity) in one unit of area, the slope represents the rate of species accumulation (beta diversity) (MacArthur 1965;Connor and McCoy 1979).
According to the neutral model, the slope of the species-area curve only depends on geographical isolation of a region (MacArthur and Wilson 1967) and is expected to approximate 0.25 (Preston 1962). However, there is large spatial variation of slopes within regions (Williamson 2003;Manne et al. 2007), which cannot be explained by purely neutral models. A number of possible factors have been studied to explain this variation. Relevant plant specific factors include dispersal ability (Patino et al. 2014), functional group (Qiao et al. 2012), plant origin (Hulme 2008) and organism size (Azovsky 2002). But also local population size (Drakare et al. 2006), spatial aggregation and distribution of species and individuals (Connor and McCoy 1979;Plotkin et al. 2000;Tjorve et al. 2008), as well as habitat type and diversity (Connor and McCoy 1979;Drakare et al. 2006), latitude (Drakare et al. 2006), spatial scale (Crawley and Harral 2001) and methodology (Williamson 2003) can affect the rate of species accumulation. However, some factors driving SARs are still controversial and heavily debated (Harte et al. 2008;Sizling et al. 2011).
In a time of global change, it is increasingly important to study how human-induced changes affect species richness patterns. Land use and climate change, which can best be 1 3 studied along elevational gradients (Körner 2021), have been identified as the major drivers of biodiversity changes (Sala et al. 2000). When studying land use and elevational effects on species richness patterns, it is crucial to understand potential area effects. If area modulates the effects of land use and elevation, area needs to be accounted for when interpreting results (Connor and McCoy 1979). However, to our knowledge it has not been investigated yet how land use affects SAR, and also the effect of elevational gradients on SAR has hardly been studied nor were such studies undertaken in multilayered tropical forests.
One of the few factors interacting with area that is usually considered by botanists is organism size. For vegetation records, botanists tend to choose larger plots in forests than in open canopy habitats as mean individual size in forest plots is larger (Mueller-Dombois and Ellenberg 1974;Dierschke 1994). However, so far only a few studies tested whether open-and closed-canopy habitats systematically differ in SAR (e.g. Dolnik and Breuer 2008), but no in complex tropical ecosystems thus justifying different treatment of such habitats.
In this study, we investigated if area modulates the effect of land use, elevation and canopy on plant species richness on Kilimanjaro using modified Whittaker plots. Kilimanjaro with its altitudinal range of 5000 m has a long history of land use and harbors a great variety of habitats (Hemp 2005(Hemp , 2006a ranging from uni-layered low-stature to complex multilayered systems with a high upper canopy harboring the tallest trees of Africa . Therefore, it can be considered a model system to study effects of land use, habitat, elevation and canopy on SARs. We expected increasing land use intensity and elevation to decrease species richness and decelerate species accumulation due to more homogeneous habitat conditions and a higher number of annual species, i.e. species with higher dispersal ability. We expected a higher species accumulation rate for closed-canopy habitats since mean individual size is larger in such habitats and thus less individuals fit into a given area than in open-canopy habitats.

Study area
Kilimanjaro is located in Tanzania, 300 km south of the equator on the border to Kenya. It is a dormant volcano with a diameter of 90 km and the highest mountain of Africa (5895 m a.s.l.). Due to its large altitudinal range of over 5000 m, Kilimanjaro harbors a great variety of habitats and vegetation belts. There are six main vegetation zones: the savanna (700-1100 m a.s.l.), the sub-and lower montane zone (1100-2100 m a.s.l.), the middle montane zone (2100-2800 m a.s.l.), the upper montane zone (2800-3200 m a.s.l.), the subalpine zone (3200-4000 m a.s.l.) and the alpine zone (4000-4600 m a.s.l.) (Hemp 2006a(Hemp , 2006b). In the lower vegetation zones, anthropogenic land use has considerably reduced the natural vegetation. Generally, natural habitats below 1800 m are much reduced due to crop and fodder production while the ones above are affected by logging and burning. In the colline zone, natural savanna woodlands consisting of a grass layer and trees such as Acacia and Terminalia have been converted into maize fields for food production. In the sub-and lower montane zone, natural forest has been removed in most places in favor of commercial coffee plantations, mown grasslands and traditional agroforestry systems (Chagga homegardens, Hemp 2006c). In the middle montane zone, natural Ocotea usambarensis forests have been selectively logged until 1984 (Agrawala et al. 2003). In the 1 3 upper montane zone, fires disturb the natural Podocarpus latifolius forest, which is then recolonized by fire-tolerant Erica excelsa trees (Hemp and Beck 2001). Likewise, the natural Erica trimera forest in the sub alpine zone is disturbed by fires. Only the natural alpine Helichrysum vegetation is still unaffected by human activity besides gentle tourism (Hemp 2008).
The climate on Kilimanjaro is typically tropical equatorial and depends on elevation. Mean annual temperature decreases linearly from 24 °C in the savanna at 700 m a.s.l. to − 7 °C at the top, while annual precipitation shows a hump-shaped pattern, peaking at 2200 m in the montane forest belt (Hemp 2006a). Due to wind exposure from the Indian Ocean, the southern and eastern slopes are wetter than the northern and western slopes (Hemp 2006a).

Site selection and study design
We used the permanent study plots of the DFG-Research Unit "Kilimanjaro ecosystems under global change: Linking biodiversity, biotic interactions and biogeochemical ecosystem processes (KiLi)," which cover all main natural and human-modified habitat types across the elevational gradient on Mt. Kilimanjaro (https:// www. kilim anjaro. bioze ntrum. uni-wuerz burg. de/). The 64 study plots of the KiLi project cover the 13 most important natural and anthropogenically influenced habitat types on Kilimanjaro occurring in the six main vegetation zones on the southern and eastern slope, where vegetation zonation is stronger and the habitat range wider than on the northern and western slope due to the wetter conditions and the more pronounced climatic and elevational gradient (Hemp 2006a).
Each habitat type is represented by five replicate plots, with the exception of fire-disturbed and undisturbed Erica forests (3 and 4 plots) and coffee plantations and alpine vegetation with 6 plots each, distributed along an east-west and an elevational gradient within each vegetation zone (Table 1). The 64 plots were distributed over five elevational transects on the southern slope of Mt. Kilimanjaro, covering an elevational gradient from 870 to 4550 m a.s.l. Minimum distance between transects is 12.2 km on the foothills and 4.6 km in the alpine zone, and minimum distance between plots is approximately 300 m. For each site, plot position (latitude, longitude and elevation) was recorded with GPS. For more details see e.g. Rutten et al. (2015). Moreover, for each plot, mean annual temperature was obtained from in-situ measures by temperature loggers (data logger DK320, Driesen and Kern GmbH, Bad Bramstedt, Germany) and mean annual precipitation was modeled using 1 3 long-term observations based on a 15-year dataset from a network rain gauges distributed across the whole mountain (Hemp 2006a;Appelhans et al. 2016).
The study was set up between February 2011 and June 2014 in a modified Whittaker plot design. At each of 64 sites, we established one main plot of 50 × 20 m (= 1000 m 2 ) and nine subplots of three sizes (1, 10 and 100 m 2) within each main plot (Fig. 1). All subplots were nested in the main plot but not in each other. This independent vegetation sampling design allows assessing plant species richness at multiple scales, which is required for investigating species-area relationships, without the disadvantages associated with other sampling designs such as nested vegetation plot series (Stohlgren et al. 1995;Dolnik and Breuer 2008). In previous studies using the identical plots we found no spatial autocorrelation between the 64 plots by calculating Moran's I-values (Ferger et al. 2014(Ferger et al. , 2017.

Measurements and data analysis
We recorded all vascular plant species in the main plots and all subplots between February 2011 and June 2014 following the method of Braun-Blanquet (1964). We separately recorded species in the following layers: herb layer (< 1 m tall), shrub layer (woody species, 1-10 m tall), tree layer (woody species, > 10 m tall). Additionally, tree and shrub lianas, defined as climbers reaching into the shrub or tree layer, epiphytes on living plants and epiphytes on deadwood were included as separate layers. If plots were newly disturbed (e.g. maize fields were ploughed or vegetation newly cut in grasslands), the vegetation not fully developed during dry seasons and species not determinable, plots were revisited up to four times.
To describe the species-area relationship, we used the power function (Arrhenius 1921;Dengler 2009). For each of the 64 SARs, we regressed the log-transformed species richness on log-transformed area. Additionally, we calculated the regression coefficients for annuals and perennials (life duration) separately.
We analyzed the effect of land use on the regression parameter estimates of the SAR (intercept and slope) using analysis of variance with land use as fixed factor. Because altitudinal zones differed in land use types, we analyzed the effect of land use in each altitudinal zone separately. To test the effect of life duration, we additionally included life duration and its interaction with land-use as fixed factors in the model. We analyzed the effect of elevation on intercept and slope of the SAR in natural habitats using multiple linear regression with backward model selection with the linear and quadratic terms of elevation as initial explanatory variables. The best model was then selected using stepwise backward selection based on AIC and further stepwise reduction based on p-values until only significant variables were left. To test whether the parameters of SAR of annuals and perennials differently changed along elevation, we included life duration and its interactions with the linear and quadratic terms of elevation as initial explanatory variables.
To study whether open and closed canopy sites differ in parameters of the SAR, we measured the leaf area index (LAI) in each main plot and all subplots at ground level using a plant canopy analyzer in combination with a remote 'above canopy' sensor (LAI 2200, LI-COR Bioscience USA, 2011). The 'above canopy' reading was measured at ground level in an open area or forest gap as close to the site as possible and LAI values were calculated using the program FV-2200 (LI-COR Bioscience USA, 2011). Data of these measurements are given in Rutten et al. (2015), the basic theory behind LAI is described in LI-COR (2012). We then calculated the mean LAI per site and classified the canopy of sites with LAI < 2 and LAI > 2 as open and closed, respectively. This threshold value separated grasslands, savanna vegetation, and alpine Helichysum scrub vegetation from all other vegetation types with a tree layer on our 64 plots on Kilimanjaro. Then, we tested the significance of plant canopy on the intercept and slope of SAR using a linear model with plant canopy as explanatory variable. Analyses with LAI as continuous variable yielded closely similar results.

Results
In the 64 sites on Kilimanjaro, we recorded 1014 vascular plant species, 553 of which occurred in natural, and 786 in disturbed habitats. Species richness increased from 336 in the savanna to 568 in the lower montane zone and then decreased to 42 in the alpine zone. The high species richness in the lower montane zone was caused by the three anthropogenic land use types, which harbored 90% of unique species, while the natural vegetation consisted of 71% of unique species. In the other vegetation zones, disturbed habitats contained 30-48% of unique species, while natural habitats harbored 15 (montane)-56% (savanna).
Values of the regression coefficient (adjusted R 2 ) of the species-area curves calculated as power function for total species richness for all Whittaker plots ranged from 0.41 (F 1,8 = 7.3, p < 0.05) to 0.96 (F 1,8 = 201.99, p < 0.001) with a mean of 0.82. The intercepts of the SAR ranged from 0.84 to 2.78, while the slopes ranged from 0.14 to 0.39.
Land use and disturbance affected intercept and slope of SAR in two vegetation zones (Fig. 2). The intercept of SAR in the subalpine zone was smaller in natural than in disturbed habitats, while in the lower montane belt it did not differ between natural and disturbed habitats, but was smaller in agroforestry systems than in grasslands. The effect of land use on the intercept was modulated by life duration in three zones (Table S1, Fig.  S1). In the savanna and the subalpine zone, the slope of SAR was steeper in natural than in disturbed habitats (Fig. 2). Perennial plants generally showed a steeper slope of SAR compared to annuals, and this effect was modulated by disturbance only in the subalpine zone (Table S1, Fig. S1).
Elevation significantly affected both intercept and slope of the SAR in natural habitats ( Table 2). The intercept showed a quadratic decrease with elevation, while the slope decreased linearly (Fig. 3a, b). Annual and perennial plants showed different change in intercept, but not in slope with elevation (Fig. 3c, d). The models with elevation, life duration and their interactions explained most of the data variance for both intercept (R 2 = 0.9) and slope of SAR as response variable (R 2 = 0.82).

Values of SAR coefficients varied widely among plots
Concurring with previous observations, we found large spatial variation in SAR coefficient values (Williamson 2003;Manne et al. 2007). Species richness per unit of area (intercept of SAR) quadratically decreased with elevation, i.e. geographic isolation, suggesting at least partial suitability of the neutral model. However, the variation in rate of species accumulation (slope of SAR) cannot be explained by a neutral model as we would expect similar slopes in all plots. Also factors that were kept constant in our study like study region, spatial scale or methodology cannot explain this variation. Therefore, plant or environment specific factors must be responsible for the large spatial variation of SAR coefficients.

Land use does not systematically affect the species-area relationship
Land use and disturbance did not affect the rate of species accumulation consistently. While in three out of five vegetation zones, land use had no effect on species accumulation rate, it had a negative effect in the savanna and the subalpine zone. Life duration of plants is unlikely to have affected this result, as only in the subalpine zone annuals and perennials differed in response to land use. A more parsimonious explanation is habitat heterogeneity (Rosenzweig 1992). Only in the savanna and subalpine zones conditions were more Table 2 Results of multiple linear regressions with backward model selection to identify variables best explaining the parameters of the species-area curve (power function) in natural habitats on Kilimanjaro calculated for all species and annual and perennial species separately Shown are t-and F-values for continuous and categorical variables, respectively, as well as levels of significance (***p < 0.001, **p < 0.01, *p < 0.05) df Intercept Slope t-/F-value t-/F-value heterogeneous in natural than disturbed habitats. An indication for this is that in natural habitats of the subalpine zone species richness was lower in smaller, but not in larger areas compared with disturbed habitats. Thus, species must have been more heterogeneously distributed in natural than disturbed habitats, resulting in a steeper SAR. The explanation is probably the leveling effect of strong disturbance, in case of the subalpine vegetation fire and in case of the savanna ploughing. Every year fires occur in the subalpine ericaceous belt (Hemp 2005) of Kilimanjaro, the last one in October 2020 destroying over 90 km 2 of bushland. This creates a vegetation of similar age and structure. The annual ploughing in the savanna has a similar effect. In contrast to our expectations, land use and disturbance did not negatively affect species richness per unit area in any vegetation zone. In fact, in the subalpine belt, species richness was even smaller in natural than in disturbed habitats, though in larger plots the difference was levelled out. There are two possible explanations for this observation. First, disturbance may not have been severe enough to negatively affect species richness. Low to moderate disturbance is known to have a neutral or even positive effect on species richness (Grime 1973;Boch et al. 2013). As in many vegetation zones on Mt. Kilimanjaro disturbance intensity is low or has occurred many years ago (Agrawala et al. 2003;Hemp Fig. 3 Relationship between elevation and the two parameters of the species-area curve (power function) calculated for all species (a, b), and annual and perennial species (c, d) in natural habitats on Kilimanjaro shown as multiple linear regression after model selection with elevation and quadratic elevation (a, b) and life duration (c, d) as initial explanatory variables. Statistical results of the multiple linear regressions after model selection are shown in Table 2 1 3 2006b), it may have led to similar species numbers in natural and disturbed habitats. This holds true in particular in the middle montane forest belt, where disturbance intensity is low or has occurred many years ago (Agrawala et al. 2003;Hemp 2006b), resulting not only in similar species numbers but also in similar structure and biomass of natural and disturbed habitats (Rutten et al. 2015). Second, disturbance may have created new habitats, resulting in change in species composition instead of richness. Indeed, the number of unique species was high in both natural and disturbed habitats of each vegetation zone. This species turnover may have outbalanced the species loss due to disturbance of natural habitats.

Elevation changes the species-area relationship in natural vegetation
Species richness per unit area decreased quadratically with increasing elevation. The peak in the montane zone was possibly caused by the strong vertical vegetation structure and thus occurrence of epiphytes and vines, and due to spillover effects from upper and lower vegetation zones (Rosenzweig 1992). In contrast to the perennial species richness, which followed the same pattern as overall species richness, annual richness was highest in the lower vegetation zones, where frequent droughts inhibit longevity of herbaceous species.
Species accumulation rate decreased linearly with elevation. Although we expected this pattern, it could not be attributed to the changing proportion of annual and perennial species with elevation. Probably, habitat heterogeneity decreased with elevation, leading to a decrease in species accumulation rate. Savanna vegetation is known to be extremely heterogeneous (Solbrig et al. 1996), while at least on Kilimanjaro, the vegetation in the alpine zone is less heterogeneous. A decrease of species accumulation rate with elevation has been found before (Qiao et al. 2012). Whether this pattern can generally be observed, remains to be studied.

Open and closed canopy plots do not differ in species-area relationship
Vegetation studies in temperate regions emphasize higher homogeneity and slower species accumulation rate of open than closed habitats (Mueller-Dombois and Ellenberg 1974;Barkman 1989;Dolnik and Breuer 2008). However, plant canopy on Kilimanjaro did not affect any of the parameters of the SAR. The most likely explanation for this is the vertical vegetational zonation. In tropical forests, tree species are often covered with epiphytes and lianas, leading to a layer structure, while in open-canopy habitats there is usually only one vegetation layer. In certain forest types of Kilimanjaro more species grow in the epiphyte layer than on the forest floor (Hemp 2006b). Therefore, the number of individuals per unit of area is probably similar despite the occurrence of large space demanding tree individuals. This is likely to lead to similar species numbers and accumulation rate in open-and closed-canopy habitats.

Conclusion
Area systematically modulates the effect of elevation on species richness. In higher elevations, a smaller plot size is needed to capture a large proportion of total species richness than in lower elevations, possibly due to lower habitat heterogeneity and lower total species number. In contrast, area does not systematically modulate the effect of land use and disturbance on species richness, although it may in special cases, and this result is equally valid for annual and perennial species. In particular, in the subalpine zone species richness of undisturbed forests was lower in smaller, but not in larger areas compared with disturbed habitats. Here, fires as the main disturbance that occur regularly in the ericaceous belt of Kilimanjaro have a strong levelling effect, reducing heterogeneity and variety of microhabitats of large areas.
There is no indication that open-and closed-canopy habitats differ in SARs, probably due to a distinct vertical vegetational zonation in tropical forests. Therefore, a 50 × 20 m plot seems to be suitable to record both forest and grassland vegetation in afromontane regions.
Our study shows the importance of taking area into account when conducting ecological experiments. Especially studies looking at elevational effects on species richness should carefully choose plot size to avoid confounding effects on their results. Furthermore, for studies of diversity patterns in highly structured tropical montane forests it is necessary to include all functional groups (in particular the epiphyte layer) and not only the trees as in most studies. Despite more than a hundred years of studying species-area relationships, the pattern and mechanism behind it is still not fully understood.