Pocket gopher (Thomomys talpoides) soil disturbance peaks at mid-elevation and is associated with air temperature, forb cover, and plant diversity

ABSTRACT

Burrowing mammals can be ecosystem engineers by increasing soil aeration and erosion and altering the structure of plant communities. Studies that characterize the constraints on the distributions of fossorial mammal disturbances to soil can help predict changes in ecosystem engineering under future climates. We quantified the density of soil disturbances caused by Thomomys talpoides (northern pocket gopher) over replicate elevation gradients spanning 2,700–4,000 m a.s.l. in the Upper Gunnison Basin, Colorado, USA. As a conceptual framework for predicting biogeographic variation in soil disturbance, we used the abundant center hypothesis (ACH), which proposes that species abundance declines monotonically away from the most abundant location in its distribution, with the assumption that ecosystem engineering scales with gopher abundance. We also evaluated the relative importance of abiotic and biotic variables as correlates of soil disturbance. Gopher disturbance peaked at mid elevations (~3,150 m), supporting the ACH. The best model for predicting gopher-caused soil disturbance contained both abiotic and biotic variables, with increased soil disturbance where mean annual temperature, forb cover, and plant diversity were greatest. Results suggest that mountain ecosystems may experience increases in gopher-caused soil disturbance as climate warms, possibly accompanied by increases in plant diversity and forb cover.

KEYWORDS:

• Alpine
• subalpine
• climate change
• Colorado Rocky Mountains
• niche/distribution modeling

Introduction

Burrowing mammals can have large effects on ecosystems that cascade to soil physical properties, plant composition, and primary productivity. However, little is known about the distribution and abundance patterns of bioturbations caused by fossorial mammals. The northern pocket gopher (Thomomys talpoides, Geomyidae; hereafter, pocket gopher or gopher) is an ecosystem engineer that creates physical soil disturbances (e.g., eskers, tunnels, mounds) by burrowing belowground and pushing subsoil to the surface. Gopher disturbances loosen the soil (Butler and Butler 2009; Turner et al. 1973), which can increase water infiltration (Hansen and Morris 1968; but see Yurkewycz et al. 2014) and soil erosion (Sherrod and Seastedt 2001). Pocket gophers dig an elaborate network of tunnels (Knight 2009), and, on average, can excavate 26.1 m³/ha a year of soil (Smallwood and Morrison, 1999). This large soil disturbance creates multiple microenvironments (Huntly and Inouye 1988). For example, pocket gopher disturbance can raise median subsurface temperatures while decreasing median surface temperatures when compared to undisturbed soils (Whitesides and Butler 2016). Furthermore, burrowing and plant foraging by pocket gophers can influence succession and plant species composition (Reichman and Seabloom 2002; Sherrod, Seastedt, and Walker 2005). Pocket gophers will come to the surface to cut plants away from their tunnel openings, and their diet includes bulbs, roots, tubers, and, occasionally, aboveground plant tissues (Rios and Álvarez-Castañeda 2007). For example, in the high elevations of central Utah (3,000–3,400 m), gophers primarily consumed and cut dandelion (Taraxacum officinale), penstemon (Penstemon rydbergii), sweet sage (Artemisia discolor), and yarrow (Achillea lanulosa) (Aldous 1951).

Despite the ecological dominance of the pocket gopher and its documented impacts on plant community structure, ecosystem processes, and geomorphology, little is known about the ecology and distribution of pocket gophers or their ecosystem engineering effects at high elevations. Although the range of pocket gophers can extend to approximately 3,800 m a.s.l., most studies have focused on low elevations (<2,500 m; e.g., Case, Halpern, and Levin 2013; Huntly and Inouye 1988). Among past observational studies, the highest elevation with reported data on pocket gopher abundance or incidence was 3,962 m (Table 1). Most prior studies examined a single altitudinal gradient (one mountain), and the lack of replication limits inference on the influence of elevation per se, because other factors could drive patterns along a single gradient. Most studies only acquired presence data or did not perform statistics on abundance patterns (Table 1). Furthermore, many prior studies included few sites along the altitudinal gradient, precluding the ability to detect nonlinear patterns, such as a hump-shaped distribution with elevation. Theory on species abundance patterns predicts the possibility of a quadratic relationship: the abundant center hypothesis (ACH; Brown 1995; Gaston 2003) proposed that a species abundance declines monotonically away from the most abundant location within its range. Thus, additional data are needed to examine biogeographic patterns in light of existing theory.

