Unravelling spatial scale effects on elevational diversity gradients: insights from montane small mammals in Kenya

Abstract

Background

Montane ecosystems play crucial roles as global biodiversity hotspots. However, climatic changes and anthropogenic pressure increasingly threaten the stability of montane community dynamics, such as diversity-elevation interactions, creating a challenge in understanding species biogeography and community ecology dynamics in these crucial conservation areas. We examined how varying sampling spatial grains influence small mammal diversity patterns within Kenya’s tallest montane ecosystems.

Methods

Employing a combination of multidimensional alpha diversity metrics and multisite beta diversity characteristics (species richness, phylogenetic and functional diversity and divergence, and multisite beta diversity) alongside spatial generalized additive multivariate regression analyses, we tested how spatial scaling influences elevational diversity gradient patterns and their associations with environmental and human activity variables.

Results

The diversity-elevation associations were generally homogeneous across spatial grains; however, idiosyncratic patterns emerged across mountains. The total (taxonomic, phylogenetic, and functional) beta diversity, nestedness, and turnover resultant components monotonically increased or decreased with varying spatial grains. The associations between the diversity patterns and the environmental and human footprint variables increased with spatial grain size but also presented variations across mountains and indices. Species richness and phylogenetic and functional richness indices were more strongly influenced by spatial scale variations than were the divergence and community structure indices in both the diversity distribution patterns and their associations with the environmental and human variables.

Conclusions

The diversity-elevation and diversity-environment (including human activity pressure) relationships across spatial grains suggest that montane small mammal diversity patterns portray subtle but systematic sensitivity to sampling spatial grain variation and underscore the importance of geographical context in shaping these elevational diversity gradients. For improved effectiveness, conservation efforts should consider these spatial effects and the unique geographical background of individual montane ecosystems.

Introduction

Achieving sustainable ecosystem management and conservation strategies that better preserve terrestrial biodiversity and promote sustainable land use [1] requires a refined understanding of the mechanisms maintaining the geographically uneven distribution of species [2,3,4,5,6,7]. Montane ecosystems are of particular interest in this regard because they are global hotspots for species richness and endemism but face more intense pressure from human activities and are more vulnerable to the accelerating impacts of climate change [8,9,10]. Mountains are thus pivotal to understanding the interplay between biodiversity distribution trends, ecological-evolutionary processes, and anthropogenic correlates. Notably, while most observations of montane species diversity are governed by robust species-area relationship hypotheses [11,12,13], focused studies spanning local scales, multiple diversity dimensions, poorly studied areas such as Afrotropical regions, and the relative influence of environmental constraints remain sparse. Such studies could bolster the robustness of ecological theories governing elevational diversity gradient patterns [14] and enhance ecosystem management and conservation strategies [9].

Understanding the spatial scale dependency of species’ ecological and evolutionary patterns and processes is crucial for interpreting observed biodiversity patterns, but it remains challenging to resolve between geographies and species groups [15,16,17]. Levin [18] demonstrated some of these challenges, noting that different spatial levels reveal different processes, all pivotal for accurate ecological modeling and effective conservation planning. The current broad consensus in species-area relationships [19,20,21] is the confluence of several guiding hypotheses, including beta diversity and spatial scale [22, 23], scale-dependence of diversity [18, 24], environmental heterogeneity [25, 26], and energy-availability [13, 27]. A common observation across these hypotheses is that ecological patterns and processes and conservation management strategies that are generalizable when sampling at a larger grain size may not apply at smaller scales [28], where species interactions and more local environmental conditions can have stronger influences on community structuring [15,16,17, 29,30,31].

Studying mountains’ high biodiversity value, evolutionary significance, and ecological uniqueness [32, 33] across a broader spatial continuum could help to better scale responses to changing climatic regimes and anthropogenic influences. Recent studies have highlighted the importance of spatial scaling in sampling montane biodiversity observations, emphasizing its role in interpreting elevational diversity gradients and the mechanisms driving these patterns [34,35,36,37,38]. To illustrate the complexities of scaling elevational diversity patterns, Rahbek [16] found that the well-documented mid-elevation peak in species richness—a pattern often observed in tropical mountains—considerably varies and can even reverse when analyzed at different spatial scales. Similarly, McCain and Grytnes [39] showed that the hump-shaped pattern in species richness along elevational gradients is common only at broader spatial scales, whereas finer, localized scales yield more variable patterns where such a peak shifts or even disappears. Others like Tello et al. [38] reported scale-dependent variation in community assembly mechanisms that determined beta diversity patterns across elevations, Graham et al. [36] demonstrated how spatial variation in species richness along elevational gradients is influenced by ecological factors and evolutionary mechanisms, and Montes et al. [37] discussed optimization options for scale-aware biodiversity sampling by accounting for changes in biodiversity across different spatial extents.

