Among‐species overlap in rodent body size distributions predicts species richness along a temperature gradient

Temperature is widely regarded as a major driver of species richness, but the mechanisms are debated. Niche theory suggests temperature may affect richness by filtering traits and species in colder habitats while promoting specialization in warmer ones. However, tests of this theory are rare because niche dimensions are challenging to quantify along broad thermal gradients. Here, we use individuallevel trait data from a long-term monitoring network spanning a large geographic extent to test niche-based theory of community assembly in small mammals. We examined variation in body size among 23 communities of North American rodents sampled across the National Ecological Observatory Network (NEON), ranging from northern hardwood forests to subtropical deserts. We quantified body size similarity among species using a metric of overlap that accounts for individual variation, and fit a structural equation model to disentangle the relationships between temperature, productivity, body size overlap, and species richness. We document a latitudinal gradient of declining similarity in body size among species towards the tropics and overall increase in the dimensions of community-wide trait space in warmer habitats. Neither environmental temperature nor net primary productivity directly affect rodent species richness. Instead, temperature determines the community-wide niche space that species can occupy, which in turn alters richness. We suggest a latitudinal gradient of trait space expansion towards the tropics may be widespread and underlie gradients in species diversity.


Introduction
An enduring challenge in ecology has been to understand the major drivers of species richness. Perhaps the most general and noted pattern is the latitudinal diversity gradient, in which species diversity increases towards the tropics (Fischer 1960, Hillebrand 2004. Temperature covaries with latitude and is thought to underlie this pattern (Hawkins et al. 2003, Tittensor et al. 2010, but the mechanisms are unclear. One classic hypothesis is that low winter temperatures filter out species that cannot subsist through periods of cold and deprivation, reducing the available niche space that can be occupied (MacArthur 1972). Conversely, in warmer environments, the reduction of environmental filters and greater stability of resources should lead to higher niche diversity and increased richness. In addition, as abiotic filters decline in importance, the relative role of biotic filters, such as competition, should increase. With stronger biotic filters, species in warmer environments are predicted to exhibit more niche partitioning and specialization, resulting in increased local richness (Schemske et al. 2009).
Evaluating how temperature drives community assembly processes has been challenging at large scales. The development of a trait-based approach to community ecology (McGill et al. 2006, Levine 2016) that incorporates intraspecific variation (Violle et al. 2012), offers a tractable approach for assessing biogeographic theories of diversity. Traits are mechanistically linked to niche dimensions (Kearney et al. 2010, Lamanna et al. 2014; the relationship between trait distributions of co-occurring species illuminates how niche-based processes structure local communities (Fig. 1). Reduced environmental filters in warmer, more stable environments may permit the survival of individuals with more extreme trait values, allowing more species to coexist (MacArthur 1972). In addition, an increase in biotic filters in warmer habitats -particularly competitive pressure -may narrow or shift the distributions of trait values within species, limiting similarity among species and reducing competitive pressure (MacArthur andLevins 1967, Schemske et al. 2009). Overall, as trait space expands and trait widths decline in warmer habitats, trait overlap is expected to diminish and opportunities for coexistence and diversification should increase.
The expansion of long-term ecological networks and community-wide trait sampling has opened new avenues for exploring spatial gradients of traits. Assessment of functional trait variation within species at community and ecosystem scales has been largely restricted to plants, reflecting the relative ease of sampling immobile taxa and historical efforts to standardize trait measurement (Díaz et al. 2004). However, the recently established National Ecological Observatory Network (NEON) uses standardized sampling protocols to measure small mammal body sizes in a range of habitats across the United States. Body size is an important trait that directly or indirectly affects many aspects of an individual's niche, including dietary preferences, habitat requirements, physiological rates, biological times, access to resources, and competitive interference (Peters 1983, Calder 1984, Gravel et al. 2013. While ecologists have found evidence that body size distributions within single communities are more even than expected, suggesting niche partitioning (Bowers and Brown 1982), thermal or latitudinal gradients in community size distributions have received comparatively little attention.
To test several hypotheses about how temperature affects niche dimensions and diversity, we examined the body mass distributions of 23 nocturnal rodent communities across the United States, ranging from subtropical scrub to northern hardwood forest. We assess spatial patterns of trait similarity and breadth and evaluate the role of environmental temperature as an abiotic filter, while also accounting for the role of primary productivity (Brown 2014). We quantify body size similarity in communities by measuring the degree of overlap between all species pairs in a community, and calculate the median pairwise overlap value to characterize each rodent community.
NEON data on community-wide intraspecific trait variation allows us to treat body size distributions as proxies for species niche dimensions (Kearney et al. 2010). To determine whether temperature directly predicts richness or acts indirectly through trait dimensions, we employ structural equation models (SEM) and null models to test the following predictions about North American small mammal community assembly (Fig. 1): 1) small mammal richness increases with temperature. 2) Temperature affects richness by modifying community-level niche dimensions (i.e. body size overlap).
3) The effect of body size overlap on richness is driven by: a) increases in the range of body sizes within Figure 1. Niche theory predictions for traits along a thermal gradient. Here we show communities in which each species has a unique trait distribution, each represented by a different color. As temperature increases, environmental filters are expected to become weaker and biotic interactions such as competition to become increasingly important. Weaker environmental filters and stronger biotic interactions may promote increased richness in three ways: 1) an increase in the range of potential trait values in the community (top panel); 2) a decrease in the width of individual species' trait distributions (middle panel); or 3) species shifting their trait distributions relative to one another to minimize overlap (bottom panel). communities; b) decreases in the average width of specieslevel body size distributions; c) species shifting their body size distributions to become more even relative to one another.

