Better living through conifer removal: A demographic analysis of sage-grouse vital rates

Sagebrush (Artemisia spp.) obligate wildlife species such as the imperiled greater sage-grouse (Centrocercus urophasianus) face numerous threats including altered ecosystem processes that have led to conifer expansion into shrub-steppe. Conifer removal is accelerating despite a lack of empirical evidence on grouse population response. Using a before-after-control-impact design at the landscape scale, we evaluated effects of conifer removal on two important demographic parameters, annual survival of females and nest survival, by monitoring 219 female sage-grouse and 225 nests in the northern Great Basin from 2010 to 2014. Estimates from the best treatment models showed positive trends in the treatment area relative to the control area resulting in an increase of 6.6% annual female survival and 18.8% nest survival relative to the control area by 2014. Using stochastic simulations of our estimates and published demographics, we estimated a 25% increase in the population growth rate in the treatment area relative to the control area. This is the first study to link sage-grouse demographics with conifer removal and supports recommendations to actively manage conifer expansion for sage-grouse conservation. Sage-grouse have become a primary catalyst for conservation funding to address conifer expansion in the West, and these findings have important implications for other ecosystem services being generated on the wings of species conservation.


Introduction
The sagebrush (Artemisia spp.) biome in western North America has experienced numerous threats since the late 1800s including conversion to agriculture, energy development, altered fire regimes, overgrazing, and invasive species [1][2][3]

Conifer management
The Bureau of Land Management (BLM) removed juniper on public land while the Natural Resources Conservation Service in association with the Oregon Department of Fish and Wildlife assisted landowners with juniper removal on private land within and surrounding the treatment area (Fig 1). Treatments generally occurred from late fall to early spring and were designed to maximize shrub retention. Conifer encroachment has been classified into three successional phases: Phase-I is early successional with conifers subdominant in the plant community; Phase-II is mid-successional with conifers codominant in the community; Phase-III is late-successional with conifers dominant in the community [8]. Most of the treated areas were Phase-I to Phase-II encroachment with generally intact understory herbaceous and shrub communities having similar composition and abundance to areas without conifer encroachment (Fig 2). Most treatments (95%) were conducted during late fall to early spring by hand-cutting with brush-and chainsaws to minimize direct negative affects to grouse breeding activities; however, 444 ha were machine cut (e.g., feller-buncher) in fall 2013 to spring 2014. Additional slash treatment of cut conifers was conducted where necessary to reduce woody fuels and vertical structure. Treatments were implemented depending on tree size and density, understory, and landowner preference [on private land] but mostly consisted of cut-leave, cut-lop, cut-burn, and cut-pile-burn. Cut-leave involved cutting trees without additional slash treatment and generally occurred in areas with trees of low size and density. Cut-lop consisted of felling trees and removing tall branches from tree boles to reduce vertical structure and avian predator perches. Cut-burn occurred with larger, denser trees to expose the understory and encourage growth. Generally, cut trees were left to dry for~1 year and then burned individually. Effort was made to burn only individual trees to reduce shrub mortality and burn scars to minimize negative affects of treatments to beneficial habitat. Cut-pile-burn involved felling trees, cutting into manageable pieces, and stacking in small piles for burning when soils were frozen. This technique was used less often due to cost but was deemed necessary in some areas of high tree density to reduce area impacted by slash-burning. Across all treatments, the objective was complete conifer removal, but an attempt was made to leave pre-settlement trees in locations that historically supported juniper, thus some areas still had standing trees after treatment [33]. BLM biologists identified pre-settlement trees using criteria including size, leader growth, crown form, bark, and habitat [8]. For our analysis we grouped all treatments together. Although some treatments occurred from 2007-2011 (<10%), most occurred from 2012-2014; slash-burning also began in 2012. We used spring of 2012 as the break between before and after treatment. Within the treatment area, 6488 ha of trees were cut and 2277 ha of trees were slash-burned, while 9443 ha and 3540 ha were cut and slash-burned, respectively, in the greater area with an average treatment size of 87 ha (Table 1; Fig 1).