Despite a general understanding of montane eco-evolutionary dynamics in mammal communities [39, 40], increasingly unpredictable climatic conditions and human encroachment into previously pristine elevation bands suggest that existing theories may no longer reliably explain species diversity and distribution dynamics across some montane landscapes. For example, upward shifts in species habitats due to changing climatic regimes and human pressures [41,42,43,44] fundamentally alter classical diversity-elevation patterns and may establish new, unknown trends [45]. Moreover, most research on elevational diversity gradients has focused on species richness as a measure of biodiversity; however, decoupled distribution patterns between species richness and various indices of phylogenetic and functional diversity [46,47,48,49,50,51,52,53,54] underscore the importance of preserving not only the number of species but also their evolutionary history and functional roles [53, 54]. Phylogenetic diversity indices quantify species evolutionary history represented within a community, emphasizing regions characterized by unique speciation trajectories, which helps prioritize the conservation of evolutionarily distinct species [46], while functional diversity assesses species ecological roles, essential for maintaining ecosystem processes and overall biodiversity health [53, 54]. For instance, while the turnover component of beta diversity is often linked to environmental filtering or species replacement, the nestedness component highlights species loss due to habitat degradation [55,56,57]. Furthermore, distinct responses of different diversity dimensions to environmental filters, especially in unique ecological contexts, such as at the extreme ends of suitability for most species [58], challenge the standard practice of using species richness as an umbrella biodiversity index [59]. Applying multifaceted biodiversity indices is a more holistic approach to identifying areas critical for conservation, ensuring that both trait and evolutionary history are preserved across landscapes [55, 60,61,62]. More research incorporating a multidimensional perspective on montane biodiversity—considering phylogenetic and functional diversity dimensions alongside species richness—would provide a comprehensive framework for conservation by identifying functionally and phylogenetically distinct species that may otherwise be overlooked when considering taxonomic diversity alone.

As typical tropical montane ecosystems [33, 63], mountains in Kenya are characterized by rich mammalian diversity driven by sharp climatic and vegetational shifts along the elevation gradients [52, 64, 65]. In these ecosystems, altitude range, geographic location, and aspect are primary drivers of the faunal diversity and community structure through intricate abiotic and biotic feedback [39, 52, 66,67,68]. For example, Mount Kenya experiences more precipitation on its southeastern slopes, fostering lush montane forests that support diverse mammal species, in contrast to the leeward side, which receives less rainfall but similarly supports rich and diverse assemblages. The distinct regional and, in many cases, local variations in temperature and precipitation between and within mountains result in differentiated habitats, shaping the composition and diversity of montane mammal communities [52, 69, 70]. In many Afrotropical regions, such as Kenya, the extent and variability of elevational diversity gradient patterns across diverse geographic areas and species groups remain largely unexplored empirically, often relying on inferential assumptions for interpretation.

Here, we investigated how spatial grain variations in sampling influence observed diversity–elevation association patterns in small mammal communities across Kenya’s montane ecosystems, using an ecologically sensitive species group—small mammals [71], a comprehensive spatially structured community dataset, and multiple biodiversity indices. We defined and adopted ten spatially varied sampling grains, ranging from microhabitats to landscape scales, across Kenya’s major montane ecosystems. The main objectives were (i) to examine how the spatial grain of field sampling affects diversity-elevation patterns and (ii) to determine whether variations in sampling grain influence the relationships between diversity patterns and environmental variables and human ecological footprint. We predicted that species, phylogenetic, and functional richness would increase with sampling spatial scale following the species-area relationship theory, which suggests that because larger areas encompass more habitats, environmental heterogeneity, and resources, they offer more opportunities for species coexistence and reduce extinction risk through stochastic processes, thereby preserving more species [19,20,21]. Similarly, we expected beta diversity to increase with spatial scale, because at smaller scales, low abiotic variability and limited habitats result in similar species assemblages, which leads to low beta diversity between communities. However, as the spatial scale increases, greater environmental heterogeneity fosters more differentiation among communities, thereby increasing beta diversity [22, 23]. We also predicted that as sampling grain size increased, stronger relationships between diversity patterns and environmental variables would emerge, in line with the environmental heterogeneity and energy-availability hypotheses, which propose that more diverse conditions in larger areas support more species due to the availability of varied habitats and resources, fostering coexistence and enhancing biodiversity [25, 26]. Finally, we anticipated a stronger correlation between diversity patterns and human ecological footprint at intermediate sampling grains, consistent with the intermediate disturbance hypothesis. The hypothesis posits that species diversity is maximized at moderate disturbance levels, where neither too rare nor too frequent disturbances allow both early and late successional species to coexist, preventing competitive exclusion. In contrast, rare disturbances can lead to dominance by competitive species, while frequent disturbances may eliminate many species—both scenarios reducing the overall diversity [72].

Material and methods

Study area and sampling design

