**The Coupled Influence of Thermal Physiology and Biotic Interactions on the Distribution and Density of Ant Species along an Elevational Gradient**

#### **Lacy D. Chick 1,2,\*, Jean-Philippe Lessard 3, Robert R. Dunn 4,5 and Nathan J. Sanders <sup>6</sup>**


Received: 31 October 2020; Accepted: 26 November 2020; Published: 30 November 2020

**Abstract:** A fundamental tenet of biogeography is that abiotic and biotic factors interact to shape the distributions of species and the organization of communities, with interactions being more important in benign environments, and environmental filtering more important in stressful environments. This pattern is often inferred using large databases or phylogenetic signal, but physiological mechanisms underlying such patterns are rarely examined. We focused on 18 ant species at 29 sites along an extensive elevational gradient, coupling experimental data on critical thermal limits, null model analyses, and observational data of density and abundance to elucidate factors governing species' elevational range limits. Thermal tolerance data showed that environmental conditions were likely to be more important in colder, more stressful environments, where physiology was the most important constraint on the distribution and density of ant species. Conversely, the evidence for species interactions was strongest in warmer, more benign conditions, as indicated by our observational data and null model analyses. Our results provide a strong test that biotic interactions drive the distributions and density of species in warm climates, but that environmental filtering predominates at colder, high-elevation sites. Such a pattern suggests that the responses of species to climate change are likely to be context-dependent and more specifically, geographically-dependent.

**Keywords:** ants; community structure; physiology; interactions; temperature

#### **1. Introduction**

One of the most striking patterns in nature is that the number of species varies, often systematically, along environmental gradients. Explaining this variation has attracted the attention of ecologists and biogeographers for decades [1], if not longer [2], and has inspired empirical studies in fields ranging from physiological ecology to macroecology and global change biology [3]. However, why does the number of species that coexist in a particular assemblage vary? One possibility is that, broadly speaking, species differ in how they respond to biotic and abiotic factors along environmental gradients, and these differences among species, in turn, influence abundance, distribution, community composition, and broad-scale patterns of diversity. For instance, temperature tends to decrease systematically with elevation and latitude [4,5] and as a result, the abiotic environment at high-elevation and high-latitude sites might be more physiologically stressful for potential colonizers than at low-elevation and low-latitude sites. In such a model, temperature acts as a filter, permitting the occurrence of only those species with traits that allow them to persist at low temperatures [6–8].

Of course, multiple factors can and do simultaneously operate to shape communities, and different factors might be more important in different locations [9–12]. Wallace (1878), Dobzhansky (1950), and Fischer (1960) all suggested that negative interspecific interactions (competition, predation, parasitism) might be more intense in benign, stable environments [13–15]. Indeed, a growing number of investigators have explored the geography of biotic interactions [16] with some studies suggesting that negative interactions might limit the distributions of species and pose a cap to the number of species that can coexist in benign environments [11,17]. Two studies along elevational gradients hint at such a scenario: in hummingbird assemblages in the Andes [18] and in ant assemblages in the U.S. and Europe [19], there is evidence that interspecific interactions shape community membership at low elevations, but that more stressful environmental conditions (e.g., cold temperatures) shape communities at high elevations. Such studies are important because they suggest a mechanism, but do so based on community phylogenetic approaches, which rely on numerous underlying assumptions and can give misleading answers about the processes that actually structure communities [9,20].

More compelling evidence for geographic variation in the relative influence of climate and biotic interactions on the species in assemblages might come from field-based measurements of physiological tolerances (e.g., [8,21–23]) and/or detailed studies of the outcomes of interactions among species, i.e., actual measurements of individual-level functional traits and observations of interactions in the field [24–27]. Such studies, however, are rare due to logistical constraints and many traits that are often measured do not directly relate to tolerance of the abiotic environment.

