Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954–2012: Moving towards recovery?

JGLR-00702; No. of pages: 12; 4C: Journal of Great Lakes Research xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Journal of Great Lakes Research journal homepage: www.elsevier.com/locate/jglr

Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954–2012: Moving towards recovery? Heather M.H. Siersma a,1, Carolyn J. Foley b,2, Carly J. Nowicki a,1, Song S. Qian c,3, Donna R. Kashian a,⁎ a b c

Department of Biological Sciences, Wayne State University, Detroit, MI 48202, USA Illinois-Indiana Sea Grant, Purdue University, West Lafayette, IN 47907, USA Department of Environmental Sciences, The University of Toledo, Toledo, OH 43606, USA

a r t i c l e

i n f o

Article history: Received 19 February 2013 Accepted 18 February 2014 Available online xxxx Communicated by Craig Stow Index words: benthic degradation biomonitoring ecosystem health Hexagenia recovery Saginaw Bay

a b s t r a c t Since at least the 1940s, multiple anthropogenic disturbances to the Laurentian Great Lakes have had detrimental effects on benthic habitats and biota including decimating the environmentally sensitive burrowing mayfly genus Hexagenia around the mid-1950s. While remediation efforts have facilitated recovery of some populations, benthic surveys in Saginaw Bay, Lake Huron in the last 50 years have only occasionally discovered Hexagenia nymphs. Recently, adult Hexagenia swarms have been reported near the bay; therefore, we corroborated the local presence of Hexagenia adults and evaluated the current status of Saginaw Bay Hexagenia nymphs. We quantified adults during mayfly emergence events in 2010 at three Tawas City, Michigan, USA area locations, and found N 17 Hexagenia/m2/site. We quantified nymphs from Ponar grab samples collected at 57 sites in Saginaw Bay between 2009 and 2012, and found 1.5 nymphs/m2 overall with nymphs present at 15.8% of sites sampled, their greatest documented distribution in Saginaw Bay since 1956. Additionally, we mapped bay sediment composition and related sampling site abiotic conditions with both Hexagenia presence and abundance using Zero-Inflated Poisson regression. Model results indicate that the probability of observed Hexagenia absence being true absence is positively related to both sediment sandiness and surficial dissolved oxygen concentration while Hexagenia abundance is greatest where surficial temperatures are ~18.6 °C and is also related to sediment sand content. The documentation of nearby adults and in bay nymphs may indicate the beginning of a Hexagenia return to Saginaw Bay, and, therefore, a possible improvement of the ecosystem's benthic health. © 2014 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved.

Introduction In accordance with the Great Lakes Water Quality Agreement (GLWQA), the United States Environmental Protection Agency (US EPA) and Environment Canada have established indicators of environmental conditions for the International Joint Commission (IJC) (Bertram and Stadler-Salt, 1999) to use to assess and report progress toward restoring and maintaining the health and integrity of the Laurentian Great Lakes (Edsall, 2001). Benthic macroinvertebrates are often used as bioindicators of water quality (e.g., Hilsenhoff, 1987; Kashian and Burton, 2000; Reynoldson et al., 1989); in particular, burrowing mayfly nymphs of the genus Hexagenia (Ephemeroptera: Ephemeridae) are considered highly sensitive to benthic degradation (Edsall et al., 1991).

⁎ Corresponding author at: Biological Sciences Building, Detroit, MI 48201. USA. Tel.: +1 313 577 8052. E-mail addresses: [email protected] (H.M.H. Siersma), [email protected] (C.J. Foley), [email protected] (C.J. Nowicki), [email protected] (S.S. Qian), [email protected] (D.R. Kashian). 1 Tel.: +1 313 577 8107. 2 Tel.: +1 765 494 3601. 3 Tel.: +1 419 530 4230.

Hexagenia are included in the original IJC environmental condition indicator suite (Bertram and Stadler-Salt, 1999), and they have been used extensively in water quality monitoring programs in the Great Lakes (e.g., Edsall, 2001; Fox, 1994). Hexagenia nymphs were historically abundant and widespread in unpolluted, mesotrophic Great Lakes habitats (e.g., Edsall et al., 2005; Schloesser et al., 2014; Surber, 1955; Williamson and Greenbank, 1939; Wood, 1963), typically occupying well-oxygenated (Hunt, 1953), shallow water (b 20 m; Mozley and LaDronka, 1988) environments. High Hexagenia densities have commonly been found in organically enriched sediments (e.g., Carlander et al., 1967; Schneider et al., 1969) composed of various combinations of fine sands, silt loams, soft clays, marls, or firm mucks (Hunt, 1953; Wright and Mattice, 1981). These nymphs spend the majority of their one to about three year lifespans (McCafferty, 1975) burrowed in the top 10 cm of sediments (Charbonneau et al., 1997) prior to emerging to molt into short lived adults (Hunt, 1953). Hexagenia nymphs provide crucial trophic links between the pelagic food web and both the benthic food web and the aquatic detrital loop, serving as a major prey item for many fish species (e.g., Duffy et al., 1987; Schaeffer et al., 2000). Beginning around the mid-1950s, Hexagenia declined in density and distribution following benthic habitat degradation (Beeton, 1961;

http://dx.doi.org/10.1016/j.jglr.2014.03.001 0380-1330/© 2014 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved.

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

2

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

Howmiller, 1971; MDEQ, 2012; Schneider et al., 1969), and they were extirpated from areas including western Lake Erie (Carr and Hiltunen, 1965), Green Bay, Lake Michigan (Howmiller and Beeton, 1971), and Saginaw Bay, Lake Huron (Schloesser et al., 2014; Shannon et al., 1967). These declines have been presumed to be in response to anthropogenic land use changes that resulted in nutrient enrichment and hypoxia (Rasmussen, 1988), increased sedimentation (MDEQ, 2012), and chemical contamination (Edsall et al., 1991; PSC, 2002). Hexagenia nymphs experience high contaminant exposure due to both their direct, prolonged contact with the contaminant enriched sediments and their ingestion of those sediments as they feed indiscriminately on detritus (Smock, 1983). Furthermore, by the 1960s western Lake Erie, Green Bay, and Saginaw Bay had transitioned from mesotrophic to eutrophic conditions (Beeton, 1965; Howmiller and Scott, 1977; Meyers and Takeuchi, 1981; Reynoldson et al., 1989), likely resulting in hypoxic events of increased frequency and duration. With these changes came shifts in the benthic macroinvertebrate community structures in all three locations to domination by organisms tolerant of pollution and low dissolved oxygen (DO) such as Oligochaeta and Chironomidae (e.g., Brandon et al., 1991; Burns, 1985; Howmiller, 1971). However, remediation efforts following GLWQA phosphorus loading and toxin input reduction mandates (e.g., Burns, 1985; Harris et al., 1987) as well as dreissenid mussel (Dreissena polymorpha and D. bugensis) invasion (Leach, 1993; Qualls et al., 2007) have likely facilitated recovery of Hexagenia in some Great Lakes regions (Edsall et al., 2005; Schloesser et al., 2000). For example, following dreissenid (e.g., DeVanna et al., 2011; Hecky et al., 2004) and GLWQA facilitated habitat improvements, Hexagenia have increased in both abundance and distribution in western Lake Erie to near historical levels (e.g., Krieger et al., 2007; Wood, 1963) and Hexagenia have also returned to Green Bay, though at an abundance and a distribution still well below historical levels (Edsall et al., 2005; Williamson and Greenbank, 1939). In Saginaw Bay, however, there has been little documented evidence of Hexagenia population recovery (Edsall et al., 2005; Nalepa et al., 2002) despite not only dreissenid invasion (Nalepa et al., 1995) and GLWQA related improvements but also the stocking of Hexagenia eggs (N600 million) and nymphs (~ 4 million) in Saginaw Bay by the Michigan Department of Natural Resources from 1988 through 1992 (Bryant, 1992). In addition, surveys conducted since 1965 have repeatedly found benthic community dominance in Saginaw Bay by Oligochaeta and Chironomidae (e.g., Barbiero and Tuchman, 2002; Nalepa et al., 2003; Schuytema and Powers, 1966). The IJC has designated Saginaw Bay and its major tributary, the Saginaw River, as an Area of Concern (AOC) under the GLWQA with 12 of 14 possible Beneficial Use Impairments (BUIs) including Degradation of Benthos (MDEQ, 2012). To delist the benthic degradation BUI for the Saginaw River/Bay AOC, its 2001 Remedial Action Plan (RAP) Update set a target restoration condition that Hexagenia be found in the AOC at a density of 30 nymphs/m2 for two consecutive years (PSC, 2002). However, the 2012 Stage 2 RAP for this AOC does not include an assessment of Hexagenia and instead states that the Degradation of Benthos BUI will be considered restored when “…all remedial actions for known contaminated sediment sites with degraded benthos are completed…” (MDEQ, 2012). At this time, this beneficial use remains impaired (MDEQ, 2012). Although the nymphal Hexagenia population in Saginaw Bay has been almost non-detectable in recent decades, there have been anecdotal reports by researchers and local residents since at least 2009 of adult Hexagenia swarms near the bay. The presence of these adults may indicate the beginning of a recovery of Hexagenia within Saginaw Bay or a source population of nymphs in one or more bay tributaries or other water bodies near the bay. Therefore, the objectives of this study were to: 1) corroborate reports of adult Hexagenia swarms near the bay; 2) document recent changes in nymphal Hexagenia distribution and abundance in Saginaw Bay utilizing methods and sampling stations replicated from 1954 to 2001 surveys; and 3) analyze water quality

