C4 grass functional traits are correlated with biotic and abiotic gradients in an African savanna

Savanna is a species-rich biome, that includes many modern mammal lineages and C4 grass (Poaceae) species. The greater productivity and grazing pressure associated with savannas is likely attributable to the foliar traits of the grasses they support. Thus, it is important to understand the complex relationships between the abiotic environment, foliar attributes and the floristic composition of savanna grasses, and the supported grazer densities. We sampled 37 grass communities in the Kruger National Park (South Africa) across three soil types and along a rainfall gradient and found that these communities lack strong phylogenetic structure. We then measured specific leaf area and leaf tensile strength for 384 individuals representing 66 species and found that both traits were strongly phylogenetically structured and associated with both rainfall and soil type. Finally, we found that grazer densities in the Park are correlated with the foliar traits of the associated grass communities, but the resolution of our data do not allow for a thorough analysis of this association. Our results demonstrate the complex interactions between climate, soils and grazer densities relative to C4 grass functional traits.


Introduction
Savanna ecosystems cover around 20% of the global land area and are characterized by both a high diversity of ungulates and C 4 grasses (Turpie and Crow 1994;Lehmann and Parr 2016). These ecosystems are of socioeconomic importance because they support a substantial portion of the world's population and most of its rangeland and herbivore biomass (Sankaran et al 2005;Scholes and Archer 1997). It is therefore important to develop a sound mechanistic understanding of the processes underpinning the complex mechanisms that drive and maintain the tremendous biological diversity within these ecosystems.
Either directly for herbivores or indirectly for carnivores, large mammals depend on plant community composition (Augustine and McNaughton 1998). Plant communities are in turn influenced by abiotic environmental gradients. In savannas, for example, there is a strong, positive, linear relationship between vegetation biomass and annual precipitation (McNaughton and Georgiadis 1986). A potential increase in vegetation biomass is, however, controlled by a combination of both herbivory and fire to around 700 mm MAP after which fire takes over as the dominant determinant (Sankaran et al. 2005;Archibald and Hempson 2016). It has, however, also been demonstrated that there is an increase in herbivore biomass as MAP increases to around 700 mm (East 1984;Grant and Scholes 2006;Cromsigt et al. 2009;Hempson et al. 2015;Archibald and Hempson 2016). At our study site in the Kruger National Park (South Africa), the highest rainfall is around 750 mm in the south of the Park (Holdo et al 2018) and as a result we would expect herbivore biomass to increase as rainfall increases. Intensive grazing is, however, a feature of African savanna ecosystems, with grass consumption positively and linearly related to plant productivity (Augustine et al. 2003). In addition, presumably because of high precipitation fosters, high grass productivity, and standing biomass (February et al. 2013), these conditions also promote high fire frequencies (Archibald and Hempson 2016).
Tropical savannas are often associated with nutrient poor soils (Scholes 1990;Medina 1993), especially in higher rainfall areas such as in West Africa and South America (Februrary et al. 2019). Where rainfall amount may be associated with soil nutrient status, a fundamental feature of most savannas is heterogenous soil fertility (Medina 1993;Venter et al. 2003).The greater productivity and grazing pressure associated with high-rainfall, high-nutrient savanna sites are likely attributable to the foliar quality of the grasses they support. Low nutrient and moisture availability are both thought to select for leaves having a higher tissue density than those from moister and/or more fertile sites (Wright et al. 2001;Wright and Westoby 2003;Ordoñez et al. 2009). High-density leaves offer clear benefits in resource-limited sites, since their greater durability and tensile strength increases leaf longevity (Wright et al. 2001) and resistance to herbivory (Franks et al. 2008). By contrast, lowdensity leaves and despite rendering plants more palatable to herbivores (Cingolani et al. 2005), are generally beneficial in resource-rich sites since they facilitate rapid growth, thereby enhancing the competitive ability of plants as well as their ability to tolerate grazing (Cingolani et al. 2005). Such changes in foliar traits between communities experiencing different environmental conditions could be due to two radically different, but non-exclusive, processes: species turnover and within-species trait variation (Lepš et al. 2011). Many plant functional traits display significant phylogenetic signal, with closely related taxa having trait values that are generally more similar than expected by chance (Liu et al. 2015;Flores et al. 2014). In this context, there is reason to expect the environmental correlates of leaf traits to translate into floristic patterning with moist/fertile and dry/infertile environments favouring contrasting suites of lineages. There is some evidence of such patterning in savanna grasses, with the distributions of C 4 grass subfamilies and tribes being predicted by contrasting environmental conditions (Visser et al. 2012) and grazing pressure (Hempson et al. 2015). In particular, where species of Andropogoneae dominate moist, infertile sites that experience limited grazing pressure but high fire frequencies, species of Chloridoideae dominate drier sites on richer soils with heavier grazing pressure and lower fire frequencies (Visser et al. 2012). Paniceae and Aristidoideae occupy intermediate positions, with mesic-adapted Paniceae responding positively to both fire and grazing, and arid-adapted Aristidoideae responding positively to grazing but not fire (Visser et al. 2012).
Against this background, we set out to examine the foliar attributes and floristic composition of savanna grass communities situated on sites of contrasting fertility and precipitation in the Kruger National Park, South Africa (KNP, hereafter). We tested the hypothesis that environmental gradients influence the phylogenetic composition and palatability of grass assemblages because low nutrient and moisture availability are both thought to select for leaves having a high tissue density and lower SLA than those from moister and/or more fertile sites (Wright et al. 2001;Wright and Westoby 2003), and that these in turn are linked to grazer densities. To address this, we measured two important leaf traits: specific leaf area (hereafter SLA), a key trait associated with leaf economics and lifespan that correlates strongly with environmental gradients (Wright et al. 2001), and leaf tensile strength, which is routinely used to determine the palatability of grass species (Martens and de V. Booysen 1968). We then assessed whether (i) there was a correlation between SLA and leaf tensile strength, (ii) phylogenetic signal in foliar traits resulted in grass lineages being non-randomly distributed locally, (iii) environmental variables influenced grass foliar attributes at the community level, and (iv) average grass foliar attributes at the community level correlated with grazer densities.