Like other ectotherms, ants exhibit thermal sensitivity, and species differ in their thermal tolerances (i.e., the ability to tolerate either extreme temperatures or a broad range of temperatures;[28–31].Thermal tolerance in ants may be related to total abundance and range size [32–34], foraging activity [35,36], and broad-scale patterns of diversity [8,32,34,37–39]. If a species occurs at all locations with suitable environmental conditions, then the environment alone would be the sole driver of its distribution. However, if the observed range of a species is smaller than its expected range based on environmental tolerance alone, then some other factor, such as competitive interactions or dispersal limitation, acts to shape the distribution of species among local communities along the gradient [40,41]. If the same suite of factors affects the distribution of many species, then such factors are expected to also influence the distribution of diversity, as diversity is simply a collective property.

Competitive interactions are widely thought to influence the structure and dynamics of some local assemblages, and might shape broad-scale patterns in the distribution of species as well. Competition likely structures local ant assemblages [42], yet its effects are mediated by temperature altering interactions between dominant and subordinate species [28,35,43] and the activities of particular species [27,35]. Here, we examine co-occurrence patterns among ant species along an environmental gradient to assess whether species occur less than expected by chance as would be expected if competition structures communities [44]. We note that considerable debate still exists about whether co-occurrence patterns are evidence of ecological interactions (e.g., [45]). However, we use them here as only one several lines of evidence to ask a series of inter-related questions about the factors that govern the distribution, abundance, and density of ant species along an extensive elevational gradient. Specifically, we ask (1) Are abundance and species density related to environmental temperature? (2) Does thermal tolerance predict elevational range size and species density, as would be the case if temperature were the sole driver of species distributions? and (3) Do species co-occur less among assemblages than would be expected if temperature alone limits membership? Based on the suspicions of early biogeographic pioneers (e.g., [13–15]), we predicted that physiological constraints would limit community membership at high-elevation sites, filtering species that have the physiological capacity to withstand more extreme temperatures, but that interspecific interactions shape assemblages at lower elevation sites that are more environmentally benign. We further predicted that thermal tolerance would be the best predictor of the occurrence and number of species in high elevation communities where environmental filtering predominates, but not at low elevation where biotic interactions are most frequent and intense.

#### **2. Materials and Methods**

#### *2.1. Sampling*

This work was done in Great Smoky Mountains National Park, TN, USA at sites that were situated in mixed hardwood forests and were located in areas away from roads, heavily visited trails, or other recent human disturbances. All sites were highly suitable for a variety of ant species, with similar vegetative cover and habitat to account for differences in habitat specificity of the ant species in the area. We systematically sampled 29 sites (from 379 to 1828 m) in June to August 2004–2007. These sites had a total temperature range of −8.0 to 29.7 ◦C (mean annual temperatures ranged from 7.7 to 13.3 ◦C) and ranged 1308–1928 mm in annual precipitation, respectively. Winkler samplers were used to extract ants from the leaf litter in 16 1-m<sup>2</sup> quadrats at each site in a haphazardly placed 50 <sup>×</sup> 50 m quadrat. At each site, species density is the observed number of species collected in the 50 × 50 m quadrat, and abundance is the number of 1-m<sup>2</sup> quadrats in which any species was detected. This estimate of abundance, which is actually "occurrence" [9,46–49] is preferable to a count of worker number because ants are social, and because counts of colonies is challenging when species have multiple nests per colony and occur in the leaf litter. We differentiate abundance from occurrence because our measure of abundance combines all species, whereas occurrence implies the presence of only a single species. At eight of the sites, ants were also collected using an array of ten pitfall traps over two years [50]. The number of species collected by pitfall traps did not differ from the number collected by the Winkler samplers (paired *t* = 1.88, *n* = 8, *p* = 0.11). Similarly, the fauna sampled by the pitfall traps was similar to the fauna sampled by the Winkler samplers [50]. At most of the sites, an asymptotic species richness estimator (Chao2 in this case) plateaued, suggesting that sampling within sites approached completeness [51]. Moreover, a Chao2 estimate of richness among all sites suggests, at least using these sampling techniques at similar sites, that there would be approximately 45 species in total, and we captured 38 species in our systematic sampling (Table 1). Therefore, we suggest that these communities are adequately sampled, while acknowledging that we undoubtedly missed some subterranean species, species that might be partially arboreal, or other exceedingly rare species in these habitats. Nevertheless, we systematically sampled each site.

In July 2011 and 2012, 31 sites were visited (from 375 to 1825 m; 17 of which were in the previously sampled sites in 2004–2007) in order to collect live individual ants for physiological tolerance estimates. At each of these 31 sites, the same Winkler extraction methods as in the previous sampling were used to extract ants from the leaf litter. However, litter was collected from only 10 1-m2 quadrats per site instead of 16, and live ants were extracted from the leaf litter by sifting through the litter in the field rather than returning them to the lab to use Winkler extractors, as is typically the case with Winkler sampling methods. These modifications were made because we were not aiming to sample the entire community and because live specimens were needed. Finally, we also baited for ants by placing laminated index cards stocked with ~5 g of tuna in oil and individuals were hand-collected at the site. For any species detected either at the bait, the Winkler extraction method, or in general hand collecting, ten live individuals were obtained and returned to the lab (1–2 h from the field site) to conduct thermal tolerance assays.

#### *2.2. Assessing Thermal Tolerance*

Critical thermal minima (CTmin) and maxima (CTmax) were used to examine the physiological constraints imposed on species across the environmental gradient. For each species collected at each site, thermal tolerance assays were performed on 5 individuals for CTmin and 5 individuals for CTmax, which were estimated by documenting the temperature at which individuals lost the ability of righting response. Loss of righting response is measured as the point in which an organism is flipped on its dorsum and can no longer independently right itself. This measure is considered an ecologically relevant endpoint for physiological tolerance because as an organism becomes incapacitated, it can no longer forage or escape predation [52,53].

We used methods described in Warren and Chick (2013) to estimate thermal tolerance for each species at each of the 31 sampled sites [33]. Individuals were transferred to 16 mm glass test tubes, plugged with cotton to reduce thermal refuges, and were placed in an Ac-150-A40 refrigerated water bath (NesLab, ThermoScientific, IL, USA). Water bath temperatures were raised or lowered at a rate of 1 ◦C min−<sup>1</sup> until thermal tolerance was reached. Thermal tolerance was characterized as the highest and lowest temperatures at which an individual could no longer retain locomotor ability, respectively. One vial contained only a copper-constantan Type-T thermocouple (Model HH200A, Omega, CT, USA) and was used to monitor temperature inside the tubes and to ensure accurate readings. All tolerance tests were performed within 5 h of field collection to reduce potential acclimation to the lab thermal environment; however, a subsequent common garden experiment indicated no effects of acclimation on thermal limits. A mean temperature of the loss of righting response served as the index for thermal tolerance for each species at each site. All ants were preserved individually in 2.0-mL vials containing 95% ethanol, and placed in NJS's private collection at the University of Tennessee.

#### *2.3. Are Abundance and Species Density Related to Environmental Temperature?*

Current mean temperatures (1950–2000) were extrapolated for each site from the WorldClim database [54] at a resolution of 30 arc-seconds (758 m distance at a latitude of 35◦ N). Previous work in this region [4] modeled climate based on empirical data collected from a 120-sensor temperature logger network. While data from the 120-sensor network are more fine-scaled, they were not used in this study because they were collected for a shorter time period (2005–2006) and because temperature measured in the data loggers is related to elevation in much the same way as WorldClim data. Similarly, data from weather stations arrayed in the region indicate that temperature declines in a manner comparable to the model used by WorldClim. For these reasons, the WorldClim dataset was used here so that our findings may be more comparable to studies along other gradients where fine-scale resolution may not exist.

Total abundance (the total number of 1-m<sup>2</sup> quadrats in which a species was detected, combined for all species) and species density (the number of species in the 50 × 50 m quadrat) were plotted against mean annual temperature (MAT; we note that MAT was strongly correlated with both January minimum and July maximum). Least squares regressions were used to examine the relationship between temperature and total abundance as well as the relationship between temperature and species density. If temperature is an important determinant of species density and abundance, and colder temperatures filter species from the regional species pool, we would expect to find a linear relationship in which both species density and abundance declined with decreasing temperature.

#### *2.4. Does Thermal Tolerance Predict Elevational Range Size and Species Density as Would Be the Case If Temperature Were the Sole Driver of Species Occurrence?*

To test whether physiological tolerance of environmental temperature influences spatial variation in species density, the relationship between the thermal ranges of species (i.e., CTmax–CTmin) and the environmental conditions across the gradient were examined. We first asked whether species with broader thermal tolerances had broader elevational ranges and higher elevational midpoints. For each species, we combined the sampling data and plotted the highest elevation at which it was collected minus the lowest elevation at which it was collected and determined the elevational range of each species. To calculate elevational midpoints, we calculated the mean of the highest elevation and lowest elevation at which each species was collected [55]. These values were then compared to the thermal range of each species. We predicted that if temperature were an important determinant of the range sizes of species, then species at higher elevations that are able to withstand colder temperatures would have broader thermal ranges than species at lower elevations that may be confined by their physiological temperature tolerances. Additionally, species with broader thermal ranges typically have broader geographic ranges and thereby also have higher elevational midpoints [56].

Many species likely overlap in the range of temperatures at which they can occur based on their physiological thermal ranges. Yet, if species do not occupy the same environmental conditions as would be predicted by their thermal tolerances alone, then some other factor accounts for at least some of the variation in species density and occurrence. To determine whether thermal tolerance influenced species density, we asked whether the species occurring in a particular community were simply the collection of those species whose thermal tolerances overlapped the annual range of temperatures of that particular place. To do this, we extracted the annual range of temperatures (maximum temperature of the warmest month–minimum temperature of the coldest month) for each of the 27 of the 29 sites for which we had estimates of species density and calculated the mean maximum and minimum thermal tolerances of each of the 18 species across the 27 sites for which we had thermal tolerance data (two sites were omitted because species found at these sites did not have thermal tolerance data and therefore we could not estimate expected densities). The extent of overlap between physiological ranges of the ants and environmental temperatures of the sites (henceforth, thermal overlap) was then calculated. For any given species × site combination, this is simply the range of shared temperatures for both the species and the site. These values were then used to estimate a probability of occurrence for each species at each site using logistic regression models.

In the logistic regression models, the probability of occurrence of one species was determined based on its thermal overlap, as well as the thermal overlaps and recorded presences of the other 17 species in the regional species pool. This approach allowed us to determine a probability of occurrence for each species at each site based on overlapping physiological and environmental conditions (thereby incorporating variation in physiological thermal ranges and environmental thermal ranges between sites), as well as actual occurrences of other species (thereby incorporating the possibility for species co-occurrences). So as not to bias the models, presence data for the focal species were not included when estimating the probability of occurrence of that species, as including the actual occurrence of a species would inherently increase its probability of occurring in a given area.

Finally, to estimate expected species density based on thermal overlap alone, the independent probabilities of occurrence for each species were simply summed at each site. This expected species density would then be the number of species that could occur at a particular place along the gradient if temperature and temperature alone limited community membership. We then plotted observed species density against the expected species density based on thermal limits alone. If the slope of the line of expected species density plotted against observed species density equals 1, then temperature would be the sole predictor of species density. Both presences and absences of species are evident in the site × species matrix. One possibility is that the absences were not true absences. So, as a test whether the potential pseudo-absences in the site × species matrix could influence the result, we filled in the matrix so that all sites between the highest and lowest elevation at which a species was recorded were counted as presences. We then compared the expected species density if each species occurred at every site within its range to the predicted range based on thermal tolerance alone.

#### *2.5. Do Species Co-Occur Less among Assemblages than Would Be Expected If Temperature Alone Limits Membership?*

Null model analyses were used to ask whether species co-occur non-randomly (some species pair combinations being less frequent than expected by chance alone) among sites, as would be predicted if competitive interactions influenced the distribution of ants. In particular, the C-score of Stone and Roberts (1990) was used to quantify co-occurrence patterns (for methods on C-score analysis, see Appendix A) [57]. C-scores that are not significantly larger than expected by chance indicate random species distributions among sites (i.e., segregation), and C-scores that are smaller than expected by chance indicate species non-random patterns of occurrence (i.e., aggregation).

This analysis was conducted for all 29 sites to determine a general pattern, and then for the 12 communities at high (>1000 m) and 17 communities at low (<1000 m) elevations separately to examine whether the signature of competition varied along the environmental gradient.

#### **3. Results**

#### *3.1. Are Abundance and Species Density Related to Environmental Temperature?*

Species density (the number of species per site) ranged from 1 to 22 (mean = 9.44), and abundance (the total number of occurrences) per site ranged from 2 to 140 (mean = 49.5). Abundance (*r*<sup>2</sup> = 0.34, *p* < 0.001; Figure 1a) and species density (*r*<sup>2</sup> = 0.47, *p* < 0.0001; Figure 1b) both declined as mean annual temperature (MAT) declined. Species had a high degree of overlap at the lower elevations compared to higher elevations in the 2004–2007 sampling (Table 1).

**Figure 1.** The relationship between (**a**) total abundance and temperature and (**b**) species density and temperature. Temperatures are current (1950–2000) mean temperatures for each site extrapolated from the WorldClim database [54] at a resolution of 30 arc-seconds. The line in each figure is the best-fit linear regression.

**Table 1.** Species occurrence matrix for the 29 elevations sampled during 2004–2007. Filled boxes indicate species were observed at the given elevation (i.e., presence)and empty boxes indicate species were not observed at that elevation [58].


#### *3.2. Does Thermal Tolerance Predict Range Size and Species Density, as Would Be the Case If Temperature Were the Sole Driver of Species Occurrence?*

We first asked whether species with broader thermal ranges had broader elevational ranges and higher elevational midpoints, as would be predicted if temperature were an important factor determining range sizes of species. There was a positive relationship between elevational ranges and thermal ranges (*r*<sup>2</sup> = 0.48, *p* = 0.001, Figure 2a) as well as a positive relationship between elevational midpoints and thermal ranges (*r*<sup>2</sup> = 0.38, *p* = 0.004, Figure 2b). Species with broader thermal ranges occurred at more elevations and tended to have higher elevational midpoints. We stress that these were lab-measured thermal tolerances and not simply the temperatures of sites at which species were collected. While not much data exist on critical thermal minima for these species, our critical thermal maxima data are similar to those found in other North American ant species [29,59].

**Figure 2.** Thermal ranges (CTmax–CTmin) show a positive relationship with (**a**) elevational ranges (*r*<sup>2</sup> = 0.48, *p* = 0.001) and (**b**) elevational midpoints (*r*<sup>2</sup> = 0.38, *p* = 0.004) of 20 species for which we obtained both physiological and distributional data. Elevational ranges were calculated as the highest elevation at which a species was recorded minus the lowest elevation at which a species was recorded.

To determine if physiological thermal tolerances alone could predict species density, we examined the relationship between the environmental conditions at each site and the composite thermal ranges of species found at that site. In comparing sites, the thermal limits of species within sites declined with mean annual temperature (CTmax = *r*<sup>2</sup> = 0.66, *p* < 0.0001, Figure 3a; CTmin = *r*<sup>2</sup> = 0.67, *p* < 0.0001, Figure 3b). So, all species occurring at the warmest sites had, on average, the highest CTmax values (Figure 3a) and all species occurring at the coldest sites had, on average, the lowest CTmin values (Figure 3b). Furthermore, CTmax (Figure 3c) and CTmin (Figure 3d) values within and among species also follow this pattern, exhibiting both intra- and interspecific variation along the elevational gradient.

It is common to interpolate the sites at which a species could occur based on its upper and lower elevations. Here, we do something similar; we interpolate the sites at which a species could occur based on its thermal tolerance. We assume a species can occur at all sites along the gradient within its physiological thermal range (where MAT of the site is higher than its CTmin but lower than its CTmax). When we performed this interpolation, we found that at lower temperatures, observed richness more closely matched expected richness based on thermal constraints alone; however, at warmer temperatures, there was more deviation in observed richness from the null expectation (Figure 4), indicating that at low elevations, temperature is not the sole driver of species density.

#### *3.3. Do Species Co-Occur Less among Assemblages than Would Be Expected If Temperature Alone Limits Membership?*

When all 29 sites along the gradient were considered, species co-occurred much less than expected by chance (i.e., they were strongly segregated; observed C-score = 12.66; simulated C-score = 10.92; SES = 7.81; *p* < 0.0001). However, when co-occurrence patterns were examined separately, species in the warm, low-elevation sites (<1000 m) were significantly segregated among assemblages (SES = 2.06; *p* = 0.02), but species in cold, high-elevation sites (>1000 m) showed no significant deviation from randomness with respect to one another (SES = 0.94; *p* = 0.17).

**Figure 3.** Species (**a**,**b**) and populations (**c**,**d**) that occur at the warmest sites have, on average, the highest CTmax values (**a**,**c**), and those that occur at the coldest sites have, on average, the lowest CTmin values (**b**,**d**). Each point is the mean of the thermal limits for all species pooled at each site (**a**,**b**) and for each population of a species (**c**,**d**) averaged for each site. Each color in (**c**,**d**) corresponds to different ant species collected across multiple populations along the elevational gradient, and each point within the color corresponds to the pooled individual CTmax (**c**) and CTmin (**d**) for an ant species at each site. Temperatures were extrapolated from WorldClim [54] at a resolution of 30 arc-seconds and represent mean temperatures from 1950–2000.

**Figure 4.** Residuals of the observed and expected species density based on thermal overlap in physiological limits and environmental temperatures. Black line indicates the 1:1 line to illustrate deviation from the null expectation.

#### **4. Discussion**

One of the most striking patterns in nature is that species density, abundance, and diversity vary, often systematically, along environmental gradients. Importantly, the influence of biotic interactions relative to abiotic factors shifts with elevation and environmental conditions. Such a finding supports the notion that the processes shaping community structure are context-dependent, and that both biotic and abiotic factors interact to determine the distribution and density of species among assemblages.

Temperature is related to total abundance and species density, and temperature (especially cold temperature) likely limits the ranges of species as well [8,34,38]. Species with broad elevational ranges also have broad thermal ranges. However, we need to elucidate why temperature matters. In this case, we can rule out some temperature-dependent mechanisms as an influence on patterns of diversity. One such influence that we can rule out is the metabolic theory of ecology, which depends on temperature-dependent activation energies [60], as previous work with this system [9] and others [61,62] have demonstrated. Similarly, the pattern of ant diversity here probably does not arise because of variation in in situ temperature-dependent speciation rates, since none of the studied species are endemic to the study region. In addition, temperature and net primary productivity (NPP) are not correlated spatially in this system, and NPP is weakly and negatively, rather than positively, related to ant diversity [9]. Finally, although it has been suggested that one effect of temperature on ant diversity is via the effects of temperature on ant foraging and access to resources, experimental manipulation in this same ant study system found that changes in temperature did not limit access to resources by ants [63].

One possibility is that low temperatures limit abundance, and might alter community evenness and species-abundance distributions and in turn, species density at high-elevation sites [64], but see [65]. Specifically, underlying physiological constraints exert a filter on community membership by allowing only certain cold-tolerant taxa to establish and persist at high elevations, as low temperatures limit both overwintering success as well as slow the rate of brood development. Previous work in this system [19] found that assemblages at high-elevation sites are characterized by the presence of fewer and clustered lineages, as might be expected if only the species of the restricted subset lineages with the ability to tolerate cold climates persist at high elevations. Our measurements

of physiological tolerance lend support to the conclusions from previous community phylogenetic approaches. On average, populations at high-elevation sites tended to have lower CTmin values than did populations at warmer low-elevation sites. In fact, thermal breadth also increased with elevation, suggesting that communities at higher elevations consist of individuals that can withstand a wider range of environmental temperatures than low-elevation species.

We found that at warmer sites, observed species density varied more from the densities predicted by physiological–environmental matching alone, but observed densities in colder and more stressful conditions approximated expected densities. Populations at higher elevations (and latitudes) often persist in areas that are colder than would be expected based on their CTmin values [55] and low-elevation (and low-latitude) populations can persist in regions that are much warmer than their CTmax would suggest. Here, we found that populations in high elevation assemblages occur at temperatures that are colder than would be expected given their thermal tolerances; this has been referred to as "overfilling" the thermal niche space [66]. In contrast, populations at low elevations do not occur in all of the places they might, based on their thermal tolerances alone, which has been dubbed "underfilling" of thermal niche space. It has been suggested that overfilling of niche space is due to winter survival mechanisms of physiological cold tolerance and behavioral avoidance strategies (e.g., diapause) [29]. This seems likely in our case, though we did not specifically examine any potential overwintering mechanisms.

Underfilling of thermal niche space in warm sites might be due to biotic interactions, as has often been suggested in the literature [66]. Here, our null model approach lends support to the idea that interspecific interactions, especially in warm sites, limit community membership. The idea that interactions structure ant assemblages is not new [25,67,68]. In fact, the strongest evidence for the effects of competition on ant assemblages comes from the collapse of many ant assemblages in the face of competitively dominant invasive species [69], non-random patterns of co-occurrence among assemblages [70,71], temporal, spatial, and resource partitioning within assemblages [24,72,73] and the influence of competitively dominant native species [25,67,68]. Here, we focused on the co-occurrence of native species, and use observed interactions as well as the C-score of Stone and Roberts [57] and null model analyses to show evidence of the role of biotic interactions within communities. When we examined the warmer sites (below <1000 m elevation) and the colder sites (above >1000 m in elevation) separately, we found evidence for the signature of interspecific interaction in low-elevation sites, but not high-elevation sites. That is, species co-occurred less than expected at low-elevation sites, as would be predicted if competitive exclusion structured communities [74] at high-elevation sites, species co-occurred randomly with respect to one another. These null models alone do not directly implicate interactions, but they agree with other independent lines of evidence. First, community phylogenetic evidence points to the role of interspecific competition in shaping low-elevation but not high-elevation sites [19]. Second, there is more deviation from the null expectation in low-elevation assemblages based on physiology-environment associations alone.

#### **5. Conclusions**

One of the fundamental tenets of biogeography is that abiotic and biotic factors interact to shape the distributions of species and the organization of communities, with interactions being more important in benign environments, and environmental filtering more important in physiologically stressful environments [75,76]. Null models and community phylogenetic studies at large spatial scales, and manipulative experiments at small spatial scales, have hinted at such a scenario [18]. Taken together, our results, using a combination of observational data, null models, and physiological measurements, provide a strong test that interspecific interactions drive the distributions and density of species in warm climates, but that physiologically-driven environmental filtering predominates at high-elevation sites, at least in temperate forests of the US.

Our results also have implications for predicting the responses of biodiversity to ongoing climate change, a topic at the forefront of ecology and diversity science [77]. Forecasts of biodiversity change in response to climate change, often rely on matching the thermal tolerances of species to thermal environments in the future, and most show that species in the warmest places are the most susceptible to ongoing warming because species are operating close to their thermal maxima, so any increase in temperature essentially pushes these species over the thermal edge [66,78]. While some studies have pointed out that organisms in the warmest places can modulate their behavior to escape stressfully high temperatures, they have generally overlooked the fact that these warm places are also where organisms are likely to face the most negative consequences of interspecific interactions. For instance, our thermal constraints models showed that diversity varied most from our expectation in the warmest places, because that is also where biotic interactions among species are the most important in limiting community membership. So, while positive interactions among species might buffer species in the face of climate change, negative interactions such as competition might exacerbate the effects of climate change on biodiversity in warm environments. Models that focus on the future of biodiversity in warm environments, where most of biodiversity is, should also examine the combined and relative effects of biotic interactions and abiotic constraints and how these processes scale up to influence patterns of diversity.

**Author Contributions:** Conceptualization, L.D.C., J.-P.L., R.R.D., and N.J.S.; methodology, L.D.C., J.-P.L., N.J.S.; validation, L.D.C. and N.J.S.; formal analysis, L.D.C., J.-P.L., R.R.D., and N.J.S.; investigation, L.D.C. and J.-P.L.; resources, L.D.C., J.-P.L. and N.J.S.; data curation, L.D.C. and N.J.S.; writing—original draft preparation, L.D.C., R.R.D. and N.J.S.; writing—review and editing, L.D.C., J.-P.L., R.R.D. and N.J.S.; visualization, L.D.C. and N.J.S.; supervision, R.R.D. and N.J.S.; project administration, L.D.C., J.-P.L., and N.J.S.; funding acquisition, N.J.S. and R.R.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the U.S. Department of Energy PER award (DEFG02-08ER64510) and a National Science Foundation Dimensions of Biodiversity grant (NSF-1136703). Additional support was granted to J.-P.L. and L.D.C. from the Department of Ecology & Evolutionary Biology at the University of Tennessee.

**Acknowledgments:** We are grateful to Quentin Read, Nicholas Gotelli, Carsten Rahbek, and Michael Borregaard for their comments and computational expertise on the null models. We thank Kathryn Stuble, Courtney Patterson, and David Fowler for field assistance. We also thank anonymous reviewers for comments on previous versions of this manuscript. We acknowledge the National Park Service for collection permits for Great Smoky Mountains National Park (GRSM-2011-SCI-0050, GRSM-2013-SCI-1071).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

The C-score quantifies the number of "checkerboard units" for each species pair, where a checkerboard unit is a 2 × 2 submatrix of the form 01 10 or 10 01. For each species pair, the number of checkerboard units is (Ri-S)(Rj-S), where Ri is the number of occurrences (equal to the row total) for species i, Rj is the number of occurrences for species j, and S is the number of sample plots in which both species occur. The C-score is the average number of checkerboard units for each unique species pair. If this index is unusually large compared with a null distribution, there is less pairwise species co-occurrence (segregation) than expected by chance. If the index is unusually small, there is more species co-occurrence (aggregation) than expected. We compared the observed C-scores to those generated from 5000 randomly constructed assemblages (using null models in EcoSim version 7.72) [79].

A fixed-fixed null model [80] was used for which both row totals and column totals are fixed within sites and among species, which maintains differences in species density among sites and total occurrences among species. Gotelli [80] suggests that SIM9 is appropriate for analyzing co-occurrence patterns of species from "island lists" and has a low probability of Type I errors.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