and sediment factors potentially affecting nymphal distribution and abundance. Such information can be used by resource managers to guide future habitat improvement and other Hexagenia restoration efforts for Saginaw Bay. Methods Study site Saginaw Bay, Lake Huron (Fig. 1) is divided along a broad shoal (Nalepa et al., 1996) into a seasonally stratified (May through October; Smith et al., 1977), meso-oligotrophic (Gardner et al., 1995) outer bay of 13.7 m mean depth (Nalepa et al., 1996) and a well-mixed and rarely thermally stratified (Nalepa et al., 2003), eutrophic (Gardner et al., 1995) to meso-eutrophic (PSC, 2007) inner bay of 5.1 m mean depth (Nalepa et al., 1996). Outer bay substrates are mostly cobble, sand, and silty sand, with only a few clay containing areas, while ~25% of inner bay sediments are silty clays mixed with fine sands and ~75% are primarily sand with gravel and cobble areas (Nalepa et al., 1996; Robbins, 1986). Sand bars lie along both the northwestern and southeastern sides of the inner bay, and between these lies the deepest bay section (~ 14 m; Nalepa et al., 2003) with the finer-grained silt and clay sediments (Robbins, 1986). Nymph sampling We collected Hexagenia nymphs via both a general benthos survey in 2009 and 2010 and an intensive Hexagenia survey in 2012 (Fig. 1). We sampled nine sites in the general benthos survey: five offshore (four in the inner and one in the outer bay) and four nearshore (two each in the inner and outer bays) (Electronic Supplementary Information (ESM) Table S1). Offshore sites were chosen to coincide with the main sampling effort for the National Oceanic and Atmospheric Administration/Great Lakes Environmental Research Laboratory's (NOAA/GLERL) research in Saginaw Bay (Nalepa et al., 2002; Stow et al., this issue), while nearshore sites included two vegetated and two non-vegetated, historically unsampled areas. We sampled the offshore sites once per month from April through November with duplicate full Ponar grabs (0.052 m2 opening) and the nearshore sites twice per month from April through August with single petite Ponar grabs (0.023 m2 opening). Samples were sieved through 500 μm mesh, preserved in 10% formalin with rose Bengal stain, and returned to the lab for additional processing. The 48 sampling sites for the intensive Hexagenia survey were selected by identifying those at which Hexagenia had previously been found and for which latitude and longitude were available; 47 such sites were identified (Edsall et al., 2005; Nalepa et al., 2002; Schneider et al., 1969; Surber, 1954, 1955). As ~ 70% of bay sites found to have Hexagenia presence in the 1950s and N80% of bay nymphal density discovered in the 1950s were located in inner Saginaw Bay, we focused our sampling effort on the inner bay by reducing these 47 sites to the 34 which were in the inner bay (Fig. 1; ESM Table S1; ESM Table S2). We then overlaid the 34 inner bay sites onto the 28 inner bay grid cells that Edsall et al. (2005) established; the sites were collectively located in 14 grid cells. To ensure complete survey coverage of the inner bay, we added to our 34 sampling sites the Edsall et al. (2005) site from each of the remaining 14 grid cells, two of which were in the outer bay near its division from the inner bay and all of which lacked documented prior Hexagenia presence (Fig. 1; ESM Table S2). Each of the 48 intensive survey sites was sampled with triplicate full Ponar grabs one time between 29 May and 7 June 2012 to precede the annual Hexagenia emergence expected in mid- to late-June (Corkum et al., 1997, 2006). Actual sampling locations were generally within 40 m of the corresponding historical sampling site as we attempted to sample in 2012 at the same sites sampled historically, and all were N50 m from each other and from general benthos survey sites such that all nymph sampling sites of the current study were unique locations within

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

3

Fig. 1. Locations of the 2009–2012 and historical (1954–2001; Edsall et al., 2005; Nalepa et al., 2002; Schneider et al., 1969; Surber, 1954, 1955) Hexagenia nymph sampling sites for Saginaw Bay, Lake Huron surveys indicating nymphal presence and absence; also, 2010 adult Hexagenia sampling site and weather monitoring station locations. The dashed black line differentiates the inner and outer bays. The 2012 sampling sites were generally within 40 m of the corresponding historical sampling site as we attempted to sample in 2012 at the same sites sampled historically; therefore the filled and empty circles that indicate historical presence and absence respectively appear as rings around the 2012 sites that overlap them. Additionally, all 2012 sampling sites were N50 m from each other and from 2009–2010 survey sites such that all nymph sampling sites of the current study were unique locations within the bay. All spatial data is projected using the GCS_WGS_1984 coordinate system.

the bay (ESM Table S1; ESM Table S2). Samples were sieved through 500 μm mesh, preserved in 90% ethanol, and returned to the lab for additional processing. In the laboratory, we picked Hexagenia nymphs from individual samples spread into gridded pans. We quantified both whole organisms and heads only, photographed all whole 2012 nymphs using a Unitron Z850 dissecting microscope (Unitron Ltd.) and a Micrometrics 318CU camera (Accu-Scope, Inc.), and, for size-frequency distribution analysis, measured the lengths of all whole 2012 nymphs from top of head to end of abdomen (cerci excluded) using the segmented line tool in the software ImageJ (Schneider et al., 2012). Hexagenia nymphs range from ~ 0.8 mm at hatching to ~ 30 mm at emergence (Edsall, 2001), and their maturation is degree-day dependent; at the mid-Michigan latitude temperatures of shallow Great Lakes waters including Lake Erie and Saginaw Bay, complete nymphal development requires ~ 2 years (Hunt, 1953). Because of this, size-frequency distribution analysis at a given site and time typically reveals two cohorts of nymphs present where Hexagenia have established (Edsall, 2001).

Adult sampling In order to corroborate adult Hexagenia presence near Saginaw Bay at the locations of reputable reports of their sighting, we sampled adult mayfly swarms at three illuminated sites along an easily accessible 2.5 km stretch of the Saginaw Bay shoreline between 15 and 135 m from the water in/near Tawas City, Michigan, USA following methods similar to Corkum et al. (2006) (Fig. 1; ESM Table S1). We sampled one site per night between 22 and 24 June 2010 by placing six circular quadrats of 1 m diameters on matte white sheets on the ground under lights, counting the adult Hexagenia present in each quadrat every 10 min from 10:30 pm until midnight, summing the Hexagenia across the quadrats for a given time-point, and calculating the average (n = 10) Hexagenia/m2/site over the 90 minute sampling period. Circular quadrats were used to minimize counting bias due to edge effects (e.g., Krebs, 1999; Wheater et al., 2011). Wind speed and direction for each sampling event were acquired from the nearest official National Weather Service (NWS) monitoring station (at Alpena County Regional

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

4

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

Airport in Alpena, Michigan, USA) and from the nearer monitoring station at Oscoda-Wurtsmith Airport in Oscoda, Michigan, USA (Fig. 1; ESM Table S1). It is noteworthy that this type of sampling unit (matte white sheets on the ground) may underestimate aerial mayfly densities due to characteristics it does not share with still water surfaces (e.g., Cronin and Marshall, 2011; Kriska et al., 1998; Luo et al., 2007). Sediment mapping During the 2012 Hexagenia survey, we collected ~ 1 L of homogenized sediment from a fourth Ponar grab sample at 47 of the sites; sediment texture analysis using the methods of Bouyoucos (1962) was performed for these by A & L Great Lakes Laboratories, Inc., Fort Wayne, Indiana, USA. Percentage sand, silt, and clay isoline maps for inner Saginaw Bay were created using the Spline (Tension type) Raster Interpolation function of the 3D Analyst extension of ESRI® ArcGIS10.1, Advanced version, Build 3035 (ESRI, 2012). The ArcGIS Spline tool is a local, exact, deterministic interpolator that subsets a study area into neighborhood units for prediction calculation and creates a smoothly contoured surface with minimal overall curvature that passes exactly through the input data points (ESRI, 2012). We used the minimum number of points for each neighborhood (8) and the default weight (0.1) parameters, and when unrealistic values of b0 or N100% were calculated, we truncated these to 0 and 100% respectively. Water quality and sediment texture analysis modeling To determine which water quality and sediment texture variables potentially affect Hexagenia in inner Saginaw Bay, we investigated relationships between the 2012 nymphal presence and abundance data and both the sediment sand, silt, and clay percentages from texture analysis as well as water pH, DO, temperature, conductivity, and depth measurements taken at each 2012 sampling site. To maintain continuity with existing Saginaw Bay data sets (e.g., Nalepa et al., 2002; Stow and Hook, 2013, this issue), water quality variables were measured within two meters of the lake's surface with meters calibrated daily to manufacturers' specifications. These measurements were considered representative of the entire depth profile as the inner bay is generally wellmixed and vertically isothermal; these characteristics can be seen in vertical profiles that reveal little difference between surface and bottom conditions including DO and temperature (e.g., Nalepa et al., 1996; Stow and Hook, 2013; ESM Fig. S1). For analysis of potential relationships between abiotic conditions and nymphal distribution and abundance, we used Zero-Inflated Poisson (ZIP) regression, a model often used in zero-inflation situations (e.g., Cunningham and Lindenmayer, 2005; Lambert, 1992; Martin et al., 2005). When more zeros are observed in a data set than can be fit to a standard distribution, a situation of zero-inflation exists (Heilbron, 1994; Tu, 2002), and zero-inflation can arise in Hexagenia survey data for multiple reasons. First, nymph detection requires tediously searching for small organisms in heterogeneous sediment samples, and it is possible even for careful analysts to miss organisms. Secondly, sediment sample collecting can be ‘hit or miss’; a large, random sample can ameliorate this difficulty, but resource constraints limit sample collection and analysis. Besides such opportunities for obtaining false zeros due to detection probability being less than one (MacKenzie et al., 2002), sampling efforts may also reveal true zeros, those representative of either that Hexagenia do not appear at a site due to the ecological process or effect under study or due to that they by chance have not saturated their entire suitable habitat (Martin et al., 2005); in cases of true zeros, there can be no other outcome but zero regardless of the sampling intensity. Thus, in the process of Hexagenia surveying, an observed zero may reflect either that nymphs were present at the site but were missed in the sampling effort (a false zero) or that they were not present at the site (a true zero); however, which type of zero is the