Telemetry and nest data
Sage-grouse females were captured during winter and spring from 2009-2014 in the treatment area and 2010-2014 in the control area using spotlighting techniques [34,35] near leks and wintering habitat with a goal of~40 radio-collared (22-g VHF radio-collars, model #A4060, Advanced Telemetry Systems, Isanti, MN, USA) females at start of nesting in each of the two areas; 93% of females were captured prior to the onset of nest initiation (1 Apr), and 100% well in advance of median nest initiation (29 Apr). The sample of monitored grouse declined throughout the year between trapping seasons due to mortality and transmitter failure, but the sample was then increased each year to our goal of 40 females per area by the beginning of nesting. Capture and handling methods were approved by the University of Idaho Animal Care and Use Committee (Protocol Number: 2012-16) and the Oregon Department of Fish and Wildlife. While transmitters were equipped with 8-hr mortality switches, encounters generally consisted of telemetry locations acquired by approaching to within 30 m of the birds without flushing. For the survival analysis, we tried to encounter each bird 1-2 times per week during February-July and 1-4 times per month during August-January. For the nest survival The greater treatment area included the treatment area as well as the immediate surrounding area (see Fig 1). analysis, we monitored radio-marked females twice per week during the potential nesting seasons from 2011-2014. When a female was observed in the same place on two consecutive locations, she was then observed visually, without flushing, to verify nesting. Nests were subsequently monitored twice per week until incubation was terminated (e.g., hatched, depredated). All nests were included as independent replicates for the analyses even though some females nested in multiple years (n = 33) or re-nested after failure during the same year (n = 19). Although some autocorrelation in these instances likely exists, we believe interannual variation, such as weather patterns, was more likely to increase independence, and thus, including all data was more beneficial than disregarding these potential pseudo-replicates. However, we assessed nest survival structures including within-year nest type (i.e., first nest or renest) to determine and control for the underlying structure.

Female survival BACI analysis
We estimated female sage-grouse survival with the nest survival model in MARK [36,37] using the RMark package [38] as an interface within R 3.1.2 [39]. Because the control area data collection started in the fall of 2010, we used treatment area data also beginning at this time. As treatments generally concluded at the start of the sage-grouse nesting season (March or April), we defined years biologically as April to March. Survival histories were relatively complete from April to July (1-2 encounters per week), but were sparse during the rest of the year (1-2 encounters per month), and we therefore, pooled samples into monthly encounter histories throughout the year. We initially compared monthly, seasonal, and annual survival models with Akaike's information criterion corrected for small sample sizes (AICc) [40] to select a within-year temporal structure to control for survival heterogeneity in the subsequent conifer removal analysis ( Table 2). We assessed 4 seasonal structures including 1) breeding where April-July was separated from the rest of the year, 2) seasonal where 2 intervals were used: April-September and October-March, 3) seasonal where 3 intervals were used: April-July, August-November, and December-March, and 4) seasonal where 4 intervals were used: April-June, July-September, October-December, and January-March.
In an environmental impact analysis, the important parameter is the impact × area (e.g., control or treatment) interaction and the additive effects are unimportant [31]. The impact can be represented as a before-after impact factor or continuous time variable [30]. Because conifer removal accumulated over time, in our analysis we assessed 2 BACI models with different measures of impact including interactions between study area and 1) before-after impact using April 2012 as the split and 2) linear time trend (Table 2). We also assessed a model where the control area had constant survival and the treatment area had an annual trend and a model with a trend among breeding periods (nesting and early brood rearing; April-July) but constant survival in nonbreeding periods. We also assessed models with just year, study area, or before-after to use as null models. We compared all models with AICc and model weights [40].
To assess the magnitude of treatment effect on survival for the best treatment effect trend model, we calculated the annual survival and the standardized difference between the control and treatment areas and calculated the 85% confidence intervals [41]. Additionally, we estimated independent survival trends in the treatment and control areas using the annual trend BACI model.

