Small mammal community composition varies among Ozark glades

Abstract

Island biogeography theory (IBT) explains and estimates large-scale ecological patterns among islands and isolated habitat patches. Specifically, IBT predicts that the number of species per habitat patch differs as a function of area and isolation as a result of local colonization and extinction. Accurate estimates of species richness are essential for testing predictions of IBT, but differences in detectability of species can lead to bias in empirical data. Hierarchical community models correct for imperfect detection by leveraging information from across the community to estimate species-specific occupancy and detection probabilities. Using the fragmented Ozark glades as our model system, we constructed a hierarchical community model to 1) estimate site-level and regional species richness of small mammals while correcting for detection error, and 2) determine environmental covariates driving occupancy. We sampled 16 glades in southwestern Missouri in summer 2016–2017 and quantified mammal community structure within the glade network. The detected species pool included eight species, and the model yielded a regional species estimate of 8.6 species, with a mean of 3.47 species per glade. Species richness increased with patch area but not isolation, and effects of patch shape varied between species in the community.

community structure, hierarchical occupancy model, Island biogeography theory, Ozark glades, small mammals

Island biogeography theory (IBT) predicts that species richness on an island increases with area and decreases with isolation due to higher rates of extinction on small islands and lower rates of colonization on distant islands (MacArthur and Wilson 1963, 1967). IBT has been extended to a variety of habitat types in the years since it was first published, with varying results (Browne 1981; Lomolino 1984; Fox and Fox 2000; Belmaker et al. 2006; Watling and Donnelly 2006). Although the species-area relationship predicted by IBT is typically present in terrestrial systems (Lomolino 1984; Fox and Fox 2000; Cook et al. 2002; Watling and Donnelly 2006), isolation effects tend to be less consistent (Watling and Donnelly 2006; Prevedello and Vieira 2009). Further, evaluations of predictions rely on completeness of an inventory (Kodric-Brown and Brown 1993; Cameron and Pokryszko 2005; Triantis et al. 2008), which is hampered by issues of detectability (Iknayan et al. 2014). For instance, analysis of small mammal communities on mountaintops and conclusions by Brown (1971) were shown to be incorrect due to incomplete samples (Lawlor 1998).

The discrepancy in isolation effects between oceanic and habitat islands largely is due to the habitat or habitats surrounding the focal patches (Ricketts 2001; Cook et al. 2002; Ewers and Didham 2005; Evans et al. 2017). While oceanic islands are surrounded by open water, habitat islands represent discrete patches of one habitat (e.g., meadows or mountaintops) surrounding by a matrix of another type (e.g., forests). Matrix types between habitat islands vary considerably in their permeability to dispersing individuals (Ewers and Didham 2005; Kupfer et al. 2006; Prevedello and Vieira 2009; Evans et al. 2017), and many habitat islands are embedded in a heterogeneous matrix consisting of a variety of habitat types (Ricketts 2001; Kupfer et al. 2006; Biswas and Wagner 2012; Deans and Chalcraft 2016). In addition to variation in matrix permeability, isolation effects can be masked by spill-over of species from the matrix to the habitat patch of interest (Brown and Dinsmore 1988; Cook et al. 2002).

The Ozark glades, open and rocky habitats surrounded by forest, can be considered habitat islands embedded in a forest matrix. Glades are remnant habitats from the Xerotheric maximum (Mondy 1970) and are characterized by shallow soils, the absence of a developed tree canopy, and aridity relative to the surrounding forest (Kucera and Martin 1957; Nelson and Ladd 1983; Heikens 2007; Nelson 2010). Glades require periodic fires to maintain ecosystem health (Guyette and McGinnes 1982; Nelson and Ladd 1983; Ladd 1991; Baskin and Baskin 2000; Nelson et al. 2013); isolation therefore has the potential to affect recolonization after each disturbance event. Glades support a variety of small mammals including habitat generalists such as the white-footed mouse (Peromyscus leucopus) and eastern woodrat (Neotoma floridana) as well as open-canopy specialists such as the Texas mouse (P. attwateri) and fulvous harvest mouse (Reithrodontomys fulvescens). The patchy distribution of the glades and the diversity of small mammal species they contain makes them a suitable system for testing predictions of IBT.

