Predictors of elevational biodiversity gradients change from single taxa to the multi-taxa community level

Abstract

The factors determining gradients of biodiversity are a fundamental yet unresolved topic in ecology. While diversity gradients have been analysed for numerous single taxa, progress towards general explanatory models has been hampered by limitations in the phylogenetic coverage of past studies. By parallel sampling of 25 major plant and animal taxa along a 3.7 km elevational gradient on Mt. Kilimanjaro, we quantify cross-taxon consensus in diversity gradients and evaluate predictors of diversity from single taxa to a multi-taxa community level. While single taxa show complex distribution patterns and respond to different environmental factors, scaling up diversity to the community level leads to an unambiguous support for temperature as the main predictor of species richness in both plants and animals. Our findings illuminate the influence of taxonomic coverage for models of diversity gradients and point to the importance of temperature for diversification and species coexistence in plant and animal communities.

ARTICLE

Received 12 Jan 2016 | Accepted 28 Oct 2016 | Published 22 Dec 2016

DOI: 10.1038/ncomms13736 OPEN

Predictors of elevational biodiversity gradients change from single taxa to the multi-taxa community level

Marcell K. Peters1, Andreas Hemp2, Tim Appelhans3, Christina Behler4, Alice Classen1, Florian Detsch3, Andreas Ensslin5, Stefan W. Ferger6, Sara B. Frederiksen7,8, Friederike Gebert1, Michael Haas7, Maria Helbig-Bonitz4, Claudia Hemp1, William J. Kindeketa1,9, Ephraim Mwangomo3,10, Christine Ngereza7,11, Insa Otte3, Juliane Rder7, Gemma Rutten5, David Schellenberger Costa12, Joseph Tardanico1, Giulia Zancolli1,13, Jrgen Deckert14, Connal D. Eardley15,16, Ralph S. Peters17, Mark-Oliver Rdel14, Matthias Schleuning6, Axel Ssymank18, Victor Kakengi19, Jie Zhang1, Katrin Bhning-Gaese6,20, Roland Brandl7, Elisabeth K.V. Kalko4,21,z, Michael Kleyer12,

Thomas Nauss3, Marco Tschapka4,21, Markus Fischer5,6 & Ingolf Steffan-Dewenter1

The factors determining gradients of biodiversity are a fundamental yet unresolved topic in ecology. While diversity gradients have been analysed for numerous single taxa, progress towards general explanatory models has been hampered by limitations in the phylogenetic coverage of past studies. By parallel sampling of 25 major plant and animal taxa along a3.7 km elevational gradient on Mt. Kilimanjaro, we quantify cross-taxon consensus in diversity gradients and evaluate predictors of diversity from single taxa to a multi-taxa community level. While single taxa show complex distribution patterns and respond to different environmental factors, scaling up diversity to the community level leads to an unambiguous support for temperature as the main predictor of species richness in both plants and animals. Our ndings illuminate the inuence of taxonomic coverage for models of diversity gradients and point to the importance of temperature for diversication and species coexistence in plant and animal communities.

1 Department of Animal Ecology and Tropical Biology, Biocenter, University of Wrzburg, Am Hubland, Wrzburg 97074, Germany. 2 Department of Plant Systematics, University of Bayreuth, Universitatsstra 30, Bayreuth 95440, Germany. 3 Environmental Informatics, Faculty of Geography, University of

Marburg, Deutschhausstra 12, Marburg 35032, Germany. 4 Institute for Evolutionary Ecology and Conservation Genomics, University of Ulm, Albert-Einstein-Allee 11, Ulm 89069, Germany. 5 Institute of Plant Sciences, University of Bern, Altenbergrain 21, Bern 3013, Switzerland. 6 Senckenberg Biodiversity and Climate Research Centre (BiK-F), Senckenberganlage 25, Frankfurt am Main 60325, Germany. 7 Department of Ecology, Animal Ecology, University of Marburg, Karl-von-Frisch-Stra 8, Marburg 35032, Germany. 8 Zoological Museum, Natural History Museum of Denmark, University of Copenhagen,

Universitetsparken 15, Copenhagen DK-2100, Denmark. 9 Tanzania Commission for Science and Technology, Department of Life Sciences, Ally Hassan Mwinyi Road, PO Box 4302, Dar es Salaam, Tanzania. 10 Mount Kilimanjaro National Park, PO Box 96, Marangu, Moshi, Tanzania. 11 National Museum of Tanzania, Shaaban Robert Street, Dar es Salaam, Tanzania. 12 Landscape Ecology Group, Institute of Biology and Environmental Sciences, University Oldenburg, Oldenburg 26111, Germany. 13 Molecular Ecology and Fisheries Genetics Lab, School of Biological Sciences, Environment Centre Wales, Bangor University, Bangor LL57 2UW, UK. 14 Museum fr Naturkunde, Leibniz Institute for Evolution and Biodiversity Science, Invalidenstra 43, Berlin 10115,

Germany. 15 Agricultural Research CouncilPlant Protection Research: Plant Health and Protection, Private Bag X134, Queenswood, Pretoria 0121, South Africa. 16 School of Life Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa. 17 Zoological Research Museum Alexander Koenig, Department Arthropoda, Adenauerallee 160, Bonn 53113, Germany. 18 Falkenweg 6, Wachtberg 53343, Germany. 19 Tanzania Wildlife Research Institute, PO Box 661, Arusha, Tanzania. 20 Institute for Ecology, Evolution and Diversity, Goethe University Frankfurt, Biologicum, Max-von-Laue-Stra 13, Frankfurt am Main 60439, Germany. 21 Smithsonian Tropical Research Institute, PO Box 0843-03092, Balboa Ancn, Republica de Panam.

Correspondence and requests for materials should be addressed to M.K.P. (email: mailto:[email protected]

Web End [email protected] ).

zDeceased 26 September 2011

NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736

The search for the primary factors that determine the distribution of biodiversity on earth has challenged naturalists for more than two centuries13. The main

hypotheses for explaining broad-scale diversity gradients are: (1) the temperature hypothesis, relating higher species richness to higher rates of biotic processes and interactions, or higher evolutionary diversication rates46; (2) the water availability hypothesis, focusing on direct or indirect (via effects on net primary productivity (NPP)) constraints of water for the maintenance of biodiversity7,8; (3) the productivity hypothesis, emphasizing the positive effect of resources on population persistence and species coexistence911; (4) the area hypothesis, assuming greater opportunities for the maintenance of richness and speciation in larger areas of land (or sea)1214; (5) the geometric constraints hypothesis, focusing on spatial constraints for the distribution of species ranges, predicting geographic gradients in species richness even in the absence of environmental drivers15,16; and (6) the plant diversity hypothesis, relating consumer richness to the species richness of plant resources17,18. While numerous studies on the diversity of individual taxa along different environmental gradients found support for one or the other hypothesis, comparative cross-taxon studies along the same environmental gradients are scarce, impeding consensus on the determinants of broad-scale diversity gradients2,19.

A persistent shortcoming in published gradient studies is that trends in species richness were mostly studied for organismic groups at a narrow taxonomic level (for example, for taxa such as ferns, ants or birds) but not for the multi-group communities of animals and plants which coexist and share resources in ecosystems. However, as species and clades tend to retain their niches and related ecological traits over evolutionary time (niche conservatism), narrowly dened taxa can be more strongly constrained to specic biotic and abiotic conditions along larger environmental gradients than communities composed of multiple, phylogenetically largely independent taxa20,21. Patterns of species richness of a specic taxon may, therefore, be idiosyncratic rather than congruent to patterns in other taxa or to those of the community as a whole. Along large environmental gradients, taxonomically diverse communities have larger opportunities for diversication and the use of niche space than a single taxon. For example, while bee species are typically restricted to warm environments and their species richness peaks in hot climates, syrphid ies have their highest density at considerably cooler temperatures, such that the thermal niche space used by pollinators in total is much larger than the thermal niches of the individual taxa22,23. Broadening the taxonomic scope may therefore lead to changes in the importance of some potential drivers of diversity over others. In particular, energy variables like temperature and NPP have been predicted to increase in importance for more broadly covered species communities6,11,24. However, to our best knowledge, the inuence of taxonomic coverage for the extrapolation of broad-scale drivers of diversity has hitherto not been tested.