Nest survival BACI analysis
We estimated daily nest survival with the nest survival model in MARK [36,37] using the RMark package [38] as an interface within R 3.1.2 [39]. For the BACI analysis, we used only nests from 2011 onward because this was the first nest season of data collection in the control area.
We initially compared weekly, biweekly, and monthly nest survival models along with female age (juvenile or adult), nest initiation, continuous time trend within year, and nest type (first or renest) with AICc [40] to select a structure to control for survival heterogeneity in the subsequent conifer removal analysis (Table 3). Similar to the female survival treatment effect models, we assessed BACI models with interactions between study area and 1) before-after impact using 2012 as the split and 2) linear time trend (Table 3). We also assessed a model where the control area had constant survival and the treatment area had an annual trend. We assessed models with just year, study area, or before-after to use as null models. We compared all models with AICc and model weights.
To assess the magnitude of treatment effect on nest survival for the best treatment effect trend model, we calculated the total nest survival (using a 37-day incubation period) [42,43] and the standardized difference between the control and treatment areas and calculated the 85% confidence intervals. Additionally, we estimated independent nest survival trends in the control and treatment areas using the annual trend BACI model.

Projected population growth rates
Nest success and female survival are two of the most important demographic vital rates among sage-grouse [29,44]. As management actions targeting demographics ultimately aim to  influence population growth, we used a model-based approach to assess the influence of treatment effects on asymptotic population growth rate (λ). Specifically, we generated two-stage (yearling and adult) female-based population projection models employing the formulation developed by Taylor et al. [29], using vital rate estimates for the control and treatment areas from the most supported models. Projection matrices were based on a pre-breeding birthpulse census with annual time steps composed of elements for stage-specific fecundity and survival. Because our interests for demographic estimation lay solely in nest and annual female survival, we did not track lower level components for fecundity including nest initiation rate, clutch size, hatching rate, and chick (hatch to 35 d) and juvenile (35 d to following breeding season) survival. For these vital rates, we used mean estimates from a recent rangewide metaanalysis [29]. Accounting for variation among vital rates is critical to analyses of population growth [29,44]. Therefore, we generated beta distributions (using stretched-beta distribution for clutch size) to represent the central tendency and variance of each vital rate. Variance components for vital rates estimated by Taylor et al. [29] include only process variance, which was estimated using mixed effects generalized linear models. Our nest and annual female survival variance components are conflated with both process and sampling variance such that estimated variation in growth rates is likely inflated. However, we believe this method still provides an informative comparison of growth rates between our control and treatment areas. We generated 10,000 population projection matrices for both the control and the treatment distributions and calculated λ as the dominant eigenvalue across projection matrices.

Female survival BACI
We captured 108 and 116 female sage-grouse from 2010-2014 in the control and treatment areas, respectively. Five birds were removed from the analysis because they left the study area or were not encountered after capture resulting in a sample of 107 and 112 in the control and treatment areas, respectively. Of the 219 birds, 141 died during the study and 78 were right censored because the transmitter failed or the study concluded in August 2014. The breeding season model (breeding: April-July, nonbreeding: August-March) was the best temporal structure to estimate survival with a weight of 0.65 (Table 2A) and was used as the base structure in the subsequent conifer removal analysis. All treatment models had higher AICc values than the base model; however, two of these models differed from the best supported model by only a fractional amount (ΔAICc = 0.35, 0.58) indicating some support for treatments improving survival. The treatment breeding trend model had the lowest AICc and therefore had the greatest support (model weight: 0.27), but the model with a trend throughout the year in the treatment area was highly competitive (model weight: 0.24; Table 2B). The two interactive models, annual trend and before-after, had low support (models weights: 0.03 and 0.02) but were within 6 AICc units and should not be discounted.
The survival estimates of the top model showed a positive trend in the breeding season through time in the treatment area (Figs 3A and 2B), but the confidence interval overlapped zero (β = 0.135 [85% CI: -0.019-0.290]). The confidence intervals were large relative to the effect, but the estimated treatment effect in annual survival was a 6.6% (85% CI: -1.6-14.7) increase in the treatment area relative to the control area by 2014 (Fig 3C).
The annual trend interactive model (Trend × Area) also showed a positive trend though time (Figs 3D and 2E), but the confidence interval of the interaction overlapped zero (β = 0.101 [85% CI: -0.121-0.323]). The confidence intervals were large relative to the effect, but the estimated treatment effect in annual survival was a 9.7% (85% CI: -5.4-24.8) increase in the treatment area relative to the control area by 2014 (Fig 3F).

