Countless wetlands have been lost and degraded globally making amphibians the most threatened vertebrate class. However, despite facing extensive threats and stressors, coastal wetlands of the Laurentian Great Lakes of North America (lakes Superior, Michigan, Huron, Erie, and Ontario) still support sizable populations of frogs (order Anura, including toads). We used data from the Great Lakes Coastal Wetland Monitoring Program to quantify the first-ever annual occurrence probabilities and trends (2011–2023) of eight marsh-breeding frog species, or groups of species, at 1,550 point count locations in 747 coastal wetlands throughout the Great Lakes, and to assess 11 potential drivers of occurrence. Across our study area, Green Frog (Lithobates clamitans) occurrence increased by 8% per year, whereas Chorus Frog (Pseudacris maculata, P. triseriata) occurrence decreased by 14% per year. We found more positive or stable trends among lakes and species (85%) than negative trends (15%). Occurrence of all species was negatively associated with one or two indicators of poor water quality: specific conductance, ammonium nitrogen, nitrate nitrogen, and urban and agricultural land cover in the surrounding watershed (median area: 12 km2). Occurrence of multiple species was positively associated with high lake levels and surrounding wetland (< 250 m) and forest (< 2.5 km) land cover and negatively associated with surrounding road density (< 2.5 km). Even though occurrence of most species was increasing or stable and was relatively high (> 50%), all will likely benefit from conservation actions. Fifty to 90% of Great Lakes coastal wetlands have been lost and converted to anthropogenic land uses leaving frog populations at a fraction of their former, original sizes. Therefore, extra precaution is critical to help ensure their growth and persistence. Improving water quality, increasing natural forest and wetland cover within 2.5 km, and reducing roads within 2.5 km of Great Lakes coastal wetlands will help conserve these important indicator species in this globally-recognized but imperiled ecosystem.

Study area and design

We surveyed frogs in coastal wetlands throughout the Great Lakes basin (Figure 1). Coastal wetlands were selected according to CWMP protocol using a stratified, random sampling procedure (Uzarski et al. 2017, 2019). Further details regarding the study design are in Burton et al. (2008). The sampling domain included all marshes greater than 4 ha in area with a permanent or periodic surface-water connection to an adjacent Great Lake or their connecting river systems (Uzarski et al. 2017). The selection of wetlands was stratified 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 one of the five 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 five years. In addition, we resampled 10% of wetlands between years according to a rotating panel design. Sampled wetlands were safely accessible and were dominated by emergent herbaceous vegetation and shallow water (< 2 m deep) containing floating and/or submerged vegetation.

Frog surveys

We conducted frog surveys at one to six 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 > 500 m apart to avoid double detections. Each point count location was surveyed for 3 min on each of three occasions at least 15 days apart between late-March and early-July, which is the main breeding season for frogs in the region. During each point count, observers recorded all species heard. Chorus frogs (Pseudacris spp.) were combined together and treefrogs (Dryophytes spp.) were combined together due to challenges with identifying species within each of these groups by sound in the field. Surveys occurred at night starting at least 0.5 hr after local sunset and only under weather conditions that were favourable for detecting all species present (no persistent or heavy precipitation; wind: Beaufort 0–3, 0–19 km/h). The first survey in the season was conducted when the air temperature was > ~5°C, the second when the air temperature was > ~10°C, and the third when the air temperature was > ~17°C. This approach maximized the odds that early, mid, and late-season breeders would be detected by timing surveys during peak calling for each group of species regardless of annual variation in timing of spring melt. We trained observers so they thoroughly understood the field protocols, and we required each observer to pass an aural frog identification test in order to collect data. For a detailed description of the sampling protocol visit greatlakeswetlands.org/‌Sampling-protocols.

Water quality surveys

We sampled water quality once between mid-June and early-September in each wetland in each year that a wetland was selected for surveys. Three water quality sampling points (replicates) were surveyed within each monodominant vegetation zone (1–4 zones per wetland). Replicates were > 15 m apart in water > 5 cm deep. Monodominant vegetation zones were visually identified as patches of macrophytes in which structurally-similar species (often a single genus) represented at least 75% of the aquatic plant community (e.g., Typha, wet meadow). Water quality was sampled at the mid-depth of the water column at each sampling point using a multi-parameter sonde for specific conductance (standardized to 25°C). For ammonium nitrogen, nitrate nitrogen, and soluble reactive phosphorus, a composite water sample from the three sampling points within each vegetation zone was filtered and frozen for later analysis in the lab following standard methods (Lipps et al. 2023). Values within each vegetation zone were averaged, and then zone-level means were averaged to yield a single measure for each wetland. For a detailed description of the sampling protocol see Uzarski et al. (2017) or visit greatlakeswetlands.org/‌Sampling-protocols. Water quality surveys were completed at fewer wetlands each year than frog surveys; therefore, data for specific conductance, ammonium nitrogen, nitrate nitrogen, and soluble reactive phosphorus were available for only 1,021 of the 1,913 (53%) wetland-year combinations with frog survey data.