cause of a given zero observation cannot readily be determined. Conventional probabilistic analyses (e.g., standard Poisson or negative binomial models) assume that true underlying mean abundance is strictly positive and, hence, all zeros in a data set are false zeros (Evans et al., 2000); thus, using conventional statistical methods may lead to biased estimates of rare species population abundance (MacKenzie et al., 2002; Martin et al., 2005). The ZIP model, in contrast, is a probabilistic model designed to accommodate the problem of zero-inflation by using covariate values to quantify the underlying probability of an observed zero being a true zero. This represents advancement in our approaches to enumerating uncommon taxa like Hexagenia and making statistical inferences about the relationships between such counts and their potential ecological causes (Martin et al., 2005). The ZIP model is a mixed model with two components: the zero model (which uses a binomial distribution with a logit link for modeling presence/absence via the likelihood of an observed zero being a true zero) and the abundance model (which uses a Poisson distribution with a log link for modeling abundance counts). Intuitively, we can understand the modeling process using the following four hypothetical steps: 1) Hexagenia count data are first grouped based on covariate values, which typically results in some positive counts and some zero counts existing for a given covariate condition set; 2) the positive counts are used to guess the mean count parameter (e.g., the Poisson model parameter λ) for the specific covariate set; 3) given the mean count parameter estimate, the probability of observing a false zero is determined by the probability model (e.g., in the Poisson model, this probability is e-λ); and 4) the difference between the observed number of zeros and the number of false zeros expected, based on the Poisson or negative binomial model, yields an approximation of the number of true zeros in the data set and thus the probability of an observed zero indicating true Hexagenia absence from a site. In reality, there is no definite information on whether an observed zero is a true zero, and even the covariates we used cannot definitively identify true zeros because whether or not a site provides suitable habitat for Hexagenia likely depends on more than our measured covariates. However, we have no definite knowledge as to the true model determining Hexagenia habitat suitability in Saginaw Bay, so we use the model to determine the probability of an observed zero being a true zero as a function of one or more of our measured covariates. Therefore, using individual water quality and sediment texture data from each sampling site (ESM Table S3), we used a process described in ESM Appendix A1 to search for the most likely model form. Within the computer software R (RCT, 2013), we first used Classification And Regression Tree (CART) analysis (Qian, 2010) to determine which of our eight potential covariates (water depth; the sediment variables %sand, %silt, and %clay; and the water column variables DO, temperature, pH, and conductivity) were important, and then we explored different ZIP model forms using combinations of those covariates with both the zero and the abundance models. The fits of these different models to our observed Hexagenia data were evaluated based on: residual deviances (− 2 * log-likelihood, where loglikelihood is the sum of the log-likelihoods from the abundance and the zero components of the ZIP model); Pr N |z| values; the standard errors of the model determined equation coefficients (α and β below); and comparisons of the former three model outputs plus the determined coefficients themselves to those produced by re-running the models using the negative binomial distribution as a check for overdispersion. The best model was identified via comparisons of these various model outputs as is described in ESM Appendix A1. We used the ZIP model implemented in the R function ‘zeroinfl’ from the package ‘pscl’ (Zeileis et al., 2008). Because the effects of the potential water quality covariates (i.e., DO, temperature, conductivity, pH, and water depth) on the abundance and distribution of Hexagenia are unlikely to be linear, we used a ‘gamma’ model (Qian and Pan, 2006) within our ZIP model; gamma models are a family of flexible models that can be used in situations ranging from monotonic increase, to unimodality, to monotonic decrease. Specifically, for modeling

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

Hexagenia abundance we use log(λi) = β0 + β1xi + β2log(xi), and for modeling the probability of absence we use logit(πi) = α0 + α1xi + α2log(xi); in both, λi is the ith mean count, πi is the probability of absence, and xi is the corresponding covariate value. In predicting Hexagenia abundance, we first used the abundance model to predict the conditional expected abundance (where Hexagenia were already present) and then calculated predicted abundance as the product of conditional expected abundance and probability of presence (1 − πi). Sediment composition covariates (%sand, %silt, and %clay) are strongly correlated and were therefore used in the model one at a time; also, because they were measured as percentages of a total, they were logit transformed for use in the model.

5

NW side of the inner bay; the seventh was located in Wildfowl Bay, the most northeastern part of the inner bay (Fig. 1). Nymphal abundance at individual sites ranged from one to four, translating to a total density of 2.0 nymphs/m2 for the entire sampling area (Table 1). The nymphs ranged in length from 10.9 to 26 mm, and this range represents a size-frequency distribution consistent with the Edsall (2001) description of multiple cohorts; also, all sites with N 1 nymph present had both a small and a large nymph present, further suggesting there are two nymphal cohorts in Saginaw Bay. While we did not observe nymphs larger than 26 mm in our samples, there were substantial amounts of large exuviae on the water in the NW side of the inner bay on 7 June that were likely Hexagenia based solely on their size. This may be evidence that the 2012 emergence occurred just prior to our sampling effort, despite our attempt to sample pre-emergence. The 2012 Michigan spring was the warmest on record, to date, with an average temperature 3.1 °C higher than the 20th century average for March through June (NOAA, 2012); these high temperatures likely accelerated nymphal maturation, prompting an earlier Hexagenia emergence.

Results Nymph sampling, 2009–2010 During the 2009–2010 general benthos survey, we found 12 Hexagenia nymphs, one in 2009 collected post-emergence (on 6 October), and 11 in 2010 collected both pre- and post-emergence (five on 3 June, one on 8 June, four on 10 August, and one on 23 September; Table 1). Two of the 2010 nymphs were from the same northwestern (NW) inner bay offshore site as was the 2009 nymph, while the others were from one outer bay nearshore site (Fig. 1; ESM Table S1); Hexagenia were therefore found at 11% of sites in 2009 and at 22% in 2010 (Table 1). Nymphal abundance at individual sites ranged from one to five, translating to a total density of 0.2 nymphs/m2 in 2009 and 2.2 nymphs/m2 in 2010 (Table 1).

Adult sampling, 2010 We found adult Hexagenia at all three sampling sites and counted a total of 2997 adults across the three sampling days. Two sites had very similar average Hexagenia densities: we counted 17.1 ± 2.7 (mean ± SE) adults/m2 at the most southwestern site on 22 June and 18.1 ± 1.0 adults/m2 at the most northeastern site on 23 June (Fig. 1). Approximately 60% more Hexagenia were observed at the most central site on 24 June: 28.3 ± 1.5 adults/m2 (Fig. 1). Winds of 5.6 km/h from the southeast on 22 June, from the north-northwest on 23 June, and from the west-northwest on 24 June were reported by the NWS monitoring station for 10:54 pm; at the same time, the Oscoda station reported winds from the southwest of 16.7 km/h on 22 June and of 11.1 km/h on 23 June, and winds from the northeast of 35.2 km/h on 24 June (Fig. 1).

Nymph sampling, 2012 During the 2012 intensive Hexagenia survey, we found 15 nymphs at seven (14.6%) of the 48 sites sampled (Fig. 1; Table 1), and five of the seven were locations at which Hexagenia had been found previously (Fig. 1; ESM Table S2). Six of the seven sites were located along the

Table 1 Hexagenia distribution and abundance data from published surveys of Saginaw Bay conducted between 1954 and 2001 for which latitude and longitude data were available, and new data collected between 2009 and 2012 as part of the current study. Reference

Year of survey

Number of locations sampled