We focused on small mammals in the orders Rodentia, Eulipotyphla, and Macroscelidea across the highest peaks in Kenya—Mount Kenya, Mount Elgon, the Aberdare Range, the Cherangani Hills, the Mathews Range, Mount Kulal, and the Chyulu Hills (Fig. 1, Additional file 2). These mountains feature unique topographies associated with Miocene and Pleistocene tectonic and volcanic activities that also created unique geomorphological features [73,74,75,76,77,78,79,80,81,82]. Mount Kenya, Mount Elgon, the Aberdare Range, Mount Kulal, and the Chyulu Hills were formed through volcanic activity, with their structures dominated by volcanic cones and craters, whereas the Cherangani Hills and the Mathews Range were formed through tectonic processes. The diverse landscapes and microclimates resulting from their formation have led to rich biodiversity and high levels of endemism, as the fertile volcanic soils and varied altitudes support diverse ecosystems, providing habitats for diverse flora and fauna [32]. These mountains also act as ecological refuges, enabling species to diverge, adapt, and persist during climatic changes [43, 83]. Collectively, they are critical biodiversity hotspots, serving essential roles in water catchment, ensuring ecological balance, and supporting wildlife conservation and the livelihoods of local communities [65]. However, they are all characterized by a lack of dedicated formal accounts of their biodiversity and the processes driving them, especially for small mammals, for which data are nonexistent for most taxa.

Fig. 1

Topographical maps of the study sites and illustrations of the sampling scheme. The main figure shows the locations and names of the studied mountains, with the color representing the elevation variations (natural breaks classification [Jenks]). The inset figures show (a) the geographical locality of the study area within Kenya and (b) an illustration of the geographic and spatial layout of the sampling scheme implemented in the study—strictly nested quadrat design

For each mountain, we designed spatially varied sampling windows/grain sizes [57] using the strictly nested quadrat method [12, 84]. Thus, we maintained the spatial extent and shape of the sampling window while varying the sampling grain/unit, guided by spatial scale concept’s two key aspects: the spatial extent (overall size of the area considered in a specific study) and spatial grain or resolution (the dimensions of the individual spatial units within that area, for which observations or predictions are made). This approach ensures more accurate estimates of expected species richness for a randomly located plot within a given area [12, 57]. We overlaid sampling grids spanning the extents of mountains for each spatial grain size—0.0001, 0.001, 0.01, 0.1, 1, 2, 4, 6, 8, and 10 km²—resulting in ten community datasets for each mountain. The spatial grain sizes were defined based on leveraging our data resolution, previous studies, and the field sampling scale at which small mammal communities demonstrate unique biogeographic diversity [52, 85]. The largest grain size was used as the sampling baseline grain size for the smaller units. Each subsequent smaller-grain grid was created by first determining the centroid of the baseline grid and then creating a buffer of varied radii around it. The radii lengths were back-calculated to match the predetermined grid sizes using the formula Pi multiplied by the squared radius (A = πr²). Due to the different geomorphological characteristics of each mountain (lateral extents and peak-base heights), we ultimately sampled variable numbers of grids and elevation limits across mountains based on the transition in environmental conditions (histograms of the climate and elevation data), vegetation, and species turnover from the surrounding lowlands (Additional file 2).

Species checklist and occurrence records

Across the seven mountains, we retrieved 124 small mammal species in the orders Eulipotyphla (30 species), Macroscelidea (4 species), and Rodentia (90 species) from local inventories and global checklists [86,87,88,89,90], ranging in body mass from 2.9 to 3,327.5 g (see Additional file 1 Table S1). Species distribution records were also determined from these inventories and checklists, in addition to the IUCN Red List [91] spatial extent estimations and species occurrence records from the Global Biodiversity Information Facility (GBIF) [92]. For species whose ranges could not be obtained from these sources, we created new polygon layers based on the literature distribution accounts. All distribution layers were merged into a single species dataset in QGIS 3.34.2 [93].

Community composition matrix

The community composition matrix (site × species matrix describing interactions between species and the studied community) was obtained by superimposing the species distribution layer with the spatially-structured grid datasets using the ‘Join Attributes by Location’ tool in QGIS. The extracted dataset listing all the species whose ranges overlapped each grid was transformed to a corresponding incidence-based site × species dataset using the ‘pivot_wider’ function from the ‘tidyr’ R package version 1.3.1 [94], resulting in a matrix-like structure, with each row a unique site (grid) and each column a unique species.

Phylogenetic reconstruction

The molecular phylogeny for estimating phylogenetic diversity indices was constructed using Cytochrome b gene sequences downloaded from GenBank [95]. Single sequences retrieved for each species in the community dataset were aligned using MAFFT v7.511 [96]. Phylogenetic analyses were performed using maximum likelihood (ML) to explore the topology of species associations and Bayesian inference (BI) to construct a time-calibrated tree as input for estimating the phylogenetic diversity indices. Both analyses were implemented with the GTR+I+G substitution model selected as the best fit under a Bayesian inference criterion in ModelFinder [97]. The ML analysis was performed in IQ-TREE v 2.3.2 [98], where branch support was estimated from 100,000 ultrafast bootstrap replicates [99]. The BI analysis was run in BEAST v2.7.6 [100] using 100 million Markov chain Monte Carlo (MCMC) rounds sampled every 10,000 intervals. For time calibration, we used log-normal priors based on the most recent common ancestors of the three orders represented in the species list, with dates derived from TimeTree [101], which included Rodentia (μ = 4.231, σ = 0.0354), Eulipotyphla (μ = 4.193, σ = 0.0345), and Macroscelidea (μ = 3.754, σ = 0.1694). The BI results were visualized in Tracer [102] to assess sampling adequacy, with sample size values greater than 200 considered acceptable. The final tree was annotated using the maximum clade credibility method in TreeAnnotator [103] with a 10% burn-in. For species that were present in the community matrix but missing sequence data in GenBank, we used the ‘add.species.to. genus’ function in ‘phytools’ version 2.1–1 [104] to assign them to the corresponding genera subclades [105]. For robustness, we compared the final phylogeny with corresponding tree subsets from the Kumar et al. [101] and Upham et al. [106] phylogenies.

