Wetlands of the Laurentian Great Lakes of North America, i.e., lakes Superior, Michigan, Huron, Erie, and Ontario, provide critical habitat for marsh birds. We used 11 years (2011–2021) of data collected by the Great Lakes Coastal Wetland Monitoring Program at 1,962 point count locations in 792 wetlands to quantify the first-ever annual abundance indices and trends of 18 marsh-breeding bird species in coastal wetlands throughout the entire Great Lakes. Nine species (50%) increased by 8–37% per year across all of the Great Lakes combined, whereas none decreased. Twelve species (67%) increased by 5–50% per year in at least 1 of the 5 Great Lakes, whereas only 3 species (17%) decreased by 2–10% per year in at least 1 of the lakes. There were more positive trends among lakes and species (n = 34, 48%) than negative trends (n = 5, 7%). These large increases are welcomed because most of the species are of conservation concern in the Great Lakes. Trends were likely caused by long-term, cyclical fluctuations in Great Lakes water levels. Lake levels increased over most of the study, which inundated vegetation and increased open water-vegetation interspersion and open water extent, all of which are known to positively influence abundance of most of the increasing species and negatively influence abundance of all of the decreasing species. Coastal wetlands may be more important for marsh birds than once thought if they provide high-lake-level-induced population pulses for species of conservation concern. Coastal wetland protection and restoration are of utmost importance to safeguard this process. Future climate projections show increases in lake levels over the coming decades, which will cause “coastal squeeze” of many wetlands if they are unable to migrate landward fast enough to keep pace. If this happens, less habitat will be available to support periodic pulses in marsh bird abundance, which appear to be important for regional population dynamics. Actions that allow landward migration of coastal wetlands during increasing water levels by removing or preventing barriers to movement, such as shoreline hardening, will be useful for maintaining marsh bird breeding habitat in the Great Lakes.

Study Area and Design

We conducted our study in coastal wetlands throughout the entire Great Lakes basin (see Figure 1 in Tozer et al.). We selected coastal wetlands using a stratified, random sampling protocol (Uzarski et al. 2017, 2019). Further details regarding the study design are in Burton et al. (2008). The sampling universe was all coastal wetlands greater than 4 ha in size with a permanent or periodic surface-water connection to an adjacent Great Lake or their connecting river systems (Uzarski et al. 2017). We stratified our selection of wetlands for the study by 1) wetland hydrogeomorphic type (riverine, lacustrine, barrier protected; Albert et al. 2005), 2) region (northern or southern; Danz et al. 2005), and 3) lake (i.e., the watershed of 1 of the 5 Great Lakes). We sampled approximately 20% of all wetlands in each stratum each year, so that nearly all coastal wetlands within the Great Lakes basin meeting the selection criteria were sampled at least once every 5 years. In addition, we resampled 10% of wetlands between years according to a rotating panel design. Sampled wetlands were dominated by emergent, herbaceous vegetation and shallow water (< 2 m deep) containing floating and/or submerged vegetation. We surveyed a mean of 369 point count locations (range: 309–423) for marsh birds in each year throughout the Great Lakes basin, and the mean annual number of point count locations in each lake was 65 points for Erie (range: 31–98 points), 99 for Huron (76–116), 56 for Michigan (36–71), 95 for Ontario (77–117), and 55 for Superior (35–78; see Figure 2 in Tozer et al.).

Bird Surveys

