Ecological Limits as the Driver of Bird Species Richness Patterns along the East Himalayan Elevational Gradient

Variation in species richness across environmental gradients results from a combination of historical nonequilibrium processes (time, speciation, extinction) and present-day differences in environmental carrying capacities (i.e., ecological limits affected by species interactions and the abundance and diversity of resources). In a study of bird richness along the subtropical east Himalayan elevational gradient, we test the prediction that species richness patterns are consistent with ecological limits using data on morphology, phylogeny, elevational distribution, and arthropod resources. Species richness peaks at midelevations. Occupied morphological volume is roughly constant from low elevations to midelevations, implying that more species are packed into the same space at midelevations compared with low elevations. However, variance in beak length and differences in beak length between close relatives decline with elevation, which is a consequence of the addition of many small insectivores at midelevations. These patterns are predicted from resource distributions: arthropod size diversity declines from low elevations to midelevations, largely because many more small insects are present at midelevations. Weak correlations of species mean morphological traits with elevation also match predictions based on resources and habitats. Elevational transects in the tropical Andes, New Guinea, and Tanzania similarly show declines in mean arthropod size and mean beak length and, in these cases, likely contribute to declining numbers of insectivorous bird species richness along these gradients. The results imply that conditions for ecological limits are met, although historical nonequilibrium processes are likely to also contribute to the pattern of species richness.


Introduction
Spatial patterns of terrestrial species richness correlate strongly with climate-warm, wet areas have more species than cold, dry areas (Hawkins et al. 2003;Field et al. 2009;Pigot et al. 2016a). Explanations for these differences fall into two general categories: nonequilibrium, historical effects of time and diversification rate (e.g., tropics as museum or cradle, time for speciation, age and area), and equilibrium, ahistorical models of carrying capacity (e.g., niche packing and resource diversity) whereby some locations can maintain more species than others, termed ecological limits (Rabosky 2009). The debate over the relative contributions of these two sets of processes is as old as the documentation of the trends themselves (Wallace 1891;Mac-Arthur 1964;Rabosky 2009). The debate persists largely because ancient, climatically stable areas, which are the most important for species accumulation in nonequilibrium models, are also those with high productivity and low seasonality, conditions held to be important for coexistence at equilibrium (Pigot et al. 2016a). It is clear that historical nonequilibrium effects, including Pleistocene disturbances (Hewitt 2000), as well as processes at deeper timescales (Fine 2015), contribute to species richness patterns. Unequivocal demonstrations of contributions of ecological limits to species richness have been harder to make, despite its theoretically predicted importance (Etienne et al. 2019) and demonstrations of associations of resource diversity with diversity of species exploiting the resource, including foliage height (MacArthur 1964), arthropods (Schoener 1971), and fruit (Kissling et al. 2007).
Morphological patterns can be brought to bear on the problem (Lamanna et al. 2014;Oliveira et al. 2016;Edie et al. 2018). For both birds and mammals, a common finding has been that once above a certain threshold, high species richness is achieved by packing species into a constrained morphological space, rather than by expanding occupied morphological volume (Ricklefs 2012;Belmaker and Jetz 2015;Oliveira et al. 2016;Pigot et al. 2016b;Pellissier et al. 2018). This observation has been interpreted as support for nonequilibrium models under which the greater ages of clades or their higher diversification rates have led to the production of different numbers of species packed into a constrained niche space at some locations or in some environments (Ricklefs 2012;Belmaker and Jetz 2015). However, packing of species into functionally similar morphological space may also reflect differences in environmental carrying capacities (Valentine et al. 2008;Price et al. 2014;Oliveira et al. 2016). First, niches may be more finely partitioned in less seasonally variable environments, which favors persistence of many specialists rather than a few generalists (Valentine et al. 2008). Second, locations with high productivity may maintain more species than those with lower productivity (Hutchinson 1959;Hurlbert and Stegen 2014;Storch et al. 2018). These more productive environments maintain more individuals, lowering the risk of population extinction through demographic stochasticity (Schwilk and Ackerly 2005;Storch et al. 2018).
In addition to revealing factors that drive the number of species present in one location, morphology can be used to assess causes of individual species distributions and thereby contribute to an evaluation of ecological limit models. If all species had unlimited ranges, species richness would not vary across space. Hence, bounded ranges are a necessary but not sufficient condition for richness to vary (Lennon et al. 2000). Range boundaries are set through limits on access to sites (movement), followed by abiotic and biotic factors that affect establishment (Boulangeat et al. 2012).
Here we study a gradient of species richness variation along a mountainside for a mobile group (birds), in which range limits surely reflect differential establishment rather than movement. Abiotic and biotic factors and their interaction should all contribute to a failure to establish in a location, but they are often correlated with each other and with elevation. Nevertheless, temperature has been invoked as an important direct cause of avian range limits along environmental gradients (Parmesan and Yohe 2003;McCain 2009;Elsen et al. 2017). If temperature alone is a primary determinant of range limits, elevational assemblages would not be close to an ecological limit, because we would infer that resources are available that could be exploited if the abiotic constraints were relaxed. This in turn would support nonequilibrium models; elevational zones with high richness would have either been occupied for longer or experienced a higher diversification rate with climatic niche conservatism limiting elevational range expansions (Kozak and Wiens 2010). Morphology can address this issue, because two well-known biogeographic rules, notably, Bergmann's rule (larger body sizes at higher latitudes; e.g., see Freeman 2016) and Allen's rule (appendages relatively small at higher latitudes; see Symonds and Tattersall 2010) are explained in terms of adaptation to the direct effects of temperature. These rules predict morphological correlates expected if physiological adaptations to temperature are important, providing an opportunity to look for evidence of direct abiotic filters using morphological data.
In this article, we use measurements of resources and morphology in bird communities along a steep climatic gradient (the east Himalaya) to ask whether richness patterns are consistent with ecological limits or climatic niche conservatism. Previously, Price et al. (2014) showed that among one suborder of birds, the oscines (order Passeriformes), insectivorous species richness correlates with arthropod abundance (both are highest at the 1,200-2,000 m elevation in the east Himalaya). We extend that study to encompass 11 bird orders and use additional measurements of resources to confirm the pattern. Then we address four predictions of how morphological patterns should vary if resources and competition for resources influence the number of coexisting species ( fig. 1). Figure 1 lays out what our expectations are for the different types of data we analyze, given alternative hypotheses (in gray). Support for these predictions would be consistent with an ecological limits model, but they would not reject nonequilibrium models, which are more difficult to falsify. For example, locations with more resources may promote higher speciation rates, and even in communities that are far below carrying capacity, competition should favor diversification to exploit all available resources. We address these issues in the discussion.