Surber, 1954 Surber, 1955 Schneider et al., 1969 Schuytema and Powers, 1966 Shannon et al., 1967 Batchelder, 1973 Shrivastava, 1974 Schaeffer et al., 2000 Nalepa et al., 2002 Nalepa et al., 2002 Nalepa et al., 2002 Nalepa et al., 2002 Nalepa et al., 2002 Nalepa et al., 2002 Nalepa et al., 2002 Nalepa et al., 2002 Nalepa et al., 2002 Edsall et al., 2005 This study This study This study This study

1954 1955 1956 1965 1965 1971 1971 1986 1987 1988 1990 1991 1992 1993 1994 1995 1996 2001 2009 2010 2012 overall

14⁎ 19 52 24 42 28 5 ~51 30 30 10 10 10 10 10 10 10 28 9 9 48 57

Sites with Hexagenia present #

%

9⁎ 8 30 n/a 0 0 0 0 0 0 1 0 0 0 1 0 0 1 1 2 7 9

64⁎ 42 58 n/a 0 0 0 0 0 0 10 0 0 0 10 0 0 4 11 22 15 16

Total Hexagenia abundance (nymphs)

Mean Hexagenia density (nymphs/m2) across sites with Hexagenia present

across all sites

n/a n/a n/a n/a 0 0 0 0 0 0 2 0 0 0 3 0 0 1 1 11 15 27

n/a 150.7 n/a n/a 0 0 0 0 0 0 7.1 0 0 0 10.6 0 0 3.7 1.2 10.4 13.7 9.0

66⁎ 63.5 7.6⁎⁎ 1 0 0 0 0 0 0 0.7 0 0 0 1.1 0 0 0.1 0.2 2.2 2.0 1.5

⁎ as reported in Schloesser et al., 2014. ⁎⁎ reported as 12 in Schloesser et al., 2014.

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

6

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

Fig. 2. Sediment and bathymetry maps representing A) sand, B) clay, and C) silt percentages in Saginaw Bay, Lake Huron, each interpolated from 2012 sediment texture analysis data collected at indicated Hexagenia sampling sites, and D) Saginaw Bay bathymetry, taken from the Michigan Geographic Data Library (www.mcgi.state.mi.us/mgdl). All panels utilize Saginaw River/Bay AOC boundary data from the US EPA Great Lakes National Program Office (www.epa.gov/greatlakes/aoc/saginaw-river/) and include 2012 Hexagenia nymph presence/absence and abundance data. The legend percentage ranges in panels A–C were selected to honor both that the overall ranges across the panels are different (i.e., each has a different maximum) and to present more narrow ranges for where the Hexagenia nymphs were found, particularly at the high sand and the low silt and low clay percentages; they are therefore not uniform across the three sediment type maps. All spatial panel data is projected using the GCS_North_American_1983 coordinate system.

Sediment mapping, 2012

Abiotic water quality and sediment texture analysis modeling, 2012

Spline interpolation calculated that the sediments with higher sand contents in inner Saginaw Bay are located in the southeastern section of the bay and along its NW margin, intermediate sand content is found in Wildfowl Bay, and the least sand is found at the Saginaw River mouth and in the center of the NW side of the bay (Fig. 2A). In contrast, the interpolation indicated the higher clay contents are at the Saginaw River mouth and in the deepest bay section, intermediate clay levels are found in the center of the bay's NW side and in Wildfowl Bay, and the lowest are along the NW margin and in its southeastern section (Figs. 2B; D). Finally, it determined that sediment silt content largely parallels that of clay such that there is a pocket of fine-grained sediments in Wildfowl Bay and also an extensive concentration of silts and clays extending from the Saginaw River mouth and covering approximately one-quarter of the inner bay on its NW side (Figs. 2B; C); much of this fine-grained sediment area overlaps with the deeper areas of inner Saginaw Bay (Fig. 2D). Of the seven sites at which Hexagenia were found in 2012, 57% had sediment textural classifications of sand, 29% of sandy loam, and 14% of loam (Fig. 2; ESM Table S3).

The 2012 sampling sites that lacked Hexagenia generally had broader ranges and higher means of the abiotic variables than did those with Hexagenia present, though mean water depth and sediment silt and clay percentages were higher at sites with Hexagenia than at those without (ESM Table S3). CART analysis suggested DO concentration and both clay and sand percentages as the most important of our eight potential covariates, and evaluation of different ZIP model forms guided by these results culminated in selection of the form with the smallest overall residual deviance (75.32; ESM Appendix A1). This best model used a gamma model for the zero-model and had the following fitted equations: [logit(πi) = -28.98 + 1.71(DO) + 3.05logit(%sand)] (zero component) and [ln(λi) = -801.85 − 22.44(temp) + 416.58ln(temp) − 0.71logit(%sand)] (abundance component), with “logit(X)” representing ln(%X/(100 − %X)), with sand percentage and DO concentration as covariates in the zero model, and with sand percentage and water temperature as covariates in the abundance model (Fig. 3; ESM Appendix A1). Specifically, the zero model showed that the probability of an observed zero being a true zero is directly related to both springtime

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

A

Hexagenia, 2009–2012 The nymphal Hexagenia population in Saginaw Bay may be showing signs of recovery. We observed higher Hexagenia densities in both 2010 and 2012, 2.2 nymphs/m2 and 2.0 nymphs/m2 respectively, than had been documented between 1965 and 2001 (≤ 1.1 nymphs/m2; Table 1). More importantly, we observed a greater nymphal distribution than had been documented between 1965 and 2001 (nymphs at ≤ 1 site). In 2010, our finding of Hexagenia at two of nine sampling sites (~ 22%; Table 1) represents a greater distribution than was found in either the Edsall et al. (b4%; 2005) or the Nalepa et al. (10%; 2002) surveys, both of which included a greater number of sampling sites (28 and 10 respectively) than did our survey (Table 1). Further, this greater distribution observed in 2010 was maintained in 2012. The Edsall et al. (2005) and Nalepa et al. (2002) surveys collectively sampled 37 distinct sites (the Edsall et al. (2005) survey duplicated the site at which Nalepa et al. (2002) found Hexagenia based on a distance of N 50 m distinguishing separate sites), and Hexagenia were found at two of these in total, indicating nymphal distribution at ~ 5.4% of sampling sites; in contrast, our 2012 survey sampled 48 distinct sites including the sites at which Edsall et al. (2005) and Nalepa et al. (2002) found nymphs, and we found nymphs at seven of these, indicating Hexagenia distribution at ~14.6% of sampling sites (Fig. 1; Table 1). The distribution we document for 2012 is therefore more than twice as great as would be expected from the ~ 20% increase in sampling site quantity in our survey as compared to the collective recent surveys in addition to being proportionally greater than the distribution expected from an upscaling of either recent survey individually. Additionally, the distribution pattern we document is consistent with historical 1950s distribution patterns that collectively indicate N70% of the inner bay's nymphs were concentrated on its NW side or in Wildfowl Bay (Schneider et al., 1969; Surber, 1954, 1955; Fig. 4); in the general benthos survey, we found nymphs only on the NW side of the inner bay, and, in the intensive Hexagenia survey, we found them there and in Wildfowl Bay (Fig. 1). The different sampling schemes used in the various historical surveys and in this study plus the large number of zero observations within all of these data sets make statistical comparisons of the densities of Hexagenia nymphs in Saginaw Bay over time challenging; regardless, the trends toward a slight increase in density and toward a substantially expanded distribution (Table 1) are suggestive that a persistent though small population of recovering Hexagenia is present in Saginaw Bay. Additional indicators of potential Hexagenia recovery include nymphs found in the gut contents of round gobies and yellow perch captured in the NW region of inner Saginaw Bay (Staton et al., this issue), increased Hexagenia tusk concentrations discovered in Saginaw Bay sediments in the late-1990s (Schloesser et al. 2014), and our size range evidence of multiple Saginaw Bay nymphal cohorts. Further, the abundances and densities we report may under-represent the bay's actual nymphal population levels due both to our sampling effort likely having followed the 2012 emergence and to the possibilities that young nymphs of b 10.9 mm may not have been retained by the 500 μm sieves used in our sampling effort or may have been missed in our laboratory processing as we did not pick nymphs under magnification. In addition to finding Hexagenia in the inner bay, we also documented nymphs in the outer bay in June and August 2010 at a nearshore site

25% sand 75% sand 90% sand 95% sand 97.5% sand

DO (mg/L)

B

Abundance Model

Abundance (nymphs)

Discussion

Zero Model

Pr(0=structural)

surficial DO concentrations and sediment sand percentage; high sand content and high DO conditions increase the late spring likelihood that observed Hexagenia absence from a site is true absence (Fig. 3A). The abundance model showed that late spring Hexagenia abundance where nymphs are already present is greatest where surficial water temperatures are near 18.6 °C and that it too is directly related to sediment sand content (Fig. 3B).

7

25% sand 75% sand 90% sand 95% sand 97.5% sand o