We conducted surveys at 1–8 fixed point count locations at the edge of, or within, each wetland in each year that a wetland was selected for surveys. Point count locations were > 250 m apart to avoid double counting individuals. We surveyed each point count location twice per year, at least 15 days apart, between 20 May and 10 July, which was the peak breeding period for marsh birds in the study area. Surveys took place either in the morning (30 min before sunrise to 4 h after sunrise) or the evening (4 h before sunset to 30 min after sunset), with 1 or both of the 2 surveys being in the morning each year (Tozer et al. 2017). We conducted surveys only when there was no precipitation and wind was < 20 km/h (Beaufort 3 or less). Each point count survey lasted 10 min, consisting of an initial 5-min passive listening period followed by a 5-min call broadcast period. The call broadcast period was intended to increase detections of secretive species by eliciting auditory responses and was composed of 30 sec of vocalizations followed by 30 sec of silence for each of the following: 1) Least Bittern, 2) Sora, 3) Virginia Rail, 4) a mixture of American Coot and Common Gallinule, and 5) Pied-billed Grebe, in that order. We trained observers so they thoroughly understood the field protocols and we required each observer to pass an aural and visual bird identification test in order to collect data. CWMP bird surveys were 15 min in duration from 2011 to 2018 but were reduced to 10 min from 2019 to 2021 (Tozer et al. 2017). To accommodate changes in survey protocol, we filtered the data to only include birds detected in the first 10 min of point counts from 2011 to 2018. For a detailed description of the sampling protocol visit greatlakeswetlands.org/Sampling-protocols.

Response Variable

The response variable for each species was the maximum number of individuals observed during either of the 2 surveys at each point count location in each year (Tozer 2020, Hohman et al. 2021). We viewed these counts as indices of true density, meaning our modeled values estimated relative abundance (e.g., Thogmartin et al. 2004). We assumed that variation in species-specific detection was uncorrelated with the predictors in our models, including year. This was sufficient in our case because our objective was to quantify relative differences and changes in abundance and not to quantify actual density. Our assumption was warranted because our data were collected using standardized methods designed to reduce heterogeneity in detection, e.g., observer training and testing, as well as restrictions on survey date, time of day, and wind (Conway 2011, Uzarski et al. 2017). It was further justified by other long-term, broad-scale studies of birds based on point counts conducted using similar standardized approaches, which found no differences in year or covariate effects based on counts that were adjusted or unadjusted for detection (Etterson et al. 2009, Zlonis et al. 2019). We note that long-term (1996–2013) marsh-breeding bird monitoring data collected throughout the developed, southern portion of the Great Lakes basin showed no systematic trends in detectability over time for 14 of 15 (93%) species (Tozer 2016). We also found no trends in detectability across years for all of the species in our dataset (see Supplemental Material Figure S1 in Tozer et al.), meaning that differences in detection did not bias our estimates of annual abundance indices or trends. Therefore, we did not adjust for detectability, which has been supported, for instance, by Hutto (2016) and Johnson (2008).

The dataset consisted of 8,120 surveys completed at 1,962 point count locations in 792 coastal wetlands in 599 watersheds (defined by Forsyth et al. [2016]) over 11 years (2011–2021; see Figure 1, 2 and Supplemental Material Table S1 in Tozer et al.). There were 2.2 ± 1.6 (mean ± SD) point count locations per wetland (range: 1–8) and 1.3 ± 0.9 wetlands per watershed (range: 1–9). In total, we analyzed 18 species: 1) American Bittern, 2) American Coot, 3) Black Tern, 4) Common Gallinule, 5) Common Grackle, 6) Common Yellowthroat, 7) Forster’s Tern, 8) Least Bittern, 9) Marsh Wren, 10) Mute Swan, 11) Pied-billed Grebe, 12) Red-winged Blackbird, 13) Sandhill Crane, 14) Sedge Wren, 15) Sora, 16) Swamp Sparrow, 17) Virginia Rail, and 18) Wilson’s Snipe. We chose these species because they were of conservation interest in the Great Lakes region (e.g., Bianchini and Tozer 2023) and regularly nested or foraged in Great Lakes coastal wetlands. We attempted to model abundance and trends for Trumpeter Swan (Cygnus buccinator) and Yellow-headed Blackbird (Xanthocephalus xanthocephalus), but data were too sparse for the models to converge. We considered some regions of our study area to be out of range for some species. We accounted for this by dividing our study area into 10 regions and dropped any of them from species-specific analyses if naive occupancy was < 5% (Supplemental Material Table S2). By excluding out-of-range point count locations, we reduced the number of zero counts and focused our analysis on point count locations where zero counts were most likely to represent legitimate absences. As a result, the number of marsh-breeding bird species for which we quantified abundance and trends varied by lake due to uneven species occurrences across the study area: Superior (n = 10), Ontario (n = 12), Erie (n = 16), Huron (n = 16), and Michigan (n = 17). The CWMP bird survey data are available by request at greatlakeswetlands.org.