Response variable

The response variable for each species was the detection (occurrence) of at least one individual during any of the three surveys at each point count location in each year (e.g., similar to Grand et al. 2020, Tozer et al. 2024), whereas environmental predictors were measured within the surrounding wetland or at the landscape scale (within 200 to 6,400 m of a point count location; Zuckerberg et al. 2012). We viewed detections/occurrences as indices of true occupancy, meaning our modeled values estimate the probability of occurrence unadjusted for detectability (MacKenzie et al. 2018). We assumed that variation in species-specific detection was uncorrelated with the predictors in our models, including year. This approach was sufficient in our case because our objective was to quantify relative differences and changes in occurrence and not to quantify true occupancy. Our assumption was warranted because our data were collected using standardized methods designed to reduce heterogeneity in detection, such as observer training and testing, as well as restrictions on temperature, time of day, and wind (Uzarski et al. 2017). We found no trends in detectability across years for any of the species in our dataset, meaning that variation in detection did not bias our estimates of annual occurrence or trends (Appendix S1:Figure S1).

The frog dataset consisted of 10,803 point count surveys completed at 1,550 point count locations in 747 coastal wetlands over 13 years (2011–2023; Figure 1, 2). In each year, we completed frog surveys at a mean of 277 point count locations (range: 216–329) in a mean of 147 wetlands (range: 111-169) throughout the Great Lakes basin. The mean annual number of point count locations in each lake was 57 points for Erie (range: 33–85 points), 81 for Huron (47–102), 43 for Michigan (24–55), 62 for Ontario (47–80), and 34 for Superior (21–56; Figure 2). We surveyed 2.1 ± 1.4 (mean ± SD) point count locations per wetland (range: 1–6). In total, we analyzed eight frog species, or groups of species: 1) American Toad; 2) American Bullfrog; 3) Chorus Frog, which included Boreal Chorus Frog (Pseudacris maculata) and Western Chorus Frog (P. triseriata); 4) Gray Treefrog, which included Cope’s Gray Treefrog (Dryophytes chrysoscelis) and Eastern Gray Treefrog (D. versicolor); 5) Green Frog; 6) Northern Leopard Frog; 7) Spring Peeper (Pseudacris crucifer); and 8) Wood Frog (Lithobates sylvaticus). We chose these species because they were widespread and bred regularly in Great Lakes coastal wetlands. Observations of other species were too few for analysis (Fowler’s Toad [Anaxyrus fowleri], Mink Frog [Lithobates septentrionalis], Pickerel Frog [L. palustris]; Appendix S1:Figure S2). To assess whether parts of our study area were out of range for some species, we divided our study area into 10 regions and dropped regions from species-specific analyses if raw (un-modeled) occurrence was < 5% of point count locations (Appendix S1:Table A1). By excluding out-of-range point count locations, we reduced the number of absences and focused our analysis on point count locations where non-detections were most likely to represent legitimate absences. As a result, Wood Frog was dropped from analyses involving Lake Erie; all of the other species were analyzed across all five of the Great Lakes.

Environmental predictors

