Silvicultural treatments, fire, and insect outbreaks are the primary disturbance events currently affecting forests in the Sierra Nevada Mountains of California, a region where plants and wildlife are highly adapted to a frequent-fire disturbance regime that has been suppressed for decades. Although the effects of both fire and silviculture on wildlife have been studied by many, there are few studies that directly compare their long-term effects on wildlife communities. We conducted avian point counts from 2010 to 2019 at 1987 in situ field survey locations across eight national forests and collected fire and silvicultural treatment data from 1987 to 2016, resulting in a 20-year post-disturbance chronosequence. We evaluated two categories of fire severity in comparison to silvicultural management (largely pre-commercial and commercial thinning treatments) as well as undisturbed locations to model their influences on abundances of 71 breeding bird species. More species (48% of the community) reached peak abundance at moderate-high-severity-fire locations than at low-severity fire (8%), silvicultural management (16%), or undisturbed (13%) locations. Total community abundance was highest in undisturbed dense forests as well as in the first few years after silvicultural management and lowest in the first few years after moderate-high-severity fire, then abundance in all types of disturbed habitats was similar by 10 years after disturbance. Even though the total community abundance was relatively low in moderate-high-severity-fire habitats, species diversity was the highest. Moderate-high-severity fire supported a unique portion of the avian community, while low-severity fire and silvicultural management were relatively similar. We conclude that a significant portion of the bird community in the Sierra Nevada region is dependent on moderate-high-severity fire and thus recommend that a prescribed and managed wildfire program that incorporates a variety of fire effects will best maintain biodiversity in this region.

https://www.mdpi.com/1999-4907/12/2/150/htm

Field data were assembled from a bioregional monitoring project designed to monitor trends in upland forest birds inhabiting actively managed national forest lands. Sample locations ranged in elevation from 1003 to 2871 m and latitudes from 35.3906° to 41.2931°. At each survey location, we established two transects in adjacent 1 km grid cells. Tran-sects were made up of four point count stations at 250 m in the cardinal directions from a fifth station in the center. This resulted in a sample of 1987 stations on 398 transects distributed as 199 spatially balanced pairs.We assigned a history of fire and management at each point count station location. To do this we overlaid the points on the Forest Service Activities Tracking System (FACTS) to identify the subsample of stations that had been mechanically treated, as well as a fire severity assessment map, to identify the subsample of stations that have burned. We recorded the year in which the silvicultural activity was reported as completed in the FACTS database, and verified that changes in vegetation occurred by comparing records to time-stamped aerial imagery in Google Earth. Similarly, we identified all stations within a fire perimeter and sampled the vegetation burn severity (assessed as canopy cover % change) at each location and grouped all stations that were burned with less than 25% canopy cover change as low-severity fire and over 25% as moderate-high-severity fire, again verifying all locations identified as burned with aerial imagery in Google Earth. All locations that were mechanically treated (salvage logging or fuels or other vegetation removals) following a fire were removed from analyses. We also removed all data from the year in which a fire or treatment occurred (i.e., time since disturbance cannot be less than 1). All locations with no FACTS or fire disturbance history were included in the undisturbed subsample. Sample sizes varied by treatment category: moderate-high-severity fire = 217 point count stations on 56 transects, 1144 location-year sampling units; low-severity fire = 197 point count stations on 70 transects, 1094 location-year sampling units; mechanical treatment: 275 point count stations on 99 transects, 1584 location-year sampling units; and undisturbed: 1542 point count stations on 350 transects, 11391 location-year sampling units. The sample was dominated by the undisturbed treatment class, a very small portion of which may have been burned or mechanically treated locations that were not captured by the treatment tracking data we used or that we were unable to verify with aerial imagery.

We used standardized five-minute unlimited-distance point count surveys to sample the avian community during the peak of the breeding season. At each survey station, we recorded all birds detected (visually or audibly) and estimated their distance to the nearest 1 m from the observer. We visited each station up to twice between 1 May and 15 July in 2010 through 2017, and in 2019. At each station, we characterized the vegetation within 50 m of the survey plot center using visual estimates of the percentage of the plot that was covered by trees and shrubs (shrub cover includes all understory woody vegetation species), counted standing snags >10 cm diameter, and measured structural characteristics, including the live-tree basal area using a 10-factor key from at least three locations within the survey plot. These relevé vegetation surveys were conducted up to three times at each point count station across the 10-year time span of bird surveys. 

