Occupancy patterns and upper range limits of lowland Bornean birds along an elevational gradient

The traditional view of species’ distributions is that they are less abundant near the edges of their ranges and more abundant towards the centre. Testing this pattern is difficult because of the complexity of distributions across wide geographical areas. An alternative strategy, however, is to measure species’ distributional patterns along elevational gradients. We applied this strategy to examine whether lowland forest birds are indeed less common near their upper range limits on a Bornean mountain, and tested co‐occurrence patterns among species for potential causes of attenuation, including signatures of habitat selection and competition at the periphery of their ranges.


| INTRODUC TI ON
Ecologists have long assumed that species are most abundant at the centre of their ranges and become rarer near the edges (Brown, 1984). This biogeographical principle, often-called the "abundant centre" hypothesis, has been thought to apply in terms of geographical distance, as well as along environmental gradients (Pironon, Villellas, Morris, Doak, & García, 2015;Van Couwenberghe, Collet, Pierrat, Verheyen, & Gégout, 2013). However, the "rule" has found only sporadic empirical support (Sagarin, Gaines, & Gaylord, 2006) and describes only one possible pattern of response to environmental gradients. Recent analyses have found the principle to be poorly supported across a variety of taxa and geographical scales (Dallas, Decker, & Hastings, 2017). Better supported is a less-constrained version of the "abundant centre" hypothesis, the "rare periphery" hypothesis, which predicts that species will become rarer near their range edges without specifying a broader unimodal abundant centre pattern (Gaston, 2009).
Testing these hypothesis over the geographical range of a species can be complicated by a lack of cross-correlation among environmental variables, differences in local adaptation and changes in demographics and species interactions across large distances (Sagarin et al., 2006). Studies have tried with mixed results to minimize these challenges by testing multiple distance metrics (Dallas et al., 2017;Soberón, Peterson, & Osorio-Olvera, 2018), employing multiple definitions of range centres and margins (Santini, Pironon, Maiorano, & Thuiller, 2018) and measuring environmental rather than spatial distances (Martínez-Meyer, Díaz-Porras, Peterson, & Yáñez-Arenas, 2013;Pironon et al., 2015;Van Couwenberghe et al., 2013). In addition to these inherent complications, comparisons of abundance estimates across geographical distances are often difficult because sampling usually depends on different observers operating at different times using different methods (Santini et al., 2018).
Consequently, most studies have provided only partial tests of these hypotheses, limited by both the difficulties listed above and a lack of data from the edges of species distributions (Santini et al., 2018).
Many species do, however, appear to be more abundant at sites predicted to be more suitable by species distribution models (Lunghi et al., 2018;Weber, Stevens, Diniz-Filho, & Grelle, 2017).
Given these challenges, a more straightforward way to test the rare periphery hypothesis, albeit on a smaller spatial scale, is to examine species occupancy along elevational gradients. Such gradients confer several advantages: habitat is usually continuous rather than patchy, and it tends to co-vary consistently with abiotic gradients in temperature, precipitation and humidity (Malhi et al., 2010).
Distances are small enough that, when combined with the relative lack of habitat patchiness, effects of dispersal limitation should be minimized (Holt & Keitt, 2000), spatial distance metrics should be straightforward and sampling protocols can be easily standardized (Santini et al., 2018).
In one such test, lowland bird species in New Guinea displayed a variety of abundance patterns across their elevational ranges; peak abundance occurred in the lower quarter of the ranges of 60% of lowland species, but several were most abundant in the middle or upper ends of their ranges (Freeman & Beehler, 2018).
Thus, the New Guinea study provided some support for the rare periphery hypothesis and was consistent with the idea that some lowland species occupy a "truncated niche" at the lower limits of their distribution (Feeley & Silman, 2010). Freeman and Beehler (2018) interpreted the pattern to mean that the current thermal niche occupied by these species does not represent the full range of temperatures that they could inhabit. Indirect effects of temperature increases predicted by climate models will likely also impact many species; in fact, lowland tropical forests are already being affected (Becek & Horwath, 2017). Lower range limits of endotherms may therefore be affected by climate change even when temperature increases per se do not exceed their range of thermal tolerance.
In contrast to lower boundaries of bird distribution, upper range limits may be set by different processes. For example ecotones (Terborgh, 1985) or a sharp increase in abundance of a competitor (Diamond, 1973) may cause a species to be common up to an abrupt range edge. Alternatively, for species that fit the "rare periphery" hypothesis (Gaston, 2009) near their upper range limits, this could be due to the gradually increasing abundance of a competitor, or a linear responses to a habitat gradient, both of which could cause occupancy to decline slowly unless abrupt response thresholds exist (Doughty & Goulden, 2008).
Co-occurrence patterns among species are often examined to provide clues to the mechanisms responsible for community structure (Colorado & Rodewald, 2015;McCreadie & Bedwell, 2013).
For example lowland species may occur near their upper range limits only in patches of suitable habitat or microclimate, which may become increasingly rare at higher elevations. In this case lowland species near their upper range limits would be clustered with one another in these patches. Alternatively, if competition with montane species limits ranges, then the upper edge of lowland species distributions should vary according to the occurrence of the montane competitor (Burner et al., In Revision).
We surveyed 50 lowland bird species along an elevational gradient of primary forest in Borneo to determine the relationship between elevation and occupancy (as a proxy for abundance) in each species. We then modelled the relationship for each species and used these models to test the "rare periphery" hypothesis and to compare relative occupancy at lower distribution limits to those along the rest of the gradient. To determine potential causes limiting the upper distributional ranges of these lowland species, we examined co-occurrence patterns of species pairs and tested two hypothesized mechanisms: (a) patchy habitat, that is that multiple lowland species will cluster near their upper range edges, presumably in pockets of suitable habitat or microclimate; and (b) interspecific competition, that is that a lowland species' upper limit is determined by the lower limit of likely montane competitors (species of the same genus or foraging guild), and that these pairs of species will be segregated from each other where ranges overlap.