Field sampling
Our study site is the 20,000 km 2 Kruger National Park located in the north east of South Africa. Rainfall is distinctly seasonal with a wet season from October to May. Mean annual precipitation is lowest in the extreme north (400 mm) and highest at Pretoriuskop in the southwest (750 mm) (Holdo et al. 2018). Since herbivore densities also generally increase from North to South in the park, we expect herbivore biomass to increase as rainfall increases across the KNP. In February 2007, we sampled 37 (10 m 9 10 m) plots distributed across the whole KNP representing a wide range of climatic and geological settings (Fig. 1). For each plot we recorded the full list of grass species and determined the mean SLA and leaf tensile strength using fully formed whole fresh leaves from three distinct individuals of each species. SLA was determined using the protocol of Garnier et al. (2001), while leaf tensile strength was determined as the breaking tension per unit of leaf cross-sectional area (Martens and de V. Booysen 1968). The breaking tension of leaves was measured in the field using a standard tensilmeter (Martens and de V. Booysen 1968), while leaf cross-sectional area was estimated as the product of leaf thickness and width, respectively measured using a digital vernier calliper precise to 0.01 mm and a metal ruler precise to 0.5 mm. Across all plots, we sampled a total of 384 individuals from 68 species, all from the subfamilies Aristidoideae, Chloridoideae and Panicoideae.
We also characterized each plot in terms of its geology, mean annual precipitation (MAP), and grazer density. Soils in the KNP reflect the underlying bedrock and can broadly be divided into basaltic soils in the east and, granitic soils in the west, with a narrow, intervening strip of shale-derived soils (Venter et al. 2003). Combined grazer densities of Buffalo (Syncerus caffer), Blue Wildebeest (Connochaetes taurinus), Waterbuck (Kobus elipsiprimnus), and White Rhino (Ceratotherium simum) were derived from aerial surveys conducted between 1987 and 1993 (Smit et al. 2007).
While zebra (Equus burchelli) is one of the most abundant herbivore in the park we did not include these in our analysis as they occurred in very low densities in our plots and we hypothesis that this is because they are distributed evenly across the landscape as they are neither short nor long grass specialists (Archibald et al. 2005). Densities were averaged over a 3 km buffer around plots to account for local anomalies and to smooth out spatial inaccuracies (Izak Smit pers. comm.). We used these data in the present study because these are the only data available based on a full aerial census of the KNP giving total counts of all animals with the precision needed for our study. From 1998 a line-transect method, using fixed wing aircraft, replaced the total count method (Kruger et al. 2008). This new census method gives population estimates for monitoring trends in animal numbers over time in each of four blocks, namely: the southern KNP (from the Crocodile to the Sabie River), the central KNP (from the Sabie to the Olifants River), the northern KNP (from the Olifants to the Shingwedzi River), and the far northern KNP (from the Shingwedzi to the Limpopo River). We examined the relative abundance of herbivores within each of the four blocks and found that these have not changed substantially between 1998 and 2016 (see Supplementary Fig. S1). Analysis of variance, treating years as independent observations, identified the proportional abundance of grazers as differing significantly amongst blocks (F 3,32 = 60.67, P \ 0.0001), with a Tukey post hoc identifying the southern and central KNP as being significantly different from the northern and far northern KNP in terms of grazer density. We use this result to justify our use of the precise aerial surveys done between 1987 and 1993.