Table 1. Summary of prior studies documenting the distribution of northern pocket gophers (Thomomys talpoides) along elevation gradients

Study	Elevation Range (m)	Number of Gradients	Number of Sites/Gradient	Altitudinal Pattern
Findley and Negus (1953)	2,865–3,612	NR	24	Presence- and absence-only data; no statistics
Hansen and Bear (1964)	2,835–3,658	1	3	Highest at mid-elevation; no statistics
Miller (1964)	1,372–3,962	NR	NR	A summary of presence in Colorado
Rickart (2001)	1,615–3,355	19	NR	Presence-only data; no statistics
Tryon and Cunningham (1968)	2,255–3,099	1	3	Highest at mid-elevation in 5 of 6 years; no statistics

Note. NR = Not Reported.

Climate factors associated with elevation could influence pocket gopher distributions in mountain ecosystems either directly by affecting their physiology or indirectly by affecting their food resources. In the Colorado Rocky Mountains, air temperature, atmospheric pressure, and soil nutrients decline with increasing elevation, but precipitation increases (Dunne, Harte, and Taylor 2003; Inouye and Wielgolaski 2013). Prior work showed that pocket gopher body size and fat content increased with higher elevation, but population densities were greatest at the middle of three sites spanning 2,255–3,099 m (Tryon and Cunningham 1968). An herbicide experiment additionally suggested that this distributional pattern was associated with the relative abundance of grasses and forbs. A greater grass-to-forb ratio at higher elevations correlated with lower gopher abundance (Keith, Hansen, and Ward 1959), because pocket gophers prefer forbs over grasses (Aldous 1951; Vaughn, 1967; Ward 1960). Metabolism and thermoregulation in pocket gophers is typically an insignificant aspect of pocket gopher survival because of the relatively constant temperatures in their burrows (Gettinger 1975). However, their metabolic rate was considerably greater than other fossorial rodents following movements outside of their burrows during critical temperatures (<25–26°C), suggesting a metabolic cost of cold weather (Gettinger 1975). Conversely, gophers facing heat stress above 37.5°C experienced absolute mortality (Gettinger 1975).

Understanding the distributions of individual species and their engineering impacts on ecosystems can improve our ability to predict range shifts along with their ecosystem consequences under climate change (Parmesan 2006; Pearson and Dawson 2003). Here, we use the ACH to set up expectations for how gopher-caused disturbance varies with elevation, with the assumption that disturbance scales with gopher abundance, as shown in prior research (Smallwood and Erickson 1995). Although the ACH has been supported (e.g., Fenberg and Rivadeneira 2011; Vaupel and Matthies 2012), its generality remains inconclusive (Dallas, Decker, and Hastings 2017; Sagarin and Gaines 2002). The ACH is predicated on spatial autocorrelation among the multiple dimensions of a species’ niche, suggesting that an analysis of the underlying abiotic and biotic factors that control, or at least correlate with, abundance would improve evaluations of the ACH (e.g., Martinez-Meyer et al. 2013).

This study adds to the general evaluation of the ACH by assessing the density of soil disturbances caused by T. talpoides along replicated elevation gradients in subalpine to alpine meadows. We focus on the biogeographic patterns of disturbance caused by gopher activity to highlight the spatial extent of this important ecosystem engineering process. First, we asked: How does pocket gopher disturbance vary across elevation gradients? To address this question, we compared models of the relationship between disturbance abundance and elevation that reflect the ACH prediction that abundance will decline in all directions away from the most abundant location (hump-shaped) versus an alternative hypotheses (random with respect to elevation or linear with elevation). Second, we also evaluated climate, soil, geography, and vegetation to ask: What are the biotic and abiotic environmental correlates of gopher disturbance? We expected gopher disturbance to decrease under low temperatures (Gettinger 1975), increase with relative cover of forbs (their primary food source; Keith, Hansen, and Ward 1959), and increase with soil depth (Miller 1964; Reichman and Seabloom 2002).