We used a novel, multi-taxa perspective to study patterns and drivers of species richness along a 3.7 km elevational gradient (8714,550 m above sea level (a.s.l.)) of natural habitat on the southern slopes of Mt. Kilimanjaro (Supplementary Fig. 1). Tropical mountains with a dry base, such as Mt. Kilimanjaro, are particularly suitable to test large-scale drivers of biodiversity because, in contrast to most terrestrial latitudinal gradients, their temperature and primary productivity gradients are largely uncorrelated (Supplementary Fig. 2). We assessed species richness of eight vascular plant taxa (that is, all vascular plants), and 16 major animal taxa synchronously on the same study sites. The studied taxa were phylogenetically highly diverse, covered

multiple trophic levels and represented signicant proportions of diversity found in terrestrial plant and animal communities25 (Fig. 1a).

We, rstly, analysed congruence in the patterns and predictors of elevational diversity among single taxonomic groups, that is, for taxonomic groups which are typically studied in macroecology. Secondly, we analysed how increasing the taxonomic coverage towards the community level modies patterns of elevational diversity and the support for hypotheses explaining them. Here we show that disparate factors drive the distributions of individual taxa. However, scaling up diversity to the community level leads to a monotonic decline of species richness with elevation and provides a strong support for temperature as the major predictor of plant and animal species richness. Our study reveals high variation in diversity gradients of narrowly dened taxa and underscores the signicance of taxonomically complex data sets for understanding the broad-scale drivers of biodiversity.

ResultsElevational species richness of single taxa. We detected a variety of elevational patterns in species richness among plant and animal taxa (Fig. 1b). Approximately half of the plant and animal taxa exhibited monotonic declines in species richness with increasing elevation (Fig. 1b and Supplementary Figs 34). The other half had hump-shaped or more complex (for example, bimodal) distribution patterns. While the plant censuses were complete and the species richness in each study site was accurate, incomplete sampling, a regular problem in studies of tropical communities18,26, could have had biased patterns of elevational diversity in animals (Supplementary Table 1). Estimates of sample coverage, a measure of sampling completeness27, of most animal taxa were high (mean sample coverage 40.80, s.d. o0.17) but values in four taxa (in spiders, parasitoid wasps, other aculeate Hymenoptera and moths) were moderate to low (mean sampling coverage: 0.540.72, s.d. 0.210.38). This suggested that further sampling could have

signicantly increased estimates of species richness (Supplementary Table 1). As sample coverage was signicantly positively correlated with elevation in some animal taxa (Collembola, true bugs, bees, other aculeate Hymenoptera; Pearsons product moment correlation, r 0.370.85, Po0.05) we checked

if incomplete sampling could have biased elevational diversity patterns. We calculated for all animal taxa estimates of asymptotic species richness with the non-parametric Chao1 index27,28 and compared diversity trends to those detected in the observed data. Patterns of elevational diversity based on estimated asymptotic species richness were highly similar to those derived from the observed data (means.d. of Pearsons r 0.960.09,

Supplementary Table 2).

In the next step we quantied the support for different hypotheses to explain species richness data at the level of single taxa using information-theory based multi-model averaging29,30. As potential predictors of diversity we used temperature (in animals: mean annual temperature, MAT; and in plants: mean minimum temperature, MMT31), mean annual precipitation (MAP, water availability hypothesis), net primary productivity (NPP; productivity hypothesis), land area within an elevation belt of 300 m above and below the level of the study sites (area hypothesis), predictions of a mid-domain effect (MDE) resulting from multiple random arrangements of species ranges within the elevation boundaries on Mt. Kilimanjaro (geometric constraints hypothesis) and in case of species richness of animals, plant species richness (PSR, plant diversity hypothesis) (Supplementary Fig. 2).

Species richness was best predicted by sets of variables which differed among the individual taxonomic groups (Table 1 and

2 NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736 ARTICLE

a

1,000 2,000 3,000 4,000 1,000 2,000 3,000 4,000 1,000 2,000 3,000 4,000 1,000 2,000 3,000 4,000

Elevation (m a.s.l.)

1,000 2,000 3,000 4,000

Figure 1 | Patterns of elevational species richness of single taxa. (a) Phylogenetic distribution of studied taxonomic groups among major terrestrial plant and animal lineages. (b) Patterns of elevational species richness for vascular plants (orange) and animals (blue). Dots represent original measurements on study sites (plants: N 30; animals: N 30, except for Collembola (N 29), ground-dwelling beetles (N 29) and amphibians (N 17)). Trend lines were

calculated using generalized additive models; all trends were signicant (Po0.05). For trend lines of additional plant groups with low numbers of species see Supplementary Fig. 3. The plant and animal images used in the gure are licensed for use in the Public Domain without copyright, except for the images used for monocots (F. & K. Starr, modied), bats (O. Peles and Y. Wong), Collembola, ground-dwelling beetles (B. Lang), Orthoptera, other aculeate wasps, bees (M. Menchetti), gastropods (G. Monger) that are licenced under a Creative Commons Attribution 3.0 Unported licence (https://creativecommons.org/licenses/by/3.0

Web End =https:// https://creativecommons.org/licenses/by/3.0

Web End =creativecommons.org/licenses/by/3.0 https://creativecommons.org/licenses/by/3.0

Web End =). Oliver Niehuis kindly provided the permission to use the image of the parasitoid wasp.

Supplementary Table 3). Temperature (MMT and MAT) was an important predictor of species richness in most taxa, having signicant positive effects in 11 plant and animal taxa and a negative effect in one taxon. NPP was signicantly correlated to species richness in eight plant and animal groups. However, positive and negative effects of primary productivity were equally common. Six taxa exhibited signicant positive, and three taxa negative relationships to MAP. Richness of phylogenetically basal plant lineages (Lycopodiopsida, ferns, conifers) tended

to increase, while richness of more modern lineages (for example, monocots) tended to decrease at elevated levels of MAP (Table 1). Signicant positive effects of land area and geometric constraints (that is, the MDE) were only observed in two taxa (millipedes and Orthoptera). In moths and gastropods, two animal groups with completely or largely herbivorous diet, signicant, positive correlations with PSR were detected. Other animal taxa did not strongly respond to the variation in plant richness.

NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736

Table 1 | Synthesis models explaining species richness of plant (a) and animal taxa (b) derived by multi-model averaging.

Conditional standardized estimates Taxon #models MMT NPP MAP Area MDE

0.32

0.63

11

64

Ferns 0.40 0.62 0.42 Magnoliids 9 8 0.34 0.21 0.72 0.28 0.11

Monocots 106 6 0.92 0.55 0.40 0.11 0.28 Other eudicots 26 4 0.15 0.80 0.95 0.08 0.19

Rosids 158 8 0.88 0.28 0.18 0.03 0.31 Asterids 189 6 0.55 0.05 0.09 0.16 0.05 All vascular plants 557 10 1.08 0.22 0.17 0.42 0.19

Conditional standardized estimates

Taxon #species*

#species*

#models MAT NPP MAP Area MDE PSR

Gastropods 45 7 0.23 0.01 0.03 0.27 0.77 0.39 Millipedes 14 6 1.37 0.57 0.57