Species trait assembly

To ensure that the traits used for estimating functional diversity indices were well-defined and measurable properties [107], we selected traits representing external morphology (in grams) and diet and activity patterns (as integer percentages of ten diet categories) and three levels of activity patterns (binary counts) from Jones et al. [108], Wilman et al. [109], and Faurby et al. [110]—see Additional file 1 Table S1. These traits are related to the species’ ecological strategies through anatomy, life history, diet, and activity adaptiveness [108, 109]. Body mass generally surrogates external morphological characterization, whereas diet and activity traits surrogate ecological strategies; all are rooted in ecosystem functioning, and are widely used to derive functional diversity indices [85, 111, 112].

Inferring diversity indices

We characterized communities using species richness and several phylogenetic and functional diversity indices to better understand how they varied across spatial grains. The selection of diversity indices was guided by the study’s aim to capture both the breadth of evolutionary history and the range of functional traits within montane small mammal communities in addition to taxonomic diversity, thus providing a comprehensive view of elevational diversity gradient patterns across the range of spatial grain sizes. For each mountain, across the ten spatial grains, we estimated five alpha diversity indices [species richness (SR), phylogenetic diversity (PD), functional diversity (FD), phylogenetic divergence/community structure—phylogenetic mean nearest taxon distance (PD^MNTD), functional divergence/community structure—functional mean nearest taxon distance (FD^MNTD)] and three multisite beta diversity indices [taxonomic, phylogenetic, and functional]. The datasets were first curated to ensure that the species represented across mountains (i.e., in the species composition matrix) matched those in the phylogenetic and trait dataset using the ‘match.phylo.comm’ and ‘match.phylo.data’ functions in ‘picante’ version 1.8.2 [113]. The SR was estimated as the presence-absence sum of unique species at a grid using the ‘specnumber’ function from ‘vegan’ version 2.6–4 [114]. The PD index was estimated as the sum of the total length of the branches of the phylogenetic tree connecting all species represented within a grid, which reflects the species aggregate evolutionary history [115], implemented with the ‘pd’ function in ‘picante.’ The FD index was estimated as the amount of functional space filled by species present within a grid, reflecting the range of traits present, with higher values indicating a greater diversity of functional traits, i.e., the ecosystem supports a wider variety of ecological roles or functions, implemented using the functional richness index of Villeger et al. [62]. The FD was estimated using the ‘dbFD’ function in ‘FD’ version 1.0–12 [116]. The PD^MNTD was estimated from the average distance between each species and its nearest relative in the community, i.e., the mean phylogenetic distance between each species in a grid and its nearest phylogenetic neighbor within the same grid, thus quantifying the average closeness or dispersion of species in the phylogenetic tree with a focus on the smallest distances among species [117, 118]. Low values suggest that species are more closely related to each other on average, indicating a more clustered community in terms of phylogenetic relationships, and higher values indicate that species are more distantly related on average, suggesting a more dispersed community. The PD^MNTD was estimated using the ‘ses.mntd’ function in ‘picante’ version 1.8.2 [113]. We estimated FD^MNTD in a similar manner to the PD^MNTD approach by replacing the input phylogeny with a dendrogram derived from the functional trait matrix. In this sense, the dendrogram represented the (dis)similarity of species based on the functional traits (hierarchical clustering of species based on their functional traits rather than their evolutionary relationships) [117, 119]. The FD^MNTD, thus, reflected the average similarity in functional traits among the nearest neighbors within a grid: a low FD^MNTD suggests that species in the studied community are more closely related in terms of functional traits, indicating greater clustering in terms of ecological functions, while higher values indicate that species are more distantly related on average, suggesting a more dispersed community [117]. We estimated beta diversity indices (taxonomic, phylogenetic, and functional) based on multiple-site dissimilarity to provide insights into how species composition varied among grids, which integrates both turnover—changes in species identities—and nestedness—the degree of subset relationships among communities [120, 121]. The analysis was implemented in R package ‘betapart’ version 1.6—taxonomic beta diversity was implemented using the ‘beta.multi’ function, phylogenetic beta diversity using the ‘phylo.beta.multi’ function, and functional beta diversity using the ‘phylo.beta.multi’ function [56].

Estimating diversity-environment associations

