DOI: 10.5061/dryad.hmgqnk9v8
title: Occurrence patterns and trends of frogs in coastal wetlands of the Great Lakes call for further habitat restoration

description: 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 (&lt; 250 m) and forest
 (&lt; 2.5 km) land cover and negatively associated with surrounding
 road density (&lt; 2.5 km). Even though occurrence of most species was
 increasing or stable and was relatively high (&gt; 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.
description: <strong>Study area and design</strong>
 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 (&lt; 2 m deep) containing floating and/or submerged
 vegetation. <strong>Frog
 surveys</strong> 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 &gt; 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
 (<em>Pseudacris</em> spp.) were combined together and
 treefrogs (<em>Dryophytes</em> 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 &gt; ~5°C, the second when the
 air temperature was &gt; ~10°C, and the third when the air temperature
 was &gt; ~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. <strong>Water quality surveys</strong> 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 &gt; 15 m apart in water &gt; 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., <em>Typha</em>, 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. <strong>Response variable</strong> 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 (<em>Pseudacris maculata</em>) and Western Chorus Frog (<em>P. triseriata</em>); 4) Gray Treefrog, which included Cope’s Gray Treefrog (<em>Dryophytes chrysoscelis</em>) and Eastern Gray Treefrog (<em>D. versicolor</em>); 5) Green Frog; 6) Northern Leopard Frog; 7) Spring Peeper (<em>Pseudacris crucifer</em>); and 8) Wood Frog (<em>Lithobates sylvaticus</em>). 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 [<em>Anaxyrus fowleri</em>], Mink Frog [<em>Lithobates septentrionalis</em>], Pickerel Frog [<em>L. palustris</em>]; 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 &lt; 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. <strong>Environmental predictors</strong> 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 &lt; 250 m from the point count location, wetland &lt; 2.5 km from the point count location, forest &lt; 2.5 km from the point count location, density of roads &lt; 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 &lt; 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 &lt; r &lt; 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| &gt; 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 &lt; 250 m was based on the <em>coastal wetland inventory </em>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 &lt; 2.5 km were based on the <em>land cover and use</em> 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 <em>population and roads</em> 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 &gt; mean + 3 times the SD. Thus, we removed 28 (1%) specific conductance observations &gt; 977 µS/cm, 42 (2%) ammonium nitrogen observations &gt; 0.15 mg/L, 48 (2%) nitrate nitrogen observations &gt; 1.2 mg/L, 57 (3%) soluble reactive phosphorus observations &gt; 0.14 mg/L, and 33 (1%) road density observations &gt; 5.8 km/ha. <strong>Statistical modeling</strong> 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 (<em>у<sub>i,j,t</sub></em>) at a point count location within a given wetland (<em>j</em>), lake (<em>i</em>), and year (<em>t</em>). The expected occurrence per lake within a given year (<em>µ<sub>i,t</sub></em>)<sub> </sub>for each of the eight species took the form: log(<em>µ<sub>it</sub></em>) = α<sub>i </sub>+ τ<sub>i</sub><em>Τ</em><sub>i</sub>,<sub>j,t</sub> + κ<sub>j</sub> + у<sub>i,t </sub>+ β<sub>1</sub>X<sub>1j</sub> ... β<sub> 11</sub>X<sub>11j</sub> where α is the random lake intercept; <em>T</em> 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<sub>1</sub> – X<sub>11</sub>) 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 (α<sub>i</sub>) had an iCAR structure, where values of α<sub>i</sub> 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 (τ<sub>i</sub>) 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 (τ<sub>i</sub>) 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 (<em>T</em>) 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 (α<sub>i</sub>) 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 (у<sub>i,t</sub>) with an idd, and combined these effects with α and τ. The coefficients of X<sub> 1</sub> to X<sub> 11</sub> 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(<em>µ<sub>it</sub></em>) = α<sub>i </sub>+ τ<sub>i</sub><em>Τ</em><sub>i</sub>,<sub>j,t</sub> + κ<sub>j</sub> + у<sub>i,t</sub> 2) log(<em>µ<sub>it</sub></em>) = α + β<sub>1</sub>X<sub>1j</sub> ... β<sub> 11</sub>X<sub>11j</sub> 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. <strong>References</strong> 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.