Soil nutrient content
Our plots were intentionally located on the three major soil types in the KNP basalt, granite, and shale soils. Differences in the nutrient content of these soils have been previously determined by summing the concentrations of base cations (Venter 1990). These data show that the cation exchange capacity of granite soils along a catenal sequence ranges between 19 and 30 me/kg at the top of the sequence, and between 36 and 155 me/kg at the bottom of the sequence. For the basaltic soils these values are higher, ranging from 241 and 301 me/kg at the top to 351 and 435 me/kg at the bottom of the sequence. Values for the shale soils are slightly higher than those of the granites ranging between 52 and 57 me/kg at the top of the sequence, to 82 and 188 me/kg at the bottom of the sequence (Venter 1990). What is important for our study is that these results show that the granite soils are generally less fertile than those derived from basalt and shale. It has also been demonstrated that nitrogen concentrations are highest in the basalt soils and lowest in the granite soils. These values range from 0.22% under trees to 0.17% away from trees on the basalt soils, and from 0.17% under trees to 0.08% away from trees on the granites (February and Higgins 2010).

Phylogeny
In order to assess phylogenetic structuring of communities, phylogenetic signal in traits and phylogenetically corrected trait-trait associations, we used the phylogenetic hypothesis of Bouchenak-Khelladi et al. (2014), pruning it to retain only the species sampled in the present study. Since Bothriochloa radicans and Digitaria eriantha, two species sampled in the KNP, were not included in the phylogenetic tree, this reduced the data set to 66 species (370 individuals). Tree manipulation was done using functions in the ape package (Paradis et al. 2004), as implemented in R version 3.3.2 (R Core Team 2016).