We selected seven environmental variables as environmental and human-activity determinants of diversity patterns, following an extensive literature review on the theories underpinning diversity-environment associations with a focus on mammalian diversity in tropical regions, coupled with initial data explorations (Figs. S2 & S8). The variables—annual temperature average (TAM), annual temperature seasonality (TAS), annual average precipitation (PAM), seasonality in annual precipitation (PAS), terrain ruggedness index (TRI), annual actual evapotranspiration (AET), normalized difference vegetation index (NDVI), and human ecological footprint variable (HFP)—reflect the key climate, topography, primary productivity, and human activity determinants of animal distribution and diversity patterns, especially in tropical ecosystems [10, 33, 122,123,124]. The PAM, PAS, TAM, and TAS were obtained from CHELSA [125], AET was extracted from Running et al. [126], NDVI was extracted from Didan [127], and TRI was derived from the elevation layer from NASA’s SRTMGL1 v003 1-arc second resolution [128] using the ‘Terrain Ruggedness Index’ tool in QGIS [93]. We represented the human footprint in wildlife areas (i.e., HFP) using the average of the global human modification/footprint values from Venter et al. [129] and Kennedy et al. [130], representing human activity impacts on wildlife biogeography patterns across global to local scales. The variables were transformed using the natural logarithm to normalize the data distribution and standardized by rescaling to a 0–1 range to enhance comparability and improve results’ interpretability. We also quantified the effects of the mid-domain effect (MDE) [131, 132], bounded by geographic coordinates and elevation limits of the studied mountains, and based on observed species richness, and computed the magnitude of deviation between the observed species richness and the richness predicted under MDE null expectations using paired t test or a Wilcoxon signed-rank test.

We estimated association dynamics between diversity indices and the predictor variables using generalized mixed additive regression models (GAMM) with controls for spatial autocorrelation [133, 134]. We preferred GAMM because of its robustness in handling the nonlinear relationships that are typical of species diversity and distribution patterns [135], especially considering the range of spatial scales spanned in our study. Spatial autocorrelation in models was first assessed using the ‘moran.test’ function from the ‘spdep’ package to calculate Moran’s I [136] and then addressed based on geographic coordinates (longitude and latitude) [137][133, 134]. To handle spatial dependency in model residuals, incorporating spatial coordinates into a Gaussian spatial covariance function (i.e., as implemented in the syntax "corSpatial (form = ~ Longitude + Latitude, type = 'gaussian')" using the gamm function of the mgcv package) simulates the correlation between observations using their spatial distances, ensuring that the residuals appropriately account for systematic effects of observations’ spatial proximity. We also introduced random effects based on geographically structured grid clusters to allow the model to further accommodate any unmeasured heterogeneity between grid clusters that might otherwise confound the fixed effects. For this, we first used Euclidean distances between the environmental data and Haversine geographic distances to identify grid clusters using the ‘partitioning around medoids’ algorithm implemented in the ‘pam’ function of R package ‘cluster’ [138, 139]. By integrating these spatial components, we adapted the inherent spatial structure of the data [133, 134, 137] to enhance the reliability of parameter estimates [140,141,142,143,144,145,146,147,148]. Models were fitted with the Gaussian distribution and identity link function due to the distribution fit to our data and the restricted maximum likelihood for its unbiased estimates of the covariance and variance parameters in mixed models [149]. The performance was appraised using the adjusted coefficient of determination, R²_a, and Akaike’s information criterion (AIC). The R²_a estimates the proportion of variation in diversity indices explained by environmental variables—the correlation between the observed and predicted outcome values—while the AIC evaluates the model fit (based on prediction errors in a regression analysis) in predicting associations between the diversity indices and environmental variables. Higher R²_a values and lower AIC values indicate that the models better predict associations. These metrics comprehensively evaluate how well the models captured the diversity patterns in response to environmental gradients. The regression analyses were implemented in the ‘mgcv’ package version 1.9–1 [150].

Results

The 124 species represented across the studied mountains spanned three orders, 14 families, and 50 genera (see Additional file 1 Table S1). Most of these species are of low conservation concern according to the IUCN Red List, with 112 listed as least concern, eight as data deficient, three as endangered—golden-rumped sengi Rhynchocyon chrysopygus, Mount Kenya thicket rat Grammomys gigas, and Barbour’s vlei rat Otomys barbouri, and one as vulnerable—the East African highland shrew Crocidura allex. We also identified evolutionarily distinct and globally endangered ‘EDGE’ species [151] on the list: three ‘borderline’ EDGE species—the endangered species (golden-rumped sengi, Mount Kenya thicket rat, and Barbour’s vlei rat) and three EDGE watch list species—the maned/crested rat Lophiomys imhausi, Rudd’s bristle-furred rat Uranomys ruddi, and the naked mole-rat Heterocephalus glaber. Sampling at larger grains always captured more species than at smaller grains, although there were only 1–2 species differences between the 10 km² and 0.0001 km² grain size variations. At the baseline size, 10 km², Mt. Kenya had the highest species representation (72), followed by Aberdare Ranges (71), Cherangani Hills (63), Mt. Elgon (62), Chyulu Hills (58), Matthews Range (53), and Mt. Kulal (27) (see Additional file 1 Table S2). As such, species composition (dis)similarity between mountains remained constant across spatial scales and was only geographically patterned; Chyulu, Elgon, and Cherangani were the most distinct, Elgon was more related to Cherangani, Mt. Kenya to Aberdares, and Mathews to Mt. Kulal (see Additional file 3 Fig. S1).