We recorded a total of 148 species over nine years of surveys on these sites, but for these analyses, we removed species such as raptors, waterfowl, nocturnal species, detections not identified to species, and non-breeding migrants for which our point count methodology generated an inappropriate sample. We further removed all species with fewer than 99 individuals detected within 100 m of observers across all survey events. The remaining 78 species were included in hierarchical distance models to estimate abundance at each location-year sampling unit using the distsamp function in the Unmarked package with statistical package R version 3.3.1. Seven of the 78 species had models that fit the data overall but had one or more parameters (described below) with very large standard errors, so we dropped those species from result summaries, leaving a final total of 71 species for which we show results. These 71 species represented 98% of all individuals detected in this dataset. To characterize the responses of a large group of species to an unbalanced sample of a wide variety of vegetation disturbances over a long time period, we took the approach of fitting abundance models to establish covariate relationships and then used simplified summaries of vegetation and environmental conditions to project abundances for each species over time in each disturbance type. To maintain a representative number of sampling units within each yearly time-since-disturbance category, we capped both burned and mechanical treatment subsamples at 20 years since disturbance for model-fitting data since sample sizes were substantially smaller beyond that time period (even though our disturbance history records included up to 30 years for fires and up to 23 years for FACTS). 

For each species, we fit a single set of model covariates (described below). In addition to field measurements of the vegetation structure, we included several other variables to account for geography, climate, and topography. We sampled elevation, aspect, and slope at each station center from the Sierra Nevada Regional Digital Elevation Model. We included three additional variables to account for the influence of a widespread drought throughout the region from 2012 to 2015. Two of these were weather-related variables, and another represented the amount of tree mortality as a result of drought-induced beetle infestations. The climatic water deficit (CWD, a measure of water stress based on evapotranspiration, solar radiation, and air temperature) and a temperature index that characterizes the average June maximum temperature compared to a 2009 baseline were sampled using the California basin characterization model (270 m resolution). Tree mortality was calculated using an index based on the normalized difference wetness index (NDWI) using freely available LANDSAT imagery in Google Earth Engine.