We included 11 environmental predictors in our models that are known to influence occurrence of frogs in the Great Lakes: six described water quality (specific conductance; ammonium nitrogen; nitrate nitrogen; soluble reactive phosphorus; urban land cover in the surrounding watershed; agricultural land cover in the surrounding watershed; e.g., Rouse et al. 1999), four described characteristics of the surrounding landscape (wetland < 250 m from the point count location, wetland < 2.5 km from the point count location, forest < 2.5 km from the point count location, density of roads < 2.5 km from the point count location; e.g., Price et al. 2004), and one described lake levels (detrended, standardized Great Lakes water levels; Gnass Giese et al. 2018). We viewed specific conductance as a proxy for road salt given specific conductance was highly correlated with chloride (r = 0.83) and we had more observations for specific conductance than for chloride. Two of the water quality predictors (urban and agricultural land cover in the surrounding watershed) and all of the landscape predictors were static covariates (i.e., they were the same for all years at a given point count location), whereas the rest of the water quality predictors (specific conductance, ammonium nitrogen, nitrate nitrogen, soluble reactive phosphorus) and Great Lakes water levels were dynamic covariates (i.e., they varied annually). Landscape and water-level data at finer temporal and spatial scales, respectively, would have been preferred, but were unavailable. Nonetheless, the landscape and water-level data we used provided useful approximations of the true values, particularly at the < 2.5 km and watershed scales (e.g., Michaud et al. 2022). We measured characteristics of the surrounding landscape within 2.5 km because studies suggested frogs respond most to land cover within 2–3 km (Carr and Fahrig 2001, Houlahan and Findlay 2003, Price et al. 2004, Eigenbrod et al. 2008, 2009). The environmental predictors were not strongly correlated with each other (-0.52 < r < 0.52; Appendix S1:Figure S3). We considered additional environmental predictors such as water temperature and clarity, pH, dissolved oxygen, total nitrogen, and total phosphorus, but we did not include them in our modeling because they were either correlated with at least one predictor already in our final set (|r| > 0.6) or were not associated or only weakly associated with frog occurrence based on existing literature.  

Derivation of dynamic water quality predictors is described in section 2.3. Percent urban and agricultural land cover in the surrounding watershed were from Host et al. (2019), with watersheds defined by Forsyth et al. (2016). Percent wetland cover < 250 m was based on the coastal wetland inventory layer built by the Great Lakes Coastal Wetland Consortium (GLCWC_CWI; Albert et al. 2004, 2005, Burton et al. 2008, Uzarski et al. 2017). Percent wetland and forest cover < 2.5 km were based on the land cover and use layer built by the Great Lakes Aquatic Habitat Framework, which combined data from the 2011 version of the National Land Cover Database in the USA (mrlc.‌gov) and the 2012 version of the Southern Ontario Land Resource Information System (geohub.‌lio.‌gov.‌on.‌ca). Road density was derived from the population and roads layer built by the Great Lakes Aquatic Habitat Framework using data from the US Census Bureau (census.‌gov) and Statistics Canada (statcan.‌gc.‌ca). The coastal wetland, land cover and use, and population and roads layers are available through the Great Lakes Aquatic Habitat Framework at hub.glahf.org. We used ArcGIS 10.8.1 (Esri 2020) to overlay point count locations onto the GIS layers and extracted the relevant predictors for each point (Figure 3). Yearly lake levels were from the Great Lakes Environmental Research Laboratory of the National Oceanic and Atmospheric Administration available at glerl.‌noaa.‌gov/‌data/‌wlevels. We used the mean yearly lake level from May to July since these months overlapped with surveys in all years and regions. We detrended lake levels from year by using the residuals from a line of best fit for each lake, given that lake levels generally increased in all lakes over the course of the study. Lake levels were also standardized across lakes by subtracting the long-term mean (2011–2023) for each lake from the annual value for each lake and dividing by the standard deviation, given the reference value is the same for all lakes (International Great Lakes Datum 1985). Our detrended, standardized lake levels therefore represent lake levels without being confounded with year (Figure 4). To facilitate model convergence, we censored extreme observations for some of the predictors by dropping all values > mean + 3 times the SD. Thus, we removed 28 (1%) specific conductance observations > 977 µS/cm, 42 (2%) ammonium nitrogen observations > 0.15 mg/L, 48 (2%) nitrate nitrogen observations > 1.2 mg/L, 57 (3%) soluble reactive phosphorus observations > 0.14 mg/L, and 33 (1%) road density observations > 5.8 km/ha.

Statistical modeling

We used Bayesian hierarchical models with spatial dependencies (e.g., Thogmartin et al. 2004, Smith et al. 2015, Ethier and Nudds 2015, Ethier et al. 2022). We fit these models using Integrated Nested Laplace Approximation (INLA) and the R-INLA package (Rue et al. 2009) for R statistical computing (version 4.3.0; R Core Team 2024). Methods generally followed Tozer et al. (2024) and Meehan et al. (2019), but given that the response variable was occurrence and not abundance, we used a binomial distribution with a logit link. For each species, we modeled the expected occurrence 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). By accounting for this spatial structure, the model allowed occurrence and trend information to be shared among adjacent lakes (as described later), 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 occurrence (у_i,j,t) at a point count location within a given wetland (j), lake (i), and year (t). The expected occurrence per lake within a given year (µ_i,t) for each of the eight species took the form:

