How Lithology Impacts Global Topography, Vegetation, and Animal Biodiversity: A Global‐Scale Analysis of Mountainous Regions

Chemical and mechanical properties of lithology exert a first‐order control on landscape evolution and biological colonization of substrate. To quantify the influence of lithology on topography, vegetation density, and animal biodiversity, I compile lithologic, topographic, climatic, and biological data sets for mountainous regions globally. I find significant variations in the topographic steepness of regions underlain by different lithologies that, accounting for tectonic uplift, likely reflect lithologic differences in erosional resistance. These relative differences in erodibility are similar across different climate zones. To isolate the effect of lithology on vegetation and animal biodiversity, I account for the heterogeneous lithologic distribution among climate zones. I show that siliciclastic, plutonic, and, for some biological variables, metamorphic rocks exhibit elevated values of Normalized Difference Vegetation Index and tetrapod and amphibian species richness relative to carbonate rocks. These results likely reflect lithology‐related variation in soil nutrients and hydrology that promote or inhibit habitat suitability.


Introduction
The Earth's surface exposes a mosaic of bedrock lithologies, which span a large range of chemical and mechanical properties. Lithologic characteristics control how landscapes respond to erosion and weathering and set boundary conditions for organisms living on the substrate. Rock resistance to erosion and weathering has a fundamental impact on the topographic form, where more resistant rocks form steeper parts of a landscape, under otherwise equal conditions (Dane, 1935;Snyder et al., 2000). It follows that topographic steepness metrics can be used as a proxy for erodibility of different lithologies in a steady-state landscape with uniform climatic and tectonic conditions (Gallen, 2018;Mills, 2003). Similarly, the composition and texture of rocks greatly affect hydrology, soil fertility, and soil chemistry (Anderson, 1988;Shaler, 1892) and thus affect the organisms colonizing the substrate (Braun-Blanquet, 1932;Hahm et al., 2014;Morford et al., 2011;Unger, 1836). For instance, carbonate rocks are more susceptible to drought as they commonly have higher surface water infiltration rates (Jiang et al., 2020), and, generally, their soils have lower nutrient content compared to other rock types (Porder & Ramachandran, 2013;Tyler, 1996). Other processes, such as dust deposition, also affect nutrient input. However, since these processes are more impactful in depositional landscapes, lithology is expected to be an important control on biological parameters in erosion-dominated landscapes, such as mountains. Dürr et al. (2005) compiled a global lithologic map and analyzed the relative frequency of different rock types in lowlands, mountains, and plateaus. The study showed that lithologies traditionally viewed as erosionally resistant, such as volcanic and plutonic rocks, are overrepresented in high relief areas, while alluvial and poorly indurated sediments are overrepresented in lowlands and sedimentary basins. However, this study was unable to make more detailed topographic comparisons, for instance, between lithology and local slope, due to the low map resolution (0.5°× 0.5°, ∼55 km at equator). Recently, a higher-resolution world lithologic map (∼500 m target scale) was published (Hartmann & Moosdorf, 2012) that now enables a worldwide examination of relationships between rock type, topography, climate, and biological parameters at unprecedented resolution.
Erodibility is a variable used in many topographic analyses (e.g., river profile analysis), yet its quantification has proven difficult due to local and regional variations in climate and tectonics (Snyder et al., 2000;Stock & Montgomery, 1999). Recent studies investigating global-scale interactions between climate, tectonics, and topography have proven useful in elucidating coherent trends for such problems (Chen et al., 2019;Harel et al., 2016). All else being equal, more competent rocks form steeper slopes compared to less competent rocks (Gabet, 2020b;Gilbert, 1877;Hack, 1960), but we still lack quantification of this effect on a global scale. Similarly, the effect of bedrock lithology on vegetation has been long been established (Unger, 1836) but has mostly been studied at local or regional scales (Hahm et al., 2014;Wentworth, 1981;Whittaker & Niering, 1968). Moreover, differences in the local vegetation related to bedrock properties (e.g., hydrology) might also influence the biodiversity of animals.
Precipitation and temperature are good predictors of tetrapod biodiversity patterns (Antonelli et al., 2018). Relief and "topographic complexity" (an inconsistently defined but commonly used term in biodiversity literature referring to regions comprised of diverse geomorphic domains) also correlate with certain metrics in biodiversity data sets (Antonelli et al., 2018;Badgley et al., 2017). However, this correlation is highly dependent on the range at which relief is measured and may arise as a consequence of incorporating several ecological niches along altitudinal transects in areas of high relief into a single data point (Kerr & Packer, 1997). Additionally, the number of soil types in a given area appears to be a weak predictor of biodiversity (Antonelli et al., 2018), but the effects of specific soil types or the underlying substrate on biodiversity have not been explored globally. Soil formation and denudation style are modulated by rock type, such that the impact of climate on these processes might be rock-type specific. Similarly, evidence suggests that climate can influence hillslope steepness (Gabet et al., 2004), and hence, some lithologies might be more or less erosionally resistant under given climatic conditions.
In this study, I compile lithologic, topographic, climatic, and biological data sets to investigate the effects of bedrock variations on topographic metrics and biological parameters. To elucidate differences in relative erodibility among main lithologic units, I examine relationships between lithology and topographic metrics, such as hillslope gradient and river channel steepness. I account for climatic differences worldwide and their influence on vegetation and species richness. I then compare nonclimatic variations in biologic parameters between lithologic groups to evaluate the impact of substrate properties on the distribution of biota.

