Report on the use of passive acoustic monitoring for analysis of bird trends in Jasper National Park
This report is dynamically generated, meaning its results may evolve with the addition of new data or further analyses. For the most recent updates, refer to the publication date and feel free to reach out to the authors.
Abstract
Since 2007, Jasper National Park has conducted passive acoustic monitoring as part of its ecological integrity monitoring program. The 18 years of data were analyzed to identify trends and extract insights that inform ongoing monitoring and strengthen future species monitoring practices. The analysis assessed whether species and guild abundances shifted by ±2.5% in the alpine and montane ecoregions. Data were managed and processed in WildTrax, combining and harmonizing legacy datasets from multiple methodologies. Sampling locations were tested for independence, and trend analyses quantified changes in species counts over time across guilds and ecoregions.
Land Acknowledgement
We respectfully acknowledge that Jasper National Park is located in Treaty 6 and 8 as well as the traditional lands of the Anishinabe, Aseniwuche Winewak, Dene-zaa, Nêhiyawak, Secwépemc, Stoney Nakoda, Mountain Métis and Métis. We acknowledge the past, present, and future generations of these nations who continue to steward the land.
Introduction
Human activities have been identified as key pressures and contributors to the global decline in forest wildlife (Allan et al. (2017)). The repercussions of habitat fragmentation (Fahrig (2003)) and loss (Hanski (2011)), climate change (Mantyka-pringle, Martin, and Rhodes (2012), Sattar et al. (2021), Abrahms et al. (2023)), and increased access to sensitive areas exert direct and indirect pressures on forest biodiversity, particularly in protected areas in Canada (Lemieux et al. (2011)). In Alberta’s Rocky Mountain Natural Region, these pressures are compounded by increasing wildfire activity, which has significantly impacted montane bird monitoring. Studies have shown large declines in bird abundance in the last several decades across large regions of North America (Rosenberg et al. (2019)), suggesting the need for parkwide monitoring of bird abundance.
In 2007, Jasper National Park initiated a program incorporating passive acoustic monitoring of the park’s vocalizing wildlife. Autonomous Recording Units (ARUs) are compact environmental sensors that are designed to passively record the environment (Shonfield and Bayne (2017)), capturing vocalizing species like birds and amphibians—a method growing in use across the globe (Sugai et al. (2018)). This technology enables biologists to conduct prolonged surveys with minimal human interference, while creating a permanent, archiveable recording of the soundscape. Given the rapidity and ease of accumulating data from these units, maintaining a high standard of data integrity is paramount to ensure future data interoperability and sharing. The data collected by these units contribute valuable information to ecological integrity metrics such as species richness, diversity, occupancy, abundance trends of species and human activities in National Parks over time.
This project analyzes Jasper’s passive acoustic monitoring data from 2007 to 2025, to assess trends while accounting for the clustering of survey points within transects and determining time-to-first-detection usability in analysis. We selected species abundance as our primary monitoring measure following recommendations in Fletcher et al. (2025). Measures that reflect abundance, reproductive performance, survival, dispersal, habitat quality, and role in ecosystem functioning provide richer insights for the achievement of conservation goals than species richness or species incidence. For example, incidence alone does not mean that a species’ presence is meaningful from a conservation perspective because a single observation is treated the same as an abundant species. As changes in abundance are more closely linked to extinction risk than species incidence, abundance-based metrics have been promoted to better assess goals in biodiversity conservation. Abundance can fluctuate more than occurrence and may therefore better track key changes in the environment, although care should be taken to not overinterpret “noisy” abundance variation. Information on abundance is essential for diagnosing decline and recovery of species and may better capture the resilience of biodiversity to ongoing environmental threats
Jasper National Park reports on the condition of three park ecosystems in its State of the Park Report: freshwater, alpine and forest (montane). Since birds are indicators in both the forest and alpine ecosystems, separate analyses were conducted for the montane ecoregion (grouping of the upper subalpine, lower subalpine, and true montane ecoregions) and alpine ecoregion. To enhance accessibility and reproducibility, the findings were presented in this online report with fully documented code, allowing future updates as data collection methods become standardized. Additionally, recommendations were developed to refine data transcription priorities, improve annual reporting methods, and evaluate species guild classifications for long-term monitoring.
This project aims to analyze Jasper’s passive acoustic monitoring data from 2007 to 2025, assessing trends in species and guild abundance while accounting for the clustering of survey points within transects. Separate analyses will be conducted for montane and sub-alpine and alpine ecoregions to align with Ecological Integrity reporting requirements, determining time-to-first-detection usability in analysis. To enhance accessibility and reproducibility, the findings will be presented in this online report with fully documented code, allowing future updates as data collection methods become standardized. Additionally, recommendations will be developed to refine data transcription priorities, improve annual reporting methods, and evaluate species guild classifications for long-term monitoring. The objectives of this report are to:
- Describe the data management and processing procedures for acoustic data collected from 2007 to 2025.
- Compare how different data processing methods influence the count of species and individuals detected on recordings.
- Report on trends in bird abundance by species and by guild in the montane (forest indicator) and alpine (alpine indicator) ecoregions, incorporating time-to-first-detection analysis
- Summarize results by ecoregion as the percentage of species and guilds declining by ≥2.5% in both forest and alpine regions. Explore methods (e.g., inclusion of covariates and summarizing by guild) to interpret these changes, understanding how birds may be responding to stressors in the park or across their ranges and how changing populations may affect ecosystems.
- Estimate species richness trends to determine if they provides insight about bird community conditions.
- Provide recommendations for:
- Prioritizing previous years’ data for re-transcription to 1SPT;
- Optimizing annual reporting methods (e.g., baseline comparisons vs. 10-year trends); and
- Refining methods for evaluating species trends against thresholds and reviewing guilds and traits used in assessments
- Facilitate data publication to the public, resource managers, academic institutions, and any other relevant agencies
Methods
Bird surveys
Monitoring for forest and alpine birds in Jasper National Park has been conducted annually since 2007 across a core set of 129 sites, and 20 sites that were added in 2025. To ensure long-term feasibility, sites were clustered along transects, and transects were randomly selected from hiking trails that were at least 500 m away from roads and were at least five kilometers apart. To sample the diversity of birds in the park, transects were stratified by ecological region (Figure 1).
Songbird data were collected using autonomous recording units (ARUs) deployed by field staff to capture one 10-minute recording per point count annually. Surveys were scheduled consistently each breeding season in June and early July, starting 30 minutes before sunrise. Technicians walked transects containing nine or ten points, each spaced at least 300 m apart to prevent duplicate detections and ensure independence of locations. At each documented location, the ARU was set up, and technicians moved 10–20 m away to minimize disturbance, allowing at least 11 minutes of recording to capture voice notes and the standardized 10-minute survey period. Four covariates thought to influence detection or abundance were recorded: start time, wind speed (km/h), cloud cover (clear, broken, scattered, overcast), temperature (°C) and recording equipment. Recording equipment varied throughout the course of the study; during the 2021 and 2022 field seasons, paired recordings were collected concurrently at each site using both Wildlife Acoustics SM4 autonomous recording units and the original recording devices (RiverForks CZM or Zoom H2N Pro, depending on site-specific deployment). This dual-recording approach was implemented to support calibration and continuity across equipment types. However, to maintain consistency in data processing, only recordings from the original devices were transcribed and included in the final analyses for those years. Beginning in 2023, recordings were collected exclusively using SM4 units.
The 2012 field season was reduced significantly and was not included in some of the subsequent analyses. Refer to Table 1 for more information.
Code
datatable(janp_locs_table,
options = list(
searching = TRUE,
paging = TRUE,
pageLength = 10
)) |>
formatStyle(columns = colnames(janp_locs_table),
backgroundColor = styleEqual(c("NA"), "lightgray")) Habitat and survey variables
We used the existing Vegetation Resource Inventory (VRI) derived from 2020 orthophoto imagery, before the 2022 Chetamon and 2024 Jasper Wildfire to describe habitat conditions. We extracted variables describing dominant land cover (including conifer, broadleaf, shrub, herbaceous, and exposed soil classes), forest structure (canopy height, age), and understory characteristics (shrub height, crown closure). To account for the severe habitat alteration caused by the 2024 Jasper Wildfire, we manually reclassified the land cover variable for all sites on the TEKARRA and VALLEY5 transects as “Burned” for the 2025 data. This classification replaced the pre-fire VRI categories (e.g., conifer, shrub) in the analysis.
Weather variables
To examine the influence of weather on species richness and abundance, we fit a negative binomial generalized linear model (glm.nb()) with richness as the response and temperature, sky (clear, overcast, scattered), wind, Julian day, and hour as predictors. Model dispersion, calculated as the ratio of residual deviance to residual degrees of freedom, was 1.02, indicating minor overdispersion and suggesting an adequate fit. Julian day (β = –0.014, p < 0.001) and hour (β = –0.024, p = 0.002) were significantly negatively associated with richness, while temperature, wind, and sky type showed weak or non-significant effects. Although temperature and wind were not significant in this model, they were retained in subsequent analyses due to their known ecological influence on bat activity and detection probability. The model residual deviance was 1527.5 on 1494 degrees of freedom, with an AIC of 7684.5 and an estimated dispersion parameter (θ) of 40.77 (SE = 8.59).
Data management, processing and quality control
Recordings were clipped and organized to only include the 10-minute count. Before adopting WildTrax in 2021, processing analysts excluded the initial 20 seconds to 1.5 minutes of recordings to reduce human impact on detection probability, then logged the first detection time per species. Recordings are now uploaded as clean 10-minute files with the voice note and observer notes removed. In WildTrax, individuals were counted by users scanning both the spectrogram and listening to the audio output (MacPhail et al. (2026)). Tags were then drawn to encompass the signal within the methods indicated in each project (see Table 2). Transcribers also had site photos available to optimize their species identification by having habitat context while processing.
To ensure comparability across the 18-year dataset, we harmonized the abundance metrics derived from the differing transcription protocols prior to analysis. For legacy data collected between 2012 and 2020 using the 3-minute interval method, we calculated abundance as the maximum count observed in any single time bin rather than summing across bins. This conservative approach prevents the double-counting of individuals that sing continuously throughout the recording. Conversely, for data processed or re-transcribed using the modern 1SPT protocol (time of first detection of each unique individual of each species; 2007–2011 and 2021–2025), abundance was defined as the total count of unique individuals detected over the full 10-minute duration. Specifically, in WildTrax maximum count of individuals was used (i.e. AMRO1, AMRO2, AMRO3 max of 3 represents 3 individuals total).
Code
transcription_table <- tibble(
Years = c("2007-2020", "2021-2022", "2023-2025", "2007-2011"),
`Transcription Method` = c(
"0-3.33, 3.33-6.66, 6.66-10 min",
"1 SPT - Species per task or recording (Every new individual is tagged within the 10 minute recording).",
"1 SPT - Species per task or recording (Every new individual is tagged within the 10 minute recording).",
"1 SPT - Species per task or recording (Every new individual is tagged within the 10 minute recording). Re-transcription."
),
`Bin Method` = c(
"Abundance re-starts for each 3.33-minute bin",
"Time of first detection over 10 minute period",
"Time of first detection over 10 minute period",
"Time of first detection over 10 minute period. Re-transcription."
),
`Method Details` = c(
"No cap on abundance; abundance re-starts for each bin, no total abundance for the 10-min recording",
"Time of first detection over 10 minute period",
"Time of first detection over 10 minute period",
"Time of first detection over 10 minute period. Re-transcription."
),
`Max # of Individuals` = c(
"No cap",
"Maximum of 3 individuals per 10-minute recording",
"No cap",
"No cap. Re-transcription."
)
)
transcription_table
## # A tibble: 4 × 5
## Years `Transcription Method` `Bin Method` `Method Details`
## <chr> <chr> <chr> <chr>
## 1 2007-2020 0-3.33, 3.33-6.66, 6.66-10 min Abundance r… No cap on abund…
## 2 2021-2022 1 SPT - Species per task or recording… Time of fir… Time of first d…
## 3 2023-2025 1 SPT - Species per task or recording… Time of fir… Time of first d…
## 4 2007-2011 1 SPT - Species per task or recording… Time of fir… Time of first d…
## # ℹ 1 more variable: `Max # of Individuals` <chr>
# Render the datatable
datatable(transcription_table,
options = list(
searching = TRUE,
paging = TRUE,
pageLength = 10
)) |>
formatStyle(columns = colnames(transcription_table),
backgroundColor = styleEqual(c("NA"), "lightgray"))
We also evaluated the performance of the HawkEars and BirdNET classifiers on the acoustic recordings (?@fig-classifiers).
Analytical methods
For the purpose of these analyses, abundance was defined as the count of individuals detected during point counts, rather than as a density x area relationship. All analyses took place in R 4.5.2 ‘[Not] Part in a Rumble’.
Our analysis proceeded in the following steps. We:
- Assessed if the new forest sites were effective controls for fire-affected sites;
- Evaluated consistency in species identification and abundance estimates between transcription methods and observers;
- Selected species and grouped species by guild for trend analysis;
Location correlation
To inform our subsequent modeling choices, we first evaluated the bird survey data for spatial autocorrelation (e.g., survey points close to one another exhibit similar bird counts). To test for this, we calculated a total abundance index for each location and year by summing the maximum number of individuals detected, representing the minimum number of individuals known to be present. Because survey points were typically spaced 300 m apart, we defined spatial neighbor relationships using a 1-nearest-neighbor approach (k = 1) based on great-circle distances. We derived spatial coordinates from geographic point data, excluding non-finite values to ensure valid estimation. Using the knearneigh() and knn2nb() functions in R, we constructed a spatial weights matrix with row-standardized weights to reflect immediate adjacency between points. We assessed global spatial autocorrelation in total abundance using Moran’s I under a randomization assumption. This statistic evaluates whether the landscape is spatially structured—specifically, whether nearby locations tend to have similar abundance values more often than expected by chance. To determine if this dependence was driven by localized clustering, we further calculated Local Indicators of Spatial Association (LISA) using local Moran’s I. This allowed for the identification of potential high–high, low–low, and spatial outlier patterns, with statistical significance evaluated at α = 0.05. While global Moran’s I indicated significant positive spatial autocorrelation across the study area, local Moran’s I revealed no statistically significant clusters. This pattern suggests that similarities in bird counts among nearby survey locations were spread broadly rather than concentrated in distinct hotspots, consistent with spatial dependence arising from gradual, landscape‑scale ecological processes rather than localized aggregations. Therefore, we proceeded by treating the locations as spatially independent.
Control site comparisons
The 2024 Jasper Wildfire severely burned two established acoustic transects (VALLEY5 and TEKKARA). While only one-third of true montane (valley-bottom) forest has burnt since 2023, two of three transects were severely burned. To continue to monitor bird communities in unburned forest, we added two control transects in 2025 in the true montane (COTTONWOOD and MUSHROOM), and we plan to continue to monitor both the burned and control transects.
To assess whether the new transects were appropriate controls for montane transects, we evaluated how similar their species assemblages were to those observed at existing true montane transects over the past three years, excluding the TEKARRA5 and VALLEY5 transects post-fire because these were expected to diverge substantially due to fire severity. Specifically, we compared species detections from the COTTONWOOD and MUSHROOM transects to the distribution of detections across all other true montane points to quantify compositional similarity and departure from the true montane baseline. Representativeness of control transects was assessed using a multivariate dispersion analysis (PERMDISP) based on Bray–Curtis dissimilarity, with distances of control transects to the true montane centroid compared against the distribution of distances among established true montane transects. Statistical significance of differences in dispersion was evaluated using a permutation test with 999 iterations.