Materials and methods

Site description and gopher disturbance sampling

We assessed gopher surface disturbance of soils for sixty-seven sites spanning six replicate mountains (elevation transects, Figure 1; Table S1) in the Upper Gunnison Basin, Gunnison County, Colorado, USA. In the region, T. talpoides is the only mammal to cause significant burrowing activity in the soil (Findley and Negus 1953). Transects were chosen to avoid pseudoreplication of the lowest sites by sampling geographically distinct watersheds that were separated by ridges. Sites were sampled at every 100 m in elevation from the peak of the mountain to approximately 2,700 m in elevation. This resulted in nine to thirteen sites per transect, depending on elevation and topography. Sites were grassland habitat and were free of trees. To our knowledge, this is the largest altitudinal survey of northern pocket gopher disturbance to date.

Figure 1. Map of the sites and mountain transects (see legend) in the Upper Gunnison Basin, Colorado, USA. Inset A is a map of Colorado State, with a black box representing the area of the main map. Inset B is a map of the United States, with a gray box representing the area of inset A

At each site, we measured pocket gopher disturbance along three 40 m belt transects with a Precimeter 305 rolatape (R. C. 2002, Spokane, Washington, USA). Each transect was visually surveyed for gopher disturbance by operating the rolatape with a 1 m wide, 2.5 cm diameter PVC pipe attached perpendicular to its main axis to create a belt transect that increased the surface area of observation (following the methods of Beck and Hansen 1966; Hansen and Beck 1968). Each time a gopher disturbance was observed along the belt transect, it was counted with a hand-held tally counter. This “transect-mound count” is the most accurate non-capture-based method for estimating pocket gopher abundance (Smallwood and Erickson 1995); although, in addition, our method counted other forms of gopher disturbance. We defined gopher disturbance to soil as standing eskers/tunnels, infilled tunnels, or soil mounds. Photos of various gopher disturbances can be found in Butler and Butler (2009). Eskers and infilled tunnels are created by gophers foraging under snowpack during winter, and infills are the disintegration of these tunnels (Knight 2009; Winchell et al. 2016). In our data collection we did not differentiate between the types of disturbance. When eskers followed the transects longitudinally, each 1 m section of the esker, in relation to the 1 m wide field of view, was counted as a single disturbance. However, this observation was rare, given the typical meandering patterns of eskers. The density of gopher disturbance was estimated as the number of gopher disturbances per m², calculated over the three parallel 40 m belt transects per site.

Biotic and abiotic environment characterization

Vegetation was sampled during July 28–September 24, 2014. At each site, a 20 m vertical transect was arranged downslope, and three 20 m horizontal transects were placed across the top (0 m), middle (10 m), and bottom (20 m) positions of the vertical transect. Starting at the 0 m position of the transects, every 2.5 m of each transect was visually analyzed for percentage vegetative cover in a 20 cm × 20 cm quadrat. The percentage cover of forbs, graminoids, and bare ground was determined such that the maximum cover in a quadrat summed to 100. This sampling method resulted in thirty-three dispersed quadrats per site (nine positions per transect, no duplicate quadrats at the transect intersection), for a total of 1.32 m² per site, representing an area of 400 m². From this, six vegetation variables were constructed at the site level: forb to graminoid ratio (percent forb cover divided by the percent graminoid cover, including sedges and the like; hereafter referred to as “forb:grass ratio”), proportional forb cover (percent forb cover divided by total vegetative percent cover), total forb cover (summed forb cover across the thirty-three quadrats), total vegetative cover (summed cover across the quadrats), plant species richness, and inverse Simpson’s diversity index (calculated using the vegan package in R; Oksanen et al. 2015).