Data and Results
To assess the influence of lithology on topography and biological parameters, I limit my analysis to erosion-dominated landscapes, here defined as areas with ≥50 m of topographic relief within a 3 km radius ( Figure 1a, for details see Text S1 in the supporting information). This relief threshold effectively excludes noneroding regions, where sedimentation dominates, and bedrock erodibility has a negligible influence on topography. Focusing on mountainous and hilly areas is also important to relate biological parameters to the properties of the underlying lithology, because soils in areas of high erosion more reliably reflect the properties of their parent rock and are less likely to be influenced by other processes, such as dust input (Arvin et al., 2017).

Lithologic Variations in Mountainous Topography Worldwide
I measure elevation and slope from a digital elevation model (DEM) with 500 m resolution and subsequently resample both bilinearly at 1 km resolution (see Text S1). A sensitivity analysis shows that the relative differences in slopes are not affected by the coarse DEM resolution ( Figure S1). Due to the sensitivity of river channels to bedrock erodibility, I utilize a global normalized river channel steepness (k sn ) map (Hilley et al., 2019), a measure of stream gradient normalized for its discharge, which is commonly approximated by upstream drainage area (Wobus et al., 2006). k sn can only be defined along stream channels; however, in the absence of lithologic variation between channels and hillslopes, the gradient of the river can be extrapolated to the drainage divide. This approach eliminates biases introduced by the inclusion of thin covers of unconsolidated sediments at valley floors in the global lithologic map (Texts S1 and S2 and Figure S2). Using interpolated k sn values introduces uncertainties from global variations in concavity (Chen et al., 2019), drainage area-discharge relationships, and drainage divide location uncertainty (see Text S2). Metrics of topographic steepness not only depend on erodibility but also on local tectonic rock uplift rates, with areas of higher tectonic uplift generally forming higher relief (Burbank et al., 1996). Since there are no global maps of tectonic uplift rates, I use the Global Strain Rate Model (GSRM) (Kreemer et al., 2014) to group the data into tectonically active and inactive areas, based on the assumption that geologic strain rate records deformation that may lead to rock uplift (Price, 1975). The GSRM describes rates of horizontal strain accumulation; hence, the values cannot directly assess vertical tectonic rock uplift. However, based on previously defined areas of deformation by Kreemer et al. (2014), I define areas with the second invariant of strain rate >0.5 × 10 −9 /a as being tectonically active, and the remainder as inactive.
Topographic variables in the data set, such as elevation, slope, and normalized river channel steepness, vary significantly among lithologies within different tectonic zones ( Figure 2). In tectonically active areas, plutonic rocks occupy the highest elevations and have the highest slopes and k sn , with median values of 8.8°and 200, respectively ( Figure 2a). Metamorphic rocks form the second steepest slopes and stream channels, followed by mixed and pure carbonates that have median slopes of~6.2°and k sn > 150, with some overlap of the 95% confidence intervals of the means. Volcanic rocks occupy higher elevations than metamorphic, pure and mixed carbonate rocks, but form slightly less steep topography. The lowest elevations, as well as slopes, and k sn occur within siliciclastic sediments and, unsurprisingly, in unconsolidated sediments.
In tectonically inactive areas, the absolute values of topographic metrics are lower than in tectonically active settings, owing to low rates of rock uplift ( Figure 2a). However, most of the relative trends between lithologies remain consistent. The only notable changes are that, in tectonically inactive regions, (1) carbonate rocks tend to be steeper than metamorphic rocks and (2) volcanic rocks occupy the highest elevations and steeper than average terrain (Figure 2a). High and steep volcanic regions in tectonically inactive areas highlight the ability of volcanoes to form topography in areas with low uplift rates.
Climate is commonly thought to influence erodibility of rocks, and this behavior might be rock-type specific. However, the relative differences in mean slope between lithologic units are similar across different Köppen- Figure 2. (a) Boxplots and comparisons of group means for topographic data in tectonically active (upper panels) and inactive (lower panels) areas. Boxplots show the full data range with quartiles and median. Boxplots for the inactive areas have the same scale as active areas to visualize the effectiveness of the global strain rate model as tectonic uplift rate proxy. Above every boxplot, a separate plot shows the deviation of the mean of every group from the mean of pure carbonates (black dashed line) (see Text S3 for details on the statistics), with error bars denoting the 95% confidence interval. (b) Same as in (a) the means of lithologic units with 95% confidence intervals are plotted for main Köppen-Geiger climate zones. Note due to small sample size polar regions were excluded.
Geiger climate zones (Figure 2b), suggesting a relatively weak or imperceptible influence of climate on lithologically controlled topographic differences.
Topographic variables differ substantially between lithologic groups, relative to their global median values. For example, there is an~7°difference in median hillslope angle between the lowest (unconsolidated sediments) and the highest (plutonic rocks) lithologic group within tectonically active areas, whereas the median hillslope angle of all tectonic areas is 5.4°. The range of medians in river channel steepness within tectonically active areas is ∼120, while the total median is ∼150. This pattern also holds for tectonically inactive regions.