Temperature ( C) Fig. 3. Panel A — The influence of water column DO concentration and sediment sand content on the probability that an observation of zero nymphs in a sample represents true nymphal absence from the sampling site as indicated by the best fit ZIP model for the 2012 Saginaw Bay, Lake Huron abiotic and Hexagenia condition data. Dots at the top depict observations that were zero; dots at the bottom depict non-zero observations. Lines show the relationship between sediment sand percentage and the probability of a zero observation. At low sand percentages, this probability is zero throughout the DO range; at the highest sand percentages, this probability is near one throughout the DO range; and for sand percentages ≥75, this probability increases with increasing DO concentration. Panel B — The relationships between nymphal abundance, water column temperature, and sediment sand percentage as indicated by the best fit ZIP model. Dots depict nymph counts; lines depict model predicted nymphal abundances at different sediment sand contents. Nymphal abundance remains low throughout the temperature range in lower sand sediments; in higher sand sediments abundance peaks at approximately 18.6 °C.

~5.5 km from the adult swarms corroborated in June 2010 (Fig. 1; ESM Table S1). We think it is unlikely that the adults we observed molted from nymphs of this outer bay site as Kovats et al. (1996) found adult Hexagenia typically travel b3 km when wind speeds are b 10 km/h as they were according to the NWS monitoring station (~5.6 km/h). Additionally, while the higher wind speeds of 11 to 36 km/h reported by the Oscoda monitoring station could have lengthened inland mayfly dispersal, insect takeoff has generally been found to be inhibited by high wind velocities (Johnson, 1969), limiting this possibility. However, it is possible that these adults emerged from elsewhere in the outer bay;

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

8

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

Fig. 4. Historical (1954–1956) Hexagenia nymph distribution and abundance (Schneider et al., 1969; Surber, 1954, 1955) in Saginaw Bay; abundance units are nymphs/m2, and displayed abundance breaks were selected as the 1st quartile, median, 3rd quartile, and maximum values of the non-zero dataset values. All spatial data is projected using the GCS_WGS_1984 coordinate system.

while the outer bay is typically considered too deep to be ideal Hexagenia habitat, nymphs were historically documented at 14 outer bay sites (Schneider et al., 1969; Surber, 1954, 1955; Figs. 1; 4; ESM Table S1). Into both the inner and outer bays, it is also possible that adults from source populations located in bay tributaries or other nearby water bodies may oviposit occasionally which may explain the low nymphal densities documented in the bay over the last 25 years; during a preliminary search of local streams in June 2010, we found three Hexagenia nymphs in Big Creek, a small tributary that enters the NW side of the inner bay, at a sampling site b2 km from the bay's shore (Nowicki unpublished data). Regardless of whether the nymphs we found in the bay originated from adults hatched in the bay or outside of it, the Hexagenia distributions and densities discovered in our surveys are suggestive of the beginning of a Saginaw Bay Hexagenia recovery. The plausibility of a Hexagenia recovery in Saginaw Bay is supported by one that occurred in western Lake Erie. After extirpation in the 1950s, Hexagenia recolonized, and nymphs increased in both abundance and distribution in the 1990s such that by 2004 they had returned to near historical levels (e.g., Krieger et al., 2007; Wood, 1963) of 195 nymphs/m2 mean density and nymphal presence at 70 (Krieger et al.,

2007) to 90% (Edsall et al., 2005) of western Lake Erie sampling sites. However, the Hexagenia recovery in western Lake Erie occurred in fewer than 15 years following the establishment of dreissenid mussels there while already more than 20 years have passed since dreissenids invaded Saginaw Bay (Leach, 1993; Nalepa et al., 1995). Hexagenia recovery in Saginaw Bay may therefore be proceeding more similarly to that in Green Bay which, along with its tributary, the Fox River, is an also an AOC with the Degradation of Benthos BUI (WDNR, 2012). As is the case for Saginaw Bay, Hexagenia numbers in Green Bay remain lower than their historical levels (Edsall et al., 2005; Williamson and Greenbank, 1939) more than 20 years post-dreissenid mussel invasion (Kraft, 1991–1995), though, at a nymphal distribution of 12.5% of sites sampled and a density of 5.65 nymphs/m2 (Edsall et al., 2005), its numbers are higher than they have been since extirpation (e.g., Harris et al., 1987; Howmiller, 1971). Water quality and sediment texture analysis: 2012 mapping and modeling Slow Hexagenia recovery rates in Saginaw Bay may be attributable to several possible factors. Cochran and Kinziger (1997), upon observing

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

that adult Hexagenia along the Fox River, the major Green Bay tributary, did not have large emergence levels even after seven years of presence, proposed that an Allee effect may contribute to inefficient Hexagenia reproduction at low population levels. However, we believe the occurrence of large Hexagenia swarms near Saginaw Bay has the potential to counteract a possible Allee effect. Alternatively, extant benthic degradation could be impeding Hexagenia recovery in the bay. Ecologically, the ZIP model's probability of an observed zero being a true zero is a probability that the sampled benthic habitat is unsuitable for Hexagenia; as such, we believe the zero component of the ZIP model is more immediately important than is the abundance component in the evaluation of abiotic factors potentially affecting Hexagenia recovery in Saginaw Bay. The zero component result of our modeling suggests that surficial DO concentrations and sediment sand percentages may be influencing Hexagenia presence in inner Saginaw Bay; therefore, nonideal levels of these covariates may be impeding Hexagenia population recovery by preventing its establishment or persistence (Fig. 3A). Specifically, our zero model results suggest that as sand percentages increase across the range we measured (from ~25% to ~97.5%), we can be more confident that Hexagenia are truly not present for locations at which we failed to find them in our samples. Further, these results suggest that that confidence is increasingly enhanced across the range of high sand percentages (75% b x b 98%) by increasingly higher surficial springtime DO concentrations in the overlying waters (Fig. 3A). We note, however, that our use of abiotic variable data, particularly that of nonstatic surficial DO concentrations, collected only once at our sampling sites in the late spring limits our modeling analysis to being only suggestive of possible explanations for observed Hexagenia conditions. As such, future research is warranted in order to confirm the ZIP model suggestion that surficial DO concentrations are useful explanatory variables for Hexagenia presence, both in the spring and at other times, prior to drawing specific conclusions about oxygen-related Hexagenia recovery inhibition in Saginaw Bay based on the ZIP modeling results. Moving forward with the above caveat noted, the ZIP model results suggest that, despite meiotrophication of Saginaw Bay as a result of GLWQA mandated phosphorus loading reductions and from high capacity dreissenid filter feeding, DO levels in the inner bay may yet be influencing the inner Saginaw Bay Hexagenia nymph distribution. The surficial DO levels we measured ranged from 11.2 to 14.5 mg/L and, based on vertical DO profiles for Saginaw Bay, were likely in that same range in the bottom waters (ESM Fig. S1). Additionally, Wang et al. (2001) found that Hexagenia are typically able to maintain their burrow DO concentrations at N 75% of that in the overlying water such that in-burrow DO levels were likely N 8 mg/L at our sampling sites, a DO threshold that previous work (Winter et al., 1996) has determined to be important for nymphal Hexagenia survival and growth. Although our in-burrow DO level estimate exceeds both threshold and hypoxic (DO b2 mg/L) conditions and although to our knowledge bottom oxygen depletion in the inner bay has only been rarely documented in the deepest section (N10 m depth) where deployed instrumentation reveals periodic, short-lived oxygen dips during relatively still periods (Stow and Hook, 2013, this issue), it is possible that DO depletion sometimes does occur at the sediment–water interface or in sediment pore spaces. It has long been known that microstratification of lake bottom waters occurs within a few centimeters of the sediment surface due to the high reducing capacity of components of the oxygen-poor sediment (Alsterberg, 1922); therefore, undetected hypoxia may exist at or in the top layers of the sediment where Hexagenia nymphs reside. Further, because the surficial DO concentrations we measured were determined under the cool temperatures of the springtime which increase water's oxygen solubility, it is possible that lower surficial DO concentrations would exist later in the summer season such that both their epibenthic counterparts and the resulting in-burrow DO concentrations would be more likely to be b8 mg/L (ESM Fig. S1). The ZIP model results also suggest that sediment sand content may be influencing the inner Saginaw Bay Hexagenia nymph distribution,

9

Fig. 5. Inner Saginaw Bay sediment sand percentage interpolations representing A) historical (1975–1978; Robbins, 1986) and B) 2012 conditions presented for a consistent area using consistent legend ranges for ease of comparison. Each panel additionally indicates nymphal presence and absence, panel A using historical data (1954–1956; Schneider et al., 1969; Surber, 1954, 1955) and panel B using 2012 data. Both panels utilize Saginaw River/Bay AOC boundary data from the US EPA Great Lakes National Program Office (www.epa.gov/greatlakes/aoc/saginaw-river/), and all spatial data is projected using the GCS_North_American_1983 coordinate system.