We also characterized several aspects of the nonclimatic abiotic environment. Soil depth was assessed by inserting a tile probe (Tile Probe 1.5 m length by 2 cm diameter, AMS Inc., American Falls, Idaho, USA) into the ground until bedrock was reached for ten independent points placed every 5 m along two 20 m transects. The slope and aspect for each site was obtained on site and also by analysis of open-source digital elevation models (DEMs) provided by the U.S. Geological Survey (one-third arc-second elevation product, Dollison 2010) within the QGIS GIS platform (QGIS Development Team 2017) by using GRASS GIS algorithms (GRASS Development Team 2017).

We collected climate data from twenty-nine regional weather stations to parameterize multiple regression models predicting the site-level climate based on elevation, aspect, and slope. Other studies modeling spatial patterns of a species have used freely available climate products, such as BIOCLIM (Hijmans et al. 2005). However, given the finer spatial scale of our data, where short linear distances have high elevation differences and are predicted to have substantially different climates (Pepin and Losleben 2002), such climate products are not useful (e.g., resolution ≥800 m²). All weather stations were located between 38.1628°–39.7649° latitude and −106.3397°–107.8741° longitude and covered various time ranges (starting points span 1981–2010; Supplementary Methods S1). From each weather station, we calculated yearly summaries of growing degree days (GDD) using a base temperature of 0°C (McMaster and Wilhelm 1997), mean annual temperature (MAT), mean annual precipitation (MAP), snow depth, and mean summer precipitation (SumP). We then constructed multiple regression models where the climate variable of interest (e.g., MAT) was modeled as a function of elevation, the sine and cosine of the aspect, and the slope, using maximum likelihood methods (lme4 package in R; Bates et al. 2015) with weather station and year as random effects. The sine of aspect in degrees creates a gradient from −1 to 1, representing west- to east-facing slopes, while the cosine of aspect creates a gradient of −1 to 1, representing south- to north-facing slopes. We used an information theoretic approach (Anderson 2008) to determine which combination of predictor variables were most predictive, given the set of candidate models. The best model for every climate parameter other than SumP (R² = 0.63) had a conditional R² (conditional on the group level effect) greater than 0.88. The estimated coefficients from these models were combined with the elevation, aspect, and slope from the survey sites to predict climate data for each gopher survey site. Models were aimed at predicting the recent climate experienced at the sites. A more detailed, step-by-step methodology of our climate modeling can be found in Supplementary Methods S1, including weather station geographic information, model details, and results.

Statistical analyses

First, we used regression models to evaluate the three hypotheses on the shape of elevation–gopher disturbance relationships. Gopher disturbance was non-normally distributed; thus, we used natural log transformation (ln(Y + 1); Sokal and Rohlf 1995). We fit three models using the lm function in the stats package, and all analyses were performed in R (R Core Team 2015): (1) A null intercept-only model, predicting that the density of gopher disturbance was a random function of elevation; (2) a linear model with gopher disturbance predicted as a linear function of elevation; and (3) a model with a quadratic term specifying a mid-elevation peak in gopher disturbance as predicted under the ACH (Brown 1984, 1995). Each model also included mountain transect identity to account for spatial autocorrelation among sites within a transect. We used model selection procedures (Anderson 2008; Burnham and Anderson 2002) to compare models using AICc values. In addition, we computed AICc weights, which can be interpreted as the conditional probability that a given model is the best, given the set of candidate models. AICc and weights were derived using the AICc and Weights functions from the MuMIn package (Barton 2015). We report p values, F statistics, and r² or R² values for each model. An alternative approach using a generalized additive model that fit splines did not improve the fit to the data, nor did it suggest a different shape to the relationship, such as an asymptote at the lowest elevations rather than a decline at low elevation (data not shown). Additionally, we fit a model with gopher disturbance expressed as density per m² to show that the pattern was invariant to log transformation (see Figure 2b).

Figure 2. Gopher disturbance peaks at mid-elevation. Each point represents a sampled site, and different symbols reflect each mountain transect. Plot (a) reflects our model selection results with the relation between elevation and the natural log of gopher disturbance described by y = −41.75 + 2.94E-2*x-4.65E-6*x² (p < 0.001; R² = 0.55). Plot (b) presents the same model, but with gopher disturbance (dist.) presented as density per m², described by y = –1.42 + 0.49E-3*x-1.50E-6*x² (p < 0.001; R² = 0.41). Parameter estimates were derived from least squares regression that included a linear and a quadratic term