Lithologic Variations Among Biological Variables Worldwide
To investigate if lithology contributes to variations in vegetation density and animal biodiversity, I compile Normalized Difference Vegetation Index (NDVI) 1 km resolution data from satellite imagery (see Text S1 and Table S1) and 10 km resolution species richness data (Jenkins et al., 2013;Pimm et al., 2014). I follow the approach of Antonelli et al. (2018) and sum the number of bird, mammal, and amphibian species per grid cell to a tetrapod species richness. Amphibians are investigated separately because they require water for reproduction (Rodríguez et al., 2005;Semlitsch, 1996) and are therefore likely sensitive to hydrological properties of the substrate.
Mean annual temperature (T) and mean annual precipitation (P) impart strong control on animal richness and plant cover (Hawkins et al., 2003). Therefore, it is important to account for the influence of these two key variables, for instance by binning the data into different climatic (T-P) groups (Figure 3a and Text S4), to isolate the influence of lithology on biological parameters. Biological parameters of NDVI and tetrapod and amphibian species richness show variations between lithologies after T-P binning (Figure 3a). The strongest differences in NDVI among lithologies are observed between carbonate areas, which tend to have low NDVI values compared to siliciclastic and plutonic areas. For several T-P bins, the difference between carbonates and plutonic or siliciclastic rocks is about one standard deviation of the total data set within this bin. Similar, but more significant, trends are observed between tetrapod and amphibian richness and rock type. (Figures 3a and S5 for significance values). Regions underlain by plutonic and metamorphic rocks exhibit higher tetrapod and amphibian species richnesses than those underlain by carbonate rocks (Figure 3a). Volcanic regions have average values for NDVI and tetrapod richness, but exhibit low amphibian richness.
As a second method of accounting for T-P variations, I apply different regression models for T and P to the three biological variables (Table S2). I find that in multiple linear regression models, T and the logarithm of P (saturation effect with high precipitation) are good predictors for biological variables. The best fit linear regression model for NDVI, using the predictors T and log(P), has an adjusted R 2 value of 0.59 (Table S2). For tetrapod and amphibian richness, the best fit regression models combine log(P) and the interaction term of T and log(P). Interactions in regression models state that the influence of one predictor (here P or log(P)), on the response variable (here biological parameter), depends on the value of another predictor variable (here T). The Tukey-Anscombe diagnostic plots of the multiple linear regressions, which show the residuals over the fitted range ( Figure S6), indicate systematic biases for all biological variables. Hence, I choose to fit Generalized Additive Models (GAMs) to allow for a more complex response of the predictors with respect to the dependent variable. These models fit a spline or 2-D tensor product surface to the predictor or predictor interaction term, respectively. Out of the tested models, the interaction term of T and P performs best in the regression and can account for 69%, 65%, and 58% of the variance in NDVI, tetrapod richness, and amphibian richness, respectively (Table S2). Surface contour plots of T-P interactions for the best fit models are displayed in Figure S6.
I utilize the residuals of my best fit GAM to examine if variations exist between the lithologic groups ( Figure 3b). An analysis of variance and Tukey Honest Significant Difference test show significant variations in the mean residuals of each lithologic group, similar to the independently derived trends in Figure 3a. This analysis also indicates that NDVI is lowest in carbonate areas and highest in siliciclastic and plutonic terrain and that tetrapod and amphibian species richnesses are low in carbonate landscapes, but high in metamorphic and plutonic areas. On average, plutonic and metamorphic terrains host~30 more tetrapod species (total mean 301, median 243) and four more amphibians (total mean 15, median 8) within a square kilometer compared to carbonate areas.
Land cover data show a similar pattern; carbonate regions host higher than average fractional areas of intermediate vegetation cover and no vegetation cover ( Figure S8). Within vegetated areas, carbonate regions also appear to have a below average percentage of broadleaf forest cover and a slightly above average percentage of needleleaf forest cover, though these variations may arise from differences in climate, anthropogenic disturbance, or other unmeasured factors that differ between lithologies.