and specifically that high sediment sand content may be a limiting factor for Hexagenia presence in the inner bay. Sediments of soft clays, silt loams, and fine but not pure sands have been previously shown to be important determinants of nymphal Hexagenia distribution and abundance (e.g., Hunt, 1953; Lyman, 1943; Wright and Mattice, 1981), and Carlander et al. (1967) noted that nymphal Hexagenia abundance decreased in the early 1960s in Mississippi River pools as sediment sandiness increased from what it had been in the late 1950s. Similarly, the historical Saginaw Bay sites that supported the most Hexagenia, collectively 45% of all nymphs found in the bay in the 1950s, had sediments with sand percentages of b50, and nearly 25% of the total 1950s nymphs inhabited sediments of b2% sand (Robbins, 1986; Schneider et al., 1969; Surber, 1954, 1955; Figs. 4; 5A). While just over half of inner Saginaw Bay had sand percentages N50 between 1975 and 1978 (Robbins, 1986; Fig. 5A) and maximum single site Hexagenia abundance in the 1950s was 646 nymphs/m2 (Surber, 1954; Fig. 4), more than 80% of the inner bay sediment in our study (75% of the same inner bay area mapped by Robbins (1986)) was

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

10

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

found to have N50% sand content (Figs. 2A; 5B), and 2012 maximum single site Hexagenia abundance was just 25.6 nymphs/m2. An increase in Saginaw Bay sand content could be the result of increased sedimentation, a consequence of the inner bay being largely influenced by the Saginaw River (Nalepa et al., 2003) which drains a formerly forested and subsequently primarily agricultural watershed of 22,556 km2, ~ 15% of Michigan's total land area (Brandon et al., 1991; Meyers and Takeuchi, 1981). Regardless of its cause, increased sediment sandiness in Saginaw Bay could be problematic for Hexagenia recovery as nymphs cannot burrow into pure sand, which lacks the adhesiveness of finer particles, and they have difficulty burrowing into coarse sands even when these are mixed with finer-grained sediments (Lyman, 1943). Also, Eriksen (1963) found that Hexagenia oxygen consumption varies with particle size such that in non-ideal substrates with low penetrability the nymphs must consume more oxygen to burrow, and Wright and Mattice (1981) found that Hexagenia have the lowest burrowing times in substrates of 6 to 7 φ median diameter (i.e., silts); although our study did not discriminate between sand particles of different diameters, based on these prior studies, an increase in sand percentages in Saginaw Bay sediments could necessitate higher DO requirements for a potentially recovering Hexagenia population. In addition to non-ideal levels of sediment sand and surficial DO as potential habitat concerns for Hexagenia, Saginaw Bay has a history of chemical contamination by organic pesticide residues (Fox, 1994; PSC, 2002), industrial petroleum compounds, and metals (Edsall et al., 1991) that has contributed to its listing as not only an AOC but also as part of a US EPA Superfund site (USEPA, 2011). Recent studies indicate that, although water quality has improved in the bay, its sediments remain contaminated and could still be negatively affecting Hexagenia recovery (e.g., Barbiero and Tuchman, 2002; MDEQ, 2012); however, recent and ongoing contaminant reduction efforts in the bay (e.g., MDEQ, 2008; USEPA, 2011) may be contributing to the larger spatial distribution of nymphs observed in this study. Importantly, this study did not evaluate the impacts to Hexagenia of contaminants or other environmental factors such as the invasive dreissenid mussels that established in the bay following Hexagenia extirpation and may be playing a role, either positive or negative, in Saginaw Bay Hexagenia recovery.

Conclusions Our findings document the presence of Hexagenia in Saginaw Bay for three recent years at an overall density of 1.5 nymphs/m2 (Table 1) and the presence of multiple adult Hexagenia swarms within 135 m of its shore (Fig. 1). To our knowledge, extensive bay sampling efforts prior to this study have not found nymphs at more than a single location or found N1.1 nymphs/m2 per sampling effort since 1956 (Table 1); the 27 nymphs found collectively in 2009, 2010, and 2012 at 15.8% of our sampling sites along both the NW side of the inner bay and in Wildfowl Bay represent the greatest abundance and density and broadest distribution of Hexagenia in Saginaw Bay documented since its extirpation (Fig. 1; Table 1). Although the nymphal abundances, densities, and distributions we report remain critically low as compared to historical levels and to the density formerly required by the IJC for delisting the Saginaw River/Bay AOC benthic degradation BUI, the proportion of sites with Hexagenia present appears to be increasing (Table 1). The confirmation of nearby adults and the presence and potentially expanding distribution of nymphs may indicate the beginning of a Hexagenia return to Saginaw Bay and, therefore, a possible improvement of bay ecosystem health particularly along the inner bay's NW side. Because a larger recovery has not been observed in the last 25 years, it is likely that bay sediments may still harbor unsuitable habitat conditions for Hexagenia. Whether a lack of suitability is related to non-ideal DO concentrations, high sand percentages, persistent chemical contamination, other factors, or a combination of factors, we conclude that additional remediative work to improve benthic conditions is likely needed to restore bay