Nest survival BACI
We located 109 and 123 nests from 2011-2014 in the control and treatment areas, respectively. Seven nests were removed because they were abandoned and we suspected observer influence. Of the remaining, 116 nests were successful and 109 were unsuccessful. The null model was the best temporal structure with a model weight of 0.34 (Table 3A) and was used as the base structure in the conifer removal analysis.
The non-treatment time trends, before-after and yearly, had the lowest AICc values (model weight 0.48 and 0.21; Table 3B). The most supported treatment effect model was the trend in the treatment area (model weight: 0.13), while the interactive treatment models also garnered some support (model weight: 0.10 and 0.06; Table 3B).
The nest survival estimates of the best treatment model showed a positive trend through time in the treatment area (Fig 4A), and the confidence interval of the slope did not overlap zero (β = 0.175 [85% CI: 0.049-0.301]). The estimated treatment effect in 37-day nest survival was a 18.8% (85% CI: 6.4-31.2) increase in the treatment area relative to the control area by 2014 (Fig 4B).
The nest survival estimates of the annual trend interactive model (Trend × Area) showed a positive trend through time in the treatment area relative to the control area (Fig 4C), and the confidence interval of the interaction slope slightly overlapped zero (β = 0.218 [85% CI: -0.055-0.490]). The estimated treatment effect in 37-day nest survival was a 22.3% (85% CI: 5.6-39.5) increase in the treatment area relative to the control area by 2014 (Fig 4D).

Population growth rates
For nest survival, we used estimates from the trend model of 0.294 (85% CI: 0.232-0.355) for the control area and 0.482 (85% CI: 0.375-0.590) for the treatment area. Because our model selection approach did not identify female age as an important source of variation, we applied the estimates across both stages (juvenile and adult). Similarly, for female survival, we applied the final year's estimates of the breeding trend model of 0.493 (85% CI: 0.445-0.541) for the Sage-grouse demographics and conifer removal control area and 0.558 (85% CI: 0.492-0.624) for the treatment area. Results from the simulation revealed a mean population growth rate of 0.85 (85% CI: 0.71-1.03) for the control area and 1.10 (85% CI: 0.89-1.39) for the treatment area.