The detection process was modeled using a multinomial function where detection probability varies by location and distance with slope, live-tree basal area, and percentage shrub cover included as covariates in the half-normal (with scale parameter σ) detection function. Abundance was modeled as a Poisson function assuming closure within each year and including an offset for the number of visits. Covariates in the abundance portion of the model included variables to account for the highly variable landscape and topography, namely elevation, quadratic of elevation, latitude, interaction between elevation and latitude, and southness (i.e., aspect represented as a proportion of south-facing. We also included the vegetation structure covariates shrub cover, tree cover, snag density, and a binary variable indicating whether each location was within a ~1980–1995 clearcut patch to account for recent profound silvicultural management that preceded the FACTS data. Three variables to account for recent drought conditions were also incorporated, climatic water deficit, temperature index, mortality, plus an interaction between climatic water deficit and latitude. Finally, the model included variables designed to quantify the influence of the three treatment types (MHSF = moderate-high-severity fire, LSF = low-severity fire, and MT = mechanical treatment) and to allow for the influence of those treatments on each species to vary over time. For these variables, we included the interaction between a treatment-type binary variable (=1 where the treatment or fire occurred and =0 elsewhere) and the number of years since fire (YSF) or mechanical treatment (YSM) occurred. Treatment was modeled as a linear variable proportional to time since treatment, as well as their quadratic effects to account for non-linear associations (e.g., to fit species with the highest or lowest abundances at intermediate values of time, or abundances that plateau at intermediate values). All continuous-scale covariates, including time since disturbance, were standardized (mean = 0.0, standard deviation = 1.0) prior to calculating the squared terms (thus the quadratic terms have high values for both high and low elevations and low values for elevations near the mean). We examined the degree of collinearity between vegetation, topography, and climatic variables using the vif function in the R package HH to calculate variance inflation factors (VIF) and found no evidence of a high degree of collinearity (all VIF <3.0). Since the live-tree basal area was moderately correlated with both tree cover (R = 0.57) and shrub cover (R = 0.38), we only included the live-tree basal area as a detection covariate. 

To evaluate the differences in abundance for each species between characteristic undisturbed, burned, or mechanically treated locations, we used the fitted models to predict abundance within each treatment category across a broad range of time since disturbance (2–19 years) with a set of characteristic covariate values derived from averages across time within each disturbance type. Covariates for these predictions included overall mean values (set to zero for standardized variables) for latitude, elevation, slope, southness, live-tree basal area, mortality, CWD, and temperature index. For the vegetation covariates, we used the mean yearly value for each of tree cover, shrub cover, and snag density within each disturbance type to represent the characteristic vegetation conditions and reflect typical changes over time, and we chose three sets of covariate values that represent three common undisturbed forest conditions: (1) undisturbed dense forest (high tree cover, low shrub cover, low snags); (2) undisturbed open forest (low tree cover, low shrub cover, low snags); and (3) undisturbed shrub-dominant montane chaparral (low tree cover, high shrub cover, low snags).

Data used in the analysis for the following publication: https://www.mdpi.com/1999-4907/12/2/150/htm

This code (FORESTS-1080928_Roberts-etal2021_abundance-prediction.r) should work with any recent versions of R 3.X.X. If there are any problems please contact the author (ljroberts@pointblue.org). 

File: "data formatted for analysis 9.csv"

Bird survey detections and covariates file. Each row in this table is a detection of a single species at a particular distance (in some cases may be more than one individual). Fields in this file include: PointYrID (a unique identifier for each point visited in a given year); Point (unique identifier for each field location); Transect (unique identifier for groups of points arranged in a 5-point transect); Year (calendar year in which the detections were recorded); Site (unique identifier for a pair of two 5-point transects in close proximity); Visit (visit identifier for sequential visit within the same calendar year, 1 and 2 are only possible values); DOY (numerical identifier of day of year on which the visit took place, possible values = 1-365); TOD (numerical format of time of day); Count (number of individuals detected); Spp (species identified, conforming to 4-letter AOS codes); Detection Cue (song=S, visual=V, call=C, drumming=D, hum=H, juvenile=J); Distance (distance in meters at first detection).

File: "Site covs9.csv"

Site covariates file. Fields include: PointYrID (unique identifier for each point visited in each year); Point (unique point identifier); Year (calendar year in which the visit occurred); MaxOfVisit (number of visits in a single calendar year in which bird survey data were gathered); latitude (geographic coordinate WGS84/NAD83); longitude (geographic coordinate WGS84/NAD83); CWHR1 (California Wildlife Habitat Relationship forest type classification); elevation (elevation in meters sampled from DEM); aspect (degrees aspect, 0-360); slope (% of vertical); SRI (Solar Radiation Index); southness (proportion of completely south-facing aspect [180 degrees]); ShrubCov (shrub cover, field estimated % cover); TreeCov (tree cover, field estimated % cover); BA (basal area, square feet per acre); g30SNAGS (count of snags greater than 30cm in diameter), totSNAGS (count of snags greater than 10cm in diameter); MaxDBH (diameter in inches of largest trees within 50m); DisturbanceYear (year in which a fire or management activity took place, "NA" for locations with no disturbance history); TRT (treatment type); Landscape (estimate of whether locations within 10km have been subjected to clearcuts between 1970 and 1995); ndwi (normalized difference wetness index); tmn (average daily minimum temperature in june); tmx (average daily maximum temperature in June); dndwi (difference in normalized difference wetness index from baseline value in 2009); dtmn (difference in average daily minimum temperature in June from baseline value in 2009); dtmx (difference in average daily maximum temperature in June from baseline value in 2009); cwd (climatic water deficit in mm H2O).

File: "species.list.6.csv"

List of species included in this analysis (4-letter AOS codes).