| Study site
Mt. Mulu (4.045° N 114.929° E) is a 2,376 m peak in Sarawak, Malaysian Borneo ( Figure 1). It is one of few locations in Borneo where unbroken primary forest stretches from near sea level to the summit. Selective logging has occurred only up to 50 m elevation at the base of the mountain in the river floodplain, and the whole mountain is now protected within Mt. Mulu National Park. Most of Borneo's montane bird species, including many of the island's endemics, are found on Mulu (Burner et al., 2016). We conducted 560 avian point counts at 175 points at 50, 300, 600, 900, 1,200, 1,500 and 1,800 m elevation from June to September of 2014 (Burner, Styring, Boer, & Sheldon, 2018). Points were surveyed an average of 3.2 times (SD = 0.94).

| Avian surveys
Bird surveys followed the standardized protocol outlined by Burner et al. (2018). Briefly, 6-min aural and visual counts were conducted every 200 m along 5 km transects at each sampled elevation. At each elevation about 25 points were sampled, and each point was surveyed three to four times over a 1-week period. At each point, sound samples were recorded using a Marantz digital recorder and Sennheiser microphone for later species identification, which improved species detection rates in this species-rich environment (Haselmayer & Quinn, 2000). Counts were conducted from 06h:00 to 10h:30, and never in the rain. Repeated visits to each site were used to increase overall detection probability. Time of day, elevation, weather, latitude and longitude were recorded for each count.
A presence/absence data matrix for each species was generated from survey data. Absolute abundance of birds is difficult to obtain, so we used naive occupancy as a proxy for abundance. Naive occupancy is a minimum estimate of site occupancy, obtained by dividing the number of sites at a given elevation where the species was detected by the total number of sites (MacKenzie et al., 2002).
It provides a comparable measure of relative abundance assuming detection probabilities and abundance within occupied patches do not differ by elevation.