Discussion
Our study is one of the first to link sage-grouse demographic rates with conifer removal treatments. At the landscape scale, we observed support for potentially increased female survival and nest survival, which are two of the most important demographic parameters affecting population growth [29]. Applying demographic estimates to population models suggested that treatments increased the population growth by 25% relative to the control from slightly a declining to an increasing population. Inference from population projection matrices relies on the assumption that unmeasured vital rates (e.g., chick survival) in our study area were similar to those from the rangewide meta-analysis. As our nest and female survival estimates were within the range of the meta-analysis estimates, we suspect that other vital rates did not vary substantially from those found in the literature. Furthermore, if conifer removal benefits the other vital rates as well, we might expect an even greater increase in λ. While further monitoring is needed, our results indicate potential positive effects on sage-grouse population growth which corresponds to the desired outcome of wildlife habitat restoration [27]. Baruch-Mordo et al. [20] found conifer expansion decreased lek occupancy, which is a proxy for population size, and recommended conifer removal near leks, while Coates et al. [19] suggested that decreases in populations are in part due to decreased female survival caused by increased conifer abundance. Sanford et al. [45] observed increased nest and brood success near conifer removal areas, but we provide the first empirical evidence that conifer removal may lead to increases in population-level demographic rates.
We speculate the mechanisms affecting survival in our study include reduced avian predator populations and hunting efficiency as well as increased availability of high quality habitat [32]. Avian predators are important sources of mortality for sage-grouse adults and eggs throughout their range [46] with corvids focusing on eggs [47] and raptors largely preying on adults [48]. Avian predator distribution and abundance are limited by nesting and perching sites [49,50], and eliminating these sites may be beneficial in reducing depredation rates. Additionally, changes in available habitat as a result of conifer removal may also impact survival. Because sage-grouse select seasonal habitats and nest sites in areas susceptible to conifer invasion [19,32], conifer removal in those areas increases the availability of higher quality habitat [32]. Compared to other sagebrush habitat types, increased nest survival has been observed in mountain big sagebrush [32], which is also more susceptible to conifer invasion than low sagebrush or Wyoming big sagebrush habitat [4] implying reductions in nest survival as conifer expansion and infill continues, but also opportunities for effective restoration.
At the landscape scale that we assessed, we observed benefits of conifer removal. However, at a smaller extent, conifer removal may affect individuals differently. Coates et al. [19] documented a potential ecological trap at low conifer canopy cover, while Severson [32] suggested that abrupt edges between woodlands and open sagebrush areas may pose a risk for some individuals. However our study, along with those of Baruch-Mordo et al. [20] and Sandford et al. [45], implies that landscape-scale conifer removal benefits sage-grouse populations as a whole. More research is needed at finer scales to fully assess the effect of conifer removal on individual grouse to better inform management planning and maximize population-level benefits.
Although management of tree invasions has long been suggested for conserving prairie and sage-grouse [28,[50][51][52], few studies have actually assessed those effects [18,53]. Many studies have documented negative effects of woody encroachment on prairie grouse habitat selection [54][55][56][57], but our study represents a major step forward in evaluating the effects of a landscapescale habitat restoration on landscape-scale demographics resulting in population growth. While our results generally indicate positive outcomes of conifer removal on sage-grouse, much remains to be learned. Further monitoring at this site, as well as other sites throughout the sage-grouse range, will be necessary to understand long-term effects of conifer removal on survival at both the landscape scale and to individual birds. Future analyses should include other demographic parameters (e.g., brood survival) and response in long-term population trends from lek counts.
Sage-grouse have become a primary catalyst for large-scale investment in sagebrush habitat conservation [18], and our findings have ramifications for other ecosystem services stemming from conifer removal as well. Holmes et al. [58] found abundance of sagebrush-obligate songbirds of conservation concern, such as Brewer's sparrow (Spizella breweri), increased in response to conifer removal primarily conducted to benefit sage-grouse. At regional scales, 85% of conifer removal conducted through the Sage Grouse Initiative coincided with Brewer's sparrow abundance centers suggesting that considerable benefits for less well-known wildlife species may be accruing through continued investment in conifer removal for sage-grouse [59]. Conifer removal preserves basic sagebrush ecosystem functions by enhancing resilience to wildlife and resistance to cheatgrass (Bromus tectorum) [60,61], extending soil water availability for plant growth [62], and benefiting water storage and delivery [63]. From a socio-economic standpoint, sustainability and profitability of livestock ranching, the predominant land use in the West, also improves with conifer removal by maintaining forage production [64].

Management implications
Conifer removal appears to benefit survival of female sage-grouse and their nests, which could yield increases in population growth through time. Landscape-scale removals may be important in project planning; in our treatment area, conifer treatments totaled 6488 ha with individual cut areas averaging 87 ha and cumulatively covering~20% of our study area. While the scale of restoration necessary to have a positive effect is still unknown, given the ecology of this landscape species, effort should be made to maximize the size of treatments. Additionally, we recommend maximizing the contiguous area of conifer removal and existing treeless habitat as opposed to patchy treatments.
Supporting information S1 Taylor, Ryan Voetsch, and Kate Yates. Comments from Katie Dugger, Eva Strand, Rebecca Taylor, and Kerri Vierling improved this manuscript. We recognize the logistical support of Oregon Department of Fish and Wildlife and Craig Foster in particular for conducting winter aerial telemetry surveys. The Oregon State Police provided pilots and planes for telemetry flights. We also thank the ranchers whose property supported many of the grouse in this study and who graciously allowed us access.