Discussion
Grouping Earth's diverse rock types into only seven lithologic categories has limitations. While some lithologic groups, such as pure carbonates, should exhibit rather homogeneous properties (chemical composition, texture, etc.), others, especially metamorphic rocks, are more heterogeneous. This heterogeneity can The data are binned by T and P into 16 (4 × 4) equally sized bins (n = ∼5e7) for each lithologic group. The bins are colored by their deviation from the mean NDVI or species richness of their entire T-P bin across different lithologies (see Text S4), scaled to the standard deviation of all the data within this bin. Hence, red colors indicate that a lithologic unit mean is half a standard deviation (sd) below the mean NDVI or species richness of all terrain in the same climatic (T-P) regime. The number of pixels per bin is plotted in the lower right panel. (b) Comparison of the means for T-P-corrected NDVI, tetrapod richness, and amphibian richness relative to the reference group (pure carbonates). The three biological variables are corrected for T-P conditions by a GAM. The plots show the difference in means between groups with 95% confidence intervals; the dashed zero-line corresponds to the respective mean in pure carbonate areas. See Text S4 for the method details. amplify differences in topographic and biological variables within a lithologic group. Thus, the heterogeneity within groups should be kept in mind when interpreting this analyses, and it suggests caution should be taken when applying the findings of this study to regional settings.