Based on IBT, we predict that species richness will be a function of patch area; specifically, there will be a positive effect of area due to lower extinction rates in larger habitat islands. Whereas IBT predicts a negative effect of isolation on richness, we predict that the presence of several generalist species and permeability of the matrix will diminish or negate the effect of isolation on richness. Because patch shape, here defined as the amount of edge relative to area, should be correlated with habitat complexity, we predict that it will affect species-specific occupancy, but the direction and magnitude of these effects will vary due to variation in life history of the species in the community. Alternatively, in the absence of habitat specialists, there should be a positive effect of shape on occupancy because increasing relative edge increases the probability of encountering the patch. Because IBT predicts that small islands have higher rates of extinction than large islands, we predict that small islands will exhibit a higher rate of compositional change between sampling years.

Materials and Methods

Data collection

We identified glades using the Missouri Natural Glades shapefile developed by the Missouri Department of Natural Resources (Nelson et al. 2013; Nelson 2015). We limited our search to glades within the Ozark Highlands ecological section (Nigh and Schroeder 2002) that were of sufficient length to contain a 250-m linear transect. From these potential sampling sites, we chose 16 glades spanning an east–west distance of approximately 125 km within the Ozark Highlands (Fig. 1). The glades were located within four state or federal management units: Roaring River State Park (glades 0–3), Caney Mountain Conservation Area (glades 4–7), Mark Twain National Forest Ava-Cassville District (glades 8–10), and Drury-Mincy Conservation Area (glades 11–15). We chose these four land units due to 1) availability of sufficiently large glades, and 2) variation in management regimes, which could affect the small mammal communities within the glades. We scouted the selected glades in person prior to sampling to ensure they met the conditions described in the “Introduction.”

Fig. 1.—

Map of sampling sites across southwestern Missouri. Each point represents a glade and each management unit is noted with a rectangular box and label.

Small mammal communities were sampled in May–July 2016–2017. Small mammals were captured using two types of Sherman live-traps, 7.62 × 8.89 × 22.86 cm (“regular”) and 10.16 × 11.43 × 38.10 cm (“large”; H. B. Sherman Co., Tallahassee, Florida), baited with sunflower seeds or a peanut-butter oat mixture. Each bait type was used in 50% of the traps at each site to limit differences in catchability between sites. To standardize trapping effort, each glade was sampled using a ~ 250-m linear transect, with trap stations set ~ 10 m apart. Two traps were placed at each station and arranged to maximize capture efficiency (e.g., along runways such as fallen logs or rock ledges). Twelve of 25 stations were randomly selected to have one large and one regular Sherman trap; the other 13 stations had two regular Sherman traps. GPS receivers were used to record positions at the beginning, midpoint, and end of each transect. Traps were opened each evening and checked the following morning; each transect was sampled for a period of four consecutive trap days.

Upon checking traps, captured mammals were transferred to a cloth bag and handled according to American Society of Mammalogists guidelines (Sikes et al. 2016). Individuals were identified to species, weighed, sex was determined, and standard external measurements (ear length, tail length, and right hind foot length) were taken.

Patch characteristics were computed by first recording the perimeter of each glade with a GPS unit and then converting this output to a polygon shapefile using ArcGIS 10.3.1. Glade area was calculated from the polygons using functions in the R packages rgdal and rgeos (Bivand et al. 2016; Bivand and Rundel 2017). Isolation between glades was estimated by calculating pairwise distances between glade centroids and the shortest pairwise distances between glade edges in R using functions in the rgdal and rgeos packages. Because the longitudinal range of the management units was much greater than the latitudinal range, pairwise distances were only calculated between glades within the same management unit to prevent biased isolation measures of glades near the edges of the longitudinal gradient. Thus, “isolation” as described here is a modified version of the two-or-more nearest neighbors metric (Moilanen and Nieminen 2002). Pairwise distance values were used to calculate standardized distance, D_ST, which is a function of pairwise distance between centroids, D_C, and pairwise distance between edges, D_E (Krasnov et al. 2010). We calculated standardized distance between 1) sampled glades only, and 2) all glades within the land unit; the latter accounts for potential stepping-stone effects. We ran separate models for each calculation of standardized distance,