log(µ_it) = α_i + τ_iΤ_i,_j,t + κ_j + у_i,t + β₁X_1j ... β ₁₁X_11j

where α is the random lake intercept; T is the year, indexed to 2023; τ is the random lake slope effect; κ is the random wetland effect; and у is the random lake-year effect. Environmental predictors (X₁ – X₁₁) characterized water quality (six predictors), characteristics of the surrounding landscape (four predictors), and lake levels (one predictor), as described in detail in section 2.5.

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 occurrence 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 earlier. 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 occurrence during the final year of the time series (Meehan et al. 2019). We accounted for differences in occurrence among wetlands (κ) with an independent and identically distributed (idd) random effect. To derive an annual occurrence probability per lake, we included a random effect per lake-year (у_i,t) with an idd, and combined these effects with α and τ. The coefficients of X ₁ to X ₁₁ were given normal priors with means of zero and precision equal to 0.001. We scaled the spatial structure of 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 and lake-year effects with parameter values UPC = 1 and PC = 0.01. 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). Following analysis, we computed posterior estimates of trends (τ), transformed into constant rates of population change using 100 × (exp(τ)-1), 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). We note that all of our models successfully converged.

We explored the predictive ability of the full model (i.e., spatially-explicit structure and environmental predictors) for each species by comparing it with two reduced models: 1) a spatially-explicit hierarchical model with no environmental predictors, and (2) a generalized linear model with environmental predictors but no spatially-explicit structure. These reduced models took the form:

1) log(µ_it) = α_i + τ_iΤ_i,_j,t + κ_j + у_i,t

2) log(µ_it) = α + β₁X_1j ... β ₁₁X_11j

The predictive ability of each of the three models for each species was measured using area under the curve (AUC), with 0.7 to 0.8 indicating acceptable predictive ability, 0.8 to 0.9 excellent ability, and greater than 0.9 outstanding ability (Mandrekar 2010). AUCs were lowest for the generalized linear model with environmental predictors but no spatially-explicit structure (0.59–0.84, n = 8 species), suggesting poor to acceptable predictive ability if only environmental predictors were used to model occupancy. The spatially-explicit hierarchical model with no environmental predictors and the full model each performed much better (AUCs: 0.85–0.98 in each case), indicating excellent predictive ability for both models. Therefore, we used the full model for inference.

Nearly half (47%) of our dataset was missing observations for specific conductance, ammonium nitrogen, nitrate nitrogen, and soluble reactive phosphorus, and the R-INLA package drops cases with missing data. Therefore, we conducted analyses in two ways. Inferences involving specific conductance, ammonium nitrogen, nitrate nitrogen, and soluble reactive phosphorus were based on a reduced dataset due to missing values, whereas all other inferences were based on models that dropped specific conductance, ammonium nitrogen, nitrate nitrogen, and soluble reactive phosphorus in order to utilize the full dataset. This approach maximized the power of our data for inference.

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.

Albert, D., J. Ingram, D. Thompson, and D. Wilcox. 2004. Great Lakes coastal wetland inventory. Developed by the Great Lakes Coastal Wetland Consortium in cooperation with Great Lakes Commission; Environment Canada, Canadian Wildlife Service - Ontario Region; U.S. Geological Service, Water Resources Discipline; Michigan Natural Features Inventory; and Ontario Ministry of Natural Resources, Ann Arbor, Michigan, USA.

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.

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. 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, Illinois, USA. Great Lakes Commission, Ann Arbor, Michigan, USA.

Carr, L. W., and L. Fahrig. 2001. Effect of road traffic on two amphibian species of differing vagility. Conservation Biology 15:1071–1078.

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.

Eigenbrod, F., S. J. Hecnar, and L. Fahrig. 2008. The relative effects of road traffic and forest cover on anuran populations. Biological Conservation 141:35–46.

Eigenbrod, F., S. J. Hecnar, and L. Fahrig. 2009. Quantifying the road-effect zone: threshold effects of a motorway on anuran populations in Ontario, Canada. Ecology and Society 14.

Esri. 2020. ArcGIS 10.8.1. Environmental Systems Resource Institute, Redland, California, USA.

Ethier, D. M., and T. D. Nudds. 2015. Scalar considerations in population trend estimates: implications for recovery strategy planning for species of conservation concern. The Condor 117:545–559.