Differences in Erodibility Impact Topography Globally
The topographic analysis shows significant variations in elevation and river and hillslope steepness between lithologic groups (Figure 2). Regions underlain by plutonic rocks have the steepest hillslopes and river channels in both tectonically active and inactive areas, presumably due to their high resistance to erosion leading to more steepening of a landscape in response to rock uplift, compared to less resistant lithologies (Howard, 1994). Unconsolidated sediments, on the other hand, have low erosional resistance and form the low-lying, less steep regions of the landscape. Similarly, the low steepness of regions underlain by siliciclastic rocks may result from their lower inferred erosional resistance; whereas the erosional resistance of mixed carbonate and pure carbonate sediments is inferred to be higher, similar to that of volcanic rocks. It is interesting to note that metamorphic rocks form steeper hillslopes and stream channels in tectonically active areas compared to carbonate rocks but are generally less steep in tectonically inactive areas. This could be related to a bias in the exposure of different lithologies in certain tectonic regions; metamorphic rocks are commonly exposed in the interior of mountain ranges with typically high uplift rates, whereas carbonates are often removed earlier during orogeny. This could suggest higher-average erosional resistance of carbonate rocks compared to metamorphic rocks, or it may indicate that carbonate mountains are steeper due to higher rates of surface water infiltration, which can reduce discharge and stream power and therefore cause channels to steepen (Ott et al., 2019). Overall, these topographically inferred differences in erodibility are consistent with the findings from previous local-and regional-scale studies (Gabet, 2020a;Garcia-Castellanos & O'Connor, 2018;Harel et al., 2016;Kühni & Pfiffner, 2001).
Climate is thought to affect bedrock erodibility and hillslope steepness (e.g., Ferrier et al., 2013;Gabet et al., 2004). Yet, dividing slope data into climate zones (Figure 2b) reveals no clear pattern in the relative erodibilities of different rock types. Most variation is seen for metamorphic rocks; however, the trends are inconsistent between tectonically active and inactive regions. Carbonates appear steeper in dry regions, which matches the expectation of carbonates being less erosionally resistant in regions of elevated chemical weathering, yet this interpretation remains ambiguous as the trend is weak in tectonically active regions. Together, these observations suggest that (1) climatic influences cannot be resolved, for example, due to oversimplification in the classification scheme or too much influence of individual mountain ranges within certain subsamples; (2) the rock properties governing differences in erodibility, such as tensile strength (see Bursztyn et al., 2015;Sklar, 2001, for a discussion of erodibility) are not subjected to major differences between climate zones; or (3) climate induces more complex interactions between variations in mechanical properties and chemical weathering.
By splitting the data set into tectonically active and inactive areas, I lowered the risk that a heterogeneous distribution of lithologies, in respect to tectonic uplift rates, would impact the analysis. Unfortunately, the record of thermochronologic exhumation rates is spatially and temporally limited, and it is currently impossible to remove this potential bias completely from the data set. However, the large differences in elevation and steepness between tectonically active and inactive areas (Figure 2) suggest that global strain rates provide a reasonable proxy for tectonic uplift.

Lithologic Influence on Vegetation and Biodiversity
To account for climate differences when examining the relationship between lithology and biological parameters, I binned the data by T-P conditions. This simple approach shows some significant differences between lithologic groups (Figure 3a). Testing of different GAM regression models shows that the interaction of T and P typically does a better job of explaining variance in the data set compared to T and P as separate predictors (Table S2). This observation is conceptually sound; T and P are positively correlated with all biological variables. However, a subtropical desert can have high T and will still have low values for the biological variables due to lack of water. The effects of T and P depend on each other and, therefore, the interaction term of both produces the best fit that explains almost 70% of the total variance in NDVI. An independent test using the residuals of a GAM fitted to the interaction of T and P yields the same principal results for the variation among lithologic groups.
Variations in tetrapod and amphibian species richnesses between lithologic groups may be slightly biased by the lower resolution at which these data are published. Moreover, differences in data availability and coverage might also bias the measured differences in biological parameters (Faurby & Svenning, 2015). Anthropogenic effects could also affect NDVI and species richness measurements, though they have been shown to disproportionately influence the diversity of large mammals and amphibians (IUCN, 2020;Schipper et al., 2008) rather than the full suite of tetrapod species examined here. However, variations in satellite-derived high-resolution NDVI between lithologic groups are consistent with differences measured from field-derived range maps of species, verifying the robustness of the analysis despite the heterogeneous nature of the data set and anthropogenic disturbance.
Although the observed variations of biological variables between lithologic groups are correlative and not necessarily causative, there are good reasons to think that they may indeed be causative. Lithology is an important control on plant nutrient supply because it partially determines weathering products, soil pH, and hydrologic conditions (clay formation, fracture density, and infiltration) (Huggett, 2002;Neff et al., 2006). Phosphorous concentrations, which limit plant growth in some environments, are high in siliciclastic rocks and low in carbonate rocks (Porder & Ramachandran, 2013), which corresponds well with my observations of NDVI variations (Figure 3). Carbonate terrains often set challenging conditions for vegetation growth due to low nutrient content, low clay content in the soil, and a high susceptibility to water infiltration into the groundwater system (Kruckeberg, 2004;Wentworth, 1981;Whittaker & Niering, 1968). Lower clay content in carbonate areas decreases the ability of the regolith to retain water for plant growth, facilitates water infiltration to groundwater, and can thereby lower primary productivity (Jiang et al., 2020). In the Arizona desert, plutonic terrains have higher vegetative species richness, more than double the net primary productivity, higher growing vegetation, and flora that is better adapted to wet conditions than regions underlain by carbonate bedrock (Wentworth, 1981;Whittaker, 1970;Whittaker & Niering, 1968). A recent study in China demonstrates that primary productivity increases with increasing silicon and decreasing carbonate bedrock content (Jiang et al., 2020). These local observations are in agreement with the findings here, which indicate particularly high NDVI in plutonic and siliciclastic areas relative to carbonate regions. Carbonate regions also exhibit below average percentages of broadleaf forest and intermediate vegetation, yet slightly above average needleleaf forest coverage ( Figure S8). This might be related to needleleaf forest occupying generally drier regions, provided these measurements are not unduly biased by an overriding influence of temperature on the distribution of needleleaf forests (Woodward et al., 2004). The local response of vegetation to substrates will differ due to various other parameters, for example, varying bedrock composition and texture, as well as local climate (Whittaker, 1970). The high NDVI on plutonic and siliciclastic rocks is likely related to their higher-average nutrient content and higher clay content, which facilitates the retention of water available to plants (Hahm et al., 2014;Porder & Ramachandran, 2013).
The susceptibility of the biosphere to bedrock properties will also differ between regions since other factors influencing soil formation will be of varying importance. For instance, dust is an important contributor to soil nutrient content in some regions (Arvin et al., 2017;Okin et al., 2004). Also, erosion influences the nutrient supply by providing fresh minerals to soils (Porder et al., 2007). I focus this analysis on hilly and mountainous areas to avoid stable, noneroding areas, like the Amazon basin, where nutrient supply may be decoupled from the underlying lithology (Porder et al., 2007).
Tetrapod and amphibian species richnesses generally mimic the variations in NDVI between lithologic groups (Figure 3). Vegetation differences between lithologies likely impact the local fauna, which rely on vegetation as their primary food source. Amphibians are a small subset of the tetrapod class but are underrepresented on carbonate and volcanic bedrock (Figure 3b). Amphibians are particularly sensitive to hydrologic conditions because they rely on water for reproduction (Rodríguez et al., 2005;Semlitsch, 1996). The high infiltration rates and low surface availability of water in karstic terrains may explain the lower biodiversity of amphibians in such areas, despite a high degree of endemism. Volcanic rocks have diverse hydrological properties (Singhal & Gupta, 2010), making it difficult to directly relate the relatively low number of amphibians observed on volcanic terrain to hydrologic bedrock properties. However, young volcanic rocks commonly have high porosity and permeability, which promote high infiltration and loss of surface water (Jefferson et al., 2010). These findings therefore suggest that substrate hydrology may indeed affect the abundance and diversity of flora and fauna.

Conclusions
Analyses of topographic, climatic, biological, and geological data in mountainous regions globally show significant differences in topographic steepness, climate-corrected NDVI, and tetrapod and amphibian species richnesses between regions underlain by different lithologies. Variations in hillslope and river steepness can be interpreted as differences in erodibility among main lithologic units and highlight, for instance, the high erosional resistance of plutonic rocks.
Similarly, variations between residuals of climate-corrected NDVI indicate that increased vegetative growth in regions underlain by plutonic and siliciclastic bedrock is likely due to higher clay content, which facilitates water retention, and higher nutrient content in the soil, relative to less vegetated carbonate regions. These differences are also reflected in variations of tetrapod richness, which likely are due to the animals relying on vegetation as primary food source. The hydrologic bedrock properties may also be the main driver of climate-corrected variations in Amphibian species richness.

Data Availability Statement
The data set compiled for this study and the code used for statistical analysis can be found in the ETH research collection (https://doi.org/10.3929/ethz-b-000411821). The interpolated global-scale k sn map can be found online (https://doi.org/10.3929/ethz-b-000411632).