Environmental Predictors

We included the following environmental predictors in our models, which were known to influence abundance of marsh-breeding birds in the Great Lakes: 1) percent local wetland cover within 250 m of point count locations (as a proxy for wetland size; e.g., Studholme et al. 2023), 2) detrended, standardized Great Lakes water levels (to avoid correlation with year; e.g., Hohman et al. 2021, Denomme-Brown et al. 2023), 3) percent urban land cover in the surrounding watershed (e.g., Rahlin et al. 2022), and 4) percent agricultural land cover in the surrounding watershed (e.g., Saunders et al. 2019). The land cover predictors were static covariates (i.e., they were the same for all years), whereas detrended, standardized Great Lakes water level was a dynamic covariate (i.e., it varied annually). Land cover and water-level information at finer spatial and temporal scales would have been preferred, but such data were unavailable. Nonetheless, it is reasonable to assume that the land cover and water-level data we used provided useful approximations of the true values, particularly at the watershed scale (e.g., Michaud et al. 2022). Percent local wetland cover was based on the coastal wetland layer built by the Great Lakes Coastal Wetland Consortium (Burton et al. 2008, Uzarski et al. 2017), and percent urban and agricultural land cover were from Host et al. (2019) with watersheds defined by Forsyth et al. (2016); all of these data are available at glahf.org/data. We used ArcGIS 10.8.1 to overlay CWMP sample points onto the land cover layers and extracted the relevant predictors for each point (see Figure 3 in Tozer et al.). Yearly water levels were from the National Oceanic and Atmospheric Administration (noaa.gov). We used the mean yearly water level from May to July since these months overlapped with our survey period. We detrended water levels from year by using the residuals from a line of best fit for each lake, given that water levels generally increased in all lakes over the course of the study. Water levels were also standardized across lakes by dividing the annual value for each lake by the long-term mean (2011–2021) for each lake, given the reference value is the same for all lakes (International Great Lakes Datum 1985). Our detrended, standardized lake levels therefore represent water levels without being confounded with year (see Figure 4 in Tozer et al.). The environmental predictors were not correlated (-0.2 < r < 0.2; see Supplemental Material Figure S2 in Tozer et al.).

Statistical Modeling

We fit models in a Bayesian framework with Integrated Nested Laplace Approximation (INLA) using the R-INLA package (Rue and Martino 2009) for R statistical computing (version 4.2.0; R Core Team 2022). For each species, we modeled the expected (predicted mean) number of individuals per point count location in each Great Lake in each year, as well as the trend in these values across years in each lake, and then pooled the lake-specific trends to obtain Great Lakes-wide estimates. We included spatial structure in the models using an intrinsic conditional autoregressive (iCAR) structure (Besag et al. 1991), which allowed for information on relative abundance to be shared across lakes sharing basin boundaries. By accounting for this spatial structure in counts, the model allowed abundance and trend information to be shared among adjacent lakes (as described below), which improved estimates for lakes with limited sample sizes (Bled et al. 2013) and reduced the amount of spatial autocorrelation in model residuals (Zuur et al. 2017).

We modeled counts уi,j,t using the maximum number of individuals observed at a point count location within a given wetland (j), lake (i), and year (t). The expected counts per lake within a given year µi,t for each of the 18 species took the form:

log(µit) = αi + τiΤi,j,t + κj + ρj + уi,t + β1Wj + β2Lj + β3Uj + β4Ai