Tekarra 9 site photos views from 2024 (A) pre-fire, and 2025 (B) post-fire.
Transcription observer comparisons
Transcription methods varied (Table 2) over time. Recordings from 2007 to 2020 were classified by a small number of repeat observers (legacy), while subsequent data were classified by a random selection from many observers (WildTrax). To evaluate consistency between the legacy dataset and the modern WildTrax workflow, we analyzed a subset of recordings processed using both approaches. This included re-processed legacy recordings that 1) enabled total abundance estimates and 2) used new observers in WildTrax, allowing for direct comparison of individual observer performance between the legacy and modern protocols. For each dual-processed recording, we derived the maximum count of individuals per species identified by each specific observer. We used two different metrics to verify data continuity.
First, to quantify consistency in abundance estimates, we calculated pairwise Pearson correlations between all individual observers. These relationships were visualized using a heatmap to identify any systematic deviations in counting between legacy contractors and current WildTrax observers. Figure 3 indicates that there is general consistency in abundance counts, but not high reliability between observers (e.g., some observers are more generous with counts than others) indicating there may be a bias - OR BECAUSE METHOD CHANGED - e.g., re-start between bins, cap on counts.
Second, to evaluate consistency in species detection (composition), we binarized the data to presence/absence and calculated Bray–Curtis dissimilarities between individual observers. We then applied hierarchical clustering to these dissimilarity scores to verify that the transition to WildTrax did not introduce observer-specific biases in species identification or community composition. Figure 4 shows that most observers agree on species, and there doesn’t appear to be a systematic bias during the transition to WildTrax. If the transition introduced a major bias, we would see two distinct large branches with one containing all legacy contractors and one containing all WildTrax observers. Data using both types of transcription observer method were grouped for analysis, however when available for a given year, the modern WildTrax workflow results were used because they included total abundance estimates (Table 2).
Selection of species and guilds for trend analysis
species_select <- janp_main |>
mutate(species_code = case_when(species_code == "VESP" ~ "FOSP", TRUE ~ species_code)) |>
filter(!(species_code %in% c("NONE","DOGG","HOMA","UROPAR","HORS","VIRA"))) |>
filter(!(grepl('^Unidentified',species_common_name))) |>
group_by(species_code) |>
tally() |>
inner_join(wt_get_species() |> dplyr::select(species_code, species_common_name, species_order, species_class) |> distinct())Joining with `by = join_by(species_code)`
The monitoring program dataset produced detections for 143 species of birds and 6 species of mammals (American Pika, Elk, Hoary Marmot, Red Squirrel, Columbian Ground Squirrel and White-tailed Deer). Following recommendations in the literature Prowse et al. (2021), we focused on native terrestrial species for which: (a) Jasper National Park is their breeding range; and (b) the total count across all point count sites exceeded 20 individuals (see Table 3).
To assess trends in species’ abundances and interpret community-level ecological responses, we assigned species to guilds. Guilds are groups of species that respond in similar ways to environmental changes as a result of similar uses of the environment (Doser et al. (2021)). Directional trends in abundant species (e.g., Yellow-rumped warbler (Setophaga coronata, YRWA)) can strongly influence the trend of the guilds of which these species are members. Given this limitation, trend analyses of ecological guilds often warrant further examination of common patterns of change among species within the guild. For example, if we find large variation in abundance trends across an ecoregion but not across bird guilds, this may suggest a consistent effect on the entire bird community (Doser et al. (2021)), which can determine if management should target a whole community, specific guilds or individual species. If all or many species within a guild show similar trends in relative abundance, then factors affecting the guild-related life history attributes deserve attention. Pacifici et al. (2014) recommend using bird guilds in post-hoc assessment. Each species was categorized within guilds and a species may be included in multiple guilds. Species were assigned to guilds using the Birds of North America, Elton Traits Database and bird guidebooks. To examine responses by different guild types we assigned birds by Elton trait (e.g., invertebrate, plant/seed, omnivore), dietary guild (e.g., aerial insectivore, bark forager), and habitat guild (e.g., wetland, forest generalist, late successional forest) (Table 3).
Ecoregion community composition analysis
To characterize bird community composition, we first aggregated species-level observations into a species-by-location matrix, populated with the maximum count of each species at each location. Survey points were then classified into two primary ecoregions: alpine and montane (high elevation forest). For the purpose of this analysis, the montane category served as a broad aggregate, grouping the upper subalpine, lower subalpine, and montane ecoregions into a single unit that refers to high elevation habitat below treeline. We quantified the variation in community composition explained by these two ecoregions using Redundancy Analysis (RDA) in the vegan package (Oksanen et al. (2025)) and visualized species–ecoregion relationships with ordination plots (C. Radhakrishna Rao (1964)). Finally, to test for statistical differences in composition between the Alpine and Montane groups, we performed a permutational multivariate analysis of variance (PERMANOVA; Anderson (2001)). This test was conducted on Bray–Curtis dissimilarities using 999 permutations under a reduced model.
Functional and community-level diversity
To evaluate community-level ecological responses, we examined temporal changes in functional diversity, species richness, and community structure. First, we quantified functional diversity using Rao’s Q (C. R. Rao (1982); Laliberté and Legendre (2010)), calculated via the dbFD() function in the FD package (Laliberté, Legendre, and Shipley (2014)). This metric estimates the average difference in functional traits between any two random individuals in the community. We also calculated species richness (number of unique species per location per year) and Shannon’s diversity index, which integrates richness and evenness to describe community structure. These metrics were modeled through time using linear, mixed-effects, and segmented regression models to detect both gradual and threshold-type changes. Results were summarized graphically by ecoregion and functional guild to highlight spatial variation in diversity trajectories. We also calculated species richness (number of unique species per location per year) and Shannon’s diversity index, which integrates richness and evenness to describe community structure. We assessed temporal changes in Shannon diversity and species richness using a linear mixed-effects model, with Shannon diversity and species count as the response variable. Year was centered at the mean to estimate temporal trends, and ecoregion (alpine vs. forested) was included as a fixed effect to account for differences between the two ecoregions. We included a random intercept for each location and the model was fit using restricted maximum likelihood (REML) with the lmer() function from the lme4 package.
Code
shannon_d <- janp_main |>
filter(!grepl('COTTONWOOD|MUSHROOM', location)) |>
mutate(ecoregion = case_when(ecoregion %in% c("Alpine") ~ "Alpine",
ecoregion %in% c("Upper Subalpine","Lower Subalpine","Montane") ~ "Montane")) |>
wt_tidy_species(remove = c("mammal","amphibian","abiotic","insect","unknown"), zerofill = F) |>
inner_join(wt_get_species() |> dplyr::select(species_code, species_class, species_order), by = "species_code") |>
dplyr::select(location, ecoregion, recording_date_time, species_code, species_common_name, individual_order, abundance) |>
distinct() |>
group_by(location, ecoregion, recording_date_time, species_code, species_common_name) |>
summarise(count = max(individual_order)) |>
ungroup() |>
pivot_wider(names_from = species_code, values_from = count, values_fill = 0) |>
pivot_longer(cols = -(location:species_common_name), names_to = "species", values_to = "count") |>
group_by(location, ecoregion, year = year(recording_date_time), species) |>
summarise(total_count = sum(count)) |>
ungroup() |>
group_by(location, ecoregion, year) |>
summarise(shannon_index = diversity(total_count, index = "shannon")) |>
ungroup() |>
filter(!(year == 2012 & ecoregion == "Montane"))
shannon_d$year_c <- shannon_d$year - mean(shannon_d$year)
shannon_d$ecoregion <- factor(shannon_d$ecoregion, levels = c("Alpine", "Montane"))
shannon_mod <- lmer(shannon_index ~ year_c + ecoregion + (1|location), data = shannon_d)
summary(shannon_mod)
spp_rich_location <- janp_main |>
filter(!grepl('COTTONWOOD|MUSHROOM', location)) |>
filter(data_type %in% c("legacy","single_visit_3_max","single_visit_0_max")) |>
filter(!(data_type == "single_visit_0_max" & year < 2023)) |>
wt_tidy_species(remove = c("mammal","amphibian","abiotic","insect","unknown"), zerofill = F) |>
dplyr::select(location, year, species_code) |>
distinct() |>
group_by(location, year) |>
summarise(species_count = n_distinct(species_code)) |>
ungroup() |>
inner_join(locs_summary, by = c("location" = "Location")) |>
mutate(ecoregion = case_when(ecoregion %in% c("Alpine") ~ "Alpine",
ecoregion %in% c("Upper Subalpine","Lower Subalpine","Montane") ~ "Montane")) |>
filter(!(year == 2012 & ecoregion == "Montane"))
spp_rich_location$year_c <- spp_rich_location$year - mean(spp_rich_location$year)
spp_rich_location$ecoregion <- factor(spp_rich_location$ecoregion, levels = c("Alpine", "Montane"))
spp_rich_mod <- lmer(species_count ~ year_c + ecoregion + (1|location), data = spp_rich_location)
summary(spp_rich_mod)Trend analysis
To quantify temporal changes in bird populations and community composition from 2007 to 2025, we analyzed trends in species-specific abundance, in montane and alpine assemblages separately, and grouping species by functional guilds. Analyses were designed to separate biological change from potential sampling and methodological effects, ensuring that observed patterns represented genuine ecological responses. This was achieved through a multi-step framework that (1) evaluated (1) modeled detection probability and methodological variability, (2) estimated detection-corrected abundance, and (3) evaluated long-term directional trends and associated shifts in functional and community-level diversity.
Location correlation
To inform modeling choices, we assessed spatial autocorrelation among bird survey locations across years. We examined the spatial relationship between survey points based on total species abundance per location and year in order to determine whether spatial clustering or spatial dependence exists in the species abundance data. Because survey points were typically spaced 300 m apart, we defined spatial neighbor relationships using a 1-nearest-neighbor approach (k = 1) based on great-circle distances. We derived spatial coordinates from geographic point data, excluding non-finite values to ensure valid estimation. Using the knearneigh() and knn2nb() in the spdep package in R, we constructed a spatial weights matrix with row-standardized weights to reflect immediate adjacency between points.
We assessed global spatial autocorrelation in total abundance using Moran’s I under a randomization assumption. This statistic evaluates whether the landscape is spatially structured; whether nearby locations tend to have similar abundance values more often than expected by chance. To determine if this dependence was driven by localized clustering, we further calculated Local Indicators of Spatial Association (LISA) using local Moran’s I. This allowed for the identification of potential high–high, low–low, and spatial outlier patterns, with statistical significance evaluated at α = 0.05.
While global Moran’s I indicated significant positive spatial autocorrelation across the study area, local Moran’s I revealed no statistically significant clusters. This pattern suggests that similarities in bird counts among nearby survey locations were spread broadly rather than concentrated in distinct hotspots, consistent with spatial dependence arising from gradual, landscape‑scale ecological processes rather than localized aggregations. Therefore, we proceeded by treating the locations as spatially independent.
Detection-corrected abundance estimation
Abundance was estimated using single-visit n-mixture abundance models implemented in the detect::svabu() function (Sólymos, Lele, and Bayne (2012)). This framework jointly models site-level abundance with observation and detection probability from single-visit counts as conditional likelihood estimates. The package and framework make it easy for users to add and interpret observation and detection covariates which supports logistical considerations for single-visit surveys. We first stratified analyses by ecoregion (alpine vs. forest), fitting separate detection-abundance models given the assumption site characteristics and ecological drivers are different in both ecoregions. Expected abundance per site visit was modelled using the latent abundance component (λ) specified as a function of interannual variation, local observation covariates and detection covariates. For the observation covariates, proportional landcover was measured within a 150-m buffer and classified into shrub, herbaceous cover, conifer cover, bedrock and other. Canopy cover, and mean canopy height were computed as observation-level predictors. Separate models were fitted using canopy cover and height as detection covariates but models failed to converge due to overdispersion and low sample size. Detectability was modelled simultaneously using julian date, hour of day, temperature (degrees C) and wind (km/h). Although we modeled for differences between observers, and structural vegetation influences on detection probability (see Yip, Knight, and Bayne (2025)), which included once again canopy height and canopy cover at the point level, models failed to converge. Abundance was modeled with a log link, and detection probability with a logit link. Conditional maximum likelihood estimates were used, with model fit assessed via log-likelihood and AIC with successful models below:
where detection (or availability) probability was modeled, and,
was used to generate standardized, detection-corrected expected abundances (λ) for each site visit, which were then averaged by year to produce annual abundance indices. These indices were subsequently used to evaluate temporal trends in species abundance.
Trend estimation
Temporal trends were quantified using the Mann–Kendall test (Mann (1945a); Hamed (2009)), which detects monotonic directional change, and Sen’s Slope (Pranab Kumar Sen (1968)) to estimate the magnitude of those trends. Both tests were implemented via the modifiedmk package (Hamed and Rao (1998)) and applied to the detection-corrected abundance estimates (λ̂).Sen’s slope provides an estimate of the median annual rate of change in the abundance index over the time series. To express this rate in a standardized and interpretable way, we converted Sen’s slope to a percent change per year by dividing the estimated slope by the mean annual abundance index across the full time series and multiplying by 100. This metric represents the average proportional change in abundance per year, relative to the long-term mean abundance of the species in that ecoregion. Positive values indicate increasing abundance, while negative values indicate declining abundance.
Results
Some of these analyses are still a work-in-progress. Check back soon for updates and additional details.
Control site comparisons
Multivariate dispersion differed significantly between control and reference montane transects (PERMDISP; F = 4.05, p = 0.05). Although control transects exhibited a slightly higher median distance (mean = 0.421), their distributions overlapped strongly with those of the true montane reference transects (Mean = 0.415), and there was no evidence of increased dispersion or outliers unique to replacement sites (Figure X). This suggests that replacement transects are not compositionally distinct from existing montane transects and fall within the natural variability of montane bird communities already sampled. Consequently, the addition of replacement sites is unlikely to introduce systematic bias, supporting their use as valid representatives of montane bird communities for pre- and post-fire comparisons.
Ecoregion community composition analysis
Figure 4 shows the relationship between species and ecoregion. The PERMANOVA test was performed using Bray-Curtis dissimilarity to assess whether community composition significantly differed between ecoregion groups. The analysis revealed a significant difference in community composition between alpine and montane groups. The ecoregion grouping explained approximately 27.85% of the variation in community composition, while residual variation accounted for 72.15%. These findings indicate a substantial divergence in species composition between ecoregion groups and helps to justify subsequent analyses looking at trend differences between these areas.
Functional and community-level diversity
Species richness per location is shown in Figure 6 and Shannon’s diversity over years is shown in Figure 7. Shannon diversity increased slightly over time (β = 0.0063 per year, SE = 0.0015, t = 4.20), after accounting for differences among locations and ecoregions. Montane communities exhibited higher diversity than Alpine (β = 0.28 ± 0.067, t = 4.15), and location-level variation remained substantial (SD = 0.35). Rao’s Q averaged between about 8.2 and 10.5 across survey locations, with a clear upward tendency over time (Figure 5). The non‑parametric Mann–Kendall test gave a Kendall’s τ of 0.32 (p ≈ 0.07), indicating a positive but marginally non‑significant monotonic increase in functional diversity. A simple linear regression of mean Rao’s Q against year yielded a slope of 0.051 units per year (p ≈ 0.06), suggesting an upward trend that narrowly misses the conventional 0.05 significance threshold. When comparing models using an AIC-based approach, a segmented model clearly outperformed both the linear model (ΔAIC = 2.97, weight = 0.18) and the mixed-effects model (ΔAIC ≈ 11,883, weight ≈ 0), with the segmented model having the lowest AIC (38.5, weight = 0.82). This supports the presence of a breakpoint around 2009, suggesting that functional diversity was relatively low and stable from 2007–2009, then increased to more variable but generally higher values from 2010 onward. The linear and mixed-effects models, while informative, were either slightly less supported or poorly suited to capture this shift in the trend.
Code
spp_rich_location |>
ggplot(aes(x=as.factor(year), y=species_count, fill=year)) +
geom_boxplot() +
geom_point(alpha = 0.7, colour = "grey") +
geom_smooth(method = "lm") +
theme_bw() +
facet_wrap(~ecoregion, ncol = 1) +
scale_fill_viridis_c(option = "cividis") +
xlab('Year') + ylab('Species richness per location') +
guides(fill = guide_legend(title = "Year"))Code
shannon_d |>
ggplot(aes(x = factor(year), y = shannon_index, fill = factor(year))) +
geom_boxplot() +
geom_point(alpha = 0.6, colour = "grey") +
labs(x = "Year",
y = "Shannon diversity index per location") +
theme_bw() +
guides(fill = guide_legend(title = "Year")) +
scale_fill_viridis_d(alpha = 0.8, option = "cividis") +
facet_wrap(~ecoregion, ncol = 1)Jasper 2024 fire effects
The pre-fire and 1-year post-fire points occupy overlapping but shifted regions indicating that the fire did not completely replace the community but it reorganized species composition in a consistent direction which is a classic disturbance-driven compositional shift and not a collapse or a reset. In Figure 9, the post-fire ellipse is significantly shifted indicating greater variability among sites after fire and more heterogeneous species assemblages. Species vectors identified taxa contributing most strongly to these differences, suggesting that the Jasper Wildfire altered species associations rather than uniformly affecting overall abundance. Dark-eyed Junco (Junco hyemalis, DEJU) shows a strong directional association with post-fire sites (they tolerate reduced canopy cover and more simplified forest structure), while Swainson’s Thrush (Catharus ustulatus, SWTH), Warbling Vireo (Vireo gilvus, WAVI), and Tennessee Warbler (Leiothlypis peregrina, TEWA) were strongly associated with pre-fire assemblages, and are generally associated with intact forest canopy and vertical foliage structure. This indicates that structurally complex forest conditions were an important component of bird community composition prior to the Jasper Fire.
Trends
Discussion and recommendations
Our analysis confirms that biological phenology, rather than immediate weather conditions, was the primary driver of species detections in Jasper National Park. The significant negative relationship between species richness and both Julian day and survey hour is consistent with well-established avian activity patterns. The decline in detections later in the season reflects the cessation of territorial signaling as breeding concludes (Catchpole & Slater 2008), while the negative effect of survey hour captures the waning of the dawn chorus (Robbins 1981).
Contrary to the strong temporal signals, we found that temperature, wind speed, and sky conditions had weak or non-significant effects on overall species richness. This lack of a uniform negative response likely reflects the complexity of bird community responses rather than a total lack of sensitivity. For instance, Morelli et al. (2022) found that weather effects on detectability are highly guild-specific: while wind negatively affected the detection of insectivores and granivores, it had no effect on omnivores and was actually associated with increased detection of granivore-insectivores. In our analysis, these heterogeneous responses likely offset one another when aggregated into a single “Species Richness” metric. Furthermore, by restricting surveys to favorable conditions—avoiding high winds and precipitation—the monitoring program likely minimized the variance attributable to these factors, ensuring that the observed long-term trends in abundance are driven by ecological changes rather than artifacts of sampling conditions.