$${\displaystyle \begin{array}{l}{\displaystyle D_{ST}=\frac{D_C}{D_C-D_E}.}\end{array}}$$

(1)

Patch shape was quantified using the corrected perimeter–area ratio, CPA, which was a function of perimeter length, L, and area, A (Farina 2008). This metric is an improvement upon the standard perimeter–area ratio, which varies with the size of the patch even when patch shape is held constant. CPA reduces this bias and provides a more accurate index for quantifying shape when there is a high degree of variability in patch area,

$${\displaystyle \begin{array}{l}{\displaystyle CPA=\frac{0.282\times L}{\surd A}.}\end{array}}$$

(2)

Though IBT assumes habitat islands are homogenous with respect to environmental filters other than area and isolation, variation in insular species composition may be a result of variation in environmental characteristics such as vegetation composition (Lomolino 2000). Glades in particular require periodic burns to maintain ecosystem health and can be invaded by trees and shrubs in the absence of periodic burns. Thus, measuring vegetation characteristics can act as a proxy for variation in burn regimes. We collected data corresponding to vegetative characteristics by placing a 1-m² quadrat at 10 randomly selected trap stations; the percent cover for grasses, forbs, shrubs, and trees was recorded, along with total cover. The dimensionality of the vegetation data was reduced using a principal components analysis (PCA) using the prcomp function in R.

Multispecies occupancy model

Area and isolation effects on species richness can be masked if species are detected imperfectly during sampling. Detection rates of species are variable due to behavioral traits and landscape characteristics (Iknayan et al. 2014), resulting in biased estimates of species occurrence for rare or hard-to-detect species (MacKenzie et al. 2005). Recent advances in statistical methodology and computing have led to community-level hierarchical models (or multispecies occupancy models, MSOMs) that leverage information from across the community to estimate species-specific occupancy and detection (Dorazio and Royle 2005; Zipkin et al. 2010).

We constructed a single-season hierarchical model in which occurrence of species i at site j, denoted Z_ij, was equal to 1 if the species is present at the site and 0 otherwise. The true occupancy state, Z_ij, was assumed to be an outcome of a Bernoulli trial, Z_ij~ Bern(Ψ_ij), where Ψ_ij is the probability that species i is present at site j. In cases of few or no observations of a species, the true occupancy state was unknown and imperfectly observed, which would lead to confounded estimates of Ψ_ij. However, because site j was sampled multiple times over a short period, observation data were incorporated into a detection model that distinguished between true species absences and nondetection (MacKenzie et al. 2002). The detection model was defined as x_ijk~ Bern(p_ijk* Z_ij), where p_ijk is the probability of detecting species i at site j during sampling period k, given that species i is present at site j. The dependency of the detection model on the occupancy state ensures that detection probability was zero when the species was not present at a site.

Because MSOMs leverage information from the entire community, rare or poorly detected species which lack sufficient data to model individually can be analyzed by “borrowing” data from the community as a whole (Iknayan et al. 2014). This was done by assuming each species-level parameter was drawn from a common community-level distribution, or “hyperparameter” (Kéry and Royle 2009; Zipkin et al. 2009). Undetected species for which there are no data were represented by including all-zero encounter histories in a process known as zero augmentation (Royle 2004; Kéry and Royle 2009). Thus, MSOMs should generate more accurate estimates compared to single-species models by accounting for the similarity in the ecology of all members (Link and Sauer 1996; Sauer and Link 2002; MacKenzie et al. 2005).

