Testing Huston's dynamic equilibrium model along fire and forest productivity gradients using avian monitoring data

Many studies investigating the response of wildlife to disturbance focus predominantly on the effects of the disturbance alone but fail to account for the influence of ecosystem productivity in moderating the response of species and thus the resulting biodiversity patterns. We use Huston's dynamic equilibrium model (DEM) to explore the relationship between avian diversity and fire across the greater Rocky Mountain region of the western United States. This model provides the theoretical foundation to understand the distinct and interactive effects disturbance and productivity have on regulating species richness.


| INTRODUC TI ON
Disturbances are mechanisms widely recognized for regulating the composition and coexistence of species within a community (Sousa, 1984). Though many definitions exist (Pickett & White, 1985), we define disturbances as discrete events in time that kill or damage organisms and alter the availability of resources (sensu Mackey & Currie, 2000). Disturbances occur in a variety of forms such as wildfire and extreme weather events, as well as anthropogenic land modification (Dornelas, 2010;Sousa, 1984).
Depending on the ecosystem and temporal and spatial scales that define them, disturbances have different effects on biodiversity. Single disturbance events, like wildfire, can reduce biodiversity through direct mortality and the extirpation of organisms. In contrast, biodiversity can increase when disturbance events occur in a spatially and temporally heterogeneous manner, altering the amount and types of resources available (Holt, 2008;Pickett & White, 1985).
Periodic disturbances can also increase biodiversity by promoting the coexistence of species adapted to different environmental conditions and resource levels, thereby reducing the effects of interspecific competition (Connell, 1978;Petraitis et al., 1989).
One of the most important disturbance processes influencing resource availability and habitat heterogeneity within forested ecosystems is wildfire (Franklin et al., 2002;Millar et al., 2007). The US Forest Service defines wildfire as an unplanned fire in a natural area that has not been significantly modified by human development (development does not refer to forest management, Stein et al., 2013).
Wildfire can alter the structural heterogeneity of forests by converting living trees to standing dead or downed trees while simultaneously consuming varying amounts of understory vegetation depending on frequency and intensity (Franklin et al., 2002;Kane et al., 2013). In North America, forest ecosystems will be subject to changes in the severity and frequency of disturbance regimes, specifically fire and drought, over the next 50-150 years (IPCC, 2014).
In particular, contemporary increases in seasonal warming and spring snowmelt (1970s-present) have already resulted in elevated wildfire frequencies, increases in burned area and extended wildfire seasons (Jolly et al., 2015;Westerling, 2006). These projected increases necessitate the need for continued examination of how fire structures biodiversity to inform fire management as well as future conservation planning.

