title: Trends in bird abundance differ among protected forests but not bird guilds
authors: Jeffrey W. Doser Aaron S. Weed Elise F. Zipkin Kathryn M. Miller Andrew O. Finley
text: 
  
      
  
          
  
          Trends in bird abundance differ among protected forests but not bird guilds
  
           
  
           Jeffrey W. Doser1, 2, Aaron S. Weed3, Elise F. Zipkin2, 4, Kathryn M. Miller5, 
   Andrew O. Finley1, 2, 6
           
  
      
  
      1Department of Forestry, Michigan State University, East Lansing, MI, 48824, USA 
  
           2Ecology, Evolution, and Behavior Program, Michigan State University, East Lansing, MI, 48824, USA 
  
           3Northeast Temperate Inventory and Monitoring Network, National Park Service, Woodstock, VT 05091, USA 
  
           4Department of Integrative Biology, Michigan State University, East Lansing, MI, 48824, USA 
  
           5Northeast Temperate Inventory and Monitoring Network, National Park Service, Bar Harbor, ME, 04609, USA 
  
           6Department of Geography, Environment, and Spatial Sciences, Michigan State University, East Lansing, MI, 48824, USA 
  
            Corresponding Author: Jeffrey W. Doser, telephone: (585) 683-4170; email: doserjef@msu.edu; ORCID ID: 0000-0002-8950-9895 
  
            Running Head: Trends in forest bird abundance
           
  
  § ABSTRACT
  
   
  Improved monitoring and associated inferential tools to efficiently identify declining bird populations, particularly of rare or sparsely distributed species, is key to informed conservation and management across large spatio-temporal regions. We assess abundance trends for 106 bird species in a network of eight forested national parks located within the northeast U.S.A. from 2006-2019 using a novel hierarchical model. We develop a multi-species, multi-region removal sampling model that shares information across species and parks to enable inference on rare species and sparsely sampled parks and to evaluate the effects of local forest structure. Trends in bird abundance over time varied widely across parks, but species showed similar trends within parks. Three parks (Acadia National Park and Marsh-Billings-Rockefeller and Morristown National Historical Parks (NHP)) decreased in bird abundance across all species, while three parks (Saratoga NHP and Roosevelt-Vanderbilt and Weir-Farm National Historic Sites) increased in abundance. Bird abundance peaked at medium levels of basal area and high levels of percent forest and forest regeneration, with percent forest having the largest effect. Variation in these effects across parks could be a result of differences in forest structural stage and diversity. By sharing information across both communities and parks, our novel hierarchical model enables uncertainty-quantified estimates of abundance across multiple geographical (i.e., network, park) and taxonomic (i.e., community, guild, species) levels over a large spatio-temporal region. We found large variation in abundance trends across parks but not across bird guilds, suggesting that local forest condition might have a broad and consistent effect on the entire bird community within a given park. Research should target the three parks with overall decreasing trends in bird abundance to further identify what specific factors are driving observed declines across the bird community. Understanding how bird communities respond to local forest structure and other stressors (e.g., pest outbreaks, climate change) is crucial for informed and lasting management. 
  
  
  Keywords: avian populations, Bayesian analysis, community dynamics, hierarchical model, imperfect detection, National Park Service, removal sampling, species richness
  
  
  
  
  
  
  § INTRODUCTION
  
  
  Developing sampling designs and efficient statistical methods to monitor trends in species (e.g., ) and communities (e.g., ) is critical to inform management of landscapes, wildlife, and other natural resources <cit.>. Species occupancy and abundance trends have been used extensively over the last decade to identify and prioritize  management tasks needed to protect wildlife (e.g., ). Species distribution models and similar methods are key tools to evaluate changes in populations and communities (e.g., ), an increasingly important task as climate change and habitat loss cause new challenges to conservation. Monitoring of bird species distributions in particular has received much attention in the ecological literature as a result of their popularity, wide development of statistical methods to account for imperfect detection (e.g., ), large-scale public science programs (e.g., Breeding Bird Survey, eBird), and bird sensitivity  to environmental stressors <cit.>. 
  
  While monitoring changes in species distributions over time is a common task, assessments of abundance trends is less common and more difficult, but critical as declines in abundance can lead to reductions in ecosystem services <cit.>. Key studies have shown large declines in bird abundance in the last several decades across extensive spatio-temporal regions in North America <cit.> and Europe <cit.>. Such research, finding large declines in both common and rare bird species, suggests the need for large-scale monitoring of bird abundance.
  
  Forest bird species in particular play an essential role in forest ecosystems primarily via their foraging ecology that results in ecological services such as scavenging carcasses, nutrient cycling, pest management, and seed dispersal <cit.>. Understanding changes in forest bird abundance over time is a key task in maintaining ecological integrity of forest ecosystems. Diversity of forest birds is strongly linked to structural forest composition and complexity <cit.>, but how that complexity and composition affects bird abundance over time is unresolved.  
  
  
  Monitoring of bird abundances occurs primarily through the use of point-count surveys, in which observers record all birds seen at a location for a set amount of time. Sampling and accompanying modeling techniques to address imperfect detection in count data have been widely developed over the last 50 years and include removal sampling <cit.>, distance sampling <cit.>, repeated counts <cit.>, capture-recapture <cit.>, spatial capture-recapture <cit.>, and useful combinations of such methods <cit.>. Many extensions of these modeling techniques exist to meet objectives of different types of monitoring programs, such as to estimate population size of multiple species at a given point in time (e.g., ) and assess trends in species abundance across larger time periods <cit.>. Recently, the multi-species occupancy model <cit.> was extended to a multi-region framework that enables modeling of occupancy and species richness across independent spatial units that comprise a broader network <cit.>. An analogous approach could model abundance across a large network to provide inference on factors driving spatio-temporal changes in the abundances of common and rare species.
  
  In 2006, the Northeast Temperate Inventory and Monitoring Network (NETN) of the National Park Service (NPS) began a volunteer monitoring program to collect data on bird community abundance and composition across upland forest habitats in eight parks. We used these data in a network-wide analysis to characterize spatio-temporal dynamics of forest bird abundance from 2006-2019. We estimated park-specific trends in bird abundance to evaluate whether patterns in the bird community change among protected areas, which may signify where and whether management efforts would be most beneficial. Within each park, we also estimated trends of functionally similar groups of birds (i.e., guilds) and of individual species to determine if management should target the whole community, specific guilds, or individual species. We also sought to understand whether bird abundance is more driven by presence of breeding habitat (i.e., percent of forest within a 1km radius) or quality of the breeding habitat characterized by the amount of live tree basal area and forest regeneration, which have been shown to influence bird species diversity <cit.>. To do this, we extend the basic removal sampling model <cit.> and the multi-species removal sampling model <cit.> to a multi-region, multi-species framework that directly accounts for imperfect detection and shares information across species within each park and across the network of parks. This novel hierarchical approach enables estimation for rare species and estimation of parameters at parks where there is a paucity of data. 
  
  We predicted high variability in abundance trends across parks and species, which we hypothesized would reveal which species and locations should be targeted for management interventions. We also predicted the amount of local forest cover (percent forest within a 1km radius) would have the largest (and positive) effect on forest bird abundance <cit.>. We expected little variation in the direction and magnitude of trends of species within a guild and larger variation among the different guilds <cit.>. For instance, we predicted forest regeneration would have a positive effect on ground and shrub nesting species, as the amount of regeneration is negatively correlated with the amount of deer browsing <cit.>, which has been shown to reduce the occupancy probability of ground-nesting and understory species <cit.>. We expected forest interior obligate and canopy nesting species to show positive relationships with the amount of basal area <cit.>, while shrub nesters and forest ground nesters will show either a negative relationship with basal area or a curvilinear relationship, peaking at intermediate values <cit.>.
  
  
  
  § MATERIALS AND METHODS
  
  
  
  
   §.§ Study Site and Sampling Methods
  
  
  
  The NETN monitored eight parks across the northeastern United States of America (USA) from 2006-2019 (Figure <ref>, ). Table <ref> includes the park acronyms used throughout the remainder of the manuscript, the total amount of forested land in each park, and a summary of data available for parks (as sampling intensity varied across years and parks). The number of point count locations used to monitor forest birds (henceforth called points) ranged across parks from 5 (WEFA) to 51 (ACAD). The NETN established point count locations based on four criteria: 1) points were 200-250m apart to avoid duplicate sampling; 2) points were located at least 50m from forest edges; 3) points were located in the dominant forest type of each park; and 4) at least one forest vegetation sampling plot was located within 50m of at least one point in each dominant forest type. To meet these criteria, the NETN used a systematic grid sampling frame to select points at all parks except ACAD, where the NETN used the Generalized Random Tesselation Stratified algorithm <cit.>. See <cit.> for further details on point selection and protocol methods.
  
  Volunteers performed 10-minute surveys at points located along permanent forested, transects, annually during the breeding season (May-July). At each point, an observer recorded the time of day and species of individual birds they saw or heard within ten, one-minute increments. Observers recorded each individual bird only the first time they detected the bird during the point count. Observers performed point count surveys in the early morning and during adequate weather conditions, when vocalization is most likely. There was large variability in the number of points surveyed in any specific year within parks as a result of volunteer availability (Table <ref>). 
  
  We used the time-interval of detections to estimate species' abundances with a model that incorporated detection probabilities and the removal sampling data <cit.>. We developed a model that shares information across parks and species within a hierarchical framework because data were sparse at some parks (i.e., WEFA, SAGA, SARA) and for rarely detected species. This approach enables detection corrected estimates of abundance at lesser sampled parks and for rare and/or elusive species, and also yields more precise estimates for all species at all parks, analogous to multi-species occupancy models (e.g., ). 
  
  
  
  
  
   §.§ Forest Covariates
  
  
  In addition to bird monitoring stations, the NETN has over 300 forest vegetation monitoring plots implemented across the same eight parks <cit.>. Some of these plots are co-located within 50m of the bird point count survey points, which enables assessment of the relationship between forest structure and bird abundance. For points that did not have an associated forest plot within 50m, we used covariate data from the closest forest vegetation monitoring plot occurring in the same general forest cover type.  
  
  We evaluated the influence of three forest covariates: amount of forest regeneration, amount of live basal area, and percent forest within a 1km radius around the survey points. We used the amount of forest regeneration (stems ha^-1) as a metric to quantify advanced tree regeneration in the forest understory, and computed regeneration as the total observed density of seedlings (diameter at breast height (DBH) < 1 cm and at least 15 cm tall) plus the total observed density of saplings (1 cm ≤ DBH ≤ 10 cm) per plot. We calculated the amount of living basal area (m^2/ha) for a given canopy tree in a plot from the diameter at breast height (DBH) and then summed over all species in the plot. Field crews measured these forest variables once every four years in each plot, resulting in 3-4 measurements of basal area and regeneration for a given point over the study period. We averaged these values so that each point was associated with a single value of basal area and forest regeneration that represented the average of these forest characteristics across the 14-year study period to provide a broad overview of how forest condition influences forest bird abundance. We obtained the percent forest within a 1km radius around the survey location in 2016 using the National Landcover Database <cit.>. We standardized all covariates to have mean zero and a standard deviation of one to facilitate comparison of effect sizes and relative importance. Preliminary analyses suggested linear effects of all three variables and a quadratic effect of basal area on abundance would yield the best-fitting model. 
  
  
  
   §.§ Modeling Framework
  
  
  We developed a hierarchical model (e.g., ) to estimate trends in community abundance and individual species abundance across multiple geographically distinct parks (; see Figure <ref> for overview of modeling framework). To share information across species and effectively model communities of species within each park, we view species-level parameters as arising from common park-level hyperparameters, which enables improved precision on species level estimates <cit.>. We extend this concept to another hierarchical layer to share information across the entire meta-community by drawing park-level hyperparameters from a common network-level distribution <cit.>. The model is thus a multi-level extension of a generalized linear model (GLM) with the following hierarchical levels: 1) detection model to accommodate imperfect detection for each species; 2) process model of the true latent abundance of each species; 3) park level that models species-specific parameters from a common park (i.e., community) level distribution to induce dependence among species in a park and share information across the community; and 4) network-level that views park-level parameters as arising from a common network-level distribution to share information across all parks in the network. While alternative approaches exist to model species jointly via direct dependency of species detections and/or abundances (e.g., ), the hierarchical approach we employ is a mechanistic alternative that enables inference at multiple ecological levels while accounting for imperfect detection. This approach has been shown to provide good model fit, ecological insights, and key information to support management decisions (e.g., ). 
  
  Following the removal sampling protocol, our model uses the time period of first observation to account for imperfect detection of individuals during the survey by estimating a detection probability that is the product of availability and detectability <cit.>.  We used intervals of two-minutes in length to summarize detections from the available data resulting in a total of B = 5 intervals in which an individual bird could be detected during the survey. 
  
  
  We implemented point count surveys using the removal sampling protocol at r = 1, …, 8 geographically distinct parks for i = 1, …, n_r species at j = 1 …, J_r points for each of t = 1, …, T_r years.  We define the vector y_r, i, j, t as the number of individuals of species i encountered during year t at point j within park r in each of the B time intervals. Subsequently, we define y^*_r, i, j, t as equivalent to y_r, i, j, t with an additional value representing the number of individuals that were not detected in any time interval. Following the basic removal sampling model, we view the data y^*_r, i, j, t as arising from a multinomial distribution with cell probabilities π^*_r, i, j, t conditioned on the latent abundance N_r, i, j, t. We then view the latent abundance N_r, i, j, t as a Poisson distributed random variable with mean λ_r, i, j, t. The multinomial observation model thus allocates the total latent number of birds of species i at point j in park r during year t (N_r, i, j , t) across the B removal time periods and the additional category for individuals not observed. Accordingly our models for both the sampling and biological processes, respectively, take the following forms
  
  
      y^*_r, i, j, t   ∼Multinomial(N_r, i, j, t, π^*_r, i, j, t) 
      
          N_r, i, j, t   ∼Poisson(λ_r, i, j, t).
  
  
  Given this formulation, we can use the relationship between the multinomial and Poisson distributions to  directly model y_r, i, j, t as a series of conditionally independent Poisson distributions by analytically deriving the marginal likelihood <cit.>. For each b = 1, … , B removal period, our marginal likelihood is defined as  
  
  
      y_r, i, j, t, b∼Poisson(λ_r, i, j, t·π_r, i, j, t, b),
  
  
  where π_r, i, j, t is equivalent to π^*_r, i, j, t but does not include the probability of failing to detect an individual that was truly present. Using this conditional likelihood, we do not directly obtain estimates of the latent abundance N_r, i, j, t, but we can subsequently draw these values from a Poisson distribution with mean λ_r, i, j, t. 
  
  
  
    §.§.§ Detection Model
  
  
  Under the removal sampling protocol, the cell probability of being observed in the bth time interval is defined as 
  
  
      π_r, i, j, t, b = p_r, i, j, t(1 - p_r, i, j, t)^b - 1,
  
  
  where p_r, i ,j ,t is the probability of an individual of species i being detected in at least one time interval during year t at point j within park r. We incorporated covariate and random effects in the detection probability models as 
  
  
      logit(p_r, i, j, t) = α_0, r, i + α_1, r, i·DAY_r, j, t + α_2, r, i·DAY^2_r, j, t + α_3, r, i·TIME_r, j, t + α_4, r, t,
  
  
  where α_0, r, i is the intercept for species i in park r, and α_1, r, i, α_2, r, i,  α_3, r, i are regression coefficients representing the linear effect of Julian date (DAY_r,j,t), the quadratic effect of Julian date (DAY^2_r, j, t), and the time since sunrise (TIME_r, j, t), respectively, on the detection probability of species i at park r. We included a random year effect stratified by park, α_4, r, t, to account for variability in detection probability across years and parks that follows a normal distribution with mean 0 and variance σ^2_4, p. The species/park specific effects α_0, r, i, α_1, r, i, α_2, r, i, and α_3, r, i follow distributions with common park-specific parameters to enable sharing of information across species within a given park. Specifically, we model α_1, r, i as
  
  
      α_1, r, i∼Normal(μ_p, 1, r, σ^2_p, REG, 1),
  
  
  where μ_p, 1, r represents the mean linear effect of day on detection probability across the entire species community in park r, and σ^2_p, REG, 1 represents the variability in this effect across the different species comprising the community in park r. Models of α_0, r, i, α_2, r, i, and α_3, r, i are defined analogously.  The μ_p, 1, r coefficients provide a straightforward way to compare how covariates influence the bird communities across different parks. Further, μ_p, 1, r follow a common distribution to enable sharing of information across parks, which allows for estimates in sparsely sampled parks (i.e., WEFA, SAGA, SARA). To do this, each μ_p, 1, r is modeled as 
  
  
      μ_p, 1, r∼Normal(μ̅_p, 1, σ^2_p, META, 1),
  
  
  where μ̅_p, 1 is the mean linear effect of day on detection probability across all species and parks (i.e., the meta-community), and σ^2_p, META, 1 represents the variability in the linear effect of day on detection probability across parks within the meta-community. Models of μ_p, 0, r, μ_p, 2, r, and μ_p, 3, r are defined analogously.
  
  
  
    §.§.§ Abundance Model
  
  
  As with the detection model, we incorporated random effects and covariates at multiple levels of the hierarchy to explain variation in λ_r, i, j, t. We model λ_r, i, j, t using the following specification
  
  
      log(λ_r, i, j, t) =    β_0, r, i + β_1, r, i·YEAR_r, t + β_2, r, i·REGEN_r, j + β_3, r, i·FOR_r, j + 
         β_4, r, i·BA_r, j + β_5, r, i·BA^2_r, j,
  
  
  where β_0, r, i is the abundance intercept for species i in park r, and β_1, r, i, β_2, r, i, β_3, r, i, β_4, r, i, and β_5, r, i are regression coefficients representing the linear effect of year (YEAR_r, t), regeneration (REGEN_r, j), local forest cover (FOR_r, j), basal area (BA_r, j), and quadratic effect of basal area (BA^2_r, j), respectively, on the abundance of species i at park r. We define all species and park-specific parameters as random effects that follow a common distribution within each park. Park-specific parameters subsequently follow a common meta-community distribution. For β_1, r, i we have  
  
  
      β_1, r, i∼Normal(μ_λ, 1, r, σ^2_λ, REG, 1) 
         μ_λ, 1, r∼Normal(μ̅_λ, 1, σ^2_λ, META, 1),
  
  
  where μ_λ, 1, r is the mean year effect across the entire community in park r, σ^2_λ, REG, 1 is the variability of the year effect across all species in the meta-community, μ̅_λ, 1 is the mean year effect across the entire meta-community, and σ^2_λ, META, 1 is the variability of the park year effects across the entire meta-community. Models of β_0, r, i, β_2, r, i, β_3, r, i, β_4, r, i, and β_5, r, i are defined analogously. 
  
  
  
   §.§ Bird Guilds
  
  
  To assess the effects of local forest structure on species' abundances, we assigned birds to behavioral and physiological response guilds following <cit.>. Bird guilds are groups of species that respond in similar ways to environmental changes as a result of similar uses of the environment <cit.>. Overall, we categorized 16 guilds in three biotic integrity elements (Functional, Compositional, Structural) as specialist or generalist guilds depending on each guild's relationship to the landscape structure and function. Specialist guilds consist of species with highly specific habitat requirements, and generalist guilds consist of species that can utilize a wide range of habitats. We selected guilds to indicate different aspects of species' life history traits, which can lead to a single species being assigned to multiple guilds (; Table <ref>). Based on the recommendations in <cit.>, we only used bird guilds in post-hoc assessment.
  
  
  
  
   §.§ Model Estimation
  
  
  We implemented the model using a Bayesian framework (e.g., ). The Bayesian approach is particularly useful in our setting because it easily accommodates missing data and facilitates inference about derived quantities such as covariate effects across different bird guilds with propagated uncertainty across all hierarchical levels (e.g., species and parks). We specified vague normal priors for regression coefficients and vague gamma priors for variance parameters. We fit our models using Markov chain Monte Carlo (MCMC) in JAGS <cit.> within the  statistical environment <cit.> using the  <cit.> package. We ran models for 50,000 iterations with a burn-in period of 45,000 iterations and a thinning rate of 2. We assessed model convergence using visual assessment of trace plots and the Gelman-Rubin R-hat diagnostic, where convergence was considered to occur for all values of R < 1.1 <cit.>. We assessed model fit using a Bayesian p-value approach <cit.>, which indicated a successful model fit (mean = 0.21). We interpret a parameter as “significant” if the 95% credible interval does not include zero. We performed all subsequent analysis in  using the  package <cit.>. All code and data are available on Zenodo <cit.>.
  
  
  
  § RESULTS
  
  
  We observed 106 bird species from 2006-2019 across the eight parks, with an average of 59 species observed at an individual park (Table <ref>). SARA had the highest observed species richness with a total of 71 observed species, while WEFA had the lowest richness, with a total of 45 observed species. Most species in each park's community were rare, as the average point-level observed abundance of a species across all parks was 0.17, with average point-level abundance highest at MORR (0.22) and lowest in ACAD (0.11). 
  
  
  
   §.§ Temporal trends in bird abundance
  
  
  There was wide variation in species abundance trends among the eight parks with three parks (MABI, ACAD, MORR) showing significant declines in overall abundance (across all species), two parks (SAGA, MIMA) showing no trend, and three parks (ROVA, SARA, WEFA) showing significant increases in abundance over the time frame of the study (Figure <ref>b). Park-level species richness showed similar trends to that of abundance (Figure  <ref>). Interestingly, the wide variation across parks resulted in no temporal trend in total bird abundance across the entire network (Figure <ref>a), revealing the benefit of the hierarchical modeling approach to evaluating park-specific dynamics that might be masked at the overall network level (i.e., Simpson's paradox). To assess if the variability across spatial units was a result of bird community composition, we computed the Bird Community Index (BCI; ), a measure of bird communities that provides inference on the ecological integrity of a site (See Appendix S1 for details of BCI computation). We performed a correlation analysis between the estimated year trend and the mean estimated BCI, which revealed a strong negative correlation (median = -0.77).
  
  Individual species abundance trends were highly variable across parks, but species within a given park tended to show similar trends (Figures <ref>c, Appendix S1: Figure S1). Further, species in the same guild showed similar trends in abundance within a given park. However, the direction and magnitude of these trends varied widely across parks (Figure <ref>). 
  
  
  
  
   §.§ Forest structure effects
  
  
  At the network level, we found moderate support for a peak in bird abundance at intermediate values of basal area (85% and 78% probability of significant linear and quadratic effects, respectively), high amounts of local forest cover (92% probability of a positive effect), and to a lesser extent high amounts of forest regeneration (70% probability of a positive effect). Most parks followed the network level pattern, although there was variation among parks, especially in the effects of percent forest and forest regeneration (Figure <ref>). Percent forest has a more prominent effect than either basal area or forest regeneration, as evidenced by comparison of the magnitudes of the standardized regression coefficients. There was little variability in the effects of the covariates across different bird guilds (Appendix S1: Table S1).  
  
  
  
  
  § DISCUSSION
  
  
  With the growing threats of habitat loss, invasive species, and climate change, large-scale monitoring networks are becoming increasingly important for monitoring trends in species distributions and abundance across space and time. Such large-scale programs are often limited by resources at certain locations, indicating a need for sophisticated modeling approaches to estimate species trends across the entire region of interest. Our novel hierarchical model reveals large variations in bird abundance trends across eight protected forests but not across bird guilds within a park, suggesting that ecological processes, biological invasions, and management activities that affect local forest condition appear to have consistent effects on local forest bird communities. An understanding of how these variables influence diverse communities of bird species is crucial to informed and lasting management solutions for both forests and birds. 
  
  Trends in bird abundance differed across space (i.e., parks) but not by species guilds (Figure <ref>), suggesting that local forest conditions might have broad and consistent effects on bird communities within parks. The consistent trends in species across guilds within a park are surprising given the variable responses bird guilds show towards insect disturbances <cit.>, silviculture <cit.>, and elevation <cit.>. However, life-history characteristics, including diet, foraging strategy, habitat preference, and nesting location, do not predict the effects of climate change on bird species distributions in northeastern North America <cit.>, suggesting responses of bird species to stressors such as climate change might be independent of life-history traits. The lack of variability we found in trends across guilds is partially attributable to species being assigned to multiple guilds. However, this cannot explain the variability in trends of guilds across the different parks (Figure <ref>). Thus, our results suggest that diverse communities of forest birds might show similar responses to variations in local forest condition and structure, as well as other spatially-dependent stressors, such as climate change <cit.>. 
  
  We identified three parks (MABI, ACAD, MORR) with significant negative declines in bird abundance, two parks (MIMA, SAGA) with no trend, and three  parks (ROVA, SARA, WEFA) with significant, positive trends (Figure <ref>). Understanding the causes of variability in bird abundance trends across parks and why species within parks behave fairly consistently is critical. Our findings correspond with previous studies showing variability in bird abundance and distribution trends across large spatial regions in the U.S. <cit.>. The strong negative correlation between the BCI and the estimated year trend (median = -0.77) suggests bird communities reflective of higher ecological integrity are showing the fastest declines over time (Appendix S1: Figure S2). This result potentially suggests communities of birds within these parks might respond differently over time as a result of differences in local environmental stressors and interactions with other species <cit.>, a phenomenon commonly addressed in species distribution models <cit.>. 
  
  While overall declines of bird abundance in MABI, ACAD, and MORR are concerning, it is important to emphasize that declines in each of these guilds should be viewed differently. Declines in specialist guilds (e.g., interior forest obligates, canopy nesters) are of highest concern, as these species are indicative of bird communities with high ecological integrity <cit.>. Declines in generalist guilds (e.g., nest predators/brood parasites, forest generalists) are of lesser concern, and such declines can actually lead to an increase in the ecological integrity of the bird community as measured by the BCI. Thus, future analyses and management efforts should focus on declining specialist guilds in these three parks that require the interior and older forest habitat these parks are designed to protect. 
  
  The amount of local forest cover has the largest effect on bird abundance across the parks (Table <ref>, Figure <ref>). Four parks (MORR, MIMA, ROVA, SARA) showed significant positive relationships between bird abundance and percent forest, which is consistent with our hypothesis and previous findings <cit.>. However, ACAD showed a significant negative relationship and SAGA showed a non-significant negative relationship, indicating higher abundance of birds at points surrounded by lower amounts of forest cover. The non-significant results at SAGA are likely a result of the low amount of variation in local forest cover (range of 49-70%) near survey points. However, local forest cover at ACAD varied considerably among points (51% to 97%). The forests in ACAD and the surrounding forest matrix show distinct characteristics compared to the other seven parks in terms of their structural stage, density of large (>30 cm DBH) trees, and tree species diversity <cit.>—underscoring that local forest cover alone does not account for all forest characteristics potentially important to bird breeding and foraging ecology. Further, it is important to note that while abundance in ACAD declines with forest cover, these sites were surrounded by relatively high forest cover, and even the points with the lowest forest cover are still surrounded by substantial forested habitat that could be enough habitat to attract and maintain many forest dwelling species <cit.>. Overall, these results suggest management should focus on limiting forest fragmentation and maintaining or increasing the amount of forested cover in the surrounding landscape matrix, followed by maintenance of forest structure and diversity. However, the inconsistent effect of local forest cover across parks suggests that forest breeding birds are likely affected by interactions between local forest structure, surrounding land use, local community interactions, climate, and local stand dynamics (e.g., pest outbreaks, disturbances, succession). 
  
  Our finding that overall abundance peaked at intermediate levels of live tree basal area in five park forests is consistent with previous research <cit.>. This effect is likely related to larger variety in vertical forest structure at intermediate levels that provides suitable habitat to a wider variety of species <cit.>. More broadly, the peaking of abundance at medium levels of live basal area could be reflective of multi-aged stands supporting a variety of species, individual tree size, and vertical layering that provides habitat to a wide range of birds.
  
  Forest regeneration showed the weakest effect of the three covariates, with only a small positive relationship on overall bird abundance and significant at only one park (ROVA). Multiple parks had highly positive-skewed distributions of forest regeneration, which could indicate that there was not enough regeneration to attract birds that utilize such habitat. Given the legacy of deer overabundance, invasive species, and regional declines in regeneration that have shaped the understory in these parks <cit.>, more attention to examining the effects of regeneration and invasive plant abundance on forest birds seems warranted. Additionally, forest management aimed at increasing regeneration abundance and diversity within these forests should also be evaluated for its benefits to the overall bird community. 
  
  While there was support for effects of tree basal area, forest regeneration, and local forest cover on bird abundance, these variables alone cannot explain the consistent trends of species within parks. ACAD, MABI, and MORR on average had the most forest within a 1km radius of each point. However, these forests vary considerably in structural stage <cit.>, type, and diversity <cit.>, suggesting that declines in bird abundance might be driven by factors other than forest structure. MABI is the only park in the network that is subject to logging, which has led to decreases in basal area, crown closure, and tree density and increases in regeneration over the study period that could contribute to the strong decreasing trends in bird abundance. MORR has experienced a decrease in crown closure as a result of intense storms and the invasive Emerald Ash Borer (Agrilus planipennis), leading to increased gaps in the forests that could contribute to the decreasing trend in overall bird abundance. Reasons for the declines in ACAD are more ambiguous, but they could partially be attributed to declines in early successional habitat throughout the park. Overall, these variations in trends across spatial units are likely a composite result of numerous effects such as agricultural intensification and land use change outside of the parks <cit.>, variations in forest age, structure, and size <cit.>, differing levels of anthropogenic mortality, loss of sensitive areas or resources, migratory behavior, and climate change effects occurring during both the breeding and nonbreeding seasons <cit.>. 
  
  Because the point counts did not contain specific information on the spatial location of each individual bird, we chose to describe population sizes in terms of point-level abundances rather than a measure of density, as the true effective sampling area of the point count is not well-defined due to movement and temporary immigration/emigration of birds <cit.>. This inability to accurately describe the sampling area is an inherent problem of any non-spatial method of counting birds (e.g., repeated counts, removal sampling) and is only remedied via approaches that incorporate spatial information of each individual, such as distance sampling <cit.> and spatial capture-recapture <cit.>. However, inference on point-level abundances, as used in the current study, provides information on population dynamics in time and space, which can yield management-relevant insights <cit.>. 
  
  We developed a novel hierarchical model that enables inference on abundance for all 106 species observed across the eight parks and provides broad inference on how forest conditions affect bird species. We extended previous work on community abundance models <cit.> to efficiently model communities of species across multiple geographically distinct regions. Our flexible modeling approach is widely applicable to monitoring programs where inference is desired on abundance of rare species and abundance of species at sparsely sampled locations. Further, our modeling approach can be easily extended to other sampling protocols commonly used to estimate abundance (e.g., distance sampling). Many of the 106 species in our data set are rare (i.e., lots of zeros in data), which leads to a trade-off between the number of species for which we can estimate trends and the number of covariates we can include in the analysis. This, in addition to the temporal misalignment of the forest covariates, limited our ability to use time-varying covariates to relate variability in forest structure over time to bird abundance trends. The addition of variables reflecting changes in species or vertical structure diversity (e.g., Shannon diversity, Gini diversity), as well as other variables related to climate change could likely provide additional insight on our results. Future work might also seek to combine additional data sources (e.g., Breeding Bird Survey, acoustic recordings, survival/productivity data) with our data using an integrated modeling approach (e.g., ) to help address this trade off, provide further inference on what is driving these trends in specialist guilds indicative of high integrity bird communities, and provide a more mechanistic approach at assessing population declines <cit.>.
  
  
  
  
  § ACKNOWLEDGEMENTS
  
  
  JWD, ASW, EFZ, and AOF conceived the modeling approach; ASW and KMM managed and organized the data collection; JWD analyzed the data and led writing of the manuscript. All authors  contributed critically to the drafts and gave final approval for publication. We declare that there are no conflicts of interest. We are grateful to Camilla Seirup and the field crews that contributed to the collection of NETN forest data in this study, to the volunteers, Steve Faccio, and Ed Sharron for their contributions to the NETN bird monitoring program, and to Adam Kozlowski and Brian Mitchell for their contributions to both programs. We thank John Marzluff, Marc Kéry, and an anonymous reviewer for insightful comments that greatly improved the manuscript. This work was supported by the National Park Service's Inventory and Monitoring Division and National Science Foundation grants DMS-1916395, EF-1253225, and DBI-1954406.
  
  
  
  § OPEN RESEARCH
  
  
  Data and code <cit.> are available on Zenodo at: <https://doi.org/10.5281/zenodo.4701477>.
  
  
  apa
  
  
  
  
  
  
  § TABLES
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § FIGURE LEGENDS
  
  
   Figure <ref>: Location of NETN parks participating in the breeding landbird monitoring program in the northeastern United States. See Table <ref> for description of park codes.
  
  Figure <ref>: Overview of hierarchical model components. Observations y_r, i, j, t are obtained for species i in park r at point j in year t using traditional point count surveys. Species level detection (α_r, i) and abundance (β_r, i) parameters follow common park-level distributions with park level coefficients (μ_p, r, μ_λ, r), which subsequently follow common network level distributions with network level coefficients μ̅_p and μ̅_λ, respectively.
  
  Figure <ref>: Trends in bird abundance from 2006-2019. Panel (a) shows the median linear trend of year across the entire network with the 95% credible interval in parentheses. The trend estimates represent the average change in abundance on the log scale of a species per point every four years. Panel (b) shows the park level trends. Red highlight indicates a significant negative year trend, white highlight indicates no significant trend, and blue highlight indicates positive significant trend. Panel (c) shows the number of species with median linear trend of year estimates that are negative (red) and positive (blue) within each park. The number of species with significant trends is shown in boldface in parentheses.
  
  Figure <ref>: Trends in estimated park-level species richness from 2006-2019 at each of the eight parks. Points are the mean species richness across all points sampled in each given park, gray regions denote the 95% credible intervals. Inset text is the median (95% credible interval) linear trend estimate for a post-hoc analysis between park-level species richness and year, representing the estimated change per year in species richness per park.
  
  Figure <ref>: Mean linear year effect at each park averaged across all species defined within 16 bird guilds. * denotes significance. Gray color indicates no species in that guild were detected.
  
  Figure <ref>: Relationship between mean expected number of birds of a single species per point and basal area (a), percent forest in a 1km radius (b), and forest regeneration (c). Abundance estimates in each graph are estimated at the median levels of all other covariates within a given park.
  
  
  
  
  
  
  § FIGURES
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  