Site and sampling description
We used small mammal data collected by the National Ecological Observatory Network (NEON) in 2015. NEON is a network of observational sites comprising 47 terrestrial locations distributed across 20 climate-defined regions in the contiguous United States, Alaska, Hawaii, and Puerto Rico. At each site, climate, ecosystem, and organismal data are collected according to a set of standardized protocols. The study sites vary widely in climate, local topography, and vegetation type (Supplementary material Appendix 1). We used climate data from PRISM (PRISM Climate Group 2004) for the years 1981-2015 to extract site-level average climate data, elevation data from the National Elevation Dataset (U. S. Geological Survey 2016) to calculate topographical heterogeneity, and remotely sensed productivity data for the years 2011-2015 from NASA's Moderate Resolution Imaging Spectroradiometer (U. S. Geological Survey 2014) to characterize the environmental conditions at each site.
NEON staff trapped small mammals at 26 of the 47 sites in 2015 (Fig. 2). At each site, NEON personnel established a variable number (3-6) of trapping grids with 100 live traps (H. B. Sherman Traps, Tallahassee, FL, USA) arrayed in a grid within a 100  100 m area. They trapped, weighed, and released small mammals at each grid approximately once per month during the growing season (4-7 months per site, typically May-September but including March-October at a few sites). As traps were set at night on the ground, the sampling design targeted nocturnal, non-volant mammals. We excluded three sites from analysis (Disney Wilderness Preserve, FL; Dead Lake, AL; and Healy, AK) due to low trapping success: the former two sites were located within wetlands, and the latter site was only sampled for a single bout. Ethical approval was obtained from IACUC. Refer to the NEON small mammal trapping protocol (NEON 2016) for more details.

Processing the NEON data
We processed and analyzed all data in R 3.3.3 (R Core Team). We excluded all non-rodents from the dataset, as well as all individuals that could not be identified to species. NEON protocols designate diurnal and volant species as unintentional bycatch; we excluded individuals of those species from our analysis. Together, non-rodents, diurnal rodents, and volant rodents were approximately 10% of all captures. We also excluded grasshopper mice due to their carnivorous diet (Onychomys spp.; approximately 3% of all captures). Captured individuals were marked, so we excluded recaptures from the final analyzed dataset (Supplementary material Appendix 2). Most rodents tend to be generalist consumers, eating plant items with limited fiber content (e.g. seeds, buds) supplemented with protein-rich invertebrates (Landry 1970, Pineda-Munoz andAlroy 2014). Consequently, we evaluate similarity in body size across all co-occurring rodents rather than within feeding guilds.
We used rodent body mass as a proxy for resource niche. We selected the earliest-recorded body mass measurement if multiple measurements were taken on the same individual within the year. Mass measurements varied relatively little over time within individuals (CV of body mass 0.25 in 96% of recaptured individuals). All mass measurements were log 10 -transformed to reduce heteroscedasticity and nonnormality. Finally, we pooled all trapping grids from a single site into a single community, since most trapping grids were separated by only 1-2 km, with the median distance between grids ~1200 m across all sites.