integrity and facilitate and sustain a large-scale Hexagenia recovery in Saginaw Bay. Acknowledgments This research was sponsored by the National Oceanic and Atmospheric Administration Center for Sponsored Coastal Ocean Research and Wayne State University. Additionally, we thank the US EPA for the use of the US EPA R/V Mudpuppy II research vessel, specifically Joe Bonem, Stacy Coullard, Rosanne Ellison, MaryBeth Giancarlo, Diane Mally, William Murray, and Marc Tuchman. We also thank Bethany Coggins, Brittanie Dabney, Jacob Dombrowski, Jared Militello, Arthur Ostaszewski, Steve Pothoven, and Charles Roswell for assistance with field collections, Jake Dombrowski, Suzy Lyttle, Cody Narlock, Brenna Stow, Michelle Walsh, and Dan Wessel for help with laboratory sample and data processing, and Zach Feiner, Larry Lemke, and F. Gianluca Sperone for assistance with ArcGIS work. We give special thanks to Charlie Bauer, Art Ostaszewski, and Craig Stow for their intellectual support, and to the anonymous reviewers whose constructive comments contributed to improving the manuscript. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.jglr.2014.03.001. References [ESRI] Environmental Systems Resource Institute, 2012. ArcGIS 10.1. ESRI Redlands, California. [MDEQ] Michigan Department of Environmental Quality, 2008. Guidance for delisting Michigan's Great Lakes Areas of Concern, revised. Report MI/DEQ/WB-06-001. MDEQ, Lansing, Michigan. [MDEQ] Michigan Department of Environmental Quality, 2012. Stage 2 Remedial Action Plan for the Saginaw River/Bay Area of Concern. MDEQ. Office of the Great Lakes, Lansing, Michigan. [NOAA] National Oceanic and Atmospheric Administration, 2012. NOAA National Climatic Data Center. State of the Climate: National Overview for Annual 2012 Published online December 2012, most recently retrieved 10 January 2013 from http://www.ncdc. noaa.gov/sotc/2012/13. [PSC] Public Sector Consultants, Inc., 2002. Targeting Environmental Restoration in the Saginaw River/Bay Area of Concern (AOC): 2001 Remedial Action Plan Update. PSC, Inc., Lansing, Michigan. [PSC] Public Sector Consultants, Inc., 2007. Phosphorus Policy Advisory Committee Final Report. PSC, Inc., Lansing, Michigan. [RCT] R Core Team, 2013. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria (URL http://www.R-project. org/). [USEPA] United States Environmental Protection Agency, 2011. Community Involvement Plan: Tittabawassee River, Saginaw River and Bay Site, Midland/Saginaw/Bay City Region, Michigan. US EPA, EPA Region 5 Chicago, Illinois. Alsterberg, G., 1922. Die respiratorischen mechanism der Tubificiden. K. Fysiogr. Sallsk. Handl. N.F. 33, 1–175. Barbiero, R.P., Tuchman, M.L., 2002. Results from GLNPO's biological open water surveillance program of the Laurentian Great lakes, 1999. EPA-905/R-02-001US EPA Great Lakes National Program Office, Chicago Illinois. Batchelder, T.L., 1973. Saginaw Bay Ecological Survey, 1971. Internal ReportDow Chemical Company, Midland, Michigan. Beeton, A.M., 1961. Environmental changes in Lake Erie. Trans. Am. Fish. Soc. 90, 153–159. Beeton, A.M., 1965. Eutrophication of the St. Lawrence Great Lakes. Limnol. Oceanogr. 10, 240–254. Bertram, P., Stadler-Salt, N., 1999. State of the Lakes Ecosystem Conference 1998: Selection of indicators for Great Lakes basin ecosystem health, Version 3, Draft for Review. US EPA, Chicago, Illinois, & Environment Canada, Burlington, Ontario. Available online, most recently retrieved 14 August 2013 from http://www. epa.gov/solec/archive/2000/Selection_of_Indicators_Version_4_(FULL).pdf. Bouyoucos, J.G., 1962. Hydrometer method improved for making particle size analyses of soils. Agron. J. 54, 464–465. Brandon, D.L., Lee, C.R., Simmers, J.W., Tatem, H.E., Skogerboe, J.G., 1991. Information Summary, Area of Concern: Saginaw River and Saginaw Bay. Miscellaneous Paper EL-91-7US Army Engineer Waterways Experiment Station, Vicksburg, Mississippi. Bryant, W.C., 1992. Accelerating the return of burrowing mayflies to Saginaw Bay, Michigan. Dingell-Johnson Annual Rept., Project No. F-53-R-8, Study 456. MDNR Fisheries Division, Ann Arbor, Michigan. Burns, N.M., 1985. Erie, the lake that survived. Rowman & Allanheld Publishers, Totawa, New Jersey. Carlander, K.D., Carlson, C.A., Gooch, V., Wenke, T.L., 1967. Populations of Hexagenia mayfly naiads in Pool 19, Mississippi River, 1959–1963. Ecology 48, 873–878.

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx Carr, J.F., Hiltunen, J.K., 1965. Changes in the bottom fauna of western Lake Erie from 1930 to 1961. Limnol. Oceanogr. 10, 551–569. Charbonneau, P., Hare, L., Carignan, R., 1997. Use of X-rays images and a contrasting agent to study the behavior of animals in soft sediments. Limnol. Oceanogr. 42, 1823–1828. Cochran, P.A., Kinziger, A.P., 1997. Hexagenia bilineata (Ephemeroptera: Ephemeridae) persists at low levels of abundance in the lower Fox River, Wisconsin. Great Lakes Entomol. 30, 89–92. Corkum, L.D., Ciborowski, J.J.H., Poulin, R., 1997. Effects of emergence date and maternal size on egg development and sizes of eggs and first-instar nymphs of a semelparous aquatic insect. Oecologia 111, 69–75. Corkum, L.D., Ciborowski, J.J.H., Dolan, D.M., 2006. Timing of Hexagenia (Ephemeridae: Ephemeroptera) mayfly swarms. Can. J. Zool. 84, 1616–1622. Cronin, T.W., Marshall, J., 2011. Patterns and properties of polarized light in air and water. Philos. Trans. R. Soc. B 366, 619–626. Cunningham, R.B., Lindenmayer, D.B., 2005. Modeling count data of rare species: Some statistical issues. Ecology 86, 1135–1142. DeVanna, K.M., Armenio, P.M., Barrett, C.A., Mayer, C.M., 2011. Invasive ecosystem engineers on soft sediment change the habitat preferences of native mayflies and their availability to predators. Freshw. Biol. 56, 2448–2458. Duffy, W.G., Batterson, T.R., McNabb, C.D., 1987. The St. Mary's River, Michigan: an ecological profile. U.S. Fish Wildl. Serv. Biol. Rep. 85(7.10). U.S. Department of the Interior Fish and Wildlife Service, Washington, D.C. Edsall, T.A., 2001. Burrowing mayflies (Hexagenia) as indicators of ecosystem health. Aquat. Ecosyst. Health 4, 283–292. Edsall, T.A., Manny, B.A., Schloesser, D.W., Nichols, S.J., Frank, A.M., 1991. Production of Hexagenia limbata nymphs in contaminated sediments in the upper Great Lakes connecting channels. Hydrobiologia 219, 353–361. Edsall, T.A., Bur, M.T., Gorman, O.T., Schaeffer, J.S., 2005. Burrowing mayflies as indicators of ecosystem health: Status of populations in western Lake Erie, Saginaw Bay, and Green Bay. Aquat. Ecosyst. Health 8, 107–116. Eriksen, C.H., 1963. The relation of oxygen consumption to substrate particle size in two burrowing mayflies. J. Exp. Biol. 40, 447–453. Evans, M., Hastings, N.A.J., Peacock, B., 2000. Statistical Distributions, third edition. John Wiley and Sons, Hoboken, New Jersey. Fox, G.A. (Ed.), 1994. Bioindicators as a measure of success for virtual elimination of persistent toxic substances. International Joint Commission, Windsor, Ontario. Gardner, W.S., Cavaletto, J.F., Johengen, T.H., Johnson, J.R., Heath, R.T., Cotner Jr., J.B., 1995. Effects of the zebra mussel, Dreissena polymorpha, on nitrogen dynamics in Saginaw Bay, Lake Huron. J. Great Lakes Res. 21, 529–544. Harris, H.J., Sager, P.E., Yarbrough, C.J., Day, H.J., 1987. Evolution of water resource management: A Laurentian Great Lakes case study. Int. J. Environ. Stud. 29, 53–70. Hecky, R.E., Smith, R.E.H., Barton, D.R., Guildford, S.J., Taylor, W.D., Charlton, M.N., Howell, T., 2004. The nearshore phosphorus shunt: a consequence of ecosystem engineering by dreissenids in the Laurentian Great Lakes. Can. J. Fish. Aquat. Sci. 61, 1285–1293. Heilbron, D.C., 1994. Zero-altered and other regression models for count data with added zeros. Biom. J. 36, 531–547. Hilsenhoff, W.L., 1987. An improved biotic index of organic stream pollution. Great Lakes Entomol. 20, 31–39. Howmiller, R.P., 1971. The benthic macrofauna of Green Bay, Lake Michigan. University Microfilms, Inc., Ann Arbor, Michigan. Howmiller, R.P., Beeton, A.M., 1971. Biological evaluation of environmental quality, Green Bay, Lake Michigan. Res. J. Water Pollut. C 43, 123–133. Howmiller, R.P., Scott, M.A., 1977. An environmental index based on relative abundance of oligochaete species. Res. J. Water Pollut. C 49, 809–815. Hunt, B.P., 1953. The life history and economic importance of a burrowing mayfly, Hexagenia limbata, in southern Michigan lakes. Bulletin Institute of Fisheries Research. No. 4Michigan Department of Conservation. Johnson, C.G., 1969. Migration and dispersal of insects by flight. Methuen, London, England. Kashian, D.R., Burton, T.M., 2000. A comparison of macroinvertebrates of two Great Lakes coastal wetlands: Testing potential metrics for an index of ecological integrity. J. Great Lakes Res. 26, 460–481. Kovats, Z.E., Ciborowski, J.J.H., Corkum, L.D., 1996. Inland dispersal of adult aquatic insects. Freshw. Biol. 36, 265–276. Kraft, C., 1991–1995. Zebra Mussel Updates. University of Wisconsin-Sea Grant Advisory Services, Green Bay. Wisconsin 8–24. Krebs, C.J., 1999. Ecological methodology, second edition. Addison-Wesley Educational Publishers, Inc., Menlo Park, California. Krieger, K.A., Bur, M.T., Ciborowski, J.J.H., Barton, D.R., Schloesser, D.W., 2007. Distribution and abundance of burrowing mayflies (Hexagenia spp.) in Lake Erie, 1997–2005. J. Great Lakes Res. 33, 20–33. Kriska, G., Horváth, G., Andrikovics, S., 1998. Why do mayflies lay their eggs en masse on dry asphalt roads? Water-imitating polarized light reflected from asphalt attracts Ephemeroptera. J. Exp. Biol. 201, 2273–2286. Lambert, D., 1992. Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics 34, 1–14. Leach, J.H., 1993. Impacts of the zebra mussel (Dreissena polymorpha) on water quality and fish spawning reefs in western Lake Erie. In: Nalepa, T.F., Schloesser, D.W. (Eds. ), Zebra Mussels: Biology, Impacts, and Control. Lewis Publishers/CRC Press, Boca Raton, Florida, pp. 381–397. Luo, Y., Zhao, Y., Li, X., Wu, T., Zhao, L., 2007. Research and application of multi-angle polarization characteristics of water body mirror reflection. Sci. China Ser. D 50, 946–952. Lyman, F.E., 1943. Swimming and burrowing activities of mayfly nymphs of the genus Hexagenia. Ann. Entomol. Soc. Am. 36, 250–256.

11