| Model fitting of elevation-occupancy relationship
To assess the shape of species response curves along the elevational gradient, we used the package 'eHOF' (Jansen & Oksanen, 2013) in R (R Core Team, 2018) to fit each of the five Huisman-Olff-Fresco (HOF) models (flat, monotonic, plateau, symmetric and skewed; labelled types I-V respectively; Figure 2). These models have been shown to characterize species' gradient responses effectively (Oksanen & Minchin, 2002), even with as few as 10 detections (Freeman & Beehler, 2018). Lowland species (based on elevational range, as designated by Burner et al., 2018) detected at 10 or more points that had their entire elevational range within the sampled elevations (Sagarin & Gaines, 2002) were selected for model fitting. Montane species were excluded because most occurred above the 1,800 m limit of our sampling, preventing their full ranges from being considered. Models were fit to the presence/absence matrix for each species including all points up to 1,800 m; additional zero counts above the lowest zero count for a species do not influence HOF model selection, and all lowland species were present at or near the bottom of the gradient. Models were ranked based on AIC c score, and we report all models within 2 AIC c units of the top model (Burnham, Anderson, & Huyvaert, 2011). Exploratory analyses revealed that the type of curve selected was not related to number of detections for a given species (ANOVA; p = .35), nor did the proportions of curve types selected among frugivores (Proportion test; p = .59) or insectivores (p = .88) differ from the full set of species.
Both monotonic and symmetric (types II & IV) curves were considered consistent with the "rare periphery" hypothesis near the upper range limit. Relative occupancy of each species at the mountain's base elevation (50 m) was calculated as (% naive occupancy at sea level)/(maximum % naive occupancy at any elevation). Species were considered likely to occupy a truncated thermal niche in the

| Species co-occurrence
We determined the co-occurrence of species pairs near their range edges using the Pairs software (Ulrich, 2008). Species pairs (Tables   S2 and S3) were chosen for analysis based on hypotheses of (a) selection of patchy habitat at range edges and (b) competition between lowland and montane species at points of contact, as outlined below.
The Pairs software uses null model simulations to compare observed and expected values of co-occurrence for each species pair in the dataset, calculating a "C-score" for each pair. This measure of co-occurrence reflects the number of "checkerboard units"-pairs of sites in which both sites are occupied by a different member of the species pair-for each pair of species (Gotelli & Ulrich, 2010). Standardized Z-scores are used to reflect the difference between an observed Cscore (calculated) and the expected (simulated) score under the null model, such that negative scores indicate aggregated species pairs, whereas positive scores indicate segregation. For our co-occurrence matrix, we pooled bird detections from all visits to each point, resulting in an overall presence/absence matrix; all visits typically occurred within a single week, and local bird occupancy has been shown to be consistent over such short time periods (Royle & Dorazio, 2008).
The choice of an appropriate null model with which to compare observed co-occurrence values can influence results. Two types of null model simulations have been shown to minimize Type I and Type II error rates (Gotelli, 2000). SIM9 constrains the simulations by fixing the number of occurrences of each species (row totals in the input matrices), as well as the species richness of each site (column totals), to that found in the original data. SIM2 also fixes the number of points at which each species is present, but allows the relative species richness of each site to vary. Fixing the relative occupancy of each species reflects the biological observation that some species are consistently more common than others. The choice of whether to allow site totals to vary (as in SIM2) depends on whether sites differ considerably in size, habitat and sampling effort (Gotelli, 2000).
We report results from both simulation types for comparison.

| Patchy habitat hypothesis
To test the hypothesis that lowland species near their upper range limits are clustered together (perhaps in sites with suitable habitat or microclimate), "core" and "peripheral" elevations were designated for each species based on the following criteria. A species' core elevation was the elevation at which it had the highest occupancy, and its peripheral elevation was designated as the highest elevation at which its relative occupancy was both <0.5 and >0.1. The first of these restrictions ensured that species that were equally abundant across their elevational range were not included, whereas the latter ensured that species were not so uncommon in their designated peripheral elevations that statistical co-occurrence inferences were impossible.
For all species pairs with the same "core" and "peripheral" elevations (Table S2), we compared the pair's co-occurrence scores between these two elevations. Comparing co-occurrence of the same pairs between their core and their peripheral elevations, rather than simply examining co-occurrence rates in peripheral elevations alone, allowed correction for the possibility that some species pairs will be segregated or aggregated in their core elevations (due to species interactions or other factors) as well. A t-test was used to determine significance of the overall combined difference in co-occurrence values between core and peripheral elevations for all pairs. Individual pairs were selected for further F I G U R E 2 Examples of Huisman-Olff-Fresco (HOF) models of species abundance using representative species from Mt. Mulu: flat (I), monotonic (II), plateau (III), symmetric (IV) and skewed (V). For each type, the number (and percentage) of lowland species best fitting each curve in our study based on AIC c values is given. When all models within 2 AIC c units of the top model are considered, the percentages are similar (Table S1) examination when differences between their core and peripheral C-scores exceeded 2.0. This threshold does not correspond exactly to a p-value of 0.05 but is a reasonable cut-off based on the trade-offs associated with identifying significant pairs (Gotelli & Ulrich, 2010). Our results were robust to a range of values (1.0-2.5) used for this threshold.