Comparison of richness estimators and removal of outliers
We calculated two species richness estimators that account for inadequate sampling: the Chao1 richness estimator (Chao et al. 2009) and the asymptotic richness estimator based on Hill numbers (Chao et al. 2014). These methods use the abundance distributions of individuals in a sample of a species assemblage to estimate the true number of species in the underlying community. We found that the estimators were qualitatively similar to one another, and that most sites were well sampled (Supplementary material Appendix 3). Therefore, in all analyses below, we use the Chao1 estimator to represent rodent richness. Our analysis identified three sites (CPER, NIWO, and MOAB) as being relatively less well sampled, but the Chao1 estimator should closely approximate the true underlying richness at those sites, and excluding those sites had no qualitative effect on our results.

Calculating overlap statistics to measure trait partitioning
While niche partitioning can be measured by assessing mean trait differences among species, this may be misleading if intraspecific trait variation is large (Violle et al. 2012, Siefert et al. 2015. A more robust approach is to measure individual traits directly and then assess the degree of trait overlap between potential competitors (Mouillot et al. 2005, Mason et al. 2011. Low trait overlap indicates high trait partitioning. We generated a statistic to measure the degree of community-level trait overlap that characterizes the level of niche partitioning in trait space by fitting nonparametric kernel density functions to each species' trait distribution and calculating their areas of overlap (Mouillot et al. 2005, Geange et al. 2011). The median pairwise overlap for a community is calculated by first determining the overlap of each species pair in trait space, and then taking the median overlap of each species pair in a community.
Median pairwise overlap in a community can range from zero, with no overlap among any species pairs, to one, with complete overlap of all species pairs. Overlap values approaching zero indicate a high degree of niche partitioning. To calculate the overlap statistic, we first generated a kernel density estimate for the trait distribution for each population. We used the 'density()' function with the 'bw.nrd0' method in the base R distribution to select the appropriate bandwidth for the kernel. The overlap statistic was not sensitive to the method used to select the optimal bandwidth (for our dataset, the minimum pairwise correlation among all bandwidth selection methods was 0.96). The density functions for each species are normalized so that the integral of the function is one. We calculated the pairwise overlap for each species pair in the community by taking the integral of the minimum of the two functions and weighted each pairwise overlap by the harmonic mean of the abundance of both species in the pair. The community-level niche overlap statistic (hereafter, simply 'overlap') is the abundance-weighted median pairwise overlap of all species pairs in the community. We calculated the pairwise overlap between species in the same way as Mouillot et al. (2005), but our community-level metric differs in that we used the median instead of the mean to reduce the influence of outliers. In addition, we weighted each value by the harmonic mean abundance of the species pair to reduce the influence of rare species on the value of the metric.

Selection of environmental predictor variables for structural equation model
To assess the prediction that temperature affects rodent richness by modifying community niche dimensions, first we examined bivariate correlations among potential predictor variables to determine which to include in the SEM. The list of potential predictors included mean annual temperature and precipitation, minimum temperature of coldest month, mean precipitation of driest month, seasonality of temperature and precipitation, interannual variation of temperature and precipitation, and MODIS-derived variables including net and gross primary productivity, vegetation indices, and leaf area index. We included two principal component axes describing heterogeneity in topography and in remotely sensed productivity variables (Supplementary material Appendix 5) to test the prediction that heterogeneous habitats would promote individual adaptation to more diverse microsites (Bolnick et al. 2003), allowing species in different microhabitats to share a body size niche.
Variables related to precipitation and climate seasonality were highly correlated with other potential predictors, and variables related to topographic and environmental heterogeneity explained very little variation in body size overlap and species richness. As a result, we did not include these variables in our final model. Variables describing temperature were also correlated, so we chose to focus on the minimum temperature of the coldest month. This variable not only provides a measure of thermally stressful period but, perhaps more importantly, the degree of resource deprivation facing animals, as cold winters are associated with a scarcity of insects, fruit and leaves (MacArthur 1972). Finally, the variables describing environmental heterogeneity and interannual climate variation had low bivariate correlations with overlap and richness, so we did not include them in the SEM.

Investigating causal relationships among environmental variables, body size overlap, and species richness
We asked whether temperature acted on trait similarity to influence richness, or whether temperature alone was a better predictor of diversity. In order to examine these possible causal pathways, we fit an SEM to estimate the strengths of direct and indirect overlap-mediated effects of temperature on richness. We hypothesized that there is a direct causal pathway between temperature (T; minimum temperature of coldest month) and richness (R; Chao1 species richness estimator). The hypothesized indirect pathway between temperature and richness consists of a causal effect of temperature on interspecific body size overlap (O), and a further causal effect of overlap on richness (Fig. 4a). The structure of this model allows us to distinguish direct and indirect overlap-mediated effects on richness. We also fit SEMs including direct and indirect effects of net primary productivity and compared them to models including only temperature effects using the Bayesian information criterion (BIC; Supplementary material Appendix 5).
We fit the model in the R package 'blavaan' (Merkle and Rosseel 2016), initializing 3 MCMC chains that each took 25 000 samples from the posterior distribution. We logittransformed the overlap statistic to meet assumptions of normality and homoscedasticity. At each iteration, we estimated each of the β parameters and the indirect effect size and net effect size derived by combining path coefficients, discarding the first 5000 samples from each chain as burnin. We assessed convergence on the target posterior distribution visually and by confirming that the Potential Scale Reduction Factor for each parameter (Gelman and Rubin 1992) was 1.001. We estimated the 95% credible interval around each parameter estimate, as well as an empirical p-value for each parameter, from the posterior samples. We compared the sizes of the direct, indirect, and net effect of temperature on richness by comparing standardized coefficients.

Determining causes of variation in body size overlap
To explore the drivers of variation in body size overlap (predictions 3a-c), we calculated three metrics: community-wide quantile range tests whether expansion of community-wide trait space underlies decreased overlap, average species quantile range tests whether contraction of species niche widths underlies decreased overlap, and effect size from a null model (z-score) tests whether trait sorting within the bounds of trait space underlies decreased overlap. First, we estimated the size of trait space in each community by pooling individual measurements across all species at a site and calculating the community-wide central 90% quantile range (Q .95 -Q .05 ) of body mass. We also calculated the 90% quantile range for each species at each site and took the abundance-weighted mean across all species at each site. Finally, we ran a null model to test whether individual species' body size distributions were more evenly spaced along the body size axis than expected by chance, finding the z-score of each observed community evaluated against the distribution of 999 null communities in which species body size means were shuffled randomly, retaining the shape and width of the distribution around the mean for each species. We regressed the logit-transformed overlap statistic against each of these predictors in a multiple regression to determine which of them underlie overlap, using BIC to determine which predictors to retain in the most parsimonious model and standardizing coefficients to determine the relative importance of predictors.

Data deposition
All data and R scripts necessary to reproduce the analyses in this manuscript are archived at  https://doi.org/10.6084/ m9.figshare.5339407  (Read et al. 2017).

Results
Across 23 rodent communities, the degree of overlap among species within a site varied 50-fold, from 0.018 at the Jornada Experimental Range in New Mexico to 0.90 at Harvard Forest in Massachusetts (example sites shown in Fig. 3; all sites shown in Supplementary material Appendix 6 Fig. A7). A latitudinal gradient of size overlap is apparent across the United States, with overlap declining toward warmer and lower-latitude sites ( Fig. 2; r 2 = 0.28, p = 0.004).
Temperature was a poor direct predictor of rodent richness in the United States (Fig. 4b), which is highest in the western portion of the country. However temperature was a strong predictor of trait overlap which, in turn, was a strong predictor of richness. The SEM indicated that winter temperature had a positive indirect effect on rodent species richness via decreased body size overlap among species in warmer environments (Fig. 4). Sites with lower minimum temperature of the coldest month harbored rodent communities with significantly higher body size overlap (standardized path coefficient -0.81; 95% credible interval [-1.48, -0.14]). Communities with high body size overlap, in turn, had significantly lower species richness (coefficient -0.44; 95% CI [-0.66, -0.23]). The direct effect of temperature on species richness that was not mediated by overlap was not distinguishable from zero (coefficient -0.19; 95% CI [-0.56, 0.18]). The two indirect pathways with negative coefficients resulted in a positive overlap-mediated indirect effect of temperature on richness (coefficient 0.36; 95% CI [0.01, 0.71]). The net effect of temperature on richness was not distinguishable from zero (coefficient 0.17; 95% CI [-0.28, 0.61]). Net primary productivity was negatively associated with rodent species richness and positively associated with overlap. However, models including net primary productivity were far less parsimonious than the temperature-only model (ΔBIC -50 in all cases) and added little explanatory power (Supplementary material Appendix 5).
Decreased trait similarity among species was associated with an expansion of community-wide trait space and a decrease in species average trait width, consistent with predictions 3a and 3b (Fig. 5). We did not, however, find evidence that species trait distributions shifted relative to one another to limit overlap; our results did not support prediction 3c. As the size of community-wide trait space, represented by central 90% quantile range, decreased, body size overlap increased (standardized coefficient -1.37; 95% confidence interval [-2.17, -0.57] ; Fig. 5a). In addition, as the average niche width (central 90% quantile range) of species increased, body size overlap increased (standardized coefficient 0.76; 95% confidence interval [0.1, 1.37]; Fig. 5b). All three predictors (community quantile range, species quantile range, and null model z-score) were retained in the best model, together explaining 82% of the variation in body size overlap, but the null model z-score did not explain a significant amount of variation. These results indicate that both expansion of trait space and contraction of niche widths, but not increased evenness of distributions along the trait axis, underlie decreased overlap.

Discussion
One of the most robust predictors of species diversity at global scales is environmental temperature, with species richness generally increasing in warmer areas. Environmental and biotic filters on species traits may underlie this pattern, but this is difficult to test due to limited availability of standardized trait data across large spatial scales. Here, we used trait data from the National Ecological Observatory Network to observe a thermal and latitudinal gradient in trait similarity across the continental United States: rodent populations were increasingly dissimilar in body size in warmer locations. Community trait space had a broader range in warmer environments, suggesting a relaxation of abiotic filters on body size. In addition, body size variance of individual rodent populations was lower in warmer habitats, indicating that niches are narrower in areas with milder winters. Overall, our results support niche theory predictions of increased specialization and relaxed filtering in warmer environments (MacArthur 1972) that have been largely untested in animal communities. We did not, however, find evidence of increased evenness of trait distributions in our communities that would indicate higher competitive pressures; more detailed study may be necessary to thoroughly test for this pattern.
Intriguingly, temperature alone was a poor predictor of small mammal diversity in the United States. For endothermic mammals with a constant body temperature, environmental temperature does not directly affect the metabolic and life history rates, such as mutation rates or generation time, that drive evolutionary rates (Gillooly et al. 2005). Instead, our results indicate that environmental temperature affects mammal richness via its influence on niche dimensions. Specifically, the number of unique niches in a community increases in milder, less biotically challenging environments as total niche space widens and average species niche width shrinks. Niche space in plants shows a similar pattern, with total niche space at local scales increasing towards the equator (Lamanna et al. 2014). If broader community-wide niche space in warmer areas proves to be a general pattern, it suggests that the latitudinal diversity gradient is driven, at least in part, by a latitudinal gradient in niche dimensions.
Rodent richness in the United States peaked at sites with lower net primary productivity, particularly the arid southwest. While this predictor variable was not included in our final model, this negative relationship demands some explanation, as it runs counter to species-energy theories of diversity and many empirical patterns (Hawkins et al. 2003, Grace et al. 2016, Liang et al. 2016). First we note that plant productivity does not always translate into more available food for vertebrates. With higher rainfall and greater plant production, trees tend to dominate terrestrial vegetation, and much of the production is locked up in hard-to-reach leaves or indigestible wood. From this perspective, terrestrial species would not be expected to increase in abundance or diversity in forest biomes relative to mesic grasslands or savannas. This pattern appears to be borne out in tropical South American rodents, where richness is relatively constant from the grassy Cerrado to the Amazonian rainforest (Maestri and Patterson 2016). Rather, it is likely that the addition of volant and arboreal species, such as bats and primates, leads to higher richness of mammals in the productive Amazonian forests. In addition, higher plant diversity in warmer sites may permit an increased number of small mammal species to coexist by specializing on consuming different plant species; such an effect could exist independent of any trend in plant productivity (MacArthur 1972, Terborgh 2015. However, this causal relationship was not supported by SEM analysis, potentially due to low statistical power (Supplementary material Appendix 5).
Continental geography and historical contingency may also explain the richness peak in the arid southwest. The highest rodent richness occurred in scrub and grassland There is no net effect of temperature on richness (b), overlap declines with increasing temperature (c), and species richness declines with increasing overlap (d). Bivariate linear model fits, back-transformed for plotting, are shown; these do not exactly correspond to the coefficients in the SEM, but highlight the negative relationships between overlap and temperature, and between richness and overlap. communities that extend along the western spine of North America. The spatial proximity and similarity of habitat may permit Mexico to act as a tropical pump (Ceballos et al. 2010), exporting tropical mammal richness to the arid north rather than the wetter east. The variable topography of the southwestern United States also allowed numerous species to persist in isolated refugia through periods of decreased precipitation during the Quaternary (IUCN/SSC Rodent Specialist Group 1998). Furthermore, historic fragmentation of arid and nonarid habitat patches associated with fluctuating precipitation and glaciation may have contributed to higher rates of allopatric speciation in the southwest relative to the rest of North America (Riddle and Honeycutt 1990). These historically and geographically contingent phenomena may underlie the negative productivity-diversity relationship and the weak temperature-diversity relationship we observed.
This study focused on the trait of body size, but as more data become available, analysis of multiple traits will permit more robust tests of niche theories (Sterck et al. 2011) and better exploration of the link between animal traits and ecosystem functioning. Nonetheless, body size is a fundamental organismal trait and is therefore an informative starting point. Researchers have long noted the central role of body size in shaping metabolism, diet, ecological interactions and coexistence (Hutchinson 1959, MacArthur and Levins 1967, Bowers and Brown 1982, Abrams 1983, Brown et al. 2004, Chesson et al. 2004). Tradeoffs associated with body size may reduce competition and promote coexistence. For instance, body size generates constraints on the maximum size of food consumed by many fruit and seed eaters (Muñoz and Bonal 2008), constraining diet selection and limiting similarity to larger guild members. Among herbivores, including rodents, larger body size is associated with the increased consumption of abundant high-fiber foods, promoting dietary breadth (Clauss et al. 2013). For desert seed-eating rodents, body size is positively correlated with seed size preference Lieberman 1973, Mares andWilliams 1977) and microhabitat preferences (M'Closkey 1980). In addition, larger rodent individuals tend to be stronger and more successful in conflicts (Ziv et al. 1993). While larger individuals may gain advantages with respect to food breadth, interference competition, and longevity, this is offset not only by higher total caloric demands but also by slower life histories, such as increased time to maturity and lower reproductive output per time (Peters 1983). Thus, tradeoffs associated with body size promote diversity at a range of sizes. For these reasons, body size is likely a fundamental trait axis driving niche differentiation and community assembly across a variety of animal communities.
With the emergence of trait-based ecology (McGill et al. 2006, Levine 2016, many seminal but hard-to-test theories of community assembly and global diversity appear more tractable. Classic notions of niche partitioning can be more readily assessed using a trait-based framework (Mouillot et al. 2005). Intraspecific variation, which has often been ignored, can be recast as trait variation to inform assembly processes and niche similarity (Violle et al. 2012). By harnessing novel trait-based approaches and large-scale standardized trait data, we can gain new insights into the mechanisms driving canonical macroecological patterns in nature. There is a negative relationship between the size of trait space (community-wide central 90% quantile range of body mass) and overlap (a) and a positive relationship between average species niche width (speciesspecific central 90% quantile range of body mass) and overlap (b). Back-transformed bivariate regression lines, with the predictor variable standardized and conditional on the value of the other predictor, are shown.