To account for potential differences in occupancy between sampling years, we used an unpaired-site model (Tingley and Beissinger 2009; Tingley et al. 2012). This framework increases model accuracy when sample sizes are small by treating different sampling periods as sets of independent sites. Differences between sampling periods—in this case, sampling years—can still be modeled by including “year” as a covariate. In addition to year, patch area, patch isolation, and patch shape were incorporated into the model using the logit-link function (Dorazio and Royle 2005):

$${\displaystyle \begin{array}{ll}{\displaystyle logit\left({\mathrm{\Psi }}_{ij}\right)=}& u_i+\alpha 1_{i}Are{a}_j+\alpha 2_{i}Isolatio{n}_j\\ & +\alpha 3_{i}Shap{e}_j+a{4}_{i}Ve{g}_j+a{5}_{i}Yea{r}_j.\end{array}}$$

(3)

Detection probability p_ijk was modeled in a similar manner using ordinal date:

$${\displaystyle \begin{array}{l}{\displaystyle logit\left(p_{ijk}\right)=v_i+{\beta }_{1i}Dat{e}_{jk}.}\end{array}}$$

(4)

With the exception of year, all covariates were scaled to have a mean of 0 and a standard deviation of 1. Model parameters were estimated using a Bayesian analysis and uninformed priors for the community-level hyperparameters in the program OpenBUGS (Spiegelhalter et al. 2007), based on code modified from Zipkin et al. (2010). Significance of covariates were evaluated using the 95% credible interval (CI); CIs that did not overlap zero were considered significant. We used three chains of length 16,000 including a burn-in of 8,500, and posterior chains were thinned by 35 to reduce autocorrelation. Convergence was evaluated using the R-hat statistic (Gelman and Rubin 1992), and we assumed the model converged when R-hat was near 1.

Statistical analysis

Effects of covariates on community and species-specific occupancy were evaluated using the hyperparameter and species-level parameters, respectively. Specifically, we examined sign and whether the 95% CI overlapped zero. Species richness, including unobserved species, was calculated for each glade by summing the estimated species in the occurrence matrix (Zipkin et al. 2010). Because this modeling framework did not explicitly include relationships between site-level richness and covariates, effects of covariates on species richness were evaluated using a general linear model (GLM) using patch area, isolation, and shape as predictor variables. Year was included in each GLM to account for temporal variation in richness. Detection error-corrected estimates of species richness derived from the model were used as the response variable.

IBT predicts that small islands have higher rates of extinction than larger islands; thus, small islands should also exhibit a higher rate of temporal turnover. To test this prediction, 22,500 potential occurrence matrices for each year were drawn from the Bayesian posteriors, and the Jaccard dissimilarity index was used to compare sites between the two years using functions in the R package betapart (Baselga and Orme 2012); a large value of the Jaccard index indicates greater difference than a smaller value (Jaccard 1901). Associations between patch area, patch isolation, and patch shape and Jaccard values were assessed using the Spearman’s correlation coefficient. Confidence intervals were calculated using a general additive model (GAM).

To determine whether differences between years were due to species replacement or species gain or loss, we used functions in betapart to determine the amount of turnover (i.e., replacement) and nestedness (i.e., gain or loss) between years. To test whether species nestedness was greater than expected by chance, we treated the contributions of the two variables as the outcome of a Bernoulli trial in which a success occurred when nestedness was greater than turnover. The probability of success was set at 0.5. All analyses were carried out in R version 3.3.1 (R Core Team 2016) unless otherwise noted.

Results

A total of 344 individuals were captured, representing eight species in the glade network (Table 1). The most frequently detected species was Neotoma floridana, which was detected on 11 glades in 2016 and eight glades in 2017; the least frequently detected species was Tamias striatus, which was detected on three glades in 2016. Overall, detections in 2017 were less than in 2016 across all mammal species.

Table 1.—

Occupancy of mammal species at each site. X = present in 2016, ▲ = present in 2017, ● = present in both years.