Spatial scale inherence in diversity-elevation associations

Correlations between diversity indices generally retained a similar topology between spatial grains, albeit variably between mountains and indices (see Additional file 3 Fig. S2). Compared with larger grains, smaller grains yielded fewer species per grid within the same mountain; however, the diversity-elevation topologies were always similar between spatial grains (Fig. 2). Notably, a unimodal peak in species richness at intermediate elevations, aligning with the mid-domain effect (MDE) null model predictions was not observed in most mountains (Fig. 2, see Additional file Figs. S3 & S4). For example, Mt. Kenya showed a sharp increase in SR at low elevations followed by a plateau, nearly similar to the Aberdare Range, which also demonstrated a steady increase in species richness with elevation, consistently across scales; Mount Elgon and the Cherangani Hills exhibited significant variation in SR across scales; Chyulu Hills and Mt. Kulal presented more uniform trends across spatial scales; and Matthews Range displayed distinct divergences in richness across spatial scales, particularly at middle and high elevations; with an intermediate peak in species richness emerging more observedly when all the mountains were combined (Fig. 2, see Additional file 3 Fig. S3 & S4). Deviations from the MDE expectations were most pronounced at smaller spatial scales (see Additional file 3 Fig. S3).

Fig. 2

The influence of spatial grain on diversity‒elevation relationships in Kenya’s montane ecosystems. The figure displays the variation in diversity‒elevation patterns across different spatial grains within montane ecosystems. The curves derived from local polynomial regression models illustrate the association between diversity and elevation without confidence intervals to facilitate easier visualization. The color gradient from yellow (0.0001 km²) to dark blue (10 km.²) represents the range of spatial grains used for sampling, following a strictly nested quadrat design. For the corresponding linear associations, see Additional file 3 Fig. S4

The phylogenetic diversity (PD) patterns were mirrored in those of SR, while the functional diversity (FD) indices’ patterns were unique, being the only index with notably unsystematic distributions across spatial grains despite smaller grains always yielding lower FD than the larger grains (Fig. 2, see Additional file 3 Fig. S4). The distribution of the divergence indices (PD^MNTD and FD^MNTD) were generally comparable across different spatial grains and were only uneven between mountains: both indices decreased as elevation increased, except in the Chyulu Hills and Mathews Range. In the combined mountains’ dataset, PD and SR portrayed hump-shaped curves, PD^MNTD inverted hump-shaped curves, and FD and PD^MNTD linearly decreased, all peaking and transitioning ca. 2,500–3,000 m elevations, but with similar patterns across spatial grains (Fig. 2). Linearly, the smaller grains yielded lower values than the larger grains, and diversity-elevation patterns were homogenous across spatial grains, except for PD^MNTD and FD^MNTD in the Cherangani Hills where they increased with elevation at the larger grains, transitioning at km² to a decreasing trend (see Additional file 3 Fig. S4). The SR and PD increased with increasing elevation, except at Mt. Kulal, where they declined, and the PD decreased at Mt. Elgon; the FD increased at Mt. Kenya, Aberdares, and Mathews but decreased at other sites; and the PD^MNTD and FD^MNTD decreased at Mt. Kenya and Aberdares, increased at Chyulu, and slightly decreased or remained unchanged in the rest of the mountains (see Additional file 3 Fig. S4). In the combined dataset, SR and PD increased as elevation increased, while FD, PD^MNTD, and FD^MNTD decreased uniformly across spatial grains (see Additional file 3 Fig. S4).

For beta diversity, communities were more heterogeneous (high turnover and low nestedness) at smaller spatial grains, with decreasing trends from small to larger grains for turnover and total beta diversity but increasing nestedness (Fig. 3). However, the diversity patterns across spatial grains were not resolved across mountains, i.e., in a combined dataset, the turnover, nestedness, and total beta diversity trends were distinct between individual mountains, with no unifying trend between mountains (Fig. 3).

Fig. 3

Influence of spatial scale on multisite beta diversity in Kenyan montane ecosystems. The figure illustrates how variations in spatial grain affect beta diversity patterns, with beta diversity decomposed into nestedness and turnover components contributing to overall beta diversity. The lines and shaded areas (95% confidence interval) represent the generalized linear regression model relationship between spatial grain size (x-axis) and changes in beta diversity (y-axis), with the size of each point indicating the corresponding spatial grain

Effects of spatial scale on diversity-environment relationships

The diversity-environment models’ adjusted coefficient of determination (R²_a) increased while the Akaike’s information criterion (AIC) decreased as spatial grain increased across mountains and diversity indices (Fig. 4, see Additional file 1 Table S3). The proportion of diversity variances explained (R²_a) was consistently higher at the larger grains; the strength of the diversity‒environment associations generally increased with increasing spatial grain size (Fig. 4). Similarly, individual predictor variables portrayed stronger correlations with diversity indices at larger spatial grains and weaker correlations at smaller grains (see Additional file 1 Table S3). When grids were portioned into elevational bands (three bands in each mountain (lower, middle, and top) except on Mt. Kulal, where only two bands were feasible) there was no discernible variation in the strength of the associations between diversity patterns and environmental variables across spatial grains. However, these associations were generally stronger (lower AIC and higher R²_a) at the top elevation bands, followed by the middle and lower bands across mountains and diversity indices (see Additional file 3 Fig. S5).