Phylogenetic signal and correlation between foliar traits
Both SLA and leaf tensile strength were log transformed to reduce skewness. Phylogenetic signal in both traits was assessed at the species level, using the mean value of traits measured for different individuals of each species. The package phytools (Revell 2012) was used to quantify phylogenetic signal using Pagel's (1999) k and by comparing the mean standardized contrast for each trait to a null distribution obtained by randomizing the tip values 999 times (Blomberg et al. 2003).
We then measured the correlation between SLA and leaf tensile strength at the individual level, using phylogenetic generalized least squares regression (PGLS) as implemented in the R package caper (Orme 2012), applying the k model to structure the residuals. Given that our trait data set contained multiple measurements per species (from multiple sites), whereas our phylogenetic tree represented each species by a single terminal, it was necessary to graft additional terminals onto the tree to accommodate replicate samples of each species. These terminals were added as polytomies whose crown nodes were positioned along the terminal branch length leading to the species in question, at an arbitrary depth of 1% of its length. The decision to represent intraspecific divergences as extremely recent was based on the assumption that conspecifics at Kruger share a common gene pool. Positioning the intraspecific polytomies deeper, at a depth of 10% of each species' terminal branch length, did not, however, qualitatively alter our results.

Phylogenetic community structure and composition
Phylogenetic structuring of the grass communities within each plot was quantified using the standardized effect size of the mean phylogenetic distance between species [SES MPD ; (Kembel et al. 2010)] relative to a null distribution based on 999 randomizations in which community richness was held constant. Values of this metric, which is -1 times the net relatedness index [NRI; (Webb et al. 2002)], and their associated significance levels were determined using the ses.mpd function in the picante package (Kembel et al. 2010). To test whether there was a general tendency towards phylogenetic clustering (SES MPD \ 0) or evenness (SES MPD [ 0) across all plots, we used a Wilcoxon signed-ranked test to assess whether the mean SES MPD differed significantly from zero. In addition, to assess whether and how phylogenetic structuring was related to environmental variables, we tested the correlation of SES MPD to mean annual precipitation and whether it differed (one-way ANOVA) between plots situated on different geologies.
Focusing on five major subfamily/tribe-level PAC-MAD lineages (Andropogoneae, Aristidoideae, Cynodonteae, Eragrostideae, Paniceae; see Table S1 for details), we then determined whether the proportion of species representing each lineage in communities significantly differed among substrate types using Pearson's v 2 test. To assess whether individual lineages favour specific environmental conditions, we tested whether variation in the relative representation of lineages within plots was associated with differences in MAP (linear regression) or geology (t test).

Relationships of leaf traits to environmental variables
We investigated how environmental variables influence traits by fitting linear regressions between average trait values in each grass community and MAP plus soil type as additive effects. We went a step further and asked whether trait turnover across communities is due to turnover in species composition or to leaf trait variation within species. To do so, we separated trait variation across grass communities into three components: variation that is due to species turnover, variation that is intraspecific, and variation that is due to the covariation between species turnover and intraspecific variation. This trait decomposition was done following the method of Lepš et al. (2011). It compares changes in community averages of traits calculated either using the mean trait of each species in all communities, which reflect species turnover only, or using changes in community averages of traits calculated from the traits actually observed in each community, which reflect both species turnover and within-species variability in traits.

Relationship of grazer density to leaf traits and grass composition
We tested how community-averaged values of SLA and leaf tensile strength correlated with the density of grazers. This was done by comparing the fit of alternative linear regression models in which grazer density was predicted by SLA, leaf tensile strength and their interaction. In a separate set of analyses, we finally tested whether grazer density was correlated with the relative proportion of species from different grass families using linear models.
To visualize relationships between environmental variables, grazer density and community-averaged values of leaf traits, we ran a principal component analysis (PCA) on values of all continuous variables for all grass communities sampled. Soil type, as a discrete variable, was not included in this ordination.

Trait phylogenetic signal and trait-trait associations
Both log-transformed traits exhibited significant phylogenetic signal, whether assessed using Pagel's k or the variance of the standardized contrasts, regardless of whether polytomies representing multiple individuals within species were placed at the last 1% or 10% of terminal branches (P \ 0.05 for both indices and both polytomy depths, see trait distribution across the phylogeny on Fig. 2). In general, subfamily Aristidoideae displayed the lowest SLA and highest leaf tensile strength, while the opposite was true of tribe Paniceae (Fig. 2). Across all individuals sampled, the two traits were negatively correlated (Fig. 3). This was true both for individual polytomies placed at 1% (log-log slope = -0.53, F 1,357 = 75.76, r 2 = 0.17, P \ 0.001) and at 10% of terminal branches (log-log slope = -0.54, F 1,357 = 76.46, r 2 = 0.17, P \ 0.001). In both cases, the residuals of their relationship showed significant phylogenetic signal (k = 0.60, P \ 0.001).
The five PACMAD clades sampled in this study were unevenly represented in the data set, with 35% of individuals belonging to Paniceae, 24% to Eragrostideae, 20% to Andropogoneae, 10% to Cynodonteae, and the remaining 10% to Aristidoideae (see Table S1). Clades were unequally distributed across soil types (Pearson's v 2 test: v 2 = 18.8, df = 8, P = 0.016; Fig. 1). When looking at the influence of environmental variables on the representation of each clade separately, we found that the proportion of Paniceae species was significantly greater in communities on basalt and sedimentary soils than on granite soils (P = 0.011 and 0.038, respectively), although  and leaf tensile strength (blue dots) are shown for each species next to the corresponding tip of the phylogeny this did not correlate with MAP (P = 0.150; overall model statistics: F 3,35 = 2.69, P = 0.061, r 2 = 12%). By contrast, the proportion of Eragrostideae species was greater in granitic communities than in either basaltic (P = 0.019) or sedimentary communities (P = 0.008) and was negatively associated with MAP (P = 0.028; overall model statistics: F 3,35 = 2.96, P = 0.046, r 2 = 13%). Aristidoideae showed a similar trend, being more prevalent on granite than on basalt (P = 0.057) and sedimentary soils (P = 0.013), its prevalence also being negatively associated with MAP (P = 0.064; overall model statistics: F 3,35 = 2.38, P = 0.087, r 2 = 10%). However, neither soil type nor MAP influenced the community-level representation of Andropogoneae (overall model statistics: F 3,35 = 1.01, P = 0.400, r 2 = 0%) and in Cynodonteae, there was only a marginally significant higher incidence on sedimentary soils (P = 0.071; overall model statistics: F 3,35 = 2.20, P = 0.11, r 2 = 9%) ( Table 1).

Trait-environment associations
Variation in the community-average values of both traits was associated with environmental variables (Fig. 3). SLA was negatively associated with MAP (P = 0.020) but did not depend significantly on soil type (P [ 0.100 for all three pairwise comparison between soil types). Overall, environmental variables explained a large proportion of the variation in mean  SD standard deviation SLA across communities (r 2 = 039, overall model statistics: F 3,36 = 9.15, P \ 0.001). Leaf tensile strength did not correlate significantly with MAP (P = 0.840) but was significantly higher on granite soils than on basalt (P = 0.003) and sedimentary soils (P = 0.019). Here again, environmental variables explained a large fraction of the variation in the mean leaf tensile strength at the community level (r 2 = 0.33, overall model statistics: F 3,36 = 7.47, P [ 0.001).
When decomposing the variation in leaf traits across hierarchical levels, we found that a large portion of variation in both traits was due to intraspecific variability. For SLA, 17.8% of the variation was attributed to species turnover among communities, 50.4% was due to intraspecific variability, and the remaining 31.8% was due to the covariation between species turnover and intraspecific variability. For leaf tensile strength 36.8% of variation across communities was due to species turnover, 30.3% to intraspecific variability, and the remaining 32.9% to their covariation.

Relationship of grazer density to leaf traits and grass composition
When comparing the statistical fit of different models, we found support to be strongest for a model in which grazer density was explained by both communityaveraged SLA and tensile strength, as well as the interaction between these two variables (DAIC c [ 3.9 with all other models; overall model statistics: F 3,33 = 7.27, P = 0.00071). In this model, grazer densities generally increased with higher SLA and lower leaf tensile strength, being highest for a combination of these extreme trait values (Fig. 4). Jointly, these variables explained a large fraction of the variation in grazer densities across communities (r 2 = 0.34). Adding MAP and soil type as covariates in this model led to a lower AIC (dAIC = 4.27) with neither of these environmental variables correlating significantly with grazer densities (P [ 0.52 in both cases).

Discussion
Our results show that leaf functional traits of 66 grass species from 37 communities distributed across the KNP are associated with both soil type and MAP. In turn, the community means of these traits are highly correlated with grazer densities in this savanna ecosystem.

Effect of rainfall and soil type on leaf functional traits
Although the negative relationship between SLA and leaf tensile strength observed in this study is not ubiquitous (Wigley et al. 2016), such a relationship is expected given the dependence of both these traits on leaf fibre density. The association of both traits with MAP and soil type is equally expected given that these traits strongly influence leaf lifespan and resilience to damage. Indeed, MAP and soil type together had a large influence on leaf traits of the KNP grasses, explaining more than 33% of their variation in both  Fig. 4 Relationship between grazer density and grass functional traits. The grey surface shows the fitted response of the density of grazers at a given location as a function of the community-average SLA and leaf tensile strength of grass species as well as their interaction cases. Consistent with earlier studies (Wright et al. 2001), grasses growing on less fertile granite soils had higher leaf tensile strengths than those occurring on the higher nutrient basalts, a pattern that is presumably attributable to the nutrient conservation benefits associated with high leaf tensile strength (Wright and Westoby 2003;Pérez-Ramos et al. 2012). Despite the relatively small MAP range in our study (400 to 750 mm year -1 ), we found a negative correlation between mean SLA at the community level and rainfall. These results are similar to those reported for Stipa-dominated grasslands in Northern China (Liu et al. 2017). In these communities, this pattern was interpreted as reflecting an adaptive response to drought (Liu et al. 2017) and we favour a similar interpretation here. These results may, however, confound the influence of soil fertility at our study site with several studies revealing strong relationships between soil fertility and SLA (Craine et al. 2001;Al Hadj Khaled et al. 2005;Ordoñez et al. 2009). While the overall relationship between substrate and average SLA at the community level was not significant, there is a weak tendency for SLA to be lower on the granite than on the basalt (P = 0.056) and sedimentary soils (P = 0.270).
The relationship between MAP and SLA may, therefore, be underpinned by the fertility of the underlying soils suggesting that soil type and MAP jointly determine grass leaf traits in the KNP.

Phylogenetic distribution of grass species
Overall, our study provides limited evidence of phylogenetic structure in the grass communities studied with only four out of 37 communities exhibiting significant phylogenetic clustering and only one significant over dispersion. We did, however, find that some of the major grass clades were unevenly distributed between soil types and along the rainfall gradient. For example, while Paniceae were more prevalent on basaltic and sedimentary soils, Aristidoideae and Eragrostideae were more prevalent on granite, the latter two also declining with increasing rainfall. Interestingly, where previous authors have found strong associations of Andropogoneae with moist, infertile environments (Visser et al. 2012;Forrestel et al. 2015), we found no evidence of such an association. This is most likely a consequence of our study covering much narrower environmental gradients than these previous studies, which covered much broader geographic areas. There is also a large range of phenotypic diversity within each of these clades, and, as such, our results should be taken with caution.
As environmental gradients influence species distributions and our results show a strong phylogenetic signal in key functional traits, the distributions of the major lineages present in the KNP are likely a consequence of traits associated with environmental filtering (Scholtz et al. 2014). Therefore, the low SLA and high tensile strength of aristidoid grasses may account for their greater prevalence in low-fertility, granitic communities. In contrast, the high SLA and low tensile strength of Paniceae likely explain the association of this tribe with more fertile, basalt settings.
We did, however, find that a large fraction of leaf trait variation across grass communities of the KNP is due to intraspecific variability while species turnover among communities accounts for less variation in leaf traits. This cautions against excessive reliance on species inventories as a proxy for the functional characteristics of grass communities, and suggests that we should improve our understanding of the links between the dynamics of individual plants and ecosystem functioning (Ackerly and Cornwell 2007).

Relationship between grass traits and grazer density
At the community level, we found that grazer densities in the KNP were highest on sites characterized by low leaf tensile strength (O'Reagain 1993), and high SLA, regardless of rainfall and soil type. These relationships were strong, with these foliar traits jointly explaining 34% of the variation in grazer densities across the KNP. Given the strong negative association of grazer density with leaf traits, the observation of low grazer densities in sites with high frequencies of Aristidoideae and Eragrostideae is unsurprising. Aristidoideae are not palatable for large herbivores because of their high leaf tensile strength, and high epidermal silica densities (Bouchenak-Khelladi et al. 2009). This observation is consistent with the hypothesis that some C 4 grass lineages, such as Aristidoideae, have evolved less palatable leaves as a defence against herbivory (McNaughton and Georgiadis 1986;Bouchenak-Khelladi et al. 2009). These findings are in line with the distinction between grazing lawns (highly palatable) and tall-grass (less palatable) communities (Hempson et al. 2015), and suggest that these two types of communities differ not only in taxonomic composition but also in functional traits.
While we have demonstrated a strong correlation between grass foliar traits and grazer densities, the causal directionality underpinning this relationship remains unresolved. On the one hand, leaf traits determine both the quality and quantity of available forage (Wright et al. 2001;Pontes et al. 2007;Duru et al. 2009) and should thus influence grazer densities. On the other hand, grazers may exert a positive feedback on leaf traits. For example, high grazer densities are likely to have a positive effect on community-level grass foliar traits through the effects of their faeces in enhancing soil fertility and of their defoliation in selecting for rapid regrowth (Wright et al. 2001;Hempson et al. 2015). Although several studies have demonstrated the associations between grazing ungulates and the moisture or nutrient content of grasses (MacNaughton 1983;Holdo et al 2009;Augustine and MacNaughton 2006), the effects of grazers on specific grass traits such as SLA and leaf tensile strength remain poorly understood (van der Plas et al. 2016).
More generally, there is also an extensive literature on the extent to which top-down (herbivory and fire) and bottom-up (substrate and rainfall) factors structure plant communities in savanna ecosystems (e.g. February et al. 2013;Archibald and Hempson 2016;Lehmann et al. 2014;Forrestel et al. 2015;February and Higgins 2016). An understanding of the interactions between top-down and bottom-up controls, and how these may vary in space and time does, however, remain elusive. Understanding how these interactions influence functional traits at the community level is, however, essential for an improved understanding of African savanna vegetation structure. Unfortunately, the resolution of our herbivore data does not allow us to assess the relative importance of top-down vs. bottom-up controls on grass communities of the KNP. Such a task would at least require fine-scale measurements of herbivore densities (in contrast to the 3 9 3 km resolution of our data), given that herbivore communities show small-scale spatial heterogeneity (Burkepile et al. 2016) and respond to very small-scale variation in plant communities (Cromsigt and Olff 2008). But this important question would be even better dealt with using time-series measurement of both grazer densities and grass functional traits and/or manipulative experiments.

Conclusion
We have demonstrated complex interactions between climate, soils, and grazer densities in relation to grass functional traits in the KNP. We found that the distribution and density of grazers is associated with the foliar traits of the associated grass communities, possibly via the influence of MAP and soil fertility on both SLA and leaf tensile strength, and consequently on leaf palatability of grass communities. This has significant implications for our understanding of the evolutionary processes that have led to the expansion of C 4 grasslands and savannas worldwide. However, the resolution of our herbivore data does not allow for a clear determination of the extent to which top-down controls influence this association, and at what levels.