where α = random lake intercept; T = year, indexed to 2021; τ = random lake slope effect; κ = random wetland effect; ρ = random wetland type effect; and у = random lake-year effect. Environmental predictors included: W = percent local wetland cover within 250 m; L = detrended, standardized water level; U = percent urban land cover in the surrounding watershed; and A = percent agricultural land cover in the surrounding watershed.

The random lake intercept (αi) had an iCAR structure, where values of αi came from a normal distribution with a mean value related to the average of adjacent lakes. The random lake intercept also had a conditional variance proportional to the variance across adjacent lakes and inversely proportional to the number of adjacent lakes. We modeled the random lake slopes (τi) as spatially structured, lake-specific, random slope coefficients for the year effect, using the iCAR structure, with conditional means and variances as described above. We incorporated spatial structure into the random lake slopes (τi) to allow for information about year effects to be shared across neighboring lakes, and to allow year effects to vary among lakes. We transformed year (T) such that the maximum year was 0, and each preceding year was a negative integer. This scaling meant that the estimates of the random lake intercepts (αi) could be interpreted as the lake-specific expected counts (i.e., index of abundance) during the final year of the time series. We accounted for differences in relative abundance among wetlands (κ) and wetland types (ρ) with an independent and identically distributed (idd) random effect. To derive an annual index of abundance per lake, we included a random effect per lake-year (у) with an idd, and combined these effects with α and τ. Β1, β2, β3, and β4 were given normal priors with mean of zero and precision equal to 0.001. We scaled the spatial structure parameters α and τ such that the geometric mean of marginal variances was equal to one (Sørbye and Rue 2014, Riebler et al. 2016, Freni-Sterrantino et al. 2018), and priors for precision parameters were penalized complexity (PC) priors, with parameter values UPC = 1 and PC = 0.01 (Simpson et al. 2017). We also assigned precision for the random wetland, wetland type, and lake-year effects with a PC prior with parameter values previously stated. In general, the weakly informed priors used here tend to shrink the structured and unstructured random effects towards zero in the absence of a strong signal (Simpson et al. 2017).

We validated distributional assumptions with simulation to ensure models could handle the large number of zero counts for some species. The abundance of most species was modeled using a zero-inflated Poisson (ZIP) distribution. Common Grackle and Red-winged Blackbird, which were more frequently detected compared to the other species, better fit a negative binomial distribution, and Common Yellowthroat better fit a Poisson distribution. We further validated models by visually inspecting 1) the fit versus raw counts; 2) residuals versus predictors; and 3) the estimate for Ф, the dispersion parameter (Zuur and Ieno 2016). Our visual inspections of fit versus raw counts suggested models were not overfit and were able to capture the variation of the raw counts. In general, residuals versus fit values behaved randomly around the zero line and residuals appeared to behave randomly with each predictor, suggesting the models fit well. The dispersion statistics were around 1 for all species, ranging lowest for Common Yellowthroat (0.72) and highest for Mute Swan (3.38), suggesting some residual under and over dispersion, respectively. Mute Swan had some high counts (outliers) which may have contributed to this. Following model analysis, we computed posterior estimates of trends (τ) and associated credible intervals for the full extent of the study area (i.e., by pooling lake-specific trends) using lake watershed size to calculate area-weighted averages (Link and Sauer 2002).

References

Albert, D. A., D. A. Wilcox, J. W. Ingram, and T. A. Thompson (2005). Hydrogeomorphic classification for Great Lakes coastal wetlands. Journal of Great Lakes Research 31:129–146.

Besag, J., J. York, and A. Mollié (1991). Bayesian image restoration, with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics 43:1–20.

Bianchini, K., and D. C. Tozer (2023). Using Breeding Bird Survey and eBird data to improve marsh bird monitoring abundance indices and trends. Avian Conservation and Ecology 18(1):4.