To address our second question, we assembled site-level data to compare factors that are hypothesized to limit gopher disturbance and distribution. We first scaled all predictor variables to mean = 0 and standard deviation = 1 to evaluate predictors on comparable scales. To avoid multicollinearity (Zuur et al. 2010), we calculated variance inflation factors for candidate predictors and proceeded with the set of variables with minimal collinearity (vif <10, calculated using the vif function from the usdm package; Naimi 2013). This returned six variables with low collinearity: forb:grass ratio, inverse Simpson’s plant species diversity, mean annual temperature (MAT), natural log of soil depth, and the sine/cosine of aspect. Predictor variables were grouped as abiotic (MAT, soil depth, sine and cosine of aspect) or biotic (forb:grass ratio, inverse Simpson’s diversity). Biotic and abiotic sets were also combined to create a full model. The candidate set of four models included an intercept-only (null), abiotic-only, biotic-only, or the full model; hypotheses were evaluated in a model selection framework, as described previously (Anderson 2008; Burnham and Anderson 2002). Models were fit using maximum likelihood generalized mixed effects models, with mountain transect as the random effect, and a Poisson distribution for nonnormal count data (Bolker 2008) with a natural log link function (glmer via the lme4 package; Bates et al. 2015). We report AICc, ΔAICc, AICc model weights, and df for each model. We did not conduct model averaging of parameter estimates, because this can be problematic (Cade 2015), and model weights were strongly divergent (ΔAICc > 75; Anderson 2008).

Results

The ACH was supported for gopher disturbance, with a predicted quadratic relationship with elevation (Figure 2, Table 2). Gopher disturbance peaked at approximately 3,150 m a.s.l. and declined with elevation in both directions, with a stronger decline toward high elevation than low elevation (Figure 2). Model selection procedures overwhelmingly supported the quadratic model against the intercept-only or linear models of gopher disturbance, as evidenced by the lowest AICc values and an AICc weight of 1 (Table 2). The quadratic model also had the largest R² (0.55) and F ratio (10.41; p < 0.001; Table 2).

Table 2. Summary of regression models. Beta and gamma refer to slope terms

Model	Intercept	Linear (β)	Quadratic (γ)	AICc	ΔAICc	AICc Weight	p Value	Adjusted r^2/R^2	F Ratio	df
Quadratic	−41.96	0.029	−4.65 × 10^−6	184.9	0	1.00	1.98 × 10^−8	0.55	10.41	59
Linear	8.83	−1.51 × 10^−3	NA	206.4	21.5	2.12 × 10^−5	1.21 × 10^−4	0.29	5.57	60
Intercept	3.95	NA	NA	218.9	34.0	4.11 × 10^−8	0.018	0.13	2.94	61

When we examined potential correlates of gopher disturbance, there was strong support for the full model, which included forb:grass ratio, inverse Simpson’s plant species diversity, as well as MAT, soil depth, and the sine/cosine of aspect (Figure 3, Table 3). The AICc weights overwhelmingly support the full model over models with only abiotic or only biotic correlates (Table 3). Because predictor variables were scaled prior to analysis, the slope estimates provide effect-size metrics. Within the full model, gopher disturbance increased with forb:grass ratio, inverse Simpson’s diversity, MAT, and soil depth, while it decreased with the sine and cosine of aspect (Figure 3). The slope estimate for MAT was nearly double that of any other predictor, followed in rank order by plant species diversity then forb:grass ratio. The sine and cosine of aspect and soil depth were the weakest predictors.

Table 3. The best model of gopher disturbance supported both abiotic and biotic predictors. The best model was determined via an AICc model selection procedure. ΔAICc reports the difference in AICc that a given model is from the lowest AICc model. AICc weights are interpreted as probabilities that the model is the closest to reality given the set of candidate models