Fig. 4

Influence of spatial grain on diversity-environment associations in Kenya’s montane small mammal communities. The figure illustrates how the strength of diversity‒environment relationships, as quantified by Akaike’s information criterion (a) and adjusted R² (b), varies with changes in spatial grain. Each point’s size represents the spatial grain size (x-axis), with the diversity indices displayed on the y-axis. The smoother lines were derived from generalized linear regression models and are used to highlight trends across spatial grains

Effect of spatial scale on human footprint correlation with diversity patterns

The strength of associations between diversity patterns and the human ecological footprint (HFP) increased as spatial grain size increased, with some deviations between mountains and indices (see Additional file Figs. S6 & S7). For example, on Mt. Elgon, the AIC increased with spatial grain size, except for the FD^MNTD dataset and the combined dataset, where the AIC increased for SR, PD, and PD^MNTD and decreased for FD and FD^MNTD. The R²_a increased with spatial grain across indices and datasets except when trends were blurred (see Additional file 3 Fig. S6). Individually, SR, PD, and FD decreased as HFP increased, except at Mt. Kulal and Mathews, and were consistently higher at larger grain sizes. Overall, the topology of indices’ associations with HFP was similar across spatial grains, except when these associations transitioned from positive to negative, or there was no association. For instance, in Mt. Kenya and Aberdares, SR, PD, and FD were negatively associated with HFP, and PD^MNTD and FD^MNTD were positively associated, while in Elgon and Chyulu, the diversity patterns were negatively associated with HFP, except for SR (see Additional file 3 Fig. S7).

Discussion

The elevational diversity gradient patterns across sampling spatial grains observed here extend previous studies’ findings to more local geographical contexts based on a more holistic sampling of mountain systems, a broader continuum of sampling spatial variations, and insights from multidimensional diversity indices.

The range of diversity-elevation patterns (Fig. 2), for instance, concurs with studies such as Rahbek [16] and Kohli et al. [40], who reported that the well-documented mid-elevation peak in richness in most studies of mountains of tropical regions may vary considerably or even reverse when analyzed at different spatial scales or grain sizes, with finer spatial resolutions revealing local biodiversity patterns that coarse-scale analyses might overlook. The deviations between the observed species richness patterns and patterns expected under a mid-domain effect (MDE)-null model were more pronounced at smaller spatial scales, likely driven by local processes such as species interactions, microhabitat diversity, and localized environmental conditions, contrasting with larger spatial scales that tended to smooth out these variations, highlighting broader biogeographic trends [40]. The significant deviations from MDE predictions, such as those observed in the Cherangani Hills and Matthews Range, suggest that while geometric constraints are important determinants of diversity patterns, other factors also significantly shape biodiversity patterns along elevation gradients [40, 152, 153]. Taken together, the absence of systematic mid-elevational peaks in species richness across mountains or consistent diversity‒elevation topologies fits various water‒energy availability hypotheses [39, 40], consistent with commonly reported patterns in small mammal communities of tropical regions [16, 34, 39, 40, 154]. The transitions at lower elevations across spatial scales likely stem from compositional transitions between lowland fauna and true mountain fauna since most of the studied mountains are surrounded by belts of relatively arid lowland savanna bushlands, grasslands, and shrublands and high pressure from human activity [51, 52, 124, 155, 156].

Although there have not been explicit studies examining spatial scale effects at local levels in the Afrotropic region, such as within a single mountain ecosystem, we can draw from other local diversity‒elevation association studies. For example, on Mt. Kenya, Musila et al. [64] and Onditi et al. [52] reported that the species richness of rodents and shrews, in addition to their phylogenetic and functional diversity, peaked around middle elevations but only on the combined dataset from leeward and windward transects. The Onditi et al. [52] study also reported that species richness in the Chyulu Hills, where only one transect was established, unimodally increased as elevation increased. These observations contrast with those of Dreiss et al. [30], who reported that rodent species richness in the Manu Biosphere Reserve, southeastern Peru, decreased with elevation, a trend similar to that of Mt. Kilimanjaro’s bat and amphibian species richness [157]. Across the several mountains included in Kohli et al. [40], the hump-shaped richness‒elevation curve was only evident in the combined dataset, with individual mountains portraying varied patterns. The only Afrotropical mountain in their study (the Rwenzori Mountains) depicted a declining richness‒elevation trend. While geographical contexts (including climate, locality, and human activity pressure) strongly underlie the differences between findings, we also note the different sampling strategies where only single mountainsides were sampled, unlike our approach, which spans three-dimensional mountain ecosystems, in addition to data resolution differences between empirical field investigations and expert-modeled records that interpolate to fill sampling gaps. Ultimately, elevation-diversity relationships are inherently tied to local geographically regionalized conditions, the lateral and vertical extents of the studied elevation gradient, and contrasting patterns of different biodiversity dimensions. Our results also highlight the intricacies of interdiversity generalizations. The phylogenetic and functional indices exhibited notable elevational deviations as the spatial grain size changed. For example, despite generally synchronous diversity-elevation responses for taxonomic and phylogenetic richness, there were notable deviations, such as on Mt. Elgon, where species richness increased while phylogenetic richness decreased as elevation increased, while in Mt. Kulal, both indices decreased toward higher elevations. Compared with the phylogenetic and functional divergence indices, these patterns indicate that high elevations were dominated by closely related species (phylogenetically and functionally), unlike lower elevations, suggesting that conditions at higher elevations favored species with specific traits that lead to relatively homogenous communities in terms of evolutionary diversity and ecological functions over time [158].