0.51

1.53 0.73

0.23 0.310.56

0.34

SpidersCollembola 23 11 1.03 1.04 0.48 1.03 0.75 0.44 Orthoptera 114 9 0.57 0.05 0.30 0.74 0.01 0.13 True bugs 44 16 0.81 0.20 0.19 0.21 0.20 ~0.00 Parasitoid wasps 223 13 0.36 0.45 0.50 0.47 0.03 0.19 Grounddwelling ants 38 20 1.00 0.33 0.21 0.47 0.17 0.24

6

52

Bees 1.08 0.42 0.09 0.05 0.18 0.11 Other aculeate Hym. 45 12 1.13 0.57 0.14 0.56 0.04 0.27 Dung beetles 56 18 1.09 0.28 0.22 0.62 0.01 0.36 Grounddwelling beetles 122 27 0.41 0.58 0.04 0.25

0.79

0.34 0.37

0.02

52

30

0.01

0.59

Moths 0.38 0.54 Hoverflies 19 13 0.02 0.04 0.63 0.37 0.34 0.02

7

122

Birds 0.87 0.78 0.07 0.66 0.68 0.29 Aerial insectivorous bats 16 15 0.93 0.27 0.13 0.25 ~0.00 0.22 All animals 1184 6 0.96 0.01 0.02 0.02 0.03 0.06

Shown are standardized parameter estimates for all predictor variables derived from weighted averaging of parameter estimates over best-t models. Colours indicate signicant (Po0.05) positive (blue) or negative (red) effects from multi-model averaging analyses. Results for two additional plant groups with low numbers of species are presented in Supplementary Table 3.

*Total number of detected species/morphospecies for each taxon.wNumber of best-t models (DAICo4) used for inference on parameter estimates and variable importance.

wwStandardized estimates (standardized beta) over all best-t models including the respective predictor variable.

Other aculeate Hym., other aculeate hymenoptera; MAT, mean annual temperature; MAP, mean annual precipitation; MDE, mid-domain effect prediction; MMT, mean minimum temperature; NPP, net primary productivity; PSR, plant species richness.

199

16

1.08

0.26

Elevational species richness of multi-taxa communities. Enhancing the taxonomic coverage of the data from one to eight (in plants) or 16 taxa (in animals) increased the explained variation of statistical models (that is, of species richness as a function of elevation) and lead to an increasing linearization of elevational species richness patterns in both plants and animals (Fig. 2). At the full multi-taxa community level (that is, species richness calculated by summing species richness values of all studied plant or animal taxa), species richness declined linearly with elevation in plants and quasi-linearly in animals (Fig. 3) with high levels of explained variation (generalized additive models; plants: explained deviance 77%, estimated degrees of freedom (edf) 1, F 94.4,

Po0.001; animals: explained deviance 92%, edf 2.34, F 105,

Po0.001). The decline of biodiversity with elevation was also evident but less pronounced when analysing richness of higher taxonomic levels (Fig. 3b,d): the number of orders and families (the latter only in plants) remained relatively constant up to mid-elevations and strongly declined at higher elevations.

Increasing the taxonomic coverage of the data from single taxa to the multi-taxa community level revealed an increasing

importance of MMT and MAT as the main predictors of diversity, while the support for NPP, area and other predictor variables declined (Fig. 4). At the highest level of taxonomic coverage, MMT and MAT were the only signicant predictors of species richness (Table 1), explaining 79 and 94% of the total variation in plant and animal species richness, respectively. Community level results were highly robust against the exclusion of single or multiple taxa (Fig. 4) and to differences in sampling intensity among taxa (Supplementary Fig. 5).

To document the relationships among predictor variables and to disentangle direct from indirect effects we calculated path models (Fig. 5). As suggested by analyses along other broad-scale gradients32, precipitation and temperature were strongly linked to the NPP of ecosystems. The direct effects of temperature on PSR were much stronger than the indirect effects via temperatures positive inuence on NPP. Similarly, in path analyses including animal species richness as the nal endogenous variable, the direct effect of temperature on animal richness was much stronger than the indirect effects mediated via its positive inuences on NPP and PSR.

4 NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736 ARTICLE

1.0

a

1.0

c

Hump-shaped or more complex

Exponential decline

Linear

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

Explained deviance

Complexity of pattern

0.0

0.0

1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8

1.0

1.0

b

d

Hump-shaped or more complex

Linear

Exponential decline

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0.0

1 13 15

3 5 7 9 11

1 13 15

3 5 7 9

11

Number of taxa included

Number of taxa included

Figure 2 | Elevational species richness patterns with increasing taxonomic coverage. (a,b) The explained deviance of generalized additive models increased with increasing taxonomic coverage of plant (orange: a,c) and animal (blue: b,d) communities. (c,d) While single taxonomic groups showed a variation of elevational species richness patterns (that is, linear decline, exponential decline or hump-shaped distributions) increasing the taxonomic coverage unambiguously led to patterns of linear decline in both plants and animals. In individual box-and-whisker-plots, bold lines indicate the median, boxes the interquartile range. Whiskers extend to the maximum and minimum values but end at 1.5 the interquartile range. More extreme data are plotted as single

dots.

1,000 2,000 3,000 4,000

a c

100

200

150

Elevation (m a.s.l.)

Species richness

80

60

100

40

20

50

0

0

1,000 2,000 3,000 4,000 1,000 2,000 3,000 4,000

b d

200

150

100

Taxon richness

80

60

100

40

20

50

0

0

1,000 2,000 3,000 4,000

Figure 3 | Elevational species richness at the community level. Species richness of vascular plants (a) and animals (c) along the elevational gradient of Mt. Kilimanjaro. Lower panels show trend lines for the number of species, families and orders of vascular plants (b) and animals (d) along the elevational gradient. All trend lines were calculated using generalized additive models (N 30 and N 29 in plants and animals, respectively; all

trends were signicant at Po0.001). Please see Fig. 1 for credits to the authors of the original plant and animal images.

DiscussionDespite decades of ecological research there is little consensus concerning the determinants of diversity gradients13,19. Our study of species richness patterns of multiple plant and animal taxa along the same environmental gradient on Mt. Kilimanjaro,

Tanzania, revealed that there is no general model for explaining diversity gradients across separately analysed taxonomic groups. Instead, patterns and predictors of elevational diversity appear to be idiosyncratic rather than uniform among taxa as they depend on taxon-specic resource requirements and adaptations to the environment. However, we found that the importance of potential drivers of diversity depends on the taxonomic coverage of analyses with taxonomic upscaling unambiguously revealing temperature as the single, well-supported predictor of species richness for taxonomically broad communities. Our results suggest multi-taxa community diversity as a new approach to develop more general models to explain trends in biodiversity: while past studies focused on congruence across different taxonomic groups, a predictor may be seen as having higher generality if it explains trends in diversity for taxonomically more inclusive species communities. Even though former studies already assessed biodiversity of several plant or animal taxa along the same elevational gradient3336, our study is unprecedented in terms of taxonomic coverage allowing generalizations and conclusions which could not be achieved with former data sets.

A monotonic decrease of species richness with elevation, as found at the highest level of taxonomic coverage, appears to contrast with former meta-analyses37,38 nding unimodal patterns of elevational diversity in a majority of studies. However, all of the individual studies on animals and many of those on plants reanalysed in meta-analyses are based on taxonomically restricted data sets, that is, these meta-analyses are based on taxonomic groups as used in our single taxa analyses. Analyses at this taxonomic level revealed also for the Kilimanjaro region a high percentage of unimodal patterns (in 44% of all tested taxa diversity patterns were unimodal). However, with increasing taxonomic coverage of plant and animal communities, unimodal patterns disappeared and the elevational diversity patterns shifted towards a monotonic decline of richness. This unication of diversity patterns with increasing

NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736

a

MMT MAP

NPP Area MDE

1.0

100

0.8

80

Variable

importance

Highest variable

importance (%)

0.6

60

0.4

40

0.2

20

0.0

0

b

1.5

1.0

100

0.5

60

Std beta

0.0

0

Highest std beta

(% of samples)

0.5

1.0

1.5

1 3 7 5

1 3 7 5

1 3 7 5

1 3 7 5

1 3 7

MAT MAP

NPP Area MDE PSR

Number of plant taxa included

5

c

1.0

100

0.8

Highest variable

importance (%)

80

0.6

Variable

importance

60

0.4

40

0.2

20

0.0

0

d

1.5

1.0

100

0.5

60

Std beta

0.0

0

0.5

1.0

1.5

1 5 15 10

1 5 15

10

1 5 15 10

1 5 15 10

1 5 15 10

1 5 15

Highest std beta

(% of samples)

10

Number of animal taxa included

Figure 4 | Statistical support for predictors of species richness in relationship to taxonomic coverage. Box-and-whisker-plots show the variation in measures of variable importance (a,c) and standardized beta (b,d) in relationship to the taxonomic coverage of plant (a,b) and animal communities (c,d) (values correspond to the left y axis). Bold curved lines give the percentage of all possible taxa combinations in which a variable had the highest variable importance or standardized beta value (extending from 0 to 100%, right y axis). Variable importance is dened as the sum of the Akaike weights of all best-t models which include the respective predictor variable. Standardized beta values are standardized parameter estimates derived from conditional weighted averaging of parameter estimates over best-t models. MMT, mean minimum temperature; MAT, mean annual temperature; NPP, net primary productivity; MAP, mean annual precipitation; MDE, mid-domain effect prediction; PSR, plant species richness.

taxonomic scale was also observed within some clades, for example, in beetles and Hymenoptera. Future studies on other mountains and along other environmental gradients will be valuable to proof the generality of our ndings in other biogeographic and climatic regions.

Taxonomic groups like birds, ants or terrestrial gastropods strongly differed in their elevational diversity distributions and their species numbers were best predicted by different sets of predictor variables. In approximately 50% of the taxa temperature was the strongest predictor of species richness, but for the other half other variables appeared to be of higher importance. This idiosyncratic response of taxa to environmental gradients is best explained by niche conservatism, that is, the tendency of lineages to retain their niches and related ecological traits over evolutionary time20,3941, leading to phylogenetic autocorre-

lation in elevational distributions (Supplementary Fig. 6) and clade-specic species distributions patterns42,43 (Fig. 1b). Strong constraints in the evolution of niche space (and of related traits and resource use) may result in a low t of diversity to energy variables in narrowly dened taxa6,11. However, the t may increase when extending the taxonomic coverage of communities (Table 1 and Fig. 4).

The integration of multiple taxonomic groups in analyses unequivocally revealed that community level diversity is mainly predicted by temperature. The increasing importance of temperature on diversity for taxonomically more broadly covered communities was driven by two effects: Even though temperature was neither the most important nor a consistently signicant predictor variable across plant and animal taxa, its effect was positive, reaching from subtle to strong, in the majority of taxa

6 NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736 ARTICLE

a

MMT

Mean minimum temperature

Area

n.s.

R2 = 0.72 R2 = 0.79

NPP

Net primary productivity

Plant

species richness

n.s.

n.s.

n.s.

MAP

Mean annual precipitation

MDE

Mid-domain effect

b

MAT

Mean annual temperature

Area

n.s.

R2 = 0.72 R2 = 0.94

n.s.

n.s.

NPP

Net primary productivity

PSR

Plant species richness

Animal

species richness

n.s.

n.s.

n.s.

n.s.

n.s.

MAP

Mean annual precipitation

MDE

Mid-domain effect

Figure 5 | Path models showing direct and indirect effects of predictor variables on species richness. For plants (a) and animals (b), the path model with the lowest Akaike information criterion (AICc) is presented as solid lines. Interrupted lines indicate potential paths used for the construction of competing models (all models with DAICco3 identied by multi-model inference) but which were excluded from the nal path model. For all paths of the nal path model, arrow width is proportional to the relative strength of standardized path coefcients. Orange and black arrows indicate negative and positive effects, respectively. For each endogenous variable the relative amount of explained variance is given. n.s., not signicant.

(Table 1). In addition, many of the taxa exhibiting strong positive temperaturediversity relationships (for example, monocots, rosids, asterids; Orthoptera, all groups of aculeate Hymenoptera, birds, moths) were more diverse, contributing more species to the local species communities, than taxa with negative or weak relationships (for example, magnoliids, millipedes, Collembola; but note that there were also exceptions, for example, the diverse group of ground-dwelling beetles). The largely positive effect of temperature on species richness (Table 1) was independent of differences in sampling intensity among taxa (Supplementary Fig. 5) and contrasts with the heterogeneous effects of other predictor variables: NPP was a well-supported predictor variable only for some single taxa but not for the multi-taxa community44. This was because positive and negative effects of NPP on species richness were equally common in single-taxa analyses. Species richness of most taxonomic groups and of the animal and plant communities as a whole was highest in the dry savannah where primary productivity and precipitation is rather low, suggesting

that even low levels of primary productivity and temporal aridity do not principally put major limitations to the diversication or species coexistence of plant or animal assemblages as a whole. High levels of statistical support for the MDE, which is often used to explain unimodal patterns of elevational diversity15,16, was rare for single taxa and the support decreased with increasing taxonomic coverage. The species richness of a large number of taxa peaked at the lowest elevation and, even in taxa exhibiting unimodal distributions of species richness, the highest richness values were not generally observed at the middle of the domain but at variable elevations (Fig. 1 and Supplementary Fig. 4). These observations suggest a limited role of this type of geometric constraints for the elevational distribution of plant or animal diversity on Mt. Kilimanjaro. However, even though the MDE has often been used to explain patterns of diversity on other mountains which do not extend to the sea level45,46, its application to these systems is debated15. Moreover, recent studies proposed to rene the MDE models by including

NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736

midpoint attractors, that is, peaks of favourable environmental conditions, in models to improve predictions of diversity gradients47. The area available to species populations is regarded as a major determinant of local and regional species richness12,13. However, land area was a poor predictor of plant and animal species richness on Mt. Kilimanjaro, where signicant positive effects were only found for millipedes and Orthoptera. PSR received high statistical support as a predictor of animal species richness only in terrestrial gastropods and moths. Particularly, moths are well known for their specialization on certain food plants which probably facilitates strong associations between herbivore and plant diversity48.

The large support for MMT and MAT as the major predictors of multi-taxa species richness suggests a high importance of temperature-driven mechanisms facilitating the origination and maintenance of biodiversity24,49. While recent studies emphasize the positive effect of temperature on speciation rates2,50,51, the young age of Mt. Kilimanjaro (o2 Mio years) and the fact that the trends were found for diversity estimates taken at local rather than regional spatial scales (study sites of 50 50 m) suggest

the additional importance of ecological mechanisms, such as positive effects of temperature on rates of negative-density dependence52,53 or positive effects on resource exploitation and other biological rates23,49.