| Species interactions hypothesis
Finally, co-occurrence values were calculated for all species pairs at a given elevation consisting of one lowland and one montane member of the same guild or genus (Table S3), following Burner et al. (2018).
This allowed us to test the hypothesis that lowland species are segregated from likely montane competitors relative to the null model, a potential signature of competitive exclusion. Guild assignments followed Johns (1989) and Myers (2009).

| RE SULTS
In total, 187 bird species were detected during 560 point counts (Burner et al., 2016), including 84 species present at 10 or more points. These latter included 50 lowland and 29 montane species, as designated by Burner et al. (2018).

| Model fitting
We fitted HOF models to the 50 common lowland species, and report all models within 2 AIC c units of the top model in Appendices 1 & 2. Fifteen (30%) were best fit by a symmetric unimodal (Type IV) pattern ( Figure 2), predicted by the abundant centre hypothesis, and one species best fit an asymmetric (Type V) pattern. An additional 34 species (68%) either declined monotonically with elevation or had a low-elevation plateau (Type II or III respectively). No species had a flat distribution. When all models within 2 AIC c units of the top model were considered, the relative prevalence of each model type remained similar (Table S1). Most species (n = 34; 68%) had a relative occupancy of >0.8 at sea level ( Figure 3). We detected 29 montane species (Appendix 3) at 10 or more points to be used for co-occurrence analyses.

| Patchy habitat hypothesis
Many of the common lowland species (n = 39; Types II & IV curves) declined gradually in occupancy near their upper range limit, in keeping with the "rare periphery" hypothesis. This decline was gradual enough that they were sufficiently abundant at a "peripheral" elevation to compare the "core" and "peripheral" co-occurrence values for pairs. We compared co-occurrence scores for all 94 species pairs (Table S2) for which both pair members had the same core and peripheral elevation. Five pairs (5.3%) were significantly more aggregated at their peripheral elevations than expected by chance, whereas nine (9.6%) were significantly more segregated from each other (using SIM9). When the difference in C-score from core to peripheral elevations was averaged for each species when paired with all other suitable species, no species differed significantly from zero.
When averaging across the community as a whole, these 94 species pairs were on average neither more clustered nor more segregated at their peripheral elevations relative to their range core (p = .295).
Results were similar when SIM2 was used, except that one additional pair was more segregated in the periphery.
Both are sallying insectivores that occasionally glean foliage for insects as well, but the body-size difference between the two makes it unlikely that there would be much prey size overlap. Segregation values for this pair were similar between SIM9 (2.37) and SIM2 (2.02).

| D ISCUSS I ON
While some lowland species had the highest occupancy in their elevational range centres, most were at peak occupancy at sea level and declined gradually with elevation. These results provide broad support for the "rare periphery" hypothesis. They also imply that most lowland tropical birds in Borneo are not currently near their thermal maximum, although indirect effects of climate change could still cause range shifts in these species. At their upper range edges, cooccurrence patterns of lowland species with each other and with potential montane competitors did not differ from the null expectation.

