Reversing tree expansion in sagebrush steppe yields population-level bene ﬁ t for imperiled grouse

. Woody plant expansion into shrub and grasslands is a global and vexing ecological problem. In the Great Basin of North America, the expansion of pinyon – juniper ( Pinus spp. – Juniperus spp.) woodlands is threatening the sagebrush ( Artemisia spp.) biome. The Greater Sage-grouse ( Centrocercus urophasianus ; sage-grouse), a sagebrush obligate species, is widespread in the Great Basin and considered an indicator for the condition of sagebrush ecosystems. To assess the population response of sage-grouse to landscape-scale juniper removal, we analyzed a long-term telemetry data set and lek counts with a Bayesian integrated population model in a before-after-control-impact design. Population growth rates ( λ ) in a treatment area (Treatment) with juniper removal and a control area (Control) without juniper removal indicated the two areas generally experienced population increase, decrease, and stability in the same years. However, the difference in λ between study areas indicated a steady increase in the Treatment relative to the Control starting in 2013 (removals initiated in 2012), with differences of 0.13 and 0.11 in 2016 and 2017, respectively. Retrospective sensitivity analysis suggested the dynamics in λ were driven by increases in juvenile, adult, ﬁ rst nest, and yearling survival in the Treatment relative to the Control. These ﬁ ndings demonstrate the effectiveness of targeted conifer removal as a management strategy for conserving sage-grouse populations in sagebrush steppe affected by conifer expansion. Examples of positive, population-level responses to habitat management are exceptionally rare for terrestrial vertebrates, and this study provides promising evidence of active management that can be implemented to aid recovery of an imperiled species and biome.