While species richness data for many different taxonomic groups were assessed and sampling effort for most taxa was high, several species-rich animal groups were not studied and additional sampling would probably yield more species in local communities. Nonetheless, several points underpin that our results are robust against the inclusion of additional taxonomic groups or increased sampling effort: (i) No taxon considered in this study was included based on an a priori expectation of a certain elevational species richness pattern, and the taxa studied were phylogenetically highly diverse and represent much of the terrestrial plant and animal phylogenetic diversity25 (Fig. 1a).(ii) The sampling resembled a stratied random sampling approach (rst strata: sampling of larger taxa, second strata: sampling of species within taxa), which provides comparative estimates of community level diversity even under incomplete sampling of taxa and species (Supplementary Fig. 7). (iii) Patterns for all taxonomic groups remained consistent when controlling for differences in sampling completeness (Supplementary Table 2) and analyses restricted to taxa sampled with high effort did not change resulting patterns (Fig. 4). (iv) Moreover, community-level analyses in which we applied a rarefaction approach to standardize sampling effort across animal taxa lead to the same conclusions regarding the importance of temperature as the major predictor of species richness (Supplementary Fig. 5).(v) Last, analysing only a fraction of the full data set of taxa unambiguously revealed temperature as the most important predictor of plant and animal species richness in terms of both variable importance and effect strength (Fig. 4).

While our study was extensive in terms of studied taxonomic diversity, it was, however, restricted to a single mountain with a specic geological, climatic and biogeographic context. For example, patterns of elevational diversity have been shown to depend on the extent of the elevational gradient38. Even though the lowest elevations studied at Mt. Kilimanjaro form the natural base of the mountain (as the Mt. Kilimanjaro region and most parts of East Africa are situated at high elevations), equatorial mountains extending to the sea level could provide even more extreme environmental gradients and could show distinct elevational patterns. In addition, our interpretation of the importance of temperature and other predictor variables is based on hypotheses-driven correlative analyses of species distribution data. It will be highly interesting to evaluate

the support for drivers of multi-taxa diversity with both correlative and experimental approaches on other mountains, differing in age, climate, degree of isolation and biogeographic context.

In conclusion, our study revealed that a broad taxonomic coverage in macroecological studies provides new insights into the drivers of broad-scale diversity gradients. Whereas upscaling from regional to global patterns signicantly enhanced the perception of the major drivers of diversity8,13,44, increasing the phylogenetic coverage of studies may be another axis to consider on the way to generalization and causal understanding. While geometric constraints, area, water and productivity variables were of importance only for predicting richness of few groups of low taxonomic level, temperature gained consistent, increasing support as a main predictor of species richness with increasing taxonomic coverage of communities. Identifying the ecological and evolutionary mechanisms by which temperature governs the distribution of biodiversity will be of fundamental importance for understanding global gradients of diversity and the long-term consequences of global warming.

Methods

Time and area of study. All data were collected from December 2011 through January 2014 on the southern and south-eastern slopes of Mt. Kilimanjaro (Tanzania, East Africa; 2450-3250S, 37000-37430E). Mt. Kilimanjaro has a northwest-southeast diameter of B90 km and rises from the savannah plains at 700 m elevation to a snow-clad summit at 5,895 m a.s.l. Precipitation is bimodal with the main rainy season occurring from March through May and the more variable short rains around November. The MAT decreases in a quasi-linear manner with elevation having an overall lapse rate of B0.56 C per 100 m starting with 25 C at the foothills and decreasing to 8 C at the top of the

mountain54. Vertical precipitation distribution shows an unimodal pattern, with 600900 mm year 1 at the base of the mountain, B2,5003,000 mm year 1 at 2,200 m a.s.l. and o500 mm year 1 in the Afroalpine zone at the highest elevations54,55. Because of a long history of human impact, natural habitats in the lowlands were largely cleared54,55. Habitats above 1,800 m a.s.l. are protected as a national park (Mt. Kilimanjaro National Park).

Phylogenetic tree. To demonstrate the distribution of studied taxa across the whole phylogenetic diversity of free-living land plant and animals we constructed a phylogenetic tree of all major plant and animal taxa. Here only the major taxa, that is, taxa with more than 1,000 described species, were considered. For identifying the plant and animal taxa with more than 1,000 described species we used Chapman25. We ltered out all predominately aquatic taxa (for example, most animal phyla nearly exclusively occur in marine environments). Taxonomic and phylogenetic data in http://en.wikipedia.org

Web End =http://en.wikipedia.org (accessed 25 July 2014), Trautweinet al.56 and Misof et al.57 was used to construct a phylogenetic tree of the major terrestrial taxa using the R package ape with branch length computation based on Grafen58. Please note that while the topology of the tree well represents current knowledge, branch lengths are, however, arbitrary and not indicative of clade ages or rates of molecular evolution. We consider this approach as sufcient because the tree is solely used for demonstrating the phylogenetic distribution of studied taxa and was not used for any statistical inference.

Study design. We established thirty 50 50 m study sites, which spanned an

elevational gradient of 8714,550 m a.s.l. Lower elevations are scarcely found in the East African highland region, where Mt. Kilimanjaro is located, and the lowest elevations incorporated in our studied elevational gradient can be considered as a typical, natural base of the mountain (within a square of 2 2 (B222 222 km)

with the Kibo peak of Mt. Kilimanjaro in the centre, 100% of the area is 4500 m a.s.l., 99% 4600 m a.s.l., 94% 4700 m a.s.l. and 83% 4800 m a.s.l.). Thestudy sites were equally distributed over the six major types of natural habitats found along the south-eastern slope of Mt. Kilimanjaro: savannah woodland (B8001,150 m a.s.l.), submontane and lower montane Newtonia forest (B1,1502,050 m a.s.l.), Ocotea forest (B1,8002,800 m a.s.l.), Podocarpus forest (B2,7003,200 m a.s.l.), Erica forest and bushland (B3,2004,000 m a.s.l.), and alpine Helichrysum shrub vegetation (B3,8504,600 m a.s.l.). Five study sites per habitat type were distributed to reect a within-habitat type elevation gradient to detect ne scale changes in biodiversity with changing elevation. Spatial distances among study sites were in all cases larger than 300 m. If possible, study sites were established in core zones of larger areas of the respective habitat type, so that effects of transition zones were minimized.

8 NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736 ARTICLE

Species richness data. For each taxonomic group, standardized approaches were used to collect species richness data on all 30 study sites. For taxonomic groups where the sampling effort per site was the same, we calculated species richness as the total (cumulative) number of species per study site. These included all groups of vascular plants, terrestrial gastropods, millipedes, spiders, true bugs, Collembola, parasitoid wasps, ants, bees, other aculeate Hymenoptera, ground-dwelling beetles (without dung beetles), dung beetles and birds. Where sampling effort per site varied we used the following approaches: For moths and bats, species richness per site was calculated by averaging the species richness values of the single surveys. For frogs, Zancolli et al.59 aimed at sampling complete assemblages and therefore sampling effort per site was adjusted to measure asymptotic species richness. The same was true for Orthoptera60,61. In both taxa, we used the total number of species per site detected by the authors. In hoveries numbers of samples varied among study sites. Initially, we reduced the number of samples on sites to a common number and calculated the cumulative species richness from the reduced data set. However, as many species were missing in this data set and the number of species per site was generally low, we run the analysis with the full (unequally sampled) data set and compared the patterns to those of the reduced data set: The pattern was very similar but the variation was much smaller. Therefore, for hoveries we used the full data set. Details on the sampling procedures for all taxa are described in the Supplementary Methods.

Temperature and precipitation. All study sites were equipped with temperature sensors that were installed B2 m above the ground. Coated plastic funnels were used for radiation shielding. The temperature sensors measured temperatures in 5 min intervals for a time period of B2 years. We subsequently calculated the MAT as the average of all measurements per study site and the MMT as the average over all monthly temperature minima per study site. On six study sites data loggers were repeatedly stolen and we had to estimate temperature data for these study sites. On the basis of the observed MAT and MMT data we calculated a linear model with MAT/MMT as the response and elevation and habitat type as additive explaining variables (R2 0.99, N 24, Po0.01 for both temperature variables).