Ethier, D., R. Torrenta, and A.-L. Kouwenberg. 2022. Spatially explicit population trend estimates of owls in the Maritime provinces of Canada and the influence of call playback. Avian Conservation and Ecology 17.

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, M. Robertson, D. M. Infante, L. Wang, J. E. McKenna, and G. Whelan. 2016. The Great Lakes hydrography dataset: consistent, binational watersheds for the Laurentian Great Lakes basin. JAWRA 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.

Gnass Giese, E. E., R. W. Howe, A. T. Wolf, and G. J. Niemi. 2018. Breeding birds and anurans of dynamic coastal wetlands in Green Bay, Lake Michigan. Journal of Great Lakes Research 44:950–959.

Grand, J., S. P. Saunders, N. L. Michel, L. Elliott, S. Beilke, A. Bracey, T. M. Gehring, E. E. Gnass Giese, R. W. Howe, B. Kasberg, N. Miller, G. J. Niemi, C. J. Norment, D. C. Tozer, J. Wu, and C. Wilsey. 2020. Prioritizing coastal wetlands for marsh bird conservation in the U.S. Great Lakes. Biological Conservation 249:108708.

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.

Houlahan, J. E., and C. S. Findlay. 2003. The effects of adjacent land use on wetland amphibian species richness and community composition. Canadian Journal of Fisheries and Aquatic Sciences 60:1078–1094.

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

Lipps, W., E. Braun-Howland, and T. Baxter, editors. 2023. Standard methods for the examination of water and wastewater. 24th edition. American Public Health Association, American Water Works Association, Water Environment Federation. APHA Press, Washington, DC.

MacKenzie, D. I., J. D. Nichols, J. A. Royle, K. H. Pollock, L. Bailey, and J. E. Hines. 2018. Occupancy estimation and modeling: inferring patterns and dynamics of species occurrence. Academic Press, London, United Kingdom.

Mandrekar, J. N. 2010. Receiver operating characteristic curve in diagnostic test assessment. Journal of Thoracic Oncology 5:1315–1316.

Meehan, T. D., N. L. Michel, and H. Rue. 2019. Spatial modeling of Audubon Christmas Bird Counts reveals fine-scale patterns and drivers of relative abundance trends. Ecosphere 10:e02707.

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

Price, S. J., D. R. Marks, R. W. Howe, J. M. Hanowski, and G. J. Niemi. 2004. The importance of spatial scale for conservation and assessment of anuran populations in coastal wetlands of the western Great Lakes, USA. Landscape Ecology 20:441–454.

R Core Team. 2024. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org.

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.

Rouse, J. D., C. A. Bishop, and J. Struger. 1999. Nitrogen pollution: an assessment of its threat to amphibian survival. Environmental Health Perspectives 107:799–803.

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.

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.

Smith, A. C., M.-A. R. Hudson, C. M. Downes, and C. M. Francis. 2015. Change points in the population trends of aerial-insectivorous birds in North America: synchronized in time across species and regions. PLOS ONE 10:e0130768.

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

Thogmartin, W. E., J. R. Sauer, 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., A. M. Bracey, G. E. Fiorino, T. M. Gehring, E. E. Gnass Giese, G. P. Grabas, R. W. Howe, G. J. Lawrence, G. J. Niemi, B. A. Wheelock, and D. M. Ethier. 2024. Increasing marsh bird abundance in coastal wetlands of the Great Lakes, 2011–2021, likely caused by increasing water levels. Ornithological Applications 126:duad062.

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, T. M. Gehring, G. P. Grabas, A. Garwood, R. W. Howe, L. B. Johnson, G. A. Lamberti, A. H. Moerke, B. A. Murry, G. J. Niemi, C. J. Norment, C. R. Ruetz, A. D. Steinman, D. C. Tozer, R. Wheeler, T. K. O’Donnell, and J. P. Schneider. 2017. Standardized measures of coastal wetland condition: implementation at a Laurentian Great Lakes basin-wide scale. Wetlands 37:15–32.

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, R. W. Howe, L. B. Johnson, G. A. Lamberti, A. H. Moerke, G. J. Niemi, T. Redder, C. R. Ruetz III, A. D. Steinman, D. C. Tozer, and T. K. O’Donnell. 2019. Leveraging a landscape-level monitoring and assessment program for developing resilient shorelines throughout the Laurentian Great Lakes. Wetlands 39:1357–1366.

Zuckerberg, B., A. Desrochers, W. M. Hochachka, D. Fink, W. D. Koenig, and J. L. Dickinson. 2012. Overlapping landscapes: a persistent, but misdirected concern when collecting and analyzing ecological data. The Journal of Wildlife Management 76:1072–1080.

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.