We observed stronger diversity‒predictor correlations at larger spatial scales, a pattern explained mainly by the environmental filtering hypothesis—greater environmental heterogeneity (including human variables) within a larger area has stronger effects on species assembly dynamics than does a smaller area [38]. Notably, however, the increase in species richness with increasing elevation suggests that habitat suitability increases with increasing elevation, with no evident deterioration at the peaks, meaning that extreme elevations (lowest and highest) may not be systematically harsher than mid-elevations, in contrast with expectations of geometric constraints and habitat productivity deterioration at these elevations [31, 34, 39, 159]. Similar to Kohli et al. [40], who invoked strong community filtering influences of water availability constraints to explain why the lowlands of wet mountains and highlands of arid mountains harbored the most functionally and phylogenetically diverse rodent communities, our results also suggest that the ecological conditions along the elevation gradients could be increasingly homogenizing from the effects of climate change and human proliferation. These conditions gradually alter the community composition dynamics by favoring some lineages over others, thus driving high richness but not correspondingly increasing the evolutionary or functional uniqueness [160]. Nonetheless, the precise nature of such configuration alterations requires more primary studies over more extended periods. We also found that human footprint (HFP) effects on diversity patterns were stronger at larger spatial grains, but variably between mountains. For example, HFP associations with various diversity indices in Mt. Kenya, Aberdares, Elgon, Chyulu, and Kulal were not uniform and even transitioned with changing spatial grains, indicating that the presence of human activity hindered or, in some cases, facilitated some aspects of the observed diversity patterns, reiterating the complex relationships between human activities and biodiversity in Afrotropical montane ecosystems [124].

Our study also presents some vital conservational insights. We cataloged 124 species across the seven mountains studied, representing more than 84% of the species in these three orders in Kenya [88, 161]. This high representation of the national mammal diversity in these mountains underscores their exceptional contribution to Kenya’s conservation endeavors. Among the recorded species, those of significant conservation concern due to their endangered or evolutionarily distinct and globally endangered ‘EDGE’ status, including Grammomys gigas, Otomys barbouri, Rhynchocyon chrysopygus, and Crocidura allex, exhibit localized distributions, each being endemic to specific mountains such as Mt. Kenya, Mt. Elgon, the Chyulu Hills, and Aberdare Ranges + Mt. Kenya, respectively. While the documentation of these taxa in field surveys often reflects a concerted endeavor to capture fauna of conservation priority, our findings highlight a systematic tendency for small-scale sampling to record less species diversity, potentially excluding vital species from local biodiversity assessments. Conservation strategies should incorporate the scale-dependent nature of interactions between human activities and biodiversity patterns in montane ecosystems, including the uncoupled nature of different metrics. More repetitive field surveys that are spatially wide-ranging and traversing extensive transects could ensure accurate delineation of local species diversity. This is particularly crucial in montane ecosystems whose global biodiversity values are dynamically threatened by anthropogenic encroachment and climatic change impacts.

Notably, because comparing results across studies that apply different spatial and temporal survey scales can lead to inconsistencies in the interpretation of diversity patterns and the influence of environmental factors, any extrapolation of our results to other regions, species groups, or ecosystems should be guided by spatial and geographical similarities.

Conclusion

Our findings suggest that varied field sampling scales generally yield consistent species richness–elevation relationships. However, systematic variations emerge between mountains and diversity indices, highlighting the significant role of geographically regionalized environmental conditions in shaping species diversity and distribution patterns in mountains. While observed montane small mammal diversity patterns seem more strongly influenced by these environmental factors than by spatial grain variations in sampling, evaluating montane biodiversity at single, arbitrary sampling grain sizes may misrepresent the true impact of predictor variables, such as human activities, on species diversity and distribution patterns. There is a need for a more detailed evaluation of elevational diversity gradients and the impact of environmental and anthropogenic filtering across spatial scales, especially when comparing mountains across different environmental contexts. This evaluation is crucial for enhancing sustainable ecosystem management and conservation strategies. Because the ranges of many species inhabiting mountains surrounded by relatively drier lowland landscapes, such as Mt. Kulal, Mathews Range, and Chyulu Hills, are projected to retreat to the peaks where they become trapped and locally extinct over time because of the effects of pressure from human activity interacting with broader climate change impacts [43, 83], future studies should incorporate multiple elevationally-banded transects and employ denser traplines to uncover rare or new taxa, as well as temporal surveys to predict the shifting community assembly dynamics with better precision.