Using this model we predicted the missing six temperature values based on the elevation and habitat type of study sites. MAP was interpolated for every study site using a co-kriging approach based on a 15-year data set from a network of about 70 rain gauges on Mt. Kilimanjaro55.

Net primary productivity. We used the normalized difference vegetation index (NDVI) as a proxy for NPP. Because of the negative biases due to sensor degradation of MODIS Terra62, NDVI estimations were exclusively based on MODIS Aqua product MYD13Q1 with a horizontal resolution of 250 m 250 m. Cloud

contamination is a very prominent feature on Mt. Kilimanjaro and problematic with regard to realistic estimations of NDVI. To address this issue we rst identied all pixels with a MYD13Q1 quality ag of three and deleted these together with the eight adjacent surrounding pixels. Afterwards, we followed the approach proposed by Atzberger and Eilers63 using the Whittaker smoother based on three iterations with a lambda of 6,000. We then calculated the overall mean NDVI for the 10-year period 20032012 and extracted the pixel values corresponding to the locations of our study sites. A comparison of NDVI data with on-site measurements of the leaf-area index conducted during the growing season revealed a high level of correlation (r 0.84, Po0.001).

Area. We calculated for each study site the available land area within a range of 300 m above and below the elevational level of study sites. This area data was derived from a digital elevation model of Mt. Kilimanjaro with a resolution of 30 m of the following extent: 37.0007437.75602E, 3.5075332.750183S.

Mid-domain effect. We determined for each species the elevational range, which is dened by its minimal and maximum elevation of occurrence. Incomplete sampling routinely underestimates range sizes64, with the most extreme case being species recorded from only a single elevation (which was the case for many species found) which thus would have an observed elevational range of 0 m. To adjust range underestimation we followed the approach used by Brehm et al.16 and added 265 m to each end of each recorded range (for all ranges), which corresponds to half the maximum elevational distance between any two adjacent sampling elevations. We then systematically reassigned the location of each of the interpolated, augmented ranges within the domain (6064,815 m) at random (sampling without replacement) and then recorded the predicted richness for the elevation of each study site. This procedure was repeated for 200 times after which the mean predicted richness for each study site was calculated by averaging over the 200 predicted richness values of the individual repeats.

Statistical analysis. Generalized additive models (gam) were used to model the relationships between species richness and elevation65,66, setting the data family to Gaussian type and the basis dimension of the smoothing function to ve.

We calculated a complexity value to characterize the shape of elevational species richness patterns with one numeric measure by comparing the explained deviance of generalized linear models (EDglm) with the explained deviance of generalized

additive models (EDgam) of species richness on elevation: complexity

(EDgam EDglm)/EDgam. Supplementary Fig. 8 exemplies calculations of

complexity for a pattern with a mid-elevational peak, a linear decline and an exponential decline. Note that in case of high error variation, EDglm may theoretically equal EDgam so that complexity values B0 even though no statistically well supported trend of linear decline is evident. However, in our data set explained deviance was consistently high (40.3) and all complexity values of o0.2 were associated with signicant negative trends (slopes of glm o0, Po0.05).