description: # Occurrence patterns and trends of frogs in coastal wetlands of the Great
 Lakes call for further habitat restoration
 [https://doi.org/10.5061/dryad.hmgqnk9v8](https://doi.org/10.5061/dryad.hmgqnk9v8) ## Description of the data and file structure We surveyed frogs in coastal wetlands throughout the Great Lakes basin. Coastal wetlands were selected according to CWMP protocol using a stratified, random sampling procedure. 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 &gt; 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. 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 &gt; 15 m apart in water &gt; 5 cm deep. ### Files and variables #### File: CWMP\_MaxOcc.csv **Description:** All of the raw data are in 1 flat table: CWMP_MaxOcc.csv. Each row in the table is the detection (occurrence) or absence of at least one individual during any of the three surveys at a point count location in a particular year for 1 of the 8 species, or species groups, analyzed, along with covariates. The variables are as follows: ##### Variables * site_id: Wetland site number assigned by the Great Lakes Coastal Wetland Monitoring Program ([greatlakeswetlands.org](https://www.greatlakeswetlands.org/Home.vbhtml)). Wetland sites based on polygons in the coastal wetland layer built by the Great Lakes Coastal Wetland Consortium. Referred to as "wetland" in Tozer et al. * point_id: Sample point number within a wetland. Referred to as "point count location" in Tozer et al. * year: Year of observation in format: XXXX. * lat: Latitude of point_id (note that coordinates vary slightly among years for some points). * long: Longitude of point_id (note that coordinates vary slightly among years for some points). * class: 1 of 3 wetland hydrogeomorphic types (i.e., riverine, lacustrine, or barrier protected) assigned by 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. Not used in Tozer et al. * basin: 1 of the 5 Great Lakes (i.e., Erie, Huron, Michigan, Ontario, Superior). Referred to as "lake" in Tozer et al. * region: 1 of 10 regions within the study area (i.e., Lake Erie North, Lake Erie South, Lake Huron East, Lake Huron South, Lake Huron West, Lake Michigan North, Lake Michigan South, Lake Ontario North, Lake Ontario South, Lake Superior). Developed by Tozer et al. because 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 these 10 regions and dropped regions from species-specific analyses if raw (un-modeled) occurrence was &lt; 5%. See Table S1 in Supplemental Material of Tozer et al. for more details. * scale.lake: Detrended, standardized Great Lakes water level (to avoid correlation with year). A dynamic covariate (i.e., it varied annually). Yearly water levels were from the National Oceanic and Atmospheric Administration ([noaa.‌gov](https://www.noaa.gov/)). We used the mean yearly water level from May to July since these months overlapped with surveys in all years and regions. 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. 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 water levels without being confounded with year. See Figure 4 of Tozer et al. for an illustration. * per_forest:** **Percent forest land cover within 2.5 km. A static covariate (i.e., it was the same for all years). Based on data from the “land cover and use” layer built by the Great Lakes Aquatic Habitat Framework ([hub.glahf.org](https://hub.glahf.org/)). See Figure 3 of Tozer et al. for an illustration. * per_wetlg: Percent wetland land cover within 2.5 km. A static covariate (i.e., it was the same for all years). Based on data from the “land cover and use” layer built by the Great Lakes Aquatic Habitat Framework ([hub.glahf.org](https://hub.glahf.org/)). See Figure 3 of Tozer et al. for an illustration. * per_wetsm: Percent wetland land cover within 250 m. A static covariate (i.e., it was the same for all years). Based on data from the “coastal wetland inventory” layer built by the Great Lakes Coastal Wetland Consortium available through the Great Lakes Aquatic Habitat Framework ([hub.glahf.org](https://hub.glahf.org/)). See Figure 3 of Tozer et al. for an illustration. * per_agri: Percent agricultural land cover in the surrounding watershed. A static covariate (i.e., it was the same for all years). Based on data from 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. See Figure 3 of Tozer et al. for an illustration. * per_urban: Percent urban land cover in the surrounding watershed. A static covariate (i.e., it was the same for all years). Based on data from 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. See Figure 3 of Tozer et al. for an illustration. * Roadkm: Road density (km/ha) within 2.5 km. A static covariate (i.e., it was the same for all years). Based on data from the “population and roads” layer built by the Great Lakes Aquatic Habitat Framework ([hub.glahf.org](https://hub.glahf.org/)). See Figure 3 of Tozer et al. for an illustration. * SolReacP: Mean soluable reactive phosphorus (mg/L) for the wetland. Sampled once between mid-June and early-September. Three water quality sampling points (replicates) were surveyed within each monodominant vegetation zone (1–4 zones per wetland). Replicates were &gt; 15 m apart in water &gt; 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). 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. Values within each vegetation zone were averaged, and then zone-level means were averaged to yield a single measure for each wetland. See “Water quality surveys” in Tozer et al. for more details. See Figure 3 of Tozer et al. for an illustration. * nitrateN: Mean nitrate nitrogen (mg/L) for the wetland. Sampled once between mid-June and early-September. Three water quality sampling points (replicates) were surveyed within each monodominant vegetation zone (1–4 zones per wetland). Replicates were &gt; 15 m apart in water &gt; 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). 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. Values within each vegetation zone were averaged, and then zone-level means were averaged to yield a single measure for each wetland. See “Water quality surveys” in Tozer et al. for more details. See Figure 3 of Tozer et al. for an illustration. * ammoniaN: Mean ammonium nitrogen (mg/L) for the wetland; note that the incorrect field name of “ammoniaN” was retained for convenience to match accompanying statistical R code. Sampled once between mid-June and early-September. Three water quality sampling points (replicates) were surveyed within each monodominant vegetation zone (1–4 zones per wetland). Replicates were &gt; 15 m apart in water &gt; 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). 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. Values within each vegetation zone were averaged, and then zone-level means were averaged to yield a single measure for each wetland. See “Water quality surveys” in Tozer et al. for more details. See Figure 3 of Tozer et al. for an illustration. * cond: Mean specific conductance (µS/cm) for the wetland. Sampled once between mid-June and early-September. Three water quality sampling points (replicates) were surveyed within each monodominant vegetation zone (1–4 zones per wetland). Replicates were &gt; 15 m apart in water &gt; 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). Sampled at the mid-depth of the water column at each sampling point using a multi-parameter sonde (standardized to 25°C). Values within each vegetation zone were averaged, and then zone-level means were averaged to yield a single measure for each wetland. See “Water quality surveys” in Tozer et al. for more details. See Figure 3 of Tozer et al. for an illustration. * maxocc: detection/occurrence (1) or absence (0) of at least one individual during any of the three surveys at each point count location in each year. * maxcode:** **a numeric variable (0-3) describing the intensity of calling activity detected during a survey: 0 = no calls heard; 1 = calls not simultaneous, number of individuals can be counted; 2 = calls distinguishable, some simultaneous calling, number of individuals can be reliably estimated; 3 = full chorus, calls continuous and overlapping, number of individuals cannot be estimated. Not used in Tozer et al. * taxa_code: 1 of 8 species analyzed: 1) AMTO: American Toad, 2) BULL: American Bullfrog, 3) CHFR: Chorus Frog, which included Boreal Chorus Frog and Western Chorus Frog; 4) GRFR: Green Frog; 5) TRFR: (Gray) Treefrog, which included Cope’s Gray Treefrog and Eastern Gray Treefrog; 6) NLFR: Northern Leopard Frog, 7) SPPE: Spring Peeper, 8) WOFR: Wood Frog. See Tozer et al. for scientific names. ## Code/software All analytical scripts used to run this analysis are openly available on Birds Canada's Github page, under the repository name CWMP_Frog_Tozer2025. ## Access information Other publicly accessible locations of the data: * The frog survey data are available by request at [greatlakeswetlands.org](https://www.greatlakeswetlands.org). If one wishes to use the frog survey data published in this Dryad dataset, then a formal data request should be submitted at [greatlakeswetlands.org](https://www.greatlakeswetlands.org), although this is not required. Please note that further details regarding the Great Lakes Coastal Wetland Monitoring Program's design, sampling protocols, and other metadata are available at [greatlakeswetlands.org/Sampling-protocols](https://www.greatlakeswetlands.org/Sampling-protocols), including a detailed Quality Assurance Project Plan, and that the program's principal investigators are available to work with researchers to help you better understand and use the dataset. - The “coastal wetland layer” built by the Great Lakes Coastal Wetland Consortium, the “land cover and use” layer built by the Great Lakes Aquatic Habitat Framework, and data from Host et al. (full citations given above) are available at [hub.glahf.org](https://hub.glahf.org/). * Yearly water levels were from the Great Lakes Environmental Research Laboratory of the National Oceanic and Atmospheric Administration available at [noaa.‌gov](https://www.noaa.gov/). - The data cleaning and analysis code for Tozer et al. are available at [github.com/‌BirdsCanada/‌CWMP_Frog_Tozer2025](https://github.com/BirdsCanada/CWMP_Frog_Tozer2025).