Site	Neotomafloridana	Sigmodonhispidus	Peromyscus leucopus	Peromyscus attwateri	Peromyscus maniculatus	Reithrodontomys fulvescens	Sylvilagus floridanus	Tamiasstriatus	Observed richness 2016	Estimated richness 2016	Observed richness 2017	Estimated richness 2017
0	X	●	●		▲				3	3.54	3	3.50
1	●	●	X	X	X		X		6	6.16	2	2.18
2	●	●	X	●	X		X		6	6.15	3	2.79
3	X	●		●			X		4	4.14	2	2.67
4	X		●	X	●			X	5	5.39	2	2.41
5	X	X	●		●			X	5	5.43	2	2.30
6	●	X	●		●	X			5	5.28	3	3.09
7		X	X	X	X				4	4.81	0	1.28
8	●	X		X					3	3.51	1	1.19
9	●			●			▲		2	2.85	3	3.03
10				●					1	1.96	1	1.39
11	▲		X	▲	X				2	2.63	2	2.12
12		X	X		X	X			4	4.51	0	1.33
13		X	●		X	X			4	4.62	1	1.38
14	●	●	X	X	X	●		X	7	7.33	3	3.40
15	●	●	X	X	X	▲			5	5.33	3	3.14

Glades ranged in size from approximately 3,500 to 92,000 m². The PCA of the vegetation data yielded two principal components with an eigenvalue greater than 1, which together explained 70.8% of the variation in the vegetation data (Fig. 2). PC1 and PC2 were correlated with patch area and shape, respectively, so vegetation was not included in the model.

Fig. 2.—

Results of the PCA for vegetation data. PC1 and PC2 explained a cumulative 70.8% of the variation in the data. Each point represents a glade; colors denote years. PC1 was negatively associated with area; PC2 was positively associated with the shape variable.

The MSOM yielded an estimate of 8.63 species across the region. Year was the only significant occupancy covariate; the community mean for the coefficient was different than 0 (Estimate, 95% CI), and all species had lower occupancy in year 2 (Fig. 3). There was moderate, albeit nonsignificant, evidence for a community-wide positive effect of area on occupancy (Estimate, 75% CI), but there were no significant species-specific effects of area (Fig. 3). Occupancy probabilities of Peromyscus attwateri and P. maniculatus were higher in glades with less edge relative to area, whereas occupancy probability of Sigmodon hispidus was higher in glades with greater edge relative to area (Fig. 3).

Fig. 3.—

Species-specific responses to A) area, B) isolation, C) shape, and D) year. Error bars represent the 95% credible interval; bars that do not overlap 0 (dashed line) were considered significant and are marked with an asterisk (*). Species abbreviations are as follows: NEFL = Neotoma floridana, PEAT = Peromyscus attwateri, PELE = P. leucopus, PEMA = P. maniculatus, REFU = Reithrodontomys fulvescens, SIHI = Sigmodon hispidus, SYFL = Sylvilagus floridanus, TAST = Tamias striatus.

The general linear model was significant (F_4,27 = 28.9, P < 0.001) and explained 63.9% of the variation in species richness. Coefficients for area (β = 0.586, t₂₇ = 3.400, P = 0.001) and year (β = −2.289, t₂₇ = −9.402, P < 0.001) were significant (Fig. 4), but coefficients for isolation and patch shape were not different from zero (β = −0.095, t₂₇ = −0.568, P = 0.572; β = 0–0.095, t₂₇ = −0.500, P = 0.619, respectively).

Fig. 4.—

Site-level mean species richness plotted against log-transformed area. Mean species richness increased with patch area for both years (β = 1.898, t₂₇ = 2.469, P = 0.020). Year 2 also had significantly lower richness than year 1 (β = −2.277, t₂₇ = −6.305, P < 0.001).