To analyse the statistical support for explanatory variables in predicting plant and animal species richness we used multi-model inference based on information-theory and ordinary least-square regression (20). Multi-model inference not only accounted for uncertainty in parameter estimates but also for uncertainty in model selection29,30,67. Moreover, multi-model inference is a way to objectively deal with correlated predictor variables as strongly correlated predictor variables with the same explanatory power will get reduced support in multi-model inference (Supplementary Fig. 9). In the multi-model inference procedure, two nal variables were calculated: rst, variable importance, which is a measure of the relative support a predictor variable receives over the full model space; second, conditional model-averaged parameter estimates and signicance levels. We z-transformed all predictor variables before analyses, so that the model-averaged parameter estimates were standardized beta values. This transformation allowed us to compare the relative inuence on species richness among the set of predictor variables which were measured at different scales. For multi-model averaging analyses the R package MuMIn (http://CRAN.R-project.org/package=MuMIn

Web End =http://CRAN.R-project.org/package=MuMIn, accessed 25 April 2014) was used.

To analyse elevational species richness patterns and the support for different predictor variables with increasing taxonomic coverage we conducted the following steps: (1) we calculated for all combinations of 18 plant (N 256) and for all

combinations of 116 animal groups (N 65,535; that is, for taxa shown individually

in Fig. 1 and Supplementary Fig. 3) the pooled (cumulative) species richness; (2) we modelled pooled species richness as a function of elevation using generalized-additive models and calculated the complexity measure of the elevational species richness pattern; (3) we calculated MDE predictions for all study sites based on the species range information of the randomly combined taxa; and (4) we ran a multi-model inference analysis with pooled species richness as the response and MAT (in animals) or MMT (in plants), MAP, NPP, land area, the predictions of MDE model and PSR (for explaining richness of animal taxa) as predictor variables. These analyses were conducted in the same way as the single taxa analyses.

We validated the inuence of incomplete sampling on diversity patterns and inference on predictor variables using a battery of analyses: rst, we calculated measures of sample coverage for all taxa using the R package iNEXT68 and analysed their correlation with elevation. To analyse the inuence of incomplete sampling on patterns of elevational diversity, we estimated asymptotic species richness for all single taxa for which abundance data was available, using a non-parametric species richness estimator for abundance-based data, that is, the Chao1 index27,28,68. We used generalized additive models to model Chao1-estimated species richness as a function of elevation and predicted for each study site the number of species from these models. We correlated the predicted estimates of asymptotic species richness to the ones predicted by the original model (based on the observed species richness data) using Pearsons product moment correlation. If variation in incomplete sampling signicantly biased elevational species richness patterns, we expected predictions not to be highly correlated.

Moreover, as we counted the species of all animal taxa equal in the pooled community-level analyses, results of statistical inference could be biased due to differences in the sampling intensity (that is, here the number of individuals which were sampled per taxon) among taxa. We therefore analysed patterns of elevational species richness and the support for predictor variables additionally for a standardized data set, in which for each animal taxon the number of individuals was rareed to the lowest number of individuals observed in any of the taxa (that is, N 83) and used for the calculation of variable importance and effect strength of

all predictor variables. This procedure was repeated for 5,000 times. From the results of these analyses of standardized data sets we calculated the mean and 95% condence intervals of variable importance and standardized beta for all predictor variables.

Path analysis69 was used to disentangle the direct and indirect effect of climate (temperature and precipitation), primary productivity, area, the MDE and plant richness (the latter only for animal species richness as the response variable) on pooled plant and animal species richness. All explanatory variables were standardized by z-transformation using the scale function in R. On the basis of prior studies, we hypothesized that temperature (MAT or MMT) and precipitation (MAP) predict species richness directly and indirectly via their combined effect on primary production (NPP) and on the species richness of plants (PSR, that is, for animal richness as the response variable). Therefore, we pre-selected possible path combinations, by analysing the response variables of our path models (animal species richness, PSR, primary productivity) with all explaining variables:

For plants:

NPP MMT MAP 1

PSR MMT MAP NPP Area MDE 2

NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736

22. Khsel, S. & Blthgen, N. High diversity stabilizes the thermal resilience of pollinator communities in intensively managed grasslands. Nat. Commun. 6, 7989 (2015).

23. Classen, A. et al. Temperature versus resource constraints: which factors determine bee diversity on Mount Kilimanjaro, Tanzania? Bee species richness on Mt Kilimanjaro. Glob. Ecol. Biogeogr 24, 642652 (2015).

24. Allen, A. P., Gillooly, J. F. & Brown, J. H. in Scaling Biodiversity. (eds Storch, D., Marquet, P. A. & Brown, J. H.) 283299 (Cambridge University Press, 2007).

25. Chapman, A. D. Numbers of Living Species in Australia and the World (Australian Biodiversity Information Services, 2009).

26. Chao, A., Colwell, R. K., Lin, C.-W. & Gotelli, N. J. Sufcient sampling for asymptotic minimum species richness estimators. Ecology 90, 11251133 (2009).

27. Chao, A. et al. Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecol. Monogr. 84, 4567 (2014).

28. Chao, A. Nonparametric estimation of the number of classes in a population. Scand. J. Stat. 11, 265270 (1984).

29. Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261304 (2004).

30. Burnham, K. P., Anderson, D. R. & Huyvaert, K. P. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behav. Ecol. Sociobiol. 65, 2335 (2011).

31. Storch, D. in Metabolic Ecology: A Scaling Approach. (eds Sibly, R. M., Brown, J. H. & Kodric-Brown, A.) 120131 (John Wiley & Sons, Ltd, 2012).

32. Grace, J. B. et al. Integrative modelling reveals mechanisms linking productivity and plant species richness. Nature 529, 390393 (2016).

33. Colwell, R. K., Brehm, G., Cardelus, C. L., Gilman, A. C. & Longino, J. T. Global warming, elevational range shifts, and lowland biotic attrition in the wet tropics. Science 322, 258261 (2008).

34. Duan, M. et al. Disentangling effects of abiotic factors and biotic interactions on cross-taxon congruence in species turnover patterns of plants, moths and beetles. Sci. Rep. 6, 23511 (2016).

35. Sanders, N. J. et al.in Data Mining for Global Trends in Mountain Biodiversity. (eds Krner, C. & Spehn, E. M.) (CRC Press, 2010).

36. Kessler, M. Elevational gradients in species richness and endemism of selected plant groups in the central Bolivian Andes. Plant Ecol. 149, 181193 (2000).

37. Rahbek, C. The elevational gradient of species richness: a uniform pattern? Ecography 18, 200205 (1995).

38. Rahbek, C. The role of spatial scale and the perception of large-scale species-richness patterns: scale and species-richness patterns. Ecol. Lett. 8, 224239 (2005).

39. Wiens, J. J. & Graham, C. H. Niche conservatism: integrating evolution, ecology, and conservation biology. Annu. Rev. Ecol. Evol. Syst. 36, 519539 (2005).40. Jablonski, D. et al. Out of the tropics, but how? Fossils, bridge species, and thermal ranges in the dynamics of the marine latitudinal diversity gradient. Proc Natl. Acad. Sci. USA 110, 1048710494 (2013).

41. Merckx, V. S. F. T. et al. Evolution of endemism on a young tropical mountain. Nature 524, 347350 (2015).

42. Wu, Y. et al. Understanding historical and current patterns of species richness of babblers along a 5000-m subtropical elevational gradient: historical and current babbler richness patterns. Glob. Ecol. Biogeogr. 23, 11671176 (2014).

43. Kozak, K. H. & Wiens, J. J. Niche conservatism drives elevational diversity patterns in Appalachian salamanders. Am. Nat. 176, 4054 (2010).

44. Tittensor, D. P. et al. Global patterns and predictors of marine biodiversity across taxa. Nature 466, 10981101 (2010).

45. Baur, B., Meier, T., Baur, A. & Schmera, D. Terrestrial gastropod diversity in an alpine region: disentangling effects of elevation, area, geometric constraints, habitat type and land-use intensity. Ecography 37, 390401 (2014).

46. McCain, C. M. Elevational gradients in diversity of small mammals. Ecology 86, 366372 (2005).

47. Colwell, R. K. et al. Midpoint attractors and species richness: modelling the interaction between environmental drivers and geometric constraints. Ecol. Lett. 19, 10091022 (2016).

48. Forister, M. L. et al. The global distribution of diet breadth in insect herbivores. Proc. Natl Acad. Sci. USA 112, 442447 (2015).

49. Zhou, J. et al. Temperature mediates continental-scale diversity of microbes in forest soils. Nat. Commun. 7, 12083 (2016).

50. Allen, A. P., Gillooly, J. F., Savage, V. M. & Brown, J. H. Kinetic effects of temperature on rates of genetic divergence and speciation. Proc. Natl Acad. Sci. 103, 91309135 (2006).

51. Rolland, J., Condamine, F. L., Jiguet, F. & Morlon, H. Faster speciation and reduced extinction in the tropics contribute to the mammalian latitudinal diversity gradient. PLoS Biol. 12, e1001775 (2014).

52. Bagchi, R. et al. Pathogens and insect herbivores drive rainforest plant diversity and composition. Nature 506, 8588 (2014).

53. Johnson, D. J., Beaulieu, W. T., Bever, J. D. & Clay, K. Conspecic negative density dependence and forest diversity. Science 336, 904907 (2012).

For animals:

NPP MAT MAP 3

PSR MAT MAP NPP Area MDE 4

Animal species richness MAT MAP NPP Area MDE PSR 5 For each model we used the dredge function of the R package MuMIn to evaluate models dened by all possible variable combinations and ranked them by their AIC-based model weight. As our sample size was relatively low compared with the number of estimated parameters, we used the AIC with a second-order bias correction (AICC) for inferring the support of individual models. From the set of best supported models (D AICCo3) we calculated paths models using the R package lavaan70 and ranked them based on their AICC. In Fig. 5 we show paths coefcients (width of arrows), their statistical signicance and multiple coefcients of determination (R2) of predictor variables of the best supported path models with the lowest AICC. In addition we indicate paths of competitive models (D AICCo3)

which were not included in the best model.

Data availability. The data that support the ndings of this study are available from the following authors upon request: corresponding author (M.K.P.):data of geographic coordinates, elevation, land area, mid-domain effect predictions for study sites, data on Heteroptera, ants, parasitoid and other aculeate wasps; A.H. (mailto:[email protected]

Web End [email protected]): plant data, mean annual precipitation; T.A. (mailto:[email protected]

Web End [email protected]): NPP, temperature data; A.C. (mailto:[email protected]

Web End [email protected]): data on bees and hoveries; S.W.F. (mailto:[email protected]

Web End [email protected]): bird data; S.B.F. (mailto:[email protected]

Web End [email protected]): millipede data; F.G. (mailto:[email protected]

Web End [email protected]): dung beetledata; M.H.-B. (mailto:[email protected]

Web End [email protected]): bat and moth data; C.H. (mailto:[email protected]

Web End [email protected]): Orthoptera data; C.N. (mailto:[email protected]

Web End [email protected]): gastropod data; J.R. (mailto:[email protected]

Web End [email protected]): data on spiders, Collembola, ground-dwelling beetles; G.Z. (mailto:[email protected]

Web End [email protected]): amphibian data.

References

1. Gaston, K. J. Global patterns in biodiversity. Nature 405, 220227 (2000).2. Brown, J. H. Why are there so many species in the tropics? J. Biogeogr. 41, 822 (2014).

3. Fine, P. V. A. Ecological and evolutionary drivers of geographic variation in species diversity. Annu. Rev. Ecol. Evol. Syst. 46, 369392 (2015).

4. Rhode, K. Latitudinal gradients in species diversity: the search for the primary cause. Oikos 65, 514527 (1992).

5. Mittelbach, G. G. et al. Evolution and the latitudinal diversity gradient: speciation, extinction and biogeography. Ecol. Lett. 10, 315331 (2007).

6. Allen, A. P., Brown, J. H. & Gillooly, J. F. Global biodiversity, biochemical kinetics, and the energetic-equivalence rule. Science 297, 15451548 (2002).

7. Hawkins, B. A. et al. Energy, water, and broad-scale geographic patterns of species richness. Ecology 84, 31053117 (2003).

8. Kreft, H. & Jetz, W. Global patterns and determinants of vascular plant diversity. Proc. Natl Acad. Sci. USA 104, 59255930 (2007).

9. Hutchinson, G. E. Homage to Santa Rosalia or why are there so many kinds of animals? Am. Nat. 93, 145159 (1959).

10. Brown, J. H. Two decades of homage to Santa Rosalia: toward a general theory of diversity. Am. Zool 21, 877888 (1981).

11. Hurlbert, A. H. & Stegen, J. C. When should species richness be energy limited, and how would we know? Ecol. Lett. 17, 401413 (2014).

12. Rosenzweig, M. L. Species Diversity in Space and Time (Cambridge University Press, 1995).

13. Belmaker, J. & Jetz, W. Relative roles of ecological and energetic constraints, diversication rates and region history on global species richness gradients. Ecol. Lett. 18, 563571 (2015).

14. Romdal, T. S. & Grytnes, J.-A. An indirect area effect on elevational species richness patterns. Ecography 30, 440448 (2007).

15. Colwell, R. K., Rahbek, C. & Gotelli, N. J. The mid-domain effect and species richness patterns: what have we learned so far? Am. Nat. 163, E1E23 (2004).

16. Brehm, G., Colwell, R. K. & Kluge, J. The role of environment and mid-domain effect on moth species richness along a tropical elevational gradient. Glob. Ecol. Biogeogr. 16, 205219 (2007).

17. Novotny, V. et al. Why are there so many species of herbivorous insects in tropical rainforests? Science 313, 11151118 (2006).

18. Basset, Y. et al. Arthropod diversity in a tropical forest. Science 338, 14811484 (2012).

19. McCain, C. M. & Grytnes, J.-A. Encyclopedia of Life Sciences. John Wiley & Sons, Ltd, 2010).