Bled, F., J. Sauer, K. Pardieck, P. Doherty, and J. A. Royle (2013). Modeling trends from North American breeding bird survey data: a spatially explicit approach. PLoS ONE 8, e81867.

Burton, T. M., J. C. Brazner, J. J. H. Ciborowski, G. P. Grabas, J. Hummer, J. Schneider, and D. G. Uzarski (Editors) (2008). Great Lakes Coastal Wetlands Monitoring Plan. Developed by the Great Lakes Coastal Wetlands Consortium, for the US EPA, Great Lakes National Program Office, Chicago, IL. Great Lakes Commission, Ann Arbor, Michigan, USA.

Conway, C. J. (2011). Standardized North American marsh bird monitoring protocol. Waterbirds 34:319–346.

Danz, N. P., R. R. Regal, G. J. Niemi, V. J. Brady, T. Hollenhorst, L. B. Johnson, G. E. Host, J. M. Hanowski, C. A. Johnston, T. Brown, J. Kingston, and J. R. Kelly (2005). Environmentally stratified sampling design for the development of Great Lakes environmental indicators. Environmental Monitoring and Assessment 102:41–65.

Denomme-Brown, S. T., G. E. Fiorino, T. M. Gehring, G. J. Lawrence, D. C. Tozer, and G. P. Grabas (2023). Marsh birds as ecological performance indicators for Lake Ontario outflow regulation. Journal of Great Lakes Research 49:479–490.

Etterson, M. A., G. J. Niemi, and N. P. Danz (2009). Estimating the effects of detection heterogeneity and overdispersion on trends estimated from avian point counts. Ecological Applications 19:2049–2066.

Forsyth, D. K., C. M. Riseng, K. E. Wehrly, L. A. Mason, J. Gaiot, T. Hollenhorst, C. M. Johnston, C. Wyrzykowski, G. Annis, C. Castiglione, K. Todd, et al. (2016) The Great Lakes hydrography dataset: consistent, binational watersheds for the Laurentian Great Lakes basin. Journal of the American Water Resources Association 52:1068–1088.

Freni-Sterrantino, A., M. Ventrucci, and H. Rue (2018). A note on intrinsic conditional autoregressive models for disconnected graphs. Spatial and Spatio-temporal Epidemiology 26:25–34.

Hohman, T. R., R. W. Howe, D. C. Tozer, E. E. Gnass Giese, A. T. Wolf, G. J. Niemi, T. M. Gehring, G. P. Grabas, and C. J. Norment (2021). Influence of lake-levels on water extent, interspersion, and marsh birds in Great Lakes coastal wetlands. Journal of Great Lakes Research 47:534–545.

Host, G. E., K. E. Kovalenko, T. N. Brown, J. J. H. Ciborowski, and L. B. Johnson (2019). Risk-based classification and interactive map of watersheds contributing anthropogenic stress to Laurentian Great Lakes coastal ecosystems. Journal of Great Lakes Research 45:609–618.

Hutto, R. L. (2016). Should scientists be required to use a model-based solution to adjust for possible distance-based detectability bias? Ecological Applications 26:1287–1294.

Johnson, D. H. (2008). In defense of indices: the case of bird surveys. Journal of Wildlife Management 72:857–868.

Link, W. A., and J. R. Sauer (2002). A hierarchical analysis of population change with application to Cerulean Warblers. Ecology 83:2832–2840.

Michaud, W., J. Telech, M. Green, B. Daneshfar, and M. Pawlowski (2022). Sub-indicator: land cover. In State of the Great Lakes 2022 Technical Report. Published by Environment and Climate Change Canada and U.S. Environmental Protection Agency.

R Core Team (2022). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Rahlin, A. A., S. P. Saunders, and S. Beilke (2022). Spatial drivers of wetland bird occupancy within an urbanized matrix in the upper midwestern United States. Ecosphere 13, e4232.