| Disturbance ecology theory
A predominant community ecology theory explaining the relationship between biodiversity and disturbance is the intermediate disturbance hypothesis (IDH) proposed by Connell (1978) (but see also Grime, 1973;Horn, 1975). The IDH predicts species richness to be greatest at intermediate levels of disturbance measured by the frequency or intensity of disturbance events at the scale of species interactions (Connell, 1978;Petraitis et al., 1989). The expectation for species richness is to be low in both infrequently and frequently disturbed systems through the dominance of competitively superior or rapidly colonizing species. Levels of intermediate disturbance promote coexistence and therefore greater species richness. However, only a few empirical examples of the Gaussian relationship predicted by the IDH exist (Mackey & Currie, 2001;Shea et al., 2004). For example, Shea et al. (2004) showed 18% of 250 studies supported the IDH, while Mackey and Currie (2001) indicated only 16% of 116 richness-disturbance relationships exhibited maximum species richness at intermediate levels of disturbance. These findings have promoted discussions to abandon the IDH (Fox, 2013). However, Huston (2014) suggests that attempts fail to find patterns predicted by the IDH because they do not account for variation in ecosystem productivity. Huston proposed the dynamic equilibrium model (DEM) that accounts for variation in productivity across sampling locations (Huston, 1979(Huston, , 1994. The DEM suggests that disturbance interacts with the productivity of an ecosystem to affect species richness. Furthermore, in ecosystems where a spatial gradient in productivity occurs across study sites, stratification along this gradient allows for the true relationship between species richness and disturbance to be realized (Huston, 2014).
Productivity is often used to describe the growth rate of plants via biomass accumulation and represents a measure of energy available to plants and the broader community of primary and secondary consumers (Putman & Wratten, 1984). Though initially used to explain patterns in plant species richness, productivity has also been used extensively in studies on the distribution of animals, typically by remotely sensed data such as the Normalized Difference Vegetation Index (NDVI, Crego et al., 2020;Pettorelli et al., 2011) or net/gross primary productivity (NPP/GPP, Cusens et al., 2012;Gebert et al., 2019;Rosenzweig, 1995).
The foundation of Huston's DEM is that competitive exclusion (Gause, 1932;Grinnell, 1904) will occur at both high and low levels of productivity. Strong or frequent disturbances are required in highly productive systems to counteract the effects of competitive dominance. Systems with low productivity have reduced diversity, which can be compounded by strong disturbances, but a weak disturbance is sufficient to counteract competitive exclusion at low productiv-

| Testing Huston's DEM using avian monitoring data
Here, we test Huston's DEM to explore the relationship between bird species richness and fire throughout the greater Rocky inference to changes in bird richness at forests categorized as low or high productivity only. We argue that focusing on the extremes of either disturbance or productivity can mask the relevant dynamics of an ecological system. Therefore, examining the relationship between bird richness and disturbance across a continuum of productivity and disturbance severity may provide a more robust test of the DEM.
Many studies have attempted to identify response patterns of birds to fire and are discussed in a variety of qualitative reviews (e.g. Kotliar et al., 2002;Saab & Powell, 2005). These reviews provide evidence that birds exhibit consistent interspecific differences in their responses to wildfire (Fontaine & Kennedy, 2012;Kotliar et al., 2002;Saab & Powell, 2005). While reviews aid in the interpretation of findings from a compilation of independent ecological studies, their results are subject to multiple sources of variation that can be unique to each set of independent data (e.g. study design or sampling protocol; Gurevitch & Hedges, 1999). For this reason, past reviews have urged for additional evaluations of bird responses to fire (Fontaine & Kennedy, 2012;Kotliar et al., 2002;Saab & Powell, 2005). Data from a single source with a consistent sampling design, such as a largescale monitoring programme, provide a more rigorous evaluation of the relationship between bird species richness and fire.

| Hypotheses and predictions
Following Huston's DEM, we hypothesized bird richness is at its maximum at intermediate levels of fire severity and productivity.
However, productivity tends to increase with latitude across the interior mountain west of the United States (Mekonnen et al., 2016) as forest composition changes from drier forests dominated by pinyonjuniper (Pinus spp. and Juniperus spp.) in Utah and western Colorado to mixed coniferous forests dominated by Douglas fir (Pseudotsuga menziesii) and lodgepole pine (Pinus contorta) as one moves northward (Hicke et al., 2007). Therefore, the relationship between richness, fire severity and productivity may change across the spatial gradient of productivity. We predict the value of productivity at which species richness is maximized increases with increasing latitude ( Figure 1a,b). Furthermore, bird richness generally increases with ecosystem productivity (e.g. Hurlbert & Haskell, 2003) but across larger spatial gradients richness can decline at very high levels of productivity resulting in a quadratic relationship (Phillips et al., 2008(Phillips et al., , 2010Verschuyl et al., 2008). There are also regionally more forest avian species in higher latitudes (Goetz et al., 2014). Accordingly, we predict species richness increases with increasing latitude, and the Gaussian species richness relationship predicted by the DEM is most peaked at higher latitudes where productivity is greatest (Figure 1c,d).

| ME THODS
To test for changes in bird species richness related to gradients in productivity, fire severity and their interaction, as well as with latitude, we employed a hierarchical multispecies occupancy modelling approach  Huston, 2014). (b) The hypothesized shift in the location of maximum richness along the productivity gradient and increase in the strength of the peaked relationship at higher latitudes where species pool is larger and average productivity is higher. Dashed lines represent the level of productivity where species richness is highest and correspond to the cross-sections (c and d) below each graph enable inference to the entire community as well as individual species by separating the ecological processes (true presence or absence of a species related to some variables of interest) from the observation process such as species detectability (Iknayan et al., 2014). To provide unbiased estimates of species richness, we also included data augmentation to incorporate the effects of rare or difficult to observe species that were present but never observed during any surveys (Royle et al., 2007;Zipkin et al., 2010). Additionally, to test for the predicted stronger, more peaked relationship between species richness and fire severity at high latitudes, we compared the variance in richness values across the fire severity gradient at low and high latitudes. We used a Bartlett test to determine whether there was a statistical difference between the variances in species richness (and consequently the shape of the richness-fire severity relationships) for low and high latitudes ("bartlett.test" R base function, Bartlett, 1937).

| Avian survey data
We used avian monitoring data collected by the Integrated We used presence/absence data from the IMBCR surveys to model species occupancy from which species richness was derived (following Dreitz et al., 2017).

| Site-level fire and productivity
We used fire data from the Monitoring Trends in Burn Severity (MTBS) program administered by the US Forest Service (MTBS, 2018). MTBS is a multi-year data set beginning in 1984 and provides 30-m resolution data of burn severity for all fires ≥4 km 2 (1,000 acres) that occurred in the United States. For each fire, we classified pixels into one of four burn severity classes: unburned (0), low (1), mixed (2) and high severity burn (3). Using severity data from the most recent fire, we calculated the average burn severity score within a 1-km radius buffer around the centroid of each study site. We also calculated the elapsed time, in decimal years, between the initiation date of each fire and the date of the annual bird survey at each study site (hereafter time since fire).
We extracted the survey period mean productivity using gross primary productivity (GPP) at each study site using the same 1-km radius buffer around each site centroid used to measure fire severity.
Within each 1-km radius buffer, we calculated the average 16-day GPP (in kg C m −2 ) from May to August to encompass the timeframe of the avian surveys and the breeding season of sampled birds. GPP data were a product created from Landsat 30-m 16-day GPP with improved handling of the spatio-temporal variability that accompanies land-based productivity data (Robinson et al., 2018). We used GPP as a measure of site-level productivity because it is a better predictor of bird species richness, especially in areas of dense vegetation such as forests, compared to other productivity measures like NDVI and NPP (Phillips et al., 2008(Phillips et al., , 2010).

| Multispecies occupancy model
We modelled the true occurrence within the avian community (w i ) of all species i observed and unobserved as a Bernoulli random variable with some probability Ω, w i~ Bern(Ω), also referred to as the superpopulation process (Iknayan et al., 2014). Variable w i = 1 if species i was observed or unobserved but available for sampling. When w i = 0 , the occurrence of species i is also zero, indicating a species unavailable for sampling. True occupancy for species i at site j in year k, denoted z i,j,k , was modelled as a Bernoulli random variable, z i,jk~ and zero otherwise. Actual observations from point count surveys, y i,j,k,l , denote detection or non-detection of species i at site j in year k F I G U R E 2 Map of study sites located across 5 western states of the United States. Species richness categories are indicated by point colour.
during point survey l. Detections were modelled as a Bernoulli random variable, y i,j,k,l~ Bern(p i,j,k,l , z i,j,k ), where p i,j,k,l represents the detection probability for species i at site j in year k during point survey l.
We modelled site-level species occupancy probabilities, Ψ i,j,k , incorporating eight site-specific covariates using the logit-link function: The species-specific intercept ( 0 i ) represents the mean occupancy for species i. Parameters 1 i and 2 i are linear effects of GPP at site j in year k and fire severity from the most recent fire at site j, respectively. To account for the nonlinear relationship between species richness and both productivity and disturbance severity suggested by Huston's DEM (Huston, 1979(Huston, , 1994, we included quadratic forms of GPP and fire severity, parameters 3 i and 4 i , respectively. Huston's DEM also hypothesizes a dependency between productivity and disturbance severity such that the response of a species to either of those variables depends on the value of the other variable (Huston, 1979(Huston, , 1994. Thus, we included the interaction, 5 i , between GPP and fire severity. Parameter 6 i is the linear effect of time, defined in decimal years since the most recent fire at site j began and when site j was surveyed. We also accounted for the known variation in occurrence among species across both geographic and elevational gradients by incorporating the effects of latitude and elevation into the model on species occupancy, parameters 7 i and 8 i , respectively. The inclusion of latitude also serves as a method to account for the hypothesized north-south spatial trend in productivity and make spatially stratified predictions. All explanatory site-level variables were standardized to have a mean of zero. Finally, j is a site-level random effect, included to account for the influence of additional site-based differences such as the presence of other disturbances or other biophysical features not measured in this study (e.g. bark beetle tree die-off).
Survey-level species detection probabilities, p i,j,k,l , were modelled as: The intercept, b0 i , represents the mean detection probability for species i across all surveys. The parameters b1 i and b2 i are linear and quadratic effects of the Julian day of year k during which survey l was conducted at site j. These effects were included to account for the temporal variation in the probability of encountering species with differing breeding phenology.
In multispecies models, single-species occupancy and detection rates are linked at the community level by treating species-specific parameters as random effects with prior distributions (Royle & Dorazio, 2008). The community response is described by the hyperparameters of those prior distributions (Royle & Dorazio, 2008). For example, we modelled 1 i as a normal random variable, 1 i ~ Normal( 1 − 2 1 ), with 1 representing the community mean effect of GPP on occurrence and 2 1 as the variance of the effect among species within the community. Hyperparameters were modelled with vague prior distributions for community means (e.g. 1 Normal( = 0, 2 = 1000) and variances (e.g.  with lower latitude study sites (<44 degrees latitude) less productive (mean = 0.06 kg C m −2 16-day −1 , SD = 0.03) compared to sites at higher latitudes (>44 degrees latitude, mean = 0.09 kg C m −2 16day −1 , SD = 0.02). Furthermore, we found GPP increased with increasing latitude across the entire study region supporting our initial predictions of a shift in regional GPP with latitude (Appendix S1).

| RE SULTS
Fires occurred at avian survey sites between 0.75 and 31.75 years prior (mean = 10.34 years, SD = 6.89 years); 30% of site surveys (108 of 389) occurred ≤five years after a fire, and 55% of site surveys (215 of 389) took place ≤10 years after a fire. The average fire severity score was 0.94 (SD = 0.74), and study sites ranged in severity scores from 0.05 to 2.68 on a 3-point maximum scale. Neither fire severity nor time since fire was correlated with latitude (r 2 < .01 and r 2 < .01, respectively

| Species richness
Across all sites, species richness increased nonlinearly as a function of increasing GPP, reaching an asymptote at 0.11 kg C m −2 (Figure 3a).

| D ISCUSS I ON
We found support for Huston's DEM based on patterns in bird species richness across combinations of fire severity and productivity gradients (Figure 4). We also found a spatial trend in productivity, with lower latitudes having lower overall productivity compared to higher latitudes within the study region (Appendix S1). As predicted, the relationship between richness, fire severity and productivity varies along the north-south productivity gradient.
Lower latitudes are less productive, and bird species richness is greatest at intermediate levels of both fire severity and productivity, as suggested by DEM with the predicted Gaussian relationship.
Furthermore, concurrent with our hypothesis, the level of productivity where species richness is maximized shifts upward along the productivity spectrum as latitude and subsequently regional F I G U R E 4 Predicted mean species richness (N) as a function of fire severity and GPP at low (a) and high (b) latitudes.
Warmer colours indicate greater species richness. Lines represent a change of one species in the richness gradient productivity increase (Figure 4, Appendix S1). For example, maximum species richness at 38 degrees latitude occurs at 0.08 kg C m −2 16-day −1 , while maximum richness at 48 degrees latitude occurs at 0.12 kg C m −2 16-day −1 , a 50% increase. Across latitudes, greatest richness occurs at a fire severity rating of 1.4-1.5, which is consistent with a mixed severity burn.
We also predicted that at higher latitudes, the maximum peak richness would be more pronounced and that richness declines at a higher rate across the fire severity and productivity gradients ( Figure 1b). While higher latitude sites generally have more species and richness was greatest at intermediate fire severity, the rate of gain and loss of species richness across the fire severity spectrum was not statistically different than that at lower latitudes. Thus, this relationship was less peaked than predicted. We also did not find evidence of a decline in richness at higher productivity levels ( Figure 1b). Rather, our findings suggest that when regional productivity is higher (such as at higher latitudes in our study region, Appendix S1), species richness asymptotes as productivity increases.
This contrasts with the expected result of very high levels of productivity eliciting a reduction in richness as predicted by the DEM.
An explanation for the lack of a strong peak in species richness at higher latitudes is not straightforward. However, one foundation of the DEM is that interspecific competition is an important regulator of community composition (Huston, 1979(Huston, , 1994(Huston, , 2014. While we did not specifically test for interspecific competition, interspecific competition in highly productive areas may be reduced when F I G U R E 5 Species plotted as a function of their mean response to latitude (x-axis) and GPP (y-axis). Species (points) are coloured by mean occupancy throughout the study region; red species are more common than black species communities contain many common species (Connell et al., 1984).
Common species encounter more conspecifics than rare species increasing intraspecific competition as a control on common species populations while allowing the occurrence of rare species to increase, thus increasing species richness (Connell et al., 1984). We found avian communities at high latitudes contained more common species and higher species richness ( Figure 5). Conversely, proportionally more rare species were present in communities at lower latitudes ( Figure 5). Furthermore, lower latitude sites were less productive than at higher latitudes, and a larger proportion of species increased in occurrence with increasing latitude and productivity ( Figure 5). These findings support reductions in interspecific competition when productivity is regionally high and provide an explanation for the lack of support for Huston's DEM at higher latitudes in our study.
Alternatively, there may be forests with higher productivity outside of the study area. We have attempted to extend the DEM to wildlife species that are highly mobile and interact with the environment at multiple spatial scales. While our study area was across a multi-state region based on the sampling extent of the monitoring programme we chose, future tests of Huston's DEM on taxa like birds may provide more insight at either broader continental-level assessments or finer scales such as individual forest stands. For example, we hypothesize stronger, more definitive species richness relationships would occur at continental scales due to the greater range of productivity across those scales. Likewise, the importance of temporal scale should also be considered complimenting the spatial scales across which ecological processes are measured (Reynolds-Hogland & Mitchell, 2007). Our study was limited to the breeding months, and analyses of bird richness-fire relationships during the nonbreeding season when productivity is lowest may provide additional insights into the structuring of forest bird communities .

| Mechanisms and management implications
Our results suggest mixed severity burns maximize species richness, importantly only after accounting for the interaction of fire severity and ecosystem productivity. Previous research has well documented the vast diversity of birds that use forests and how they respond to fire in the inter-mountain west of North America (Fontaine & Kennedy, 2012;Kotliar et al., 2002;Saab & Powell, 2005). This variation across forest bird species is embodied by two primary life history components: nesting and foraging requirements. For example, the black-backed woodpecker (Picoides arcticus) nests in cavities of dead trees and forages for insects on the bark of recently burned trees (Tremblay et al., 2016). In contrast, the western tanager (Piranga ludoviciana) nests in intact tree canopy and forages by gleaning insects off of live foliage not typically present following burns of medium to high severity (Hudon, 1999). Mixed severity burns are comprised of both low and high levels of disturbance creating the potential for greater resource diversity thereby supporting a broader range of species.
Forest managers are faced with maintaining multiple resources, including timber and biodiversity, under new disturbance regimes. By 2060, warmer and drier climatic conditions are expected to increase seasonal fire severity ratings by 110%-300% (Flannigan et al., 2013) and total area burned has the potential to increase by 25%-170% (Dale et al., 2001;Yue et al., 2013). Contemporary increases in fire season severity and an expectation for continued increases due to climate change have prompted land managers to move towards more active forest management (Stephens & Ruth, 2005). Strictly through the lens of maximizing forest bird richness, managing for mixed severity fire may improve species richness; however, doing so across large spatial scales (e.g. >1,000,000 km 2 in this study) would be challenging.
One limitation of this study and a challenge for any study attempting to test ecological theory at broad spatial scales using monitoring data is the lack of fine-scale vegetation information, specifically measures of vegetation structure and classifications of the plant community. We incorporated the time since disturbance in our modelling to account for the successional changes which occur following fire in forest ecosystems. This study was conducted over a broad spatial scale, and our results suggest that at the spatial scale of our study, time since fire is a rudimentary measure and one that explained a small amount of the variation in species richness

| CON CLUS IONS
Our results provide an example of how relationships between species richness and disturbance may be missed if viewed outside the context of other environmental variables and biophysical characteristics. When examining the relationship between bird species richness and changes in fire severity alone, richness changed little across the fire severity gradient (Figure 3.2b). However, we found species richness can be highly variable when viewing species responses to fire severity in the context of changing productivity (Figure 4). These results demonstrate that studies across broader spatial scales are still subject to the same sources of variation that cause disagreement between studies conducted at local spatial scales. While challenging, accounting for the ecological context a given community of organisms experience is integral to accurately assess the response of those organisms to changing conditions (e.g. via disturbance).