| Abundance distributions
The variety of abundance distributions we observed indicates that no single pattern holds universally. Just under one third (30%) of lowland species we surveyed fit a HOF model Type IV symmetric unimodal abundance pattern ( Figure 2). This finding is similar to that for lowland birds in New Guinea, where 33% showed the same pattern (Freeman & Beehler, 2018). Many surveys of montane fauna stop below the highest elevations, but comparisons across the entire elevational gradient would allow evaluation of these patterns in montane as well as lowland species, especially in the Andes and Himalayas, where mountains are especially high and occupied by several bird community layers.
Nonetheless, most species (75%) in our study had low occupancy near their upper range limits. This is consistent with a "rare periphery" hypothesis (Gaston, 2009), which is a less-constrained version of the abundant centre hypothesis in that it predicts low abundance near range edges without specifying the broader pattern. Gradients in many abiotic factors (e.g. temperature, rainfall) are gradual and correlated with elevation (Lomolino, 2001;Lunghi et al., 2018). This means that elevational distance should be a better proxy for environmental distance than is spatial distance over large geographical areas, making "rare periphery" patterns theoretically more likely when measured elevationally (Brown, 1984;Van Couwenberghe et al., 2013) along elevational gradients. However, biotic factors can still change abruptly along smooth abiotic gradients if tolerance thresholds (particularly of plants) are crossed (Doughty & Goulden, 2008).
Abundance of animal species can also change rapidly in the absence of sharp changes in environmental factors (Sagarin et al., 2006).
A majority of the lowland bird species we tested (68%) had the highest occupancy at sea level (Figure 3 shift their ranges upwards, most lowland protected areas are not sufficiently connected to intact higher-elevation habitats (Scriven, Hodgson, McClean, & Hill, 2015). Furthermore, occupancy and abundance levels can be more sensitive to temperature than are range limits; as a result, some species may become less (or more) abundant even in the absence of obvious habitat changes (Dulle et al., 2016). Specialists, for example may be more sensitive to changes than generalists (Jiguet et al., 2010) and may respond differently.
Counterintuitive responses to increasing temperatures can also result from changing competitive interactions (Tingley, Koo, Moritz, Rush, & Beissinger, 2012) as distributions and abundance of some species change. Temperature does appear to affect the range limits of montane forest birds in Borneo (Burner et al., In Revision), many of which appear to have shifted over time (Harris et al., 2012).
Tropical montane species, whose habitat is slightly more secure from logging and development because of the difficulty and expense of working in mountains, are probably in greater danger from climate change because of limited space for upward movement.
Additional work on the physiology and competitive interactions of both lowland and montane tropical species would help to refine these conclusions.

| Causes of occupancy declines at range edges
Most species (75%) had low occupancy near their upper range limits.
These "cold-edge" range limits could be set by a variety of factors, including habitat patchiness, biotic interactions, or abiotic factors (Louthan, Doak, & Angert, 2015). We hypothesized that if lowland species were limited by suitable habitat and microclimate near their upper range margins (Holt & Keitt, 2000), then they would be more aggregated with each other in those sites (perhaps where patches of suitable habitat occurred) than in the core of their ranges. However, most such species pairs were no more segregated in the periphery of their elevational distribution than in the core, and for the minority of pairs that differed from random expectations, slightly more were segregated (n = 9) than were aggregated (n = 5). It thus appears either that patchy habitat does not play a primary role in limiting distributions of these species at this scale within forest ecosystems, or that lowland species differ from one another enough in habitat requirements that they are not restricted to the same sites in their range peripheries. It is possible, for example that each species extends upwards into patches of habitat that contain some key feature common on lower slopes, but that these features vary by species and do not necessarily co-occur spatially. Or, lowland species may be less abundant near their range limits not because suitable habitat is patchy, but because overall quality of habitat is lower for them (Holt & Keitt, 2000).
Another possible explanation for range limits is negative species interactions. Evidence that competition can limit species ranges on mountain slopes is available from other regions (Elsen, Tingley, Kalyanaraman, Ramesh, & Wilcove, 2017;Jankowski et al., 2013;Romdal & Rahbek, 2009). Similarly, based on comparisons of elevational ranges of bird species and their potential competitors (from the same genus or foraging guild) on Bornean mountains, there is evidence that competitive interactions likely play a role in determining range limits of both lowland and montane species (Burner et al., In Revision). However, we found that few lowland-montane species pairs differed in co-occurrence scores from the null predictions, even when both pair-members were in the same genus or foraging guild.
Moreover, recent simulations have shown that competition could, under different conditions, produce either segregated or aggregated patterns (Ulrich, Jabot, & Gotelli, 2017). For example the several significantly aggregated pairs of barbets in our study are likely to be at least occasional competitors and were probably concentrated around a food source, given that they are wide-ranging species that move according to unpredictable fruiting events (Lambert & Marshall, 1991).
Only one pair of insectivores was significantly segregated, a drongo and a stenostirid flycatcher. This could reflect competition that causes these species to avoid each other, but this seems unlikely because of their large difference in body size. Such segregation may be effected by direct asymmetric aggression rather than through resource competition (Terborgh, 2015), with the lower-elevation pair member often thought to be more aggressive (Freeman & Montgomery, 2015;Jankowski, Robinson, & Levey, 2010;Martin & Martin, 2001). Of course, non-random co-occurrence patterns can also be caused by unmeasured variables rather than by direct interactions (D'Amen, Mod, Gotelli, & Guisan, 2018).
The lack of significant positively or negatively co-occurring pairs in this bird community may indicate that strong pairwise competition is either uncommon or is difficult to detect due to the mixture of aggregating and segregating forces that can result from competition.
Alternatively diffuse competition within this assemblage of birds may be more common than intense pairwise competition (Jankowski et al., 2012), making it difficult to detect using most common survey and analysis methods. Our lack of detailed ecological knowledge of most species makes even hypothesizing likely networks of competitors difficult (Jankowski et al., 2012), beyond the relatively coarse metric of foraging guild or at the limited scale of genus. Diffuse competition may often result in niche differentiation based on foraging manoeuvres and strata rather than area occupancy, even among morphologically similar species with similar diets, as has been found for babblers (Timaliidae) in Malaysia (Mansor & Ramli, 2017;Styring, Ragai, Zakaria, & Sheldon, 2016). Competition may also be difficult to detect because elevational range overlap of strongly competing species may be effectively eliminated.
The results of this study are largely congruent with those of Freeman and Beehler (2018), but provide important information on species ranges in Borneo, an island with very different characteristics than New Guinea where they worked. Due to extensive alteration of lowland forest habitat, opportunities to investigate species elevational ranges in a natural setting in Southeast Asia are now extremely limited. Gaining a clearer understanding of habitat requirements, physiological limits and species interactions is key to gaining a more in-depth understanding of abundance distributions and critical to accurately predicting species responses to habitat alteration and climate change.

Sarawak Forestry Department, Forestry Corporation and
Biodiversity Centre and Parks kindly provided research permits.
Andrew Siani identified many difficult species in our recordings. Best fit HOF models for 50 lowland bird species on Mt. Mulu, Borneo. All species with detections from at least 10 points were included for modelling. Package 'eHOF' in R (Jansen & Oksanen, 2013)  List of 50 common lowland bird species (Burner et al., 2018) from Mt. Mulu, Borneo for which HOF models were fit. Each species was detected at 10 or more points. Package 'eHOF' in R (Jansen & Oksanen, 2013) was used to fit each of the five Huisman-Olff-Fresco (HOF) models (flat, monotonic, plateau, symmetric and skewed; I-V respectively). Guild designation include Arboreal Foliage-Gleaning Insectivore (AFGI), Arboreal Foliage-Gleaning Insectivore Frugivore (AFGIF), Arboreal Frugivore (AF), Nectivorous Insectivore/Frugivore (NI/F), Sallying Insectivore/Sallying Substrate-Gleaning Insectivore (SI/SSGI) and Terrestrial Insectivore/Frugivore (TI/F) based on Johns (1989) and Myers (2009). A lowland species' "core" elevation was the elevation at which it had the highest relative occupancy, and its "peripheral" elevation was the highest elevation at which its relative occupancy was both less than 0.5 and greater than 0.1. Species without "core" or "peripheral" elevations listed did not have a "peripheral" elevation that fit the above criterion. Each pair of species with similar "core" and "peripheral" elevations is listed, with co-occurrence scores, in Table S2