20. Wiens, J. J. et al. Niche conservatism as an emerging principle in ecology and conservation biology: niche conservatism, ecology, and conservation. Ecol. Lett. 13, 13101324 (2010).

21. Crisp, M. D. et al. Phylogenetic biome conservatism on a global scale. Nature 458, 754756 (2009).

10 NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms13736 ARTICLE

54. Hemp, A. Continuum or zonation? Altitudinal gradients in the forest vegetation of Mt. Kilimanjaro. Plant Ecol. 184, 2742 (2006).

55. Hemp, A. Vegetation of Kilimanjaro: hidden endemics and missing bamboo. Afr. J. Ecol. 44, 305328 (2006).

56. Trautwein, M. D., Wiegmann, B. M., Beutel, R., Kjer, K. M. & Yeates, D. K. Advances in insect phylogeny at the dawn of the postgenomic era. Annu. Rev. Entomol. 57, 449468 (2012).

57. Misof, B. et al. Phylogenomics resolves the timing and pattern of insect evolution. Science 346, 763767 (2014).

58. Grafen, A. The phylogenetic regression. Philos. Trans. R. Soc. B Biol. Sci. 326, 119157 (1989).

59. Zancolli, G., Steffan-Dewenter, I. & Rdel, M.-O. Amphibian diversity on the roof of Africa: unveiling the effects of habitat degradation, altitude and biogeography. Divers. Distrib. 20, 297308 (2014).

60. Hemp, C. Annotated list of Caelifera (Orthoptera) of Mt. Kilimanjaro, Tanzania. J. Orthoptera Res. 18, 183214 (2009).

61. Hemp, C. Annotated list of Ensifera (Orthoptera) and further records on Caelifera (Orthoptera) of Mt Kilimanjaro, Tanzania. Zootaxa 3613, 301342 (2013).

62. Wang, D. et al. Impact of sensor degradation on the MODIS NDVI time series. Remote Sens. Environ. 119, 5561 (2012).

63. Atzberger, C. & Eilers, P. H. C. A time series for monitoring vegetation activity and phenology at 10-daily time steps covering large parts of South America. Int.J. Digit. Earth 4, 365386 (2011).64. Colwell, R. K. & Hurtt, G. C. Nonbiological gradients in species richness and a spurious rapoport effect. Am. Nat. 144, 570595 (1994).

65. Crawley, M. J. The R Book (Wiley, 2013).66. Wood, S. N. Generalized Additive Models: An Introduction with R (Chapman & Hall/CRC, 2006).

67. Johnson, J. B. & Omland, K. S. Model selection in ecology and evolution. Trends Ecol. Evol. 19, 101108 (2004).

68. Hsieh, T. C., Ma, K. H. & Chao, A. iNEXT: an R package for rarefaction and extrapolation of species diversity (Hill numbers). Methods Ecol. Evol. doi:http://dx.doi.org/10.1111/2041-210X.12613

Web End =10.1111/2041-210X.12613 (2016).

69. Shipley, B. Cause and Correlation in Biology: A Users Guide to Path Analysis, Structural Equations and Causal Inference (Cambridge Univ. Press, 2004).

70. Rosseel, Y. lavaan: an R package for structural equation modelling. J. Stat. Softw. 48, 136 (2012).

Acknowledgements

We thank the Tanzanian Commission for Science and Technology, the Tanzania Wildlife

Research Institute and the Kilimanjaro National Park authority for their support and for

granting us access to the Kilimanjaro National Park area. We are grateful to all the

companies and private farmers who allowed us to work on their land. We thank the KiLi

eld staff for helping to collect data at Mt. Kilimanjaro. This study was conducted within

the framework of the Research Unit FOR1246 (Kilimanjaro ecosystems under global

change: linking biodiversity, biotic interactions and biogeochemical ecosystem processes,

https://www.kilimanjaro.biozentrum.uni-wuerzburg.de) funded by the Deutsche For

schungsgemeinschaft (DFG).

Author contributions

I.S.-D., A.H., M.F. designed the concept for ecological research at Mt. Kilimanjaro. A.H.,

M.-O.R., M.S., K.B.-G., R.B., E.K.V.K., M.K., T.N., M.T., M.F. and I.S.-D. conceptualized

and supervised single taxa studies. A.H. implemented study sites. M.K.P. and I.S.-D.

conceived the study. M.K.P., A.H., T.A., C.B., A.C., F.D., A.E., S.B.F., S.W.F., F.G., M.H.,

M.H.-B., C.H., W.J.K., E.M., C.N., I.O., J.R., G.R., D.S.C., J.T. and G.Z. collected the data.

J.D., C.D.E., R.S.P. and A.S. identied signicant quantities of specimens. A.H., C.H.,

V.K. and J.Z. organized and maintained logistic infrastructure. M.K.P. processed and

analysed the data and wrote the rst version of the manuscript with input from I.S.-D.

All authors contributed to the nal version of the manuscript.

Additional information

Supplementary Information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/

http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: The authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions/

Web End =http://npg.nature.com/

http://npg.nature.com/reprintsandpermissions/

Web End =reprintsandpermissions/

How to cite this article: Peters, M. K et al. Predictors of elevational biodiversity gradients

change from single taxa to the multi-taxa community level. Nat. Commun. 7, 13736

doi: 10.1038/ncomms13736 (2016).

Publishers note: Springer Nature remains neutral with regard to jurisdictional claims in

published maps and institutional afliations.

This work is licensed under a Creative Commons Attribution 4.0

International License. The images or other third party material in this

article are included in the articles Creative Commons license, unless indicated otherwise

in the credit line; if the material is not included under the Creative Commons license,

users will need to obtain permission from the license holder to reproduce the material.

To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/

Web End =http://creativecommons.org/licenses/by/4.0/

r The Author(s) 2016

NATURE COMMUNICATIONS | 7:13736 | DOI: 10.1038/ncomms13736 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 11