The mean Jaccard similarity index for each site ranged from small (J = 0.096) to large (J = 0.961). There were no significant associations between Jaccard dissimilarity values and patch area (ρ = −0.038, P = 0.891) or patch shape (ρ = −0.45, P = 0.082). There was a significant and positive association between the Jaccard dissimilarity values and isolation (ρ = 0.515, P = 0.044; Fig. 5). Between years, there was significantly greater nestedness than turnover in communities than predicted by chance (number of successes = 12, trials = 16, P = 0.038).

Fig. 5.—

Jaccard index values plotted with respect to scaled isolation. A GAM estimated a curvilinear relationship (black line), with the 95% confidence interval shown in gray.

Discussion

We found mixed evidence for IBT when testing against observations of the small mammal fauna in Ozark glades. Consistent with predictions of IBT, patch area had a significant effect on mean estimated species richness, yet area did not significantly influence the occupancy probability of any species, and the community-wide area effect was not significant. However, each species in the glade network exhibited a positive, nonsignificant response to area (Fig. 3A), and the significant relationship between area and species richness may result from the combined responses of the individual species. As expected for our habitat islands, predictions associated with isolation effects on occupancy and species richness were not supported. Likewise, patch shape did not affect the community estimates as a whole, supporting that the communities had few habitat or edge specialists. Finally, compositional changes did not necessarily track predictions of IBT, in part due to potential violations of assumptions of equilibrium as well as non-neutral species dynamics.

Isolation effects on occupancy and species richness were not significant, a common finding in studies of habitat islands (Brown 1971; Cook et al. 2002; Watling and Donnelly 2006). In the glade network, colonization of generalists such as N. floridana and P. leucopus from the matrix likely added noise to the data set (Cook et al. 2002), effectively masking isolation effects. Additionally, nearest-neighbor isolation metrics, such as the one used here, represent a simplified landscape, which might fail to capture the realities of a heterogeneous matrix (Ricketts 2001; Moilanen and Nieminen 2002; Kindlmann and Burel 2008; Biswas and Wagner 2012). Isolation metrics that account for heterogeneity in the matrix, such as landscape friction (McRae et al. 2008), could provide a more accurate picture of isolation effects, but they also require assumptions of metrics of conductance and resistance.

Small mammals are known to respond to environmental shifts, particularly with respect to their population dynamics (Schnurr et al. 2002; Jacob 2003; Kelt et al. 2008; Thibault and Brown 2008). Between 2016 and 2017, we observed a decline in species richness within glades. It is unlikely this decline was due to trap mortality, as we only observed one mortality event in over 600 captures. Declines in species richness were possibly a response to reduced resources related to unusually high rainfall prior to the second trapping season (Guinan 2017). Despite the regional decline in species richness, larger glades maintained higher mean estimated species richness than smaller glades. The observed differences in composition between years was more often due to species nestedness than turnover, identifying species loss as the underlying mechanism. Essentially, the population dynamics of some of these species, specifically S. hispidus and P. leucopus, may create extinction debt scenarios during resource-rich periods. We also noted a significant, positive association between the Jaccard dissimilarity index and patch isolation. This finding was unexpected given the lack of significant isolation effects on species richness and individual species occupancy, and confounding factors may be at work. However, this pattern may be explained in the context of IBT, as the low colonization rates predicted to occur on more isolated glades may be hindering recolonization following local extinction events.

Patch shape can influence species richness and composition because of the associated degree of edge available. Species richness can increase with increasing edge due to co-occurrence of focal habitat and matrix species (Ingham and Samways 1996; Magura 2002; Guirado et al. 2006; Batáry et al. 2014) or decrease due to behaviors such as edge avoidance by certain species (Orrock and Danielson 2005; Bieringer et al. 2013; Besnard et al. 2016). Indeed, several species encountered here (S. hispidus, P. leucopus, and N. floridana) have been shown to have positive relationships with edge elsewhere (Cummings and Vessey 1994; Manson et al. 1999; Beckmann et al. 2001). Additionally, increasing edge increases the probability that the patch will be encountered by a dispersing individual (Collinge and Palmer 2002) and reduces population persistence (Bevers and Flather 1999; van Kirk and Lewis 1999), thereby increasing the amount of compositional change in a habitat patch. Our observations for individual species responses to patch shape differed in significance and direction (Fig. 3), suggesting that the varying ecology of individual species in a community may be contributing to the inconsistent effects of shape in the literature (Ewers and Didham 2005).