MacKenzie, D.I., Nichols, J.D., Lachman, G.B., Droege, S., Royle, J.A., Langtimm, C., 2002. Estimating site occupancy rates when detection probabilities are less then one. Ecology 83, 2248–2255. Martin, T.G., Wintle, B.A., Rhodes, J.R., Kuhnert, P.M., Field, S.A., Low-Choy, S.J., Tyre, A.J., Possingham, H.P., 2005. Zero tolerance ecology: improving ecological inference by modelling the source of zero observations. Ecol. Lett. 8, 1235–1246. McCafferty, W.P., 1975. The burrowing mayflies of the United States (Ephemeroptera: Ephemeroidea). Trans. Am. Entomol. Soc. 101, 447–504. Meyers, P.A., Takeuchi, N., 1981. Environmental changes in Saginaw Bay, Lake Huron recorded by geolipid contents of sediments deposited since 1800. Environ. Geol. 3, 257–266. Mozley, S.C., LaDronka, R.M., 1988. Ephemera and Hexagenia (Ephemeridae: Ephemeroptera) in the Straits of Mackinac, 1955–56. J. Great Lakes Res. 14, 171–177. Nalepa, T.F., Wojcik, J.A., Fanslow, D.L., Lang, G.A., 1995. Initial colonization of the zebra mussel (Dreissena polymorpha) in Saginaw Bay, Lake Huron: population recruitment, density, and size structure. J. Great Lakes Res. 21, 417–434. Nalepa, T.F., Fahnenstiel, G.L., McCormick, M.J., Johengen, T.H., Lang, G.A., Cavaletto, J.F., Goudy, G., 1996. Physical and chemical variables of Saginaw Bay, Lake Huron in 1991–1993. NOAA Technical Memorandum GLERL-91. Great Lakes Environmental Research Laboratory, Ann Arbor, Michigan. Nalepa, T.F., Fanslow, D.L., Lansing, M.B., Lang, G.A., Ford, M., Gostenik, G., Hartson, D. J., 2002. Abundance, biomass, and species composition of benthic macroinvertebrate populations in Saginaw Bay, Lake Huron, 1987–96. NOAA Technical Memorandum GLERL-122. Great Lakes Environmental Research Laboratory, Ann Arbor, Michigan. Nalepa, T.F., Fanslow, D.L., Lansing, M.B., Lang, G.A., 2003. Trends in the benthic macroinvertebrate community of Saginaw Bay, Lake Huron, 1987–1996: responses to phosphorous abatement and the zebra mussel, Dreissena polymorpha. J. Great Lakes Res. 29, 14–33. Qian, S.S., 2010. Environmental and ecological statistics with R. Taylor and Francis Group LLC, Boca Raton, Florida. Qian, S.S., Pan, Y., 2006. Historical soil total phosphorus concentration in the Everglades. In: Burk, A.R. (Ed.), Focus on Ecological Research. Nova Science Publishers Inc., Hauppauge, New York, pp. 131–150. Qualls, T.M., Dolan, D.M., Reed, T., Zorn, M.E., Kennedy, J., 2007. Analysis of the impacts of the zebra mussel, Dreissena polymorpha, on nutrients, water clarity, and the chlorophyll–phosphorous relationship in lower Green Bay. J. Great Lakes Res. 33, 617–626. Rasmussen, J.B., 1988. Habitat requirements of burrowing mayflies (Ephemeridae: Hexagenia) in lakes, with special reference to the effects of eutrophication. J. N. Am. Benthol. Soc. 7, 51–64. Reynoldson, T.B., Schloesser, D.W., Manny, B.A., 1989. Development of benthic invertebrate objective for mesotrophic Great Lakes waters. J. Great Lakes Res. 15, 669–686. Robbins, J.A., 1986. Sediments of Saginaw Bay, Lake Huron: elemental composition and accumulation rates. Special Report No. 102 of the Great Lakes Research Division. Great Lakes Environmental Research Laboratory, Ann Arbor, Michigan. Schaeffer, J.S., Diana, J.S., Haas, R.C., 2000. Effects of long-term changes in the benthic community on yellow perch in Saginaw Bay, Lake Huron. J. Great Lakes Res. 26, 340–351. Schloesser, D.W., Krieger, K.A., Ciborowski, J.J.H., Corkum, L.D., 2000. Recolonization and possible recovery of burrowing mayflies (Ephemeroptera: Ephemeridae: Hexagenia spp.) in Lake Erie of the Laurentian Great Lakes. J. Aquat. Ecosyst. Stress. Recover. 8, 125–141. Schloesser, D.W., Robbins, J.A., Matisoff, G., Nalepa, T.F., Morehead, N.R., 2014. A 200 year chronology of burrowing mayflies (Hexagenia spp.) in Saginaw Bay. J. Great Lakes Res 40, 80–91. Schneider, J.C., Hooper, F.F., Beeton, A.M., 1969. The distribution and abundance of benthic fauna in Saginaw Bay, Lake Huron. Proc. 12th Conf. Great Lakes Res., Internat. Assoc. Great Lakes Res. The University of Michigan, Ann Arbor, Michigan, pp. 80–90. Schneider, C.A., Rasband, W.S., Eliceiri, K.W., 2012. NIH Image to ImageJ: 25 years of image analysis. Nat. Methods 9, 671–675. Schuytema, G.S., Powers, R.E., 1966. The distribution of benthic fauna in Lake Huron. Proc. 9th Conf. Great Lakes Res., Great Lakes Res. Div. Pub. No. 15. The University of Michigan, Ann Arbor, Michigan, pp. 155–163. Shannon, E.S., Alexander, H.C., Fetterolf, C.M., Newton, M.E., 1967. Saginaw Bay Biological Survey, July–August, 1965. Internal Report. Dow Chemical Company, Midland, MI. Shrivastava, H., 1974. Macrobenthos of Lake Huron. Technical Report No. 449 of the Fisheries Research Board of Canada. Canada Centre for Inland Waters, Great Lakes Biolimnology Lab, Burlington, Ontario. Smith, V.E., Lee, K.W., Filkins, J.C., Hartwell, K.W., Rygwelski, K.R., Townsend, J.M., 1977. Survey of chemical factors in Saginaw Bay (Lake Huron). US EPA Ecological Research Series, EPS-600/3-77-125US EPA Environmental Research Laboratory Office of Research and Development, Duluth, Minnesota. Smock, L.A., 1983. The influence of feeding habits on whole-body metal concentrations in aquatic insects. Freshw. Biol. 13, 301–311. Staton, J.M., Roswell, C.R., Fielder, D.G., Thomas, M.V., Pothoven, S.A., Höök, T.O., 2014. Evaluating differences in condition of yellow perch in Saginaw Bay, Lake Huron (1970–2011). J. Great Lakes Res. (in this issue). Stow, C.A., Hook, T., Eds., with Beletsky, D., Beletsky, R., Burtner, A., Bredin, J.H., Cavaletto, J., Cha, Y., De Marchi, C., DePinto, J.V., Dyble, J., Fanslow, D., Francoeur, S., Gossiaux, D., Hawley, N., He, C., Johengen, T.H., Kashian, D.R., Kaplowitz, M.D., Keiper, W., Lavrentyev, P.J., Lupi, F., Millie, D.F., Morehead, N., Nalepa, T.F., Newcomb, T.J., Palladino, D., Peacor, S.D., Pothoven, S.A., Redder, T., Selzer, M., Vanderploeg, H.A., Verhamme, E., Winslow, K., 2013. Saginaw Bay Management Report. NOAA Technical

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001

12

H.M.H. Siersma et al. / Journal of Great Lakes Research xxx (2014) xxx–xxx

Memorandum GLERL-160. NOAA, Great Lakes Environmental Research Laboratory. Ann Arbor MI. 48108. Stow, C.A., Dyble, J., Kashian, D.R., Johengen, T.H., Winslow, K.P., Peacor, S.D., Francoeur, S. N., Burtner, A.M., Palladino, D., Morehead, N., Gossiaux, D., Cha, Y., Qian, S.S., Miller, D., 2014. Phosphorus Targets and Eutrophication Objectives in Saginaw Bay: A 35 Year Assessment. J. Great Lakes Res. (in this issue). Surber, E.W., 1954. Results of a biological survey of Saginaw Bay, July 16, 17, 1954. Michigan Water Resources Commission, Lansing, Michigan. Surber, E.W., 1955. Results of a biological survey of Saginaw Bay, June 1955. Michigan Water Resources Commission, Lansing, Michigan. Tu, W., 2002. Zero-inflated data. In: El-Shaarawi, A.H., Piegorsch, W.W. (Eds.), Encyclopedia of Environmetrics. John Wiley and Sons, Chichester, West Sussex, England, pp. 2387–2391. Wang, F.Y., Tessier, A., Hare, L., 2001. Oxygen measurements in the burrows of freshwater insects. Freshw. Biol. 46, 317–327. Wheater, C.P., Bell, J.R., Cook, P.A., 2011. Practical field ecology: a project guide. John Wiley and Sons, Chichester, West Sussex, England.

Williamson, B.L., Greenbank, J., 1939. Investigations of the Pollution of the Fox and East Rivers and of Green Bay in the Vicinity of the City of Green Bay. Wisconsin State Committee on Water Pollution and State Board of Health, Madison, Wisconsin. Winter, A., Ciborowski, J.H., Reynoldson, T.B., 1996. Effects of chronic hypoxia and reduced temperature on survival and growth of burrowing mayflies, Hexagenia limbata (Ephemeroptera: Ephemeridae). Can. J. Fish. Aquat. Sci. 53, 1565–1571. Wisconsin Department of Natural Resources, 2012. Stage 2 Remedial Action Plan Update for the Lower Green Bay and Fox River Area of Concern. WDNR Office of the Great Lakes, Madison, Wisconsin. Wood, K.G., 1963. The bottom fauna of western Lake Erie, 1951–1952. Proc. 6th Conf. Great Lakes Res., Great Lakes Res. Div. Pub. No. 10. The University of Michigan, Ann Arbor, Michigan, pp. 258–265. Wright, L.L., Mattice, J.S., 1981. Substrate selection as a factor in Hexagenia distribution. Aquat. Insects 3, 13–24. Zeileis, A., Kleiber, C., Jackman, S., 2008. Regression Models for Count Data in R. J. Stat. Softw. 27 (URL http://www.jstatsoft.org/v27/i08/).

Please cite this article as: Siersma, H.M.H., et al., Trends in the distribution and abundance of Hexagenia spp. in Saginaw Bay, Lake Huron, 1954– 2012: Moving towards recovery? J Great Lakes Res (2014), http://dx.doi.org/10.1016/j.jglr.2014.03.001