Riebler, A., S. H. Sørbye, D. Simpson, and H. Rue (2016). An intuitive Bayesian spatial model for disease mapping that accounts for scaling. Statistical Methods in Medical Research 25:1145–1165.

Rue, H., S. Martino, and N. Chopin (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society Series B (Statistical Methodology) 71:319–392.

Saunders, S. P., K. A. L. Hall, N. Hill, and N. L. Michel (2019). Multiscale effects of wetland availability and matrix composition on wetland breeding birds in Minnesota, USA. Condor 121:duz024.

Simpson, D., H. Rue, A. Riebler, T. G. Martins, and S. H. Sørbye (2017). Penalising model component complexity: a principled, practical approach to constructing priors. Statistical Science 32:1–28.

Sørbye, S. H., and H. Rue (2014). Scaling intrinsic Gaussian Markov random field priors in spatial modeling. Spatial Statistics 8:39–51.

Studholme, K. R., G. E. Fiorino, G. P. Grabas, and D. C. Tozer (2023). Influence of surrounding land cover on marsh-breeding birds: implications for wetland restoration and conservation planning. Journal of Great Lakes Research 49:318–331.

Thogmartin, W. E., J. R. Sauer JR, and M. G. Knutson (2004). A hierarchical spatial model of avian abundance with application to Cerulean Warblers. Ecological Applications 14:1766–1779.

Tozer, D. C. (2016). Marsh bird occupancy dynamics, trends, and conservation in the southern Great Lakes basin: 1996 to 2013. Journal of Great Lakes Research 42:136–145.

Tozer, D. C. (2020). Great Lakes Marsh Monitoring Program: 25 years of conserving birds and frogs. Birds Canada, Port Rowan, Ontario, Canada.

Tozer, D. C., C. M. Falconer, A. M. Bracey, E. E. Gnass Giese, G. J. Niemi, R. W. Howe, T. M. Gerhing, and C. J. Norment (2017). Influence of call broadcast timing within point counts and survey duration on detection probability of marsh breeding birds. Avian Conservation and Ecology 12(2):8.

[Tozer et al.] Tozer DC, Bracey AM, Fiorino GE, Gehring TM, Gnass Giese EE, Grabas GP, Howe RW, Lawrence GJ, Niemi GJ, Wheelock BA, Ethier DM. Increasing marsh bird abundance in coastal wetlands of the Great Lakes, 2011–2021, likely caused by increasing water levels. Ornithological Applications.

Uzarski, D. G., D. A. Wilcox, V. J. Brady, M. J. Cooper, D. A. Albert, J. J. H. Ciborowski, N. P. Danz, A. Garwood, J. P. Gathman, T. M. Gehring, G. P. Grabas, et al. (2019). Leveraging a landscape-level monitoring and assessment program for developing resilient shorelines throughout the Laurentian Great Lakes. Wetlands 39:1357–1366.

Uzarski, D. G., V. J. Brady, M. J. Cooper, D. A. Wilcox, D. A. Albert, R. P. Axler, P. Bostwick, T. N. Brown, J. J. H. Ciborowski, N. P. Danz, J. P. Gathman, et al. (2017). Standardized measures of coastal wetland condition: implementation at a Laurentian Great Lakes basin-wide scale. Wetlands 37:15–32.

Zlonis, E. J., N. G. Walton, B. R. Sturtevant, P. T. Wolter, and G. J. Niemi (2019). Burn severity and heterogeneity mediate avian response to wildfire in a hemiboreal forest. Forest Ecology and Management 439:70–80.

Zuur, A. F., and E. I. Ieno (2016). A protocol for conducting and presenting results of regression-type analyses. Methods in Ecology and Evolution 7:636–645.

Zuur, A. F., E. I. Ieno, and A. A. Saveliev (2017). Beginner’s guide to spatial, temporal and spatial-temporal ecological data analysis with R-INLA. Volume I: Using GLM and GLMM. Highland Statistics, Newburgh, United Kingdom.