Previous attempts to study compositional change are subject to many criticisms. Artificially inflated turnover rates, termed pseudoturnover (Nilsson and Nilsson 1983), are usually the result of sampling error, and this error can be a result of poor study design, such as when sampling methods are inconsistent (Gilbert 1980; Nilsson and Nilsson 1982; Abbott 1983) or because of variation in detectability between species (Nilsson and Nilsson 1983). By using MSOMs to generate presence-absence matrices that have been corrected for detection error, we accounted for potential pseudoturnover in our data set and yielded estimates that should be robust.

Multispecies occupancy models can be a powerful tool for conservation. By assuming that occupancy and detection probabilities of each individual species are drawn from a common distribution, rare species for which data are sparse can be modeled by deriving information from the entire community (Zipkin et al. 2009, 2010; Iknayan et al. 2014). The ability to evaluate rare species in the context of the community, as well as the use of occupancy rather than abundance data, allows for more efficient use of data than single-species models and can reduce the time and cost necessary for effective biological monitoring (Manley et al. 2005; Zipkin et al. 2009; Rich et al. 2017).

The hierarchical framework used in MSOMs is highly flexible and can be readily applied to a variety of situations (Dorazio and Royle 2005). MSOMs have been extended to quantifying range shifts in response to climate change (Tingley et al. 2012), accounting for open communities by allowing for periods of immigration and emigration (Kéry et al. 2009), predicting species occurrences at unsampled sites (Zipkin et al. 2012), evaluating the effects of land use changes on occupancy (Goijman et al. 2015), and integrating occupancy data from multiple spatial scales (Rich et al. 2017). The efficiency and adaptability of MSOMs provides a solid framework for investigating patterns in species occupancy and informing management decisions.

IBT in its most basic form assumes ecological equivalence of all species in the island network. This assumption limits the scope of the theory to the number of species on an island, when the composition of species may be of equal (or even greater) interest (Lomolino 2000). This is particularly true in the context of conservation, where decisions based on multiple criteria are more likely to be effective than those based solely on counts of species richness (Fleishman et al. 2006). Our use of a MSOM to evaluate the applicability of IBT in the glade network added a layer of species-based information, revealing how individual species may be driving the patterns we saw in the estimates of species richness. Using a MSOM in this manner provided a more complete picture of the mammal community when combined with the species richness estimates of traditional IBT.

While IBT is useful for explaining broad patterns of species richness in islands, the paradigm fails to capture the complex processes structuring glade communities. Consistent with this study, various taxa in the glade network appear to be affected by environmental factors other than area and isolation, including vegetation structure (Bergmann and Chaplin 1992), fire regimes (Templeton et al. 2001; Brisson et al. 2003), and top-down regulatory processes (Van Zandt et al. 2005). Additionally, many species found in the glade network are affected by the intervening forest matrix (Nelson et al. 2013), which is not accounted for in traditional IBT. A landscape-level approach that accounts for matrix effects and species-specific responses to habitat variation may be more suitable for application to the communities of the Ozark glades.

Acknowledgments

We thank C. Davis, H. Williams, C. Adkins, M. Sieron, and S. Meilink for field assistance, and D. Finn, J. Greene, H. York, and two anonymous reviewers for helpful comments on the manuscript. This study was funded by the Missouri State University Graduate College and Missouri State University Biology Department. Permission for trapping small mammals was given by the Missouri Department of Conservation (Wildlife Collector’s Permits #16141 and #17155), the US Forest Service (Authorization ID AVA1609), and the Missouri Department of Natural Resources. All methods were approved by the Missouri State University IACUC (IACUC ID #16-020.0).