Data
We studied elevational breeding distributions of 479 species of birds from 11 orders (Passeriformes, Piciformes, Bucerotiformes, Upupiformes, Trogoniformes, Coraciiformes, Cuculiformes, Psittaciformes, Apodiformes, Strigiformes, and Columbiformes) that are found within a 10,000 km 2 area in the east Himalaya that includes parts of Sikkim, West Bengal, and Arunachal Pradesh in India, as well as Bhutan (Price et al. 2011;. The Passeriformes are the most speciose order, with 374 species present in the east Himalaya, including two suborders, the dominant oscines (or songbirds, 369 species) and the suboscines (five species). We omit the four species of hornbills (Bucerotiformes) from all morphological analyses because their unusual morphology makes them outliers that influence several measures we employ. We estimated species richness at each elevation by overlapping elevational ranges. Ranges were taken from the published literature referenced in Price et al. (2011) but subsequently modified by field experience. First, if we encountered a breeding population outside the published range, we expanded the range to accommodate this observation. Second, if different publications gave different ranges, we sometimes rejected the wider range based on our knowledge of the species. The original data set was compiled in Price et al. (2011). It is being continually updated and, along with other data in this article, has been deposited in the Dryad Digital Repository at https://doi.org/10.5061 /dryad.hmgqnk9bx (Schumm 2019). We also conducted breeding bird censuses in 18 different 5-ha plots covering all elevations White et al. 2019). Estimation of elevational assemblages of species by range overlap (e.g., used by Pigot et al. [2016b]) includes all species, but does not account for the likelihood that some species do not co-occur locally, and is subject to sampling artifacts given that range limits are uncertain (e.g., Price et al. 2014), whereas censuses at select locations (e.g., used by Ricklefs [2012]) sample only a fraction of species present. Here, we use both methods to assess robustness of the results.
Ecology, Morphology, Phylogeny. We compiled specieslevel ecological attributes of all species evaluated from the literature (primarily Ali and Ripley [1972] and Rasmussen and Anderton [2005]) and our own field notes. Habitat type (forest/open habitat) is a qualitative assessment of whether a bird looking up will typically see leaf or sky (Price et al. 2011). Additional categories we recorded were primary diet (insectivore, frugivore, nectarivore, omnivore, carnivore), secondary diet (e.g., a species may be primarily insectivorous but also eat fruit), and primary foraging substrate (ground, bush, tree, aerial); here we focus on primary diet. We measured 475 species in total in museum collections, using male specimens collected as near to the study location as possible White 2016). Apart from 18 rare species, we measured at least two individuals of each species. We measured beak length, beak depth, beak width, wing length, and tarsus length. These traits are reliably measured from museum specimens and Figure 1: Sources of data and hypotheses evaluated to assess drivers of east Himalayan bird species richness. Data presented in this article are categorized into four types. Four hypotheses of ecological drivers of species' distributions are shown in the gray shaded boxes. Beneath each box, we outline some predictions emerging from each hypothesis and the relevant data types that we use to evaluate each. The image for elevational ranges (left) is from White et al. (2019); it is a graphic showing the study area (up to 3,000 m), with the green/gray boundary approximately demarcating the freezing line. The image for beaks is of three species of Phylloscopus warbler (center) from Price et al. (2000). The photo of Phylloscopus affinis (right) was taken by Suresh Kumar and is reproduced with permission. are generally correlated with ecological features: the beak is the main instrument for dealing with different foods, the tarsus is associated with foraging mode, and overall size is associated with prey size (e.g., Miles and Ricklefs 1984;Price 1991;Ghosh-Harihar and Price 2014;Pigot et al. 2016b). We extracted principal components from the correlation matrix of natural-log-transformed species mean morphology. We did not control for phylogenetic nonindependence using principal component analysis (PCA; cf. Revell 2009), as we wished to summarize variation among extant species (in our use, PCA is a dimensionality reduction method rather than a test of significance, where independence becomes an issue). The first four principal components summarize 99% of the variation in morphology, and we interpret these components as biologically meaningful. We standardized these axes to have unit variance, that is, give them equal weighting because dimensions of morphological space that explain relatively little variation in morphology are important in differentiating species ecologically (Miles and Ricklefs 1984;Price 1991;Cooney et al. 2017). Morphology correlates with broad ecological categories, such as those between aerial versus nonaerial foragers and noninsectivorous versus insectivorous species (fig. S1; figs. S1-S8 are available online). We used the time-dated fivegene molecular phylogeny from Price et al. (2014) to incorporate phylogenetic corrections into regression models, to determine species closest relatives, and to estimate diversification statistics. Eight species are missing from that phylogeny, which we added assuming the placement and relative branch lengths given by Jetz et al. (2012). Omitting these species makes no qualitative difference to the results.
Resources. Because the richness patterns in our study are driven by species that are insectivorous (68% of all species in the study are primarily insectivorous), we focused on arthropods as the main food resource and arthropod size distributions as a more easily measured resource axis. We evaluated two data sets. First, Ghosh-Harihar and Price (2014) and Price et al. (2014) assessed arthropod abundance at 10 sites in the east Himalaya by enclosing a small branch in a garbage bag, chloroforming and later sorting through the bag contents, and classifying arthropods by size (30-50 bags at a site). Mousumi Ghosh-Harihar (personal communication) provided the size data, from which we divided arthropods into two size classes with body lengths !6 mm and 16 mm. Second, at our main two bird census sites at 200 m and 2,000 m in west Bengal, K. S. collected new samples in 2014 and 2015, using four methods: Malaise traps, beating, branch clipping, and pitfall traps (Ausden and Drake 2006). We used these methods to survey across a wider spectrum of insects than are picked up in the bags. At both elevations in May-June 2014, we deployed the Malaise trap at seven different locations in each place for a total of 66 hours. In May 2015, we used branch clipping and beating at 10 trees at each elevation (Ozanne 2005;Piñol et al. 2012). In the branch-clipping method, we used a tree pruner to clip a short branch from a tree and let it fall on an upturned 1 m diameter umbrella, and we then collected all the arthropods from the surface of the umbrella and from the clipped leaves. The beating method involved hitting the foliage of a tree with a stick and collecting all the arthropods that fell on to the umbrella after three beats of the stick. We repeated the process once more on another part of the same tree and added the two collections. For pitfall traps, we used plastic bowls (diameter#depth p 114#76 mm; Bioquip catalog no. 2838A) filled to onequarter of their depth with propylene glycol and left out for 24 h. We put up 18 ground pitfall traps each at low elevations and midelevations. Across all methods, we collected and measured the length of 1,349 arthropods at the 200 m elevation and 1,951 arthropods at the 2,000 m elevation.

Analysis
Diversification Rate. Using the phylogeny, we assessed diversification rate and time using statistics shown to be generally robust to species richness (Oliveira et al. 2016). First, we computed mean pairwise phylogenetic distance among species in each 200-m elevational band (using R package Picante; Kembel et al. 2010). Second, we asked whether each species belongs to an actively diversifying clade (i.e., has many close relatives) or lies on a long independent branch. To do this we computed the DR (also known as DivRate) statistic, following Jetz et al. (2012), and took the average DR for all species in the elevational band. Finally, we computed tip net diversification rates in the Bayesian analysis of macroevolutionary mixtures (BAMM; Title and Rabosky 2019) program and took the average for each assemblage. BAMM models switch between evolutionary regimes as a compound Poisson process, allowing timedependent and diversification-dependent patterns, although all rate shifts are assumed to be discrete, an approximation of true shifts in evolutionary rates.
Morphology. Following Pigot et al. (2016b), we used a greedy-search algorithm to ask how packing in morphological space changes as one moves from elevations of lower species richness to those of higher richness. The method removes from the higher-richness assemblage those species not present in the lower-richness assemblage, which when removed result in the largest decrease in the volume of the convex hull, until the volume is at or below that of the lower-richness assemblage (R code is in the data file, available online; 1 for a graphical illustration, see fig. 2).
The number of species removed is the expansion contribution to richness increase, and any remaining species new to the higher-richness assemblage are the packing contribution. Note that the method is agnostic as to the mean position and shape of the convex hull and is concerned only with its volume-we address positional changes in occupation of morphological volume using a different method described below. From the range overlap data, we partitioned the elevational gradient from 100 to 5,700 m into 200-m slices and determined species turnover between adjacent slices. For the breeding bird census data, we ordered census sites by elevation and compared adjacent sites. Packing should also be reflected in smaller mean nearest neighbor differences in morphology in more species-rich communities (Pigot et al. 2016b;Pellissier et al. 2018). We calculated mean nearest neighbor distances using the R package nndist (Baddeley and Turner 2005). For both convex volume and mean nearest neighbor distance, we simulated (for 200 times each) elevational gradients with species randomly drawn without replacement from the total species pool to fill the true 200-m assemblage species richness, and determined at each elevation the range of volume and neighbor distance values generated.
Species may be packed in some regions of morphological volume, yet more loosely positioned in other regions of the same volume. This is not reflected by measures of the convex hull or mean nearest neighbor distance. The sparsely occupied or empty regions of a convex volume are sometimes called holes, but they may lie on the edge of the volume (Blonder 2016). A single species with extreme trait values can thus drastically increase the overall volume occupied by introducing an empty space into the convex volume ( fig. 2). To address the potential effects of holes in our conclusions about morphological packing, we repeated analyses of variation in morphological volume with elevation using the hypervolume, which is the total volume of an arbitrarily complex shape or set of shapes in trait space. The hypervolume quantifies both packing and positional shifts (Lamanna et al. 2014; see the graphical illustration in fig. 2), as volume in the merged pair of communities is compared with volume of each separately. Hypervolumes and hypervolume overlap are calculated by randomly dropping points on to a Gaussian kernel centered on each species with a fixed bandwidth. Here, following Blonder (2016), we used the Silverman method to set the bandwidth to 0.309 in all analyses (the value of 0.309 arises because all traits are standardized to unit variance). We compared adjacent assemblages along the elevational gradient, considering both the volume of the unique morphological space vacated and the volume of unique morphological space added. For all analyses, we used principal component scores computed on the entire data set of species trait means. Results were similar if we used raw morphology rather than principal components ( fig. S2) or individual measurements rather than species means (N p 1,050 individuals, from one to six measurements per species; fig. S3).
We also calculated the variance in species mean morphology. The variance is a measure of the relative evenness

Change in elevation Green species drops out, orange added
Empty space contributes to convex volume , but not hypervolume

Hypervolume analysis
Trait 2 Trait 1 Figure 2: Methods of the convex hull and hypervolume. Each row compares morphology of species in two communities, with the less speciose community on the left. Blue indicates species present in both communities, green indicates species lost from the first community, and orange indicates species gained by the second community. Top, convex hull. The greedy search algorithm (Pigot et al. 2016b) used in figure 4 sequentially drops from the more speciose community the species that when removed reduce the convex hull the most, until the hull has the same or less volume as the less speciose community. In this example, that would be achieved by dropping the top right species. We conclude that the addition of this species was associated with expansion of morphological space and that the addition of the other species represents packing. Bottom, hypervolume representation, where bandwidths circle each point. The size of the hypervolume is the union of the bandwidths, that is, the gray shaded area. The dashed green circle at right represents the species lost with change in elevation; area lost is the white portion within that circle. Area gained by the hypervolume is given by the union of the gray areas that do not intersect either lost or shared species. Gain and loss of particular areas are mutually exclusive as morphological space varies in both size and position. of species mean values in different communities, which may be related to patterns of resource dispersion. If species are added in a restricted portion of the morphological space, the variance will decrease, as species means cluster closer to the grand mean. Therefore, if morphological volume remains constant as species richness increases, the variance will often decrease. Finally, we assessed the morphological difference between sister pairs on the east Himalaya bird phylogeny, focusing on insectivorous species. We plotted the absolute difference between the species pairs' (log-transformed) beak lengths against mean elevational midpoint of the pair. We used the sister pair method to maintain phylogenetic independence and because close relatives are expected to compete most strongly. This sometimes resulted in comparisons between quite distantly related species, and we restricted our analysis to congener pairs or species pairs with a common ancestor dating to !10 Ma. We correlated the magnitude of trait differences with elevation, including age of the species pair as a covariate. Finally, we regressed mean species morphology on midpoint of elevational range and other covariates controlling for phylogeny using phylogenetic generalized least squares (R package caper; Orme et al. 2012).
Generalizing beyond the Himalaya. Data on arthropod sizes along tropical elevation gradients are now available for the Ecuadorian Andes (Guevara and Aviles 2007), Mount Kilimanjaro, Tanzania (Classen et al. 2016;Schellenberg Costa et al. 2017;Tiede et al. 2018), and Mount Wilhelm, New Guinea (Sam et al. 2017). We assessed variation in beak length of insectivorous species to compare against these data. For the Andes, we obtained bird species lists for the two sites (Estación Biológica Jatun Sacha at 400 m and Yanayacu Biological Station at 2,200 m) studied by Guevara and Aviles (2007) from eBird (https://ebird.org /home), with beak morphological measurements provided by Joseph Tobias (personal communication), and insectivory defined as 150% insect diet in the EltonTraits database compiled by Wilman et al. (2014). For Mount Kilimanjaro, following the approach of Quintero and Jetz (2018), we used species elevational distributions and diet information taken from the Handbook of the Birds of the World to determine insectivorous species present in relatively low (1,000-m) and high (2,500-m) bands and obtained morphological data from Hanz et al. (2018). For New Guinea, we used the compilation of Beehler (1981) for both elevational distributions and functional categories to compare sea-level insectivores with those at 2,500 m. Ben Freeman (personal communication) provided morphological data from mist net surveys for 26 species, and we measured an additional 33 species from the Field Museum in Chicago using the same methods (traits measured from nares to tip of the culmen).

Results
In the overlapping range data, species richness reaches a maximum at ∼1,800 m, a peak that is entirely attributable to insectivores ( fig. 3) and matches the peak in arthropod abundance ( fig. S4; see Price et al. 2014). Apart from omnivores, species richness for other diet categories declines monotonically with elevation, but all diet guilds have at least one representative at 3,000 m ( fig. 3). Recent diversification rate does not correlate with the richness peak but instead increases monotonically with elevation ( fig. 3; BAMM results reveal the same pattern; fig. S7). The product of age and recent diversification rate should approximately equal richness, with the only deviation a consequence of rate varying through time. Hence, as expected, mean phylogenetic distance between species in an assemblage decreases monotonically with elevation ( fig. 3).

Morphological Volume
Results based on overlapping ranges are shown in figure 4. Almost 80% of the hypervolume calculated from all 475 species is occupied by species present at the lowest elevation assemblage, where 38% of these 475 species are present. The hypervolume shows little change in total size from the lowest elevations to the richness peak ( fig. 4), but it shifts in position as species drop out and others are added. Some species dropping out have unusual morphologies. They include both oscines (e.g., the little spiderhunter, Arachnothera longirostra) and nonpasserines, such as a bee-eater (Merops orientalis), two nightjars (Caprimulgus spp.), a woodpecker (Chrysocolaptes lucidus), and a parakeet (Psittacula alexandri). Species added also have unusual but different morphologies, including several scimitar babblers (most notably, Pomatorhinus superciliaris), the magpie, Urocissa flavirostris, and several species of cuckoos. A 3D plot of the three beak measurements illustrates the general pattern of addition and loss of species across elevations ( fig. S5). As one moves from the base elevation to midelevations, the net increase of 100 species is correlated with an increase in the convex hull (r p 0:78, across the nine elevational points), but the slope is shallow. Hence, most of the species addition results in increased packing, as can be seen from the mean nearest neighbor distance (which at the richness peak lies outside a null model of community assembly; fig. 4b) and from the relatively minor contributions of expansion or added volume to the convex volume and hypervolume at each 200-m step ( fig. 4a, 4c). Added species are not distributed uniformly within the morphological space. Instead there is a core position occupied by many small insectivorous species that is greatly increased in species richness at the midelevation peak (e.g., fig. 5b).
Coming down to the richness peak from the highest elevations, morphological volume increases rapidly with richness all the way down to the peak at 1,800 m ( fig. 4d). The transition is most striking in association with the increase in richness from 2,200 to 2,000 m, where there is a large increase in volume with minimal space vacated. We can see this transition in three-dimensional beak trait space: many species are added at the fringes of the trait space ( fig. S5). However, packing still contributes substantially.
The results described above are calculated from the range overlap data. We found similar patterns when we used censuses on 5-ha plots as a more direct measurement of a-diversity, with the exception that morphological volume and, to a lesser extent, richness peaked at 1,200 m ( fig. S6). This is at least in part because the range overlap data at 1,800 m include a number of rare and unusual low-elevation species at their upper range limit, and individuals of these species do not appear in the censuses at that elevation but several do at 1,200 m. In summary, as species richness increases from low elevations to midelevations, morphological space stays roughly constant in total occupied volume, despite moving around, as unusual morphologies are lost and gained. On the other hand, as species richness increases from the highest elevation down to the midelevations, expansion of morphological space is much more apparent. We now place these observations in the context of resource measurements.

Measurements of Resources and Morphological Diversity
Previous collections across the entire elevational gradient based on bagging of branches showed that arthropods are most abundant at midelevations, that is, 2,000 m Figure 4: a, Contribution to the convex hull in 200-m elevational slices from species not present at the adjacent site that necessarily cause expansion (number of new taxa added that necessarily increase the volume occupied), with any additional new species driving packing. Elevational slices are analyzed with reference to the adjacent lower slice between 100 and 1,800 m (as richness increases from the lowest elevation to 1,800 m) and with reference to the adjacent higher band between 3,300 and 1,800 m (as richness increases from the highest elevation down to 1,800 m). Note that the convex hull may shift position; here we are concerned only with total volume occupied. b, Convex hull for all species (as the principal component axes are standardized after computation, before their use in convex hull computation, the convex volume units [standard units] are to the nth [5th] power) and mean nearest neighbor distance (MNND) of insectivorous bird species. Transparent bands show the range of these two metrics from the 2.5 to 97.5 percentiles of values calculated from random species draws from the gradient to match richness at each elevation. c, Size of unique volume expanded into and vacated as species numbers increase toward the richness peak at 1,800 m, as in panel a. Volume added accounts for novel trait space occupation at that elevation relative to the lower richness site. Volume vacated accounts for the trait space that is lost as some species are lost from the lower richness site. d, The fraction of the hypervolume occupied by all 475 species that is occupied by species in each 200-m elevational slice. Standard deviations on the points result from 100 replicate runs. Note that the algorithm proceeds by populating Gaussian kernels randomly; hence, the variation.  (2014) is about twice as high at low elevations than at midelevations (10% vs. 5%; fig. 5c). Collections based on beats and clips show a greater disparity: 19% of all arthropods were 16 mm at 200 m, contrasting with 3% at 2,000 m ( fig. 5a). Categorizing size classes as small and large (separated at 6 mm), the binomial variance in arthropod size is two times greater at 200 m than at 2,000 m for the bags and 13 times greater for beats and clips. The reason for the lower variance at 2,000 m is the addition of many small arthropods (large arthropods are similar in number across the gradient; fig. 5a). This skews the distribution, producing clustering close to the mean value. Aerial arthropods show a similar pattern with a nonsignificant tendency toward smaller arthropods at lower elevations ( fig. S8). By contrast, ground arthropod body size frequency distributions, as assessed by pitfall traps, tend to be larger at 2,000 m than at 200 m, with the exception of the largest arthropods (19 mm) we captured ( fig. S8).
Variance in bird beak sizes matches foliage and aerial arthropod size distributions. In the census plots, beak length variance is 1.96 times lower at 2,000 m than at 200 m ( fig. 5b) despite the higher richness, because many small-billed insectivores are added. In the range overlap data, pairs of closely related forest insectivore species differ more in beak length measurements at low than at high elevations ( fig. 5d). Both of these observations are consistent with a morphological response to the more even distribution of the available resource base at low elevations.

Morphological Correlates of Elevational Distribution
Body mass is significantly positively correlated with elevational range midpoint, but elevation explains only a small fraction of the variance (in a nonphylogenetically corrected model, R 2 p 2%; fig. 6a). Open habitat species are generally larger than high-elevation forest species, which drives a change in the trend in body size from weakly negative within forests to weakly positive across the tree line (∼4,000 m), well fit by a breakpoint regression ( fig. 6a). Body size varies little with elevation within the forest oscines and forest nonoscines separately, and the weak negative trend results because the larger nonoscines drop out ( fig. 6b). A measure of shape, the ratio of beak to tarsus (PC2; see table S1; tables S1, S2 are available online), correlates with elevation ( fig. 6c). Birds with smaller beaks and longer tarsi tend to occupy higher elevations. This relationship is retained when forest/open habitat or oscine/nonoscine dummy variables are included in the model (table S2). Not only is shape different, beaks are absolutely smaller and tarsi are absolutely longer at higher elevations (table S2).

Generalization to Tropical Mountains
As in the Himalaya, average arthropod size declines with elevation in Ecuador and New Guinea ( fig. 7). On Mount Kilimanjaro, where the sampled gradient starts at 1,000 m, orthopteran size declines (Tiede et al. 2018), as does bee size (Schellenberg Costa et al. 2017; fig. 7), but moth size increases (Schellenberg Costa et al. 2017; fig. 7). Median beak length of insectivorous bird species declines with elevation across all gradients. The decline in beak length results in a higher variance in beak length at lower elevations than at middle elevations, matching measurements of resource size distributions in each location ( fig. 7). However, richness patterns do not necessarily correspond across gradients. In the east Himalaya richness peaks at midelevations, and the relatively low variance in beak lengths at midelevations is produced because so many similar short-beaked species are found there. By contrast, richness in the Andes peaks at lower elevations (Pigot et al. 2016b). In this case, the number of short-beaked species appear similar at both elevations and instead long-beaked species are relatively more frequent at lower elevations. Because short-beaked species are more abundant than long-beaked ones everywhere, this results in a more even distribution at the mountain base.

Discussion
Along the east Himalayan elevational gradient, differences in species richness are largely associated with increased morphological packing, especially as one moves up from the lowest elevations to the richness peak. However, areas of morphological space that are occupied vary along the gradient, as do the mean and variance in morphology of species, in ways that are consistent with distributions and abundances of available resources; all predictions in figure 1 are met. We argue that these results support resource abundance, diversity, and competition for resources as important controls on breeding bird distributions in the east Himalaya, and hence play an important role in setting richness patterns. In the following sections, we separately consider the roles of packing, resource diversity, determinants of range limits, and competition in driving richness patterns.

Packing
Variation in species richness is not strongly associated with morphological volume: it appears that as species accumulate they pack into the same space, at least between low elevations and midelevations. This pattern has been taken to indicate a role for diversification dynamics in setting richness patterns (Ricklefs 2012;Belmaker and Jetz 2015). Under this idea, some climatic regimes have been associated with relatively high diversification rates, or clades persisting for a longer time, and consequently more species have been produced in these regimes. The alternative ecological limits model is that locations with many species can accommodate more species than other elevations, even without overall niche expansion. While not rejecting a role for diversification dynamics, the east Himalayan data are consistent with predictions of accommodation. Notably, small arthropods are especially abundant at the richness peak. Higher resource abundance, by supporting larger populations, should lower the risk of population extinction and enable finer niche partitioning (Hutchinson 1959;Schwilk and Ackerly 2005). Furthermore, small insectivorous species dominate the bird richness peak, and such species are thought to be able to partition niches more finely than frugivores and other diet guilds, owing to more axes for specialization (Terborgh 1980). A few other studies have searched for and found correlations between total arthropod abundance and insectivorous bird species richness (Rabenold 1978;Ferger et al. 2014;Ghosh-Harihar and Price 2014). A recent review of arthropod abundance patterns across elevational gradients shows them to be quite variable, with some tropical transects having more arthropods at the base than at intermediate elevations (Supriya et al. 2019), and it would be of interest to ask whether these patterns correlate The numbers are from 10 beats and 10 branch clips at each location, sampling for foliage arthropods. Arrows indicate median size classes for the two sites. Standard errors are shown around the size class mean log counts; size classes with significantly different counts of insects between the two elevations as determined by the Wilcoxon test (P ! :05) are indicated with asterisks above the X-axis. b, Histograms of (untransformed) beak length (mm) for insectivore species recorded as present in censuses at these sites. Median lengths are 8 and 13.3 mm for midelevations and low elevations, respectively, as indicated by arrows (Wilcoxon test, P ! :01). Colors correspond to elevations in panel a. c, Proportion of insects in bag surveys across the east Himalaya recorded as being large (16 mm) against elevation. A least squares regression line is fit (P ! :005; R 2 p 0:71; N p 10 locations). d, Natural-log-transformed beak length differences between pairs of closely related insectivorous forest bird species plotted against the mean of each species midpoint elevation range. The linear regression slope is b 1 p 25:7#10 25 (df p 64; R 2 p 0:1; P p :007). Effect size and significance of the elevation predictor remain very similar in a multiple regression with pair age (b elevation p 25:8#10 25 and P p :008, with pair age insignificant, P p :51). In a univariate regression, age remains an insignificant predictor (b 1 p 4:2#10 23 ; R 2 p 0:01; P p :42).
with species richness. For example, arthropod abundance along the Mount Wilhelm New Guinea transect shows a slight decline from low elevations up to 2,000 m (Supriya et al. 2019), which appears to be matched quite closely by insectivorous species richness (Sam et al. 2017).

Resource Diversity
The decline in insectivorous bird species richness from the peak at 2,000 m to the highest elevations is matched by a decline in foliage arthropod diversity, resource abundance, and total morphological volume occupied. Hence, the decline above 2,000 m is consistent with a decrease in niche availability. From 200 to 2,000 m, however, resource abundance and species richness increases, resource diversity and morphological variance declines, and morphological volume (convex hull and hypervolume) remains more or less constant. Across the lower part of the gradient, we suggest that loss of some resources drives the loss of some species (e.g., more species of drongos and bee-eaters, which feed on large flying insects, are present at the low elevations), the gain in small insectivores results in an overall increase in richness and decrease in variance, and the morphological volume is not greatly altered because of an increase in other large insectivores, notably, ground foragers, such as laughing thrushes and scimitar babblers (consistent with the larger number of large [16 mm] arthropods in pitfall traps at 2,000 m than at 200 m; fig. S8).
A decrease in the median and variance of beak length of insectivorous bird species appears quite general along tropical gradients, but in other locations such as the Andes it is correlated with a decline of species richness along the entire gradient. A latitudinal decline in the number of large insectivorous bird species from the New World tropics to temperate North America (Schoener 1971;Greenberg 1981) has also been correlated with a decline in large insects and thus a narrowing of resource diversity (Schoener and Janzen 1968). Other groups in the east Himalaya, notably, frugivores, follow the Andean pattern for insectivores and show a monotonic decline in species richness with elevation ( fig. 3). Kissling et al. (2007) showed that frugivore richness across Africa correlates with resource diversity. This seems likely in the east Himalaya. While we did not measure fruit diversity, tree species richness in 5-ha plots is about twice as high, at 200 m rather than at 2,000 m (Rana et al. 2019).
The greater divergence in size between close relatives at the lowest elevations in the Himalaya ( fig. 5d) is similar to the pattern along the New World latitudinal gradient (Bothwell et al. 2015). Bothwell et al. (2015) suggested that physiological constraints at higher latitudes may slow evolutionary divergence among sympatric species, but the actual mechanism that would cause this is not clear. However, greater divergence along a given resource dimension is theoretically predicted to result from a broader resource distribution (Slatkin 1980). We thus attribute . The solid gray line is a significant segmented regression for all species (one-sided Davies test: P p :029, using the R package segmented [Muggeo 2008]). b, Body mass against elevational range midpoint for forest species only, divided into oscines (N p 274, orange points) and suboscines plus nonpasserines (N p 92, blue points). Lines are least square regressions fit to the two groups separately (both P 1 :05). c, Relationship of shape (principal component 2, correlated with shorter tarsi and longer beaks; table S1) with elevational midpoint. The line is the linear regression (R 2 p 0:05; P ! :0001). See table S2 for phylogenetically corrected statistics, including multiple predictors. the greater differences between related species, at both low elevations and low latitudes, to be a result of a broad distribution of resources (notably, proportionally more large arthropods).

Range Limits and Ecological Limits
Limited elevational ranges are a necessary condition for species richness to vary across elevations. Elsen et al. (2017) evaluated roles of temperature and competition as determinants of the elevational ranges of breeding birds across temperate elevations (above 2,000 m) at two sites in the west Himalaya. Using field censuses, Elsen et al. (2017) examined the importance of competition by testing for reciprocal effects of close congeners on each other's abundance, the role of habitat (considering the distribution of three forest types), and temperature, as measured in the field. All variables correlate with elevation, but Elsen et al. (2017) a Beak length (mm) No. species  Figure 7: Patterns of variation of bird beak morphology with prey insect size disparity. a, Bird beak length distributions in the Andes, using data from J. Tobias (unpublished) and species lists from Cornell University's eBird database, for species present at the Yanayacu Biological Station (2,200-m cloud forest site) and Jatun Sacha Biological Station (elevation, 400 m), the two sites surveyed by Guevara et al. (2007). Beak length distributions differ significantly between elevations (P ! :01) by both t-and Wilcoxon (nonparametric) tests; median lengths are 15.9 and 20.0 mm for the two sites (indicated approximately by arrows). b, c, Range data from Beehler (1981;  reasoned that if temperature were important, a particular species' abundance should be related to temperature in a hump-shaped manner, because maximum abundance would be centered at a thermal optimum, whereas there was no such expectation for habitat or competitive interactions. Perhaps because of this quadratic expectation, they found more support for species being controlled by temperature than by other factors. White et al. (2019) analyzed the bird communities surveyed by Price et al. (2014; as used for fig. S6 of this article) and found a transition in community composition associated with the line of regular freezing (at about 1,700 m). Both these findings imply links to temperature, but we suggest that they operate through biotic features. For example, sharp gradients in habitat across the freezing line may result because plants that invest in adaptations to freezing are poor competitors in nonfreezing environments, so few species traverse that line (Koehler et al. 2011). Furthermore, multiple abiotic and biotic features, including habitat, competitors, predators, and parasites, are all likely to contribute a competitive edge to a species that regularly inhabits a particular elevational range. Other abiotic features such as physiological response to oxygen limitation may limit elevation range in the Himalaya, and this is also likely to contribute to a competitive advantage over ecologically similar species at other elevations (Barve et al. 2016;Ishtiaq and Barve 2018). One way to assess roles of the direct effects of climate, habitat, and competition is to search for morphological correlates with elevation (and, hence, its strong correlate, temperature). Bergmann's and Allen's rules make specific predictions based on physiology and heat conservation of minimized body and appendage surface area in cold climates. In birds, Bergmann's rule has received little support along elevational gradients generally (Freeman 2016), but we found that open-country species above tree line are on average larger than forest species, similar to the pattern across latitudes (Olson et al. 2009). We suggest that the relatively small size of forest species is best explained by resources, because foraging in relatively dense foliage selects for dexterity and small size (Price and Gross 2005). A correlation of the beak-to-tarsus ratio with elevation ( fig. 6c) has also been demonstrated within Himalayan corvids and Phylloscopus warblers (Price 1991;Kennedy et al. 2012;Ghosh-Harihar and Price 2014). Hence, expectations of Allen's rule for smaller extremities at colder temperatures is met for beaks but not for tarsi. However, observed trends in both are consistent with resource distributions. First, beak size correlates with prey size (e.g., Price 1991;Price andGross 2005, Sam et al. 2017), and average arthropod size declines toward higher elevations. Second, relative length of the tarsus increases toward high elevations, where more open habitats should favor hopping over flying movements (Price and Gross 2005).

Competitive Release in Elevational Range on Mountains
Morphological correlates described in the previous section match predictions based on resource distributions, suggesting a role for competition for these resources in setting elevational range limits (Price 1991). As one moves from east to west in the Himalaya, species numbers decline. This leads to the prediction that species common to the west and east should have larger elevational ranges in the west, resulting from competitive release (Terborgh and Weske 1975;Elsen et al. 2017). In fact, species contract elevational range toward the west (Price et al. 2011;Ghosh-Harihar and Price 2014;Elsen et al. 2017). Range contractions appear contrary to the idea that competition constrains a species' elevational range (Terborgh and Weske 1975;Elsen et al. 2017). However, competition operates in the context of available resources. According to our measurements, resources are fewer in the drier west (Ghosh-Harihar and Price 2014) and we therefore predict range contractions even in the absence of competitors. Resources, especially arthropods, are thought to be measured with more error than climate (although quality of climate models remains uncertain; Soria-Auza et al. 2010), or coarse habitat variables, which form the basis for the majority of macroecological studies (Gotelli et al. 2010;Pigot et al. 2016a). Our findings imply that they are essential to a comprehensive understanding of species distributions.
More generally, evidence of a role for competition in setting elevational distributions is quite strong: (1) closely related species often partition elevational gradients (Price et al. 2011;Jankowski et al. 2012), (2) comparative studies of sister species implicate sympatry as the driver of divergence in elevational differences (Freeman 2015), (3) song playbacks demonstrate aggressive interspecific responses (Jankowski et al. 2010;Freeman et al. 2016;Boyce and Martin 2019), and (4) increases in elevational range in the absence of congeners on islands (Lack 1976) and isolated mountaintops (Diamond 1973;Terborgh and Weske 1975) have been regularly observed.

Conclusions
Our analyses of morphological diversity in birds and arthropod abundance and diversity along an elevational gradient in the east Himalaya support necessary conditions for ecological limits acting through resources to shape species diversity patterns. While these patterns are expected if communities are close to saturation and species are limited by resources, we acknowledge that all of the ecological evidence presented could be consistent with nonequilibrium diversification dynamics contributing to differences in richness. For example, more species may be produced at midelevations (i.e., ∼2,000 m) because this is the position of most dispersal barriers, and this is indeed the belt where young allopatric species, separated between the east and west Himalaya, are the most common (White 2016). Once produced, if such species come into sympatry, they could then subdivide niche space more finely, replacing generalists (Ricklefs 2012). However, in the east Himalaya, co-occurring species that belong to more actively radiating clades are not concentrated at the species richness peak ( fig. 3b), but rather occur at higher elevations. Quintero and Jetz (2018) also show the monotonic gradient in diversification rate using a phylogeny based on different molecular data and taxonomic inference. Even at high elevations, species are long separated from each other, suggesting a limit on range expansions set by niche filling . This, in addition to the positive evidence from resource measurements and morphological correlates described here, adds credence to an ecological limits model as a major contributor to patterns of richness across the east Himalayan elevational gradient and potentially elsewhere ( fig. 7).