INTRODUCTION
Native, woody plant expansion is a global phenomenon with impacts often indistinguishable from exotic invasions (Nackley et al. 2017). Expansion of trees into shrub and grassland ecosystems is driven by changes in climate, land use, and subsequent alterations to historic fire regimes (Romme et al. 2009, Staver et al. 2011, Miller et al. 2019. While these expansions occur over many decades, they have impacted biodiversity (Ratajczak et al. 2012), distributions of flora and fauna (Baruch-Mordo et al. 2013, Boggie et al. 2017, and ecosystem function and services (Huxman et al. 2005, Anadón et al. 2014. Despite the consequences and global extent of human-accelerated woody plant expansions, their role in a rapidly changing planet during the Anthropocene remains underappreciated (Nackley et al. 2017).
The Great Basin is considered one of the most imperiled ecoregions in the United States (Noss et al. 1995, Center for Science, Economics and Environment 2002, Chambers and Wisdom 2009). Threats to Great Basin sagebrush (Artemisia spp.) ecosystems primarily include altered fire regimes, invasion of exotic plants such as cheatgrass (Bromus tectorum), and expansion of pinyon-juniper (Pinus spp.-Juniperus spp.) woodlands (Wisdom et al. 2005, Chambers andWisdom 2009, Fig. 1A). No species has come to symbolize the plight of sagebrush ecosystems of the Great Basin more than the Greater Sagegrouse (Centrocercus urophasianus; sage-grouse; Fig. 1B). The sage-grouse is widespread in the Great Basin and a sagebrush obligate, using sagebrush during every phase of its lifecycle (Connelly et al. 2011). Additionally, sage-grouse often occupy distinct seasonal ranges that span a variety of sagebrush communities (Connelly et al. 2011). These characteristics make sagegrouse a useful indicator of broader sagebrush ecosystem quantity and quality. Tracking the loss and fragmentation of sagebrush ecosystems, the sage-grouse occupies 56% of its pre-European settlement distribution (Schroeder et al. 2004). Counts of sage-grouse males indicate declines in 8 of the 11 states where they occur and estimated range-wide declines of 0.83% per year rangewide since 1965 (WAFWA 2015).
Western juniper (J. occidentalis; juniper) is a common, coniferous tree native to California, Idaho, Nevada, Oregon, and Washington that occupies approximately 3.6 mil ha (Azuma et al. 2005. The range of juniper changed significantly during the late Pleistocene and into the Holocene . These changes are largely attributed to long-term variability in temperature, precipitation, and fire on the landscape (Davis 1982, Van Devender et al. 1987, Wigand et al. 1995. Range expansion of juniper woodlands since European settlement has occurred faster than at any other time during the Holocene Wigand 1994, Miller andTausch 2001). Tree ring data and evidence from pollen in lake sediment cores indicate a rapid increase in juniper establishment since European settlement of the Great Basin in the 1870s, coinciding with favorable climate for tree establishment, introduction of unregulated livestock grazing, and, later, active fire suppression (Miller et al. 2019). These factors reduced the frequency of wildfires that historically limited the expansion of juniper (Mehringer 1987. Juniper woodlands with ≥10% canopy cover in eastern Oregon increased over 79% from 1936to 1988(Cowlin et al. 1942, Gedney et al. 1999. Recent remote-sensing analyses show pinyon-juniper woodlands continue to replace shrublands across the Great Basin with 4600 km 2 of new forest (≥10% canopy cover) establishing between 2000 and 2016 at a rate of 0.46% per year . As much as 90% of juniper expansion is occurring in sagebrush ecosystems upon which sage-grouse depend (Davies et al. 2011, Miller et al. 2011. Negative effects on sage-grouse manifest rapidly in the early stages of tree invasion into sagebrush shrublands. Probability of lek (i.e., breeding arena) activity decreases significantly when juniper canopy cover within 1 km of the lek exceeds 4% (Baruch-Mordo et al. 2013). Similar thresholds (≤4% canopy cover) to probability of use have also been documented in recent studies using marked sage-grouse (Severson 2016, Severson et al. 2017a). At these low levels of tree cover, the sagebrush understory is predominantly intact , suggesting mechanisms unrelated to understory vegetation drive changes in sage-grouse space use in areas with limited juniper cover. Encroaching juniper may provide more perch sites and increase the presence and success of sage-grouse nest predators such as the common raven (Corvus corax) and predators of adult sage-grouse such as the golden eagle (Aquila chryseatos; Wolff et al. 1999, Andersson et al. 2009).
Our study built upon Severson et al. (2017b) who assessed the short-term response of female annual survival and nest survival to juniper removals in a treatment area (Treatment) relative to an area without juniper removal (Control) with treatments initiated in 2012. Severson et al. (2017b) observed increased adult survival and nest success during the first three years after juniper removals began. These estimates were incorporated with rangewide vital rates from Taylor et al. (2012) in matrix projection models to estimate annual population growth rates (λ) in the Treatment and Control. While this analysis suggested a promising increase in λ in the Treatment relative to the Control, its inference was limited due to its lack of study area-specific vital rate estimates or incorporation of lek attendance data.
Other studies have also estimated the response of individual vital rates to pinyon-juniper and its removal , Sandford et al. 2017), but a robust analysis of population-level effects of pinyon-juniper removal is currently lacking. Given the ongoing, widespread pinyon-juniper removal across the Great Basin, much of which is implemented to benefit sage-grouse, a rigorous analysis of population-level response of sagegrouse to these habitat management actions is needed. Sage-grouse population growth rates are sensitive to adult female survival, nest success, and chick survival (Taylor et al. 2012, Dahlgren et al. 2016, and management actions that affect these vital rates are most likely to yield population-level benefits (Wisdom and Mills 1997).
Leveraging a landscape-scale experiment in the northern Great Basin, we sought to refine estimates of demographic rates of sage-grouse in treated (Treatment) vs. untreated (Control) landscapes, estimate population growth, and assess contributions of individual vital rates to observed changes in population growth resulting Fig. 2. Treatment and control study areas used in the "before-after-control-impact" framework to assess the effect of juniper removal on sage-grouse. v www.esajournals.org 4 June 2021 v Volume 12(6) v Article e03551 from juniper removal. We predicted that λ would increase in the Treatment relative to the Control driven by increases in adult, nest, and chick survival following juniper removal.

Experimental design
We used a before-after-control-impact (BACI) experimental design to assess the impact of juniper removal on sage-grouse. This study design is commonly used for impact assessments and is considered a "quasi-experiment" because it lacks replication and random allocation of treatments (Stewart-Oaten et al. 1992, Block et al. 2001. Although the inference of quasi-experiments is limited, true experiments at large scales are often logistically and financially unfeasible (Michener 1997). The BACI design of this study followed the description of Block et al. (2001). There was a Control where no juniper removals occurred and a Treatment where juniper was removed (Fig. 2). Geographic features including a large valley dominated by agriculture and human development separate the two study areas and discourage sage-grouse movement between them. Movement between the two areas in a given year was rare (mean = 2.25 movements per year;~3% of marked individuals in a given year), and only two permanent movements from one area to the other were documented from 2010 to 2017. The two areas were similar to each other in habitat characteristics prior to juniper removal, and data were collected at both areas before and after juniper removal. Approximately 13,851 ha of juniper was removed in the~40,000 ha Treatment. Hand cutting was the primary juniper removal technique employed (BLM 2011). When trees were dense, felled individual trees were burned or slash was piled and burned to minimize litter and structures that may serve as perch sites for avian predators (BLM 2011). We considered the sage-grouse populations at the two areas collectively as the population of inference.

Count data
Counts of males at leks were conducted following established protocols (Connelly et al. 2003, Hagen 2011a) collaboratively by Bureau of Land Management, Oregon Department of Fish and Wildlife, University of Idaho, and Oregon State University within the study area boundaries. Counts were completed two to three times at each lek at 7-to 10-d intervals between 15 March and 30 April. Individual counts occurred within the first 2 h after sunrise under calm, mostly clear conditions.

Vital rate data
Female sage-grouse were captured using spotlighting techniques (Giesen et al. 1982, Wakkinen et al. 1992 in winter habitat from 2009 to 2017 and fitted with radio collars (22 g VHF, Advanced Telemetry Systems; Isanti, Minnesota, USA) or rump-mounted GPS backpacks (22 g PTT-100 Solar Argos/GPS PTT; Microwave Telemetry, Inc., Columbia, Maryland, USA; 22 g Solar GPS PTT with 3.5 g Holohil PD-2 VHF transmitter attached; GeoTrak, Inc., Apex, North Carolina, USA) with the goal of marking 40 females in both the Treatment and Control prior to the start of the breeding season (~1 April). Marking with GPS transmitters began in 2015 and by the breeding season of 2017, all females were marked with GPS transmitters. Females marked with VHF transmitters were located twice per week March-July and approximately once per month using aerial telemetry August-February. Locations were collected from females marked with GPS transmitters four to five times per day year round.
Nests of females marked with VHF transmitters were identified by two consecutive locations in the same spot and then visual confirmation of the female on a nest at <30 m without flushing. Nests of females marked with GPS transmitters were identified by ≥3 consecutive points at the same location and subsequent visual confirmation as above. Nests of females marked with VHF transmitters were monitored twice weekly until nest fate was determined [successful (≥1 egg hatched), failed (depredated or abandoned)]. Clutch size was determined by inspection of eggshells after hatch of successful nests. Data from depredated nests were not used for clutch size analysis due to difficulty in counting eggshells and missing eggshells.
Prior to 2015, brood flush counts were conducted for hens with successful nests at 28 and 50 d post-hatch. Marked females were located using radio telemetry near sunrise, and chicks were counted by gridding the area within~30 m of the female's location and flushing all chicks v www.esajournals.org present (Gregg and Crawford 2009). If no chicks were counted, another count was conducted within~2 d to confirm brood fate. We repeated counts occasionally when a female exhibited behavior indicative of having a brood, such as clucking or reluctance to flush, but no chicks were detected. During 2015-2017, flush counts were conducted at 14, 24, 34, 44, and 54 d post-hatch. We truncated the chick survival data set to 34 d due to the reduced detection of older chicks and because the majority of sage-grouse chick mortality occurs prior to 34 d (Johnson and Boyce 1990). The VHF transmitters had 8-h mortality switches, and the fate of mortalities was typically determined within~5 d during the breeding season. Mortalities of individuals marked with GPS transmitters were determined by remote inspection of GPS locations indicating >18 h (~3 locations) without movement and subsequent recovery and identification of fate within~5 d.

N-mixture model
We modeled lek data using open N-mixture models prior to inclusion in an integrated population model (IPM, see Integrated Population Model) to account for imperfect detection and provide improved estimates of male lek attendance (Kéry et al. 2009, McCaffery et al. 2016. The N-mixture model comprised a state and observation process (i.e., state-space model). In this case, the state process was the true, and unobserved male attendance and the observation process accounted for the imperfect detection inherent in count data. We modeled attendance (N) at lek (i) in year (t) with a Poisson distribution: We modeled the log-transformed Poisson intensity similar to Kéry et al. (2009) and McCaffery et al. (2016) and included a random effect of year (α t ) and a random effect of lek (ɛ i ): Given the above state process, we modeled repeated surveys (j) at leks (y i,j,t ) with a binomial distribution, where N i,t is the number of males on a lek in a given year and p i,j,t is the probability of detecting an individual male during the j th survey: We modeled logit-transformed detection probability with an intercept (β 0 ), effect of day of year (β date ), and quadratic effect of day of the year (β date 2 ): The quadratic effect of day of year in the observation process accounted for the rise, peak, and decline in male attendance at leks over the course of the breeding season. Given the limited sample size of leks in the project area, there was significant uncertainty around estimates of detection probability, which subsequently reduced the certainty around abundance estimates. To remedy this, we fit the N-mixture model to data from across the state of Oregon (n = 1162 leks) to more accurately estimate detection probability, and only used estimates of N from leks within the sampling frame in the IPM.
We estimated N-mixture model coefficients in JAGS (Plummer 2003) using R version 3.2.2 (R Core Team 2015) and the R2Jags package (Su and Yajima 2009). We ran two parallel chains for 10,000 iterations, and the first 5000 iterations were discarded. We determined convergence bŷ R < 1.1 and visual inspection of chains (Gelman et al. 2004, Link andBarker 2010).

Integrated population model
We used an IPM to estimate female sagegrouse vital rates and λ (Fig. 3). These models use information from multiple data sources for more precise demographic parameter estimates and can be used to estimate parameters for which no data were collected under some circumstances (Besbeas et al. 2002, Kéry and Schaub 2012. Additionally, IPMs are robust to violations of the assumption of independence between data sets and the improvement of the precision and accuracy of parameter estimates provided by these models is v www.esajournals.org often more pronounced at small sample sizes (Besbeas et al. 2003, Abadi et al. 2010. We combined count data from lek surveys and data from marked females to estimate mean annual vital rates, annual abundance, and λ for the Treatment and Control. We modeled the likelihood of lek attendance estimates from the N-mixture model with a statespace model (Besbeas et al. 2002; see Count data and N-mixture model). The state process can be described using a pre-breeding matrix projection model (Caswell 2001) as follows: where N 1 is the number of yearlings, N 2þ is the number of adults, f 1 is yearling fecundity, f 2þ is adult fecundity, ϕ juv is juvenile survival from hatch to first breeding, ϕ 1 is yearling annual survival, and ϕ 2þ is adult annual survival. Abundance and vital rates were estimated for Treatment or Control areas (h) separately for each year (t), while 0.5 is the proportion of female offspring assuming equal sex ratios at hatch (Atamian and Sedinger 2010). We modeled the state process for yearling (N 1 ) and adult females (N 2 ) as follows: Fecundity was comprised of first nest propensity (π 1 ), or the probability of nest initiation, renest propensity (π 2 ), first nest survival (ϕ nest1 ), renest survival (ϕ nest2 ), first nest clutch size (CS 1 ), renest clutch size (CS 2 ), egg hatchability (ϕ egg ), and chick survival (ϕ chick ) such that: f h,k,t ¼ π 1,k,t ϕ nest1,h,t CS 1,t ϕ egg,h,t ϕ chick,h,t 0:5 þ 1 À ϕ nest1,h,t À Á π 2,k,t ϕ nest2,h,t CS 2,t ϕ egg,h,t ϕ chick,h,t 0:5 where the first and second terms of the equation represent the component fecundities from first nests and renests, respectively, and 1 À ϕ nest,h,t À Á represents the conditional probability of a first nest failing (depredated or abandoned). Subscript h denotes Treatment or Control, and k denotes yearling or adult.
We modeled yearling and adult first (π 1,k,t ) and renest propensities (π 2,k,t ) using binomial distributions as follows: y j,k,t ∼ Binomial N j,k,t , π j,k,t À Á Fig. 3. Directed acyclic graph of a hierarchical, integrated population model for sage-grouse. Estimated parameters are within circles, and data sources are within rectangles. Arrows represent dependencies between graph components. Estimated parameters include fecundity of adults (ϕ 2þ ), fecundity of yearlings ( f 1 ), juvenile survival (ϕ juv ), yearling survival (ϕ 1 ), adult survival (ϕ 2þ ), and lambda as derived from annual population size (λ ⊢ N). Annual estimates of population size from N-mixture models were incorporated as population count data (ŷ).
v www.esajournals.org where y j,k,t is the number of first nests (j = 1) or renests (j = 2) belonging to females of age class k observed in year t, N j,k,t is the number of females of age class k followed through nesting season t in the case of first nest propensity (j = 1), and number of females with failed first nests followed in nesting season t in the case of renest propensity, and π j,k,t is the probability of initiating a first nest or renesting following a failed first nest for females in age class k in year t. We used informative prior distributions for nest propensity, which varied by age and nesting attempt, using estimates of mean and process variance from Taylor  We modeled nest survival to 37 d (combined laying and incubation periods) with a Bernoulli distribution as follows: The nest data consisted of daily encounter histories where y k,j is the nest state (1 for alive, 0 for dead, NA for not checked) for a given k on day j. Daily nest state arises from a Bernoulli distribution of the product of the previous nest state (y k,jÀ1 ) and the probability of surviving that interval (ϕ k,j ). The probability of surviving an interval is a linear function of overall mean survival over a given interval (β 0,ϕ k,j ), the time-varying effects of being in the Treatment (β treatment;t ; Treatment coded as 1) and renests (β renest,t ; renest coded as 1), and the random effect of year (α t ). We derived annual estimates of first and renest survival by area by exponentiating the logit of the linear function with the appropriate variable values to the power of 37.
We modeled first and renest clutch sizes by year (CS nest,t ) with Poisson distributions with the general format: where y i,t is the number of eggs in nest i in year t. Informative priors for annual mean clutch size are averaged across both age classes reported in Taylor et al. (2012). We modeled egg hatchability by area and year (ϕ egg,h,t ) with a binomial distribution as follows: The egg data consisted of number eggs in a given clutch from a successful nest (n k,j ) and number of eggs that hatched (y k,j ). The logit of the probability of an egg hatching in a successful nest is a linear function of overall mean egg hatchability (β 0,ϕ egg ), the time-varying effect of being in the Treatment (β treatment,t ; Treatment coded as 1), and the random effect of year (α t ). Similar to nest survival, we derived annual egg hatchability by computing the logit of the linear function using the appropriate variables.
The data for chick survival consisted of number of hatched eggs in successful nests (n k,j Þ and number of chicks counted at 34 d post-hatch (y k,j ). The model for chick survival by area and year (ϕ chick,h,t ) and the derived annual estimates v www.esajournals.org followed the same structure as the egg hatchability model above: Adult and yearling annual survival had similar model structure to the nest survival model: The survival data consisted of monthly encounter history where y k,j is the state of the marked female (1 for alive, 0 for dead, NA for not checked) for individual k in month j. Monthly survival state arises from a Bernoulli distribution of the product of the state in the previous month (y k,jÀ1 ) and the probability of surviving that interval (ϕ k,j ). The probability of surviving an interval is a linear function of overall mean survival over a given interval (β 0,ϕ k,j ), the time-varying effects of being in the Treatment (β treatment,t ; Treatment coded as 1) and being a yearling (β yearling,t ; yearling coded as 1), and the random effect of year (α, t ). We defined years as April-March to coincide with the reproductive cycle of sage-grouse (Severson et al. 2017b). We derived annual estimates of adult (ϕ 2þ,h,t ) and yearling survival by area (ϕ 1,h,t ) by raising the logit of the linear function with the appropriate variable values to the power of 12.
No data were collected on juvenile survival (ϕ juv ) from 34 d post-hatch to first breeding (~April). However, we assumed a relationship between juvenile and adult survival and derived it annually as follows: Adult annual survival is raised to the 2/3 power to reduce it to the same temporal scale as juvenile survival (~8 months). Adult survival during this shortened period is multiplied by an annual adjustment ratio (γ, h,t ), which has an informative prior distribution from Apa et al. (2017) who reported a mean ratio of juvenile to adult September-March survival for radiomarked adult and juvenile sage-grouse of 0.70 over 3 years at two study sites in Colorado.
The observation component of the state-space model follows Ross et al. (2018): where y, h,t is the population estimate for a given area in a given year from the N-mixture model, N 1,h,t and N 2þ,h,t are the estimates of yearling and adult abundance, respectively, at breeding, and τ h,t is the precision of the estimates from the Nmixture model. We used age ratio-adjusted abundance estimates for 2008 from the N-mixture model as initial values of yearling (N 1,h,1 ) and adult abundance (N 2þ,h,1 ) for the state process. The age ratio is based on the mean productivity estimate of 1.58 from hunter-donated sage-grouse wings (n = 7986) in Oregon from 1993 to 2005 reported by Hagen and Loughin (2008). We implemented the IPM in the Bayesian framework using JAGS (Plummer 2003) in R version 3.2.2 (R Core Team 2015) and the R2Jags package (Su and Yajima 2009). Three parallel chains ran for 10,000 iterations, and the first 5000 iterations were discarded resulting in 15,000 saved iterations. We determined convergence byR < 1.1 and visual inspection of history plots (Gelman et al. 2004, Link andBarker 2010).

Retrospective sensitivity analysis
We estimated the contributions of given vital rates to the observed dynamics in λ by correlating the differences (Treatment-Control) in the annual estimates of a given vital rate with the corresponding differences in λ (Treatment--Control) following Kéry and Schaub, (2012). The strength of these correlations is indicative of the strength of the contribution of annual variation in the difference in vital rates to the annual variation in λ . We used the mode of the posterior distributions of the correlation coefficients to describe these relationships because most of the posterior distributions were skewed . Additionally, we derived the probability that the correlation was positive from the output of the analysis following Kéry and Schaub (2012).

Integrated population model
Over the course of the project, 417 hens were captured and marked, 378 nests were observed, and 223 broods were monitored. Six leks were counted in the Treatment with a mean maximum male count of 16. Three leks were counted in the Control with a mean maximum male count of 34. Annual estimates from the IPM for the Treatment and Control of adult survival ranged 0.50-0.86, yearling survival ranged 0.54-0.91, first nest survival ranged 0.28-0.70, chick survival ranged 0.13-0.50, and juvenile survival ranged 0.49-0.85 (Appendix S1: Table S1). Lambda in the Treatment and Control generally followed the same pattern through time with both areas experiencing population increase (λ > 1), population decrease (λ < 1), and stable populations) (λ = 1) during the same years (Fig. 4A). This pattern indicated that both areas experienced the same climate and the Control dynamics were representative of the Treatment, as is required of a BACI study design (Block et al. 2001). However, the magnitude of differences in λ between the Treatment and Control increased steadily from 2013 to 2016 after juniper removal began resulting in λ that was 0.13 (95% CI: −0.05 to 0.30) higher in the Treatment relative to the Control by 2016 (Fig. 4B). Lambda remained higher in the Treatment in 2017 with a 0.11 difference (95% CI: −0.03 to 0.27; Fig. 4B). Difference in λ exhibited greater stochasticity during the two pre-removal years and the first year of removals Fig. 4B). Differences in λ during this period indicated that the Treatment λ was greater than the Control in 2010 (0.20, 95% CI: 0.04-0.37) and 2012 (0.19, 95% CI: 0.07-0.32) and lower in 2011 (−0.21, 95% CI: −0.33 to −0.10; Fig. 4B).

Retrospective sensitivity analysis
Of the vital rates estimated separately for the Treatment and Control, higher juvenile survival in the Treatment relative to the Control had the greatest contribution (correlation coefficient) to Δλ Based on their posterior distributions, the probability that the contributions were positive was ≥0.95 for juvenile, adult, first nest, yearling, and renest survival (Fig. 5). Given the negative and nonsignificant contribution of difference in chick survival, there was only a probability of 0.16 that its contribution was positive (Fig. 5).

DISCUSSION
Integrating lek and vital rate data in a landscape-scale experiment, we present evidence that sage-grouse population growth rates were up to 13% higher in the Treatment area where juniper removal occurred relative to the Control (Fig. 4). To our knowledge, our study is the first empirical evidence of a population-level benefit for sage-grouse from habitat management. The estimated population response adds to the body of evidence supporting targeted pinyon-juniper removal as an effective tool to restore sagegrouse habitats and increase populations.
While treatment effects on λ were not as high as previously predicted (25%; Severson et al. 2017b), a sustained difference in λ of 11-13%, as estimated in our study, could have a tremendous impact on sage-grouse abundance in the Treatment relative to the Control. These differences in λ occurred 5-6 yr after removals were initiated. Given the lower reproductive output and long v www.esajournals.org lifespan of adult sage-grouse relative to other gallinaceous birds, it may be biologically unfeasible for sage-grouse populations to respond to habitat management as rapidly as other gallinaceous birds. McKenzie (2017) modeled lek data from 1980 to 2015 in relation to juniper removals in a study area that encompassed predominantly eastern Oregon. The findings suggested a 5-to 10-yr lag before an effect of juniper removals is realized in measurable, increased lek counts (McKenzie 2017). Similarly, analysis of time-lag effects from lek data collected at 704 leks over 12 yr in Wyoming indicated a delay of 2-10 yr between energy development and measurable negative effects on lek attendance (Harju et al. 2010). Short-term (2-3 yr post) population assessments may not capture population-level effects of pinyon-juniper removals or other conservation actions for sage-grouse or may underestimate their effects. The lag effects observed in this study and others highlight the importance of long-term data sets to best assess the response of sage-grouse populations to management actions or disturbance.
The sensitivity analysis suggested that juvenile, adult, first nest, and yearling survival had the greatest and similar contributions to the observed changes in population dynamics in the Treatment relative to the Control. This finding supported our initial prediction that increases in adult and nest survival would contribute to increases in λ in the Treatment relative to the Control. Similar findings have been reported from larger data sets and Leslie-matrix style sensitivity analyses (Taylor et al. 2012, Dahlgren et al. 2016. However, these studies also reported chick survival as one of the top contributors to population dynamics among sage-grouse vital rates. Additionally, Sandford et al. (2017) documented increased brood success (≥1 chick at 50 d post-hatch) for female sage-grouse that selected brooding locations that were closer to juniper removals and in areas with minimal juniper cover in Utah. However, contrary to our predictions, chick survival did not have a significant contribution to Δλ and the juniper removals may not have benefited chick survival as much as other vital rates. The primary predators of sagegrouse chicks are largely unknown (Hagen 2011b). If reduction in avian predator perch sites and hunting efficiency is a mechanism for increases in sage-grouse vital rates following pinyon-juniper removal, we speculate that the primary predators of chicks in the project area may not be avian and/or were not affected by the removal. We also speculate that the deleterious effects on chick survival of wet, cold weather events in June 2016 and 2017 may have masked any benefits of juniper removal on this vital rate. The sensitivity of Δλ to differences in juvenile survival was not predicted, but complemented the findings of Prochazka et al. (2017) who documented increases in age-specific risk of daily mortality while moving through pinyon-juniper, with juveniles having the greatest increase in mortality risk (56% increase) compared with yearlings (42%) and adults (16%). However, juvenile survival was estimated as a latent vital rate Fig. 5. Correlation between treatment (juniper removal) effects on individual vital rates and treatment effects on population growth rates. Differences (Treatment-Control) in annual estimates of vital rates for sage-grouse are plotted against annual estimates of differences (Treatment-Control) in population growth rates (λ) between the Treatment and Control study areas, [2010][2011][2012][2013][2014][2015][2016][2017]. Vertical and horizontal lines around points are 95% credible intervals associated with these differences. The mode and 95% credible of posterior distributions for coefficients of correlation (r) are reported in addition to the probability that the correlation is positive [P(r > 0)]. in our analysis and was specified to be closely tied to adult survival, so this result should be interpreted cautiously. Other retrospective analyses of the sensitivity of λ in the project area to the estimated vital rates such as life table response experiments (Caswell 1996(Caswell , 2001 or transient population dynamics (Yearsley 2004) may better elucidate the mechanism behind observed increases in λ.
Although understory dynamics vary across ecological gradients after juniper removal, understory vegetation cover can generally be increased by removing juniper (Roundy et al. 2014). Removal of juniper by hand cutting at a longterm study site in Oregon increased understory herbaceous biomass and cover and generally increased the productivity of these plant communities post-removal (Bates et al. , 2007. Additionally, burning individual cut trees during the winter increased herbaceous and perennial grass cover 150% and 200%, respectively, while increasing the density of perennial grasses and reducing annual grass cover compared with areas that were cut but not burned (Bates and Svejcar 2009). Given the role of sage-grouse as indicators of sagebrush ecosystem function and the positive understory vegetation response often associated with juniper removal, projects directed toward sage-grouse may benefit other species, particularly sagebrush obligates and near obligates , Copeland et al. 2014, Donnelly et al. 2017, Holmes et al. 2017.
Accelerated expansion of native woody plants into adjacent plant communities is a global problem with cascading ecological consequences (Nackley et al. 2017, Stevens et al. 2017. Slowing or reversing ecosystem transitions caused by native woody plant expansion is challenging, complex, and often controversial (Van Auken 2009, Nackley et al. 2017). However, we provide evidence of the potential for spatially targeted tree removal to benefit an imperiled species and shrubland ecosystem in North America, encouraging findings during an epoch defined by unprecedented extirpations and loss of biodiversity (Dirzo et al. 2014). Despite recent efforts to accelerate pinyon-juniper removal in the Great Basin, the scale of current conifer reductions remains relatively small (1.6% of treed area; Reinhardt et al. 2020) and may just be keeping pace with the rate of expansion . Large-scale, cross-boundary efforts to reduce conifer expansion in priority watersheds, similar in scale to that implemented in our study area, provide one of few options managers have for addressing persistent, nonregulatory threats (Chambers et al. 2017) and reversing or slowing sage-grouse population declines. Scaling up this type of active management of native tree expansion may be increasingly required to save conservation-reliant species of grass and shrublands on a changing planet (Scott et al. 2010, Nackley et al. 2017.