Model	Predictors	AICc	ΔAICc	AICc Weights	df
Full model	Forb:grass, inverse Simpson’s diversity, MAT, soil depth, sine and cosine of aspect	2,273.5	0	~1	59
Abiotic	MAT, soil depth, sine and cosine of aspect	2,352.5	79.0	6.98 × 10^−18	61
Biotic	Forb:grass, inverse Simpson’s diversity	2,618.2	344.7	1.46 × 10^−75	63
Null model	Intercept only	2,678.2	404.7	1.36 × 10^−88	65

Figure 3. Gopher disturbance is predicted to vary with both abiotic and biotic factors (a). All variables were scaled (mean of 0 and standard deviation of 1) prior to analysis, making them unit-less and comparable as effect sizes. Estimates come from maximum likelihood mixed effects models with peak as a random-group level effect. Asterisks denote significance at the p < 0.05 level. Points show means ± s.e. Plots b–g present the relationship of each parameter to predicted disturbance, accounting for the other parameters and group effects. Gopher disturbance was modeled with a log link function (log disturbance counts)

Discussion

Gopher disturbance over elevation gradients supports the ACH

Gopher disturbance to soils was greatest at middle elevations (~3,150 m a.s.l.; Figure 2), supported over six replicated gradients. Given that soil disturbance is highly correlated with abundance for T. talpoides (Beck and Hansen 1966; Ellison 1946; Reid, Hansen, and Ward 1966; Smallwood and Erickson 1995), this result supported the ACH for the regional distribution of the northern pocket gopher. It is common to find support for the ACH (e.g., Fenberg and Rivadeneira 2011; Vaupel and Matthies 2012) despite recent skepticism (Dallas, Decker, and Hastings 2017; Sagarin and Gaines 2002; Sagarin, Gaines, and Gaylord 2006). It is important to note that Brown (1984, 1995) predicated the ACH on three assumptions: (1) species abundance depends on the response of a local population to local conditions, (2) species abundance is determined by how closely local conditions meet the species’ Hutchinsonian (1957) n-dimensional niche requirements, and (3) there is spatial autocorrelation in the multidimensional niche requirements that control species’ abundance and distribution. The last assumption is most central. Several tests of the ACH have failed to first assess whether these assumptions hold, thus causing premature rejection of the hypothesis based solely on observed patterns of abundance. If the Hutchinsonian niche dimensions are not spatially autocorrelated, then applying the ACH would be inappropriate. Thus, the ACH can be rejected via two avenues: (1) if there is evidence that the assumptions are poor and/or unrealistic or (2) if the assumptions hold but there is no abundant center. Spatial autocorrelation in environmental variables is widespread (Legendre 1993; Koenig 1999; but see Helmuth et al. 2002; Sagarin, Gaines, and Gaylord 2006). Our dataset showed significant spatial autocorrelation in four of five gopher disturbance niche dimensions, including the best predictor: mean annual temperature (evaluated using Moran’s I; Table S2). We suggest that a rigorous analysis of Brown’s (1984, 1995) assumptions is key to determining the usefulness and generality of the ACH.

Our results are constrained by lack of data across the full geographical distribution of T. talpoides, and this is particularly evident in the stronger resolution on elevational pattern for high elevations than for low elevations in our dataset (Figure 2). We expect that sampling additional sites at lower elevations than included in our study would increase understanding of potential factors that define the low-elevation range limits of T. talpoides. Despite this limitation, the large environmental diversity over a manageable spatial scale permitted detailed data collection on potential niche dimensions, such as plant composition. In fact, Brown (1984, 1995) initially drew from work on abundance distributions over elevation gradients (Whittaker 1956, 1960, 1965) to formalize the ACH by characterizing the spatial distributions of environmental variables over elevation gradients. The ACH may be upheld more commonly over local elevation gradients (e.g., Herrera and Bazaga 2008; Vaupel and Matthies 2012) than for a species’ full geographic range, but to our knowledge no prior study has made such a comparison.

Niche dimensions of the pocket gopher

When parsing the niche dimensions that could underlie the abundant-center pattern, we found that aspects of climate, soil, geography, and plant communities were important correlates of gopher disturbance. Mean annual temperature was predicted to have a large influence on the distribution of gopher activity, with greater disturbance at sites with warmer mean annual temperature (MAT; Figure 3). Mean annual temperature declines with elevation, suggesting that variation in MAT may drive much of the elevation pattern (Figure 2). Prior work suggested that pocket gophers experienced a high cost of low temperatures because of their relatively fast metabolic rates (Gettinger 1975). At low MAT, gophers were expected to be less abundant because a short growing season limits the vegetation production needed for their survival (Thorn 1978). In contrast, at low elevations warmer MAT may expose gophers to physiological stress related to thermal mortality limits (~37.5°C; Gettinger 1975). The metabolic stress associated with both ends of the gopher’s distribution over mean annual temperature gradients deserves further experimental attention and may drive the hump-shaped gradient in gopher disturbance with elevation.

Vegetation patterns were also strongly correlated with the abundance and distribution of the soil disturbances caused by T. talpoides (Figure 3). Past work has suggested that gophers prefer forbs over grasses (Keith, Hansen, and Ward 1959; Stuebe and Andersen 1985; Vaughan 1967), and our survey supported the expectation that gopher disturbance increases with greater proportions of their preferred food. Future work may benefit from investigating the nutritional quality of the plant species associated with gopher disturbance (e.g., rodent preference for N-rich legumes; Sirotnak and Huntly 2000). We also found that gopher disturbance increased with plant-species diversity, suggesting that pocket gopher abundance increases with potential vegetative diet breadth or that pocket gophers promote plant diversity. The direction of causality for our biotic predictors is ambiguous: Does plant diversity and greater forb cover lead to increased gopher disturbance? Or does gopher disturbance promote plant diversity and increase forb cover? There is strong evidence for the latter from Sherrod, Seastedt, and Walker (2005), who found that experimental disturbance designed to reflect gopher burial increased forb cover, and that plant-species richness increased in recently disturbed soils in alpine environments (Colorado Rocky Mountains, ~3,500 m a.s.l.). Additionally, northern pocket gophers selectively consumed forbs over graminoids and other functional groups in feeding trials (Aldous 1951; Keith, Hansen, and Ward 1959; Stuebe and Andersen 1985; Vaughan 1967). That evidence exists for both pathways of causality may indicate positive feedback loops between gopher abundance and forb abundance or plant-species diversity (Cortinas and Seastedt 1996; Sherrod, Seastedt, and Walker 2005), where gopher disturbance increases forb cover, leading to a greater abundance of the gopher’s primary food source, which in turn increases their population size and level of soil disturbance. Gopher activity could increase plant diversity by processes related to the intermediate disturbance hypothesis (Grime 1973a, 1973b; Huston 2014), given that the disturbance related to gopher burrows does not reach such extremes that it mimics landslides, fires, or other catastrophic events.

Soil depth and aspect were additional, reliable predictors of gopher disturbance across elevation gradients (Figure 3). However, compared to climate and plant biotic factors, their influence on gopher disturbance was minor. Gopher disturbance increased with soil depth in our system (Figure 3; see also Marcy et al. 2013). Pocket gophers tend to occupy the upper soil layer where plant roots are most dense, approximately 6–20 cm (Reichman and Seabloom 2002), and although the lowest average soil depth across our sites was approximately 80 cm, sites with deeper soil profiles had more gopher disturbance. Additionally, soil depth is known to have large effects on vegetation cover, productivity, and composition (e.g., Fridley et al. 2011, 2016).

Our results were partially consistent with past work that found the greatest densities of T. talpoides disturbance on southeast-facing slopes (Thorn 1982). The negative relationship between the sine of aspect and gopher disturbance suggested that gopher disturbance was greatest on western-facing slopes in the study region, while the negative relationship with the cosine of aspect suggested that gopher disturbance also increased toward southern-facing slopes (Figure 3). The wind in the Colorado Rockies generally prevails from the west, leading to snow accumulation on east-facing slopes (Thorn 1982). Deeper snowpack during the winter better insulates gopher dwellings and could reduce reliance on metabolic heat (Andersen and MacMahon, 1981); however, deep snow also may reduce the length of the plant growing season and limit food availability.

Although our models investigated many of the niche dimensions expected to limit T. talpoides disturbance and abundance, there are other limiting factors that were not addressed. For instance, variable snow depth and drifts caused by islands of trees or geological features may provide winter refugia from freezing soils (Winchell et al. 2016), but our study was unable to characterize snow-depth variation at this large number of sites. Soil texture and rock features of the soil substrate are also known to influence gopher abundance (Sherrod and Seastedt 2001; Winchell et al. 2016), but were not available for our sites. Future modeling of fossorial mammals would benefit from considering these important, although difficult to quantify, niche dimensions.

Consequences of climate change for ecosystem engineering by pocket gophers

Pocket gophers are ecosystem engineers; thus, as climate change continues, understanding pocket gopher range shifts—and, importantly, shifts in the distribution of their soil disturbances—may be important for predicting shifts in local plant communities and ecosystem processes (Berteaux et al. 2004; Parmesan and Yohe 2003; Root et al. 2003). The usual response of animals to warming is to migrate toward cooler, higher elevations or latitudes (Berteaux et al. 2004; Parmesan and Yohe 2003; Root et al. 2003). Given that climate change may be dramatic in mountainous ecosystems (Rangwala and Miller 2012), and vegetation responses to altered climate have already been observed (Zorio, Williams, and Aho 2016), we expect that climate change could shift the existing T. talpoides populations up mountainsides (Gettinger 1975). If we assume that MAT limits the upper elevation distribution of the pocket gopher, as predicted by theory (Brown 1995; Dobzhansky 1950; MacArthur 1972) and supported by our results, what will happen to the distribution of T. talpoides in response to climate change? High elevation environments are expected to experience approximately 2–3°C increase in MAT with climate change (Rangwala and Miller 2012), which could increase snowmelt date by two weeks in both spring and fall (experimentally determined by Harte and Shaw 1995, and Harte, Saleska, and Levy 2015). Gopher activity is greatest during the fall months (Bandoli 1981). Extended fall growing seasons under warmer climates with more growing degree days will prolong resource capture periods and may increase T. talpoides population sizes at high elevations. Given our prediction that the high elevational limits of pocket gophers are temperature and growing-season limited, we suggest that gopher disturbance and abundance at high elevation will increase in future climates. If a 3°C increase in the coming decades for the region reflects climatic change approximate to a 430 m increase in elevation (Pepin and Losleben 2002), we predict that a 3°C temperature increase will shift the maximum gopher disturbance from approximately 3,150 m to 3,580 m in elevation. Interestingly, throughout the climate changes of the past sixty-five years, Zorio, Williams, and Aho (2016) found that high elevation herbaceous plant communities generally decreased in plant diversity and forb representation. This result highlights the necessity of long-term experiments and observations to elucidate the interplay of various causes of gopher distributions, as changes in food resources (an indirect effect of climate change) could counteract direct climate-warming effects for pocket gophers.

Conclusion

Our data provide new support for the abundant center hypothesis across elevation gradients for the northern pocket gopher. Our study provides a more comprehensive understanding of pocket gopher disturbance at high elevations compared to previous studies that have focused primarily on low elevation grasslands or had low replication of elevation (Table 1). Our results suggest that mean annual temperature and aspects of plant resources are key niche dimensions for T. talpoides. If habitat suitability for this ecosystem engineer moves upslope with climate warming, our data predict that increases in forb cover and plant species diversity will track upslope increases in the density of T. talpoides disturbances to soil.

Acknowledgments

We thank the Rocky Mountain Biological Laboratory for logistical support and their summer REU program, which hosted SC. We also thank M. Kazenel, C. Forrester, I. Latella, and W. Brown for assistance in the field. M. Kazenel and S. Kivlin assisted with the climate interpolations.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was supported by the National Science Foundation DEB-1354972 awarded to JAR. Work was additionally supported by UNM’s Biology Department Springfield Scholarship awarded to JSL. RMBL GIS equipment and weather station infrastructure was supported by NSF DBI-0420910 and DBI-0821369, respectively.