Introduction

Estuaries are dynamic ecosystems at the intersection of marine, freshwater, and terrestrial habitats that experience large fluctuations in environmental variables over tidal, seasonal, and interannual scales (McLusky and Elliott 2004). These dynamic habitat conditions affect fish physiology (Fry 1971), growth rates (Sogard 1992), swimming performance (Webb 1978; Johnson et al. 1993; Yetsko and Sancho 2015), and trophic ecology (Nelson et al. 2015; Whitney et al. 2018), which may translate to shifts in ecological community structure (Thiel et al. 1995, Davis et al. 2019). The high frequency and magnitude of fluctuations in salinity and temperature limit the total number of species that can inhabit estuaries (Thiel et al. 1995); yet, due to an influx of nutrients from multiple ecosystems, estuaries are also among the most productive habitats on the planet (Correll 1978), supporting a large biomass of fish species that are tolerant of environmental variability (Houde and Rutherford 1993; Whitfield 2016). Estuaries provide crucial staging and migration routes for anadromous fishes (e.g., Pacific salmon, Oncorhynchus spp.; Neilson et al. 1985, Thorpe 1994, Simmons et al. 2013) and are ideal nursery grounds for many ecologically, commercially, and culturally important marine fish and invertebrate species (Beck et al. 2001, Gillanders et al. 2003; Minello et al. 2003). Approximately 8% of annual global commercial marine fish catch is estuarine in origin and just over 50% of fishery harvests in the USA are estuarine for at least part of their life history (Houde and Rutherford 1993).

While estuaries provide critical habitat for diverse species, they are simultaneously vulnerable to anthropogenic disturbance and rapid environmental change. Estuaries rank high among aquatic habitats most affected by coastal development, pollution (McLusky and Elliott 2004), and global climate change (Elliott et al. 2019). High-latitude estuaries are projected to become increasingly affected by watershed changes resulting from global warming (Arendt et al. 2009), including increased precipitation falling as rain (McAfee et al. 2013) and significant loss of glacial coverage (Radić et al. 2013) that translates to a higher volume of glacier runoff into streams (O’Neel et al. 2014). Early to intermediate stages of glacial recession can generate new stream and estuarine habitats (Pitman et al. 2020) and initiate changes in physical and chemical properties of existing watersheds (Hood and Berner 2009, Moore et al. 2009). Runoff from glacial rivers also mediates the delivery of bioavailable carbon, nitrogen, and phosphorous (Hood and Scott 2008) to delta and estuarine habitats, with potential effects on nutrient sources for nearshore marine consumers (e.g., Arimitsu et al. 2016, Whitney et al. 2018). Freshening and stratification of coastal shelf waters has been linked to changes in abundance, distribution, and seasonal cycles of phytoplankton and zooplankton, which affect energy transfer to higher trophic levels (Greene et al. 2008).

Along the Gulf of Alaska, changes in the transport of freshwater, sediments, and nutrients to estuary habitats have been driven by rapid glacier retreat (O’Neel et al. 2015), with an estimated annual glacier mass loss of ~ 75 Gt over the period 1994–2013 (Larsen et al. 2015). While only 9% of the individual watersheds bordering the Gulf of Alaska are classified as glacial runoff systems, they make up 96% of the watershed land area (Sergeant et al. 2020) and contribute an estimated 38% of the mean annual freshwater runoff into the Gulf of Alaska (Beamer et al. 2016). As glaciers continue to melt and recede, watersheds along the Gulf of Alaska will transform (Sergeant et al. 2020), with consequences for downstream ecosystems. The coastal margin along the Gulf of Alaska has served as a flagship location for investigating aquatic community dynamics in response to changes in glacierized (ice-covered) landscapes (e.g., Chapin et al. 1994, Fastie 1995, Milner et al. 2000, Arimitsu et al. 2016). For instance, studies of ecological succession in newly formed glacial streams have tracked the colonization of aquatic invertebrates and fishes over decades of glacier retreat in Glacier Bay, Alaska (Milner et al. 2008, Milner et al. 2011). In the nearshore marine environment, gradients in temperature and turbidity with distance from tidewater glaciers affect zooplankton and forage fish community dynamics (Arimitsu et al. 2012, 2016). Yet, there remain significant gaps in foundational knowledge of ecological community structure in delta and estuarine habitats, including the temporal and spatial dynamics of estuarine fish communities and how they are affected by meltwater from high elevation, land-terminating glaciers (O’Neel et al. 2015, Pitman et al. 2020). Over the coming decades, deglaciation has potential to lead to a temporary increase, followed by a long-term decline of glacial influence on estuaries (Pitman et al. 2020). Therefore, quantifying the current variability in estuarine fish community structure in regions predicted to be most affected by glacial volume reduction is crucial for anticipating future impacts.

This study contributes new knowledge of fish community structure and spatiotemporal dynamics in glacially influenced estuaries. Here, we synthesize fish community data collected monthly (April through September) over 7 years (2013–2019) along 65 km of coastline in the southeastern Gulf of Alaska. We selected intertidal river delta sites along a natural gradient of glacial to non-glacial watersheds near Juneau, Alaska, to characterize variation in fish communities that are exposed to different degrees of glacial influence. Our objectives were to quantify interannual and seasonal variation in fish community structure (taxa richness and relative abundance) along the glacial to non-glacial watershed gradient and to parse the relative importance of temporal and spatial variables in explaining estuarine fish community patterns. Findings from several studies support the hypothesis that spatial variation plays a stronger role than temporal variation in shaping estuarine communities in coastal Alaska (Robards et al. 1999; Abookire and Piatt 2005; Speckman et al. 2005; Arimitsu et al. 2012); however, these studies involve relatively larger geographic distances and/or narrower windows of temporal sampling than the current study, which complements previous work by increasing the temporal resolution of sampling events both within and among years. This large dataset allows for powerful tests of whether seasonal and annual variation in estuarine fish communities differ among sites that vary in their glacial influence.

Methods

Field sampling

Sampling was conducted in estuaries near Juneau, Alaska, along Lynn Canal, Favorite Channel, and Gastineau Channel (Fig. 1). Sites (n = 7) are located in intertidal habitats at the mouths of rivers (deltas) in watersheds that vary in their size and degree of glacier cover (Table 1, Figure S1). Total watershed area and the ice-covered area within the watershed (expressed as percent glacier cover) have been published for a majority of these sites (Hood and Berner 2009). Updated values were provided by Croix Fylpaa (University of Southeast, personal communication, 4/1/2020), who calculated watershed areas (delineated at the HUC12 boundary designation) using the National Hydrography Dataset and Watershed Boundary Layer dataset (USGS 2020) and glacier coverage using the Randolf glacier inventory (RGI Consortium 2017). The vast majority of sampling was conducted at four delta sites—Sheep Creek, Cowee Creek, Eagle River, and Mendenhall River—along a watershed gradient from 0 to 54% glacierized (i.e., land currently overlaid by glacier; Table 1). Three less frequently sampled sites—Lemon Creek, Berners Bay, and Sawmill Creek—were included in catch summaries so that our full dataset is comprehensively reported, but excluded from statistical analyses. Lemon Creek was sampled intensively during only 1 year (2019) and required a modified sampling protocol due to high current at the site, while Berners Bay and Sawmill Creek were each sampled on 2 days in 2014. Watershed metrics were also unavailable for Berners Bay and Sawmill Creek.

Fig. 1
figure 1

Map of study sites (numbered points) near Juneau, Alaska. The City of Juneau is in closest proximity to symbol #6 on the map. Dark blue lines show the rivers draining into each site (Table 1) and white polygons with blue stippling denote the glacierized portion of watersheds outlined in black. Percent glacier cover is reported for all sites in Table 1

Full size image
Table 1 Sites sampled at river mouths near Juneau, Alaska, from 2013 to 2019. Watershed area and glacier coverage data are from C. Fylpaa (University of Alaska Southeast, personal communication, 4/1/2020). U indicates that watershed metrics are undetermined for infrequently sampled locations
Full size table

Fish were collected by beach seine in shallow sloping, low-relief, intertidal areas along the beach near the outflow of each river. Substrate type was consistent across sites, characterized by sand, mud, and silt with interspersed cobble (Whitney et al. 2017; Figure S1). Data were collected for multiple studies from 2013 to 2017 and 2019 (Whitney et al. 2017; Whitney et al. 2018; Duncan and Beaudreau 2019; Lundstrom et al. 2022; Bergstrom unpublished data), so there was some variation in frequency of sampling and net type used. Seining methodology was consistent across studies and is described in detail by Whitney et al. (2017). For a given sampling event (unique combination of site and date, n = 137), between one and seven seine sets were performed [4.69 ± 1.69 (mean ± SD) sets per event] within 2 h of negative low tide during morning daylight hours between April and September (Table S1). Seine nets were towed parallel to shore by two researchers for a duration of 2–12 min per set [n = 643 sets, 5.59 ± 1.38 (mean ± SD) minutes per set]. While we strove for consistency in the number of sets per sampling event and the duration of each set, weather and tidal conditions on a given day led to some variation. A target number of 5–6 sets were established based on the maximum number of seine sets that were possible during each low tide sampling window. A target set duration of 5–7 min was drawn from our experience operating the beach seine under local conditions, allowing sufficient time to herd fish into the net while limiting the extent of fouling the net with algae or debris. Two net types were consistently used from 2013 to 2019 (black mesh dimensions: 15.2 m long × 2.4 m deep with 0.95 cm stretched knotless black mesh; white mesh dimensions: 15.2 m long × 2.4 m deep with 1.27 cm stretched knotless white mesh). In 2016 and 2017, the white mesh and black mesh nets were both used, allowing for a comparison of catch composition between them (see the “Univariate analysis of taxa richness and relative abundance” section). For 62 sampling events in 2016 and 2017, alternating seine sets were performed with the white mesh and black mesh nets (approximately 3 sets each, for a total of 6 sets/event).

Captured fish were placed in buckets with aerators and processed following completion of all sets. Fishes were identified to species or lowest taxonomic resolution possible (Johnson et al. 2015) and quantified separately for each set. Fish that were not retained for laboratory analyses were returned alive to the site of capture following quantification. Only fishes > 20 mm TL were recorded due to smaller fish escaping through seine net mesh. Sampling was performed under collection permits from the Alaska Department of Fish and Game and the Alaska Department of Natural Resources, and with approval from the University of Alaska Institutional Animal Care and Use Committee (protocols 465729, 880562, 479533, 1238650).

Data preparation and visualization

To prepare data for analysis, we checked for errors and consistency in species identification. For some taxa, fishes were identified to different levels of taxonomic resolution depending on the particular study year and objectives (e.g., juvenile salmonids, Oncorhynchus spp.). Therefore, we reclassified fish to the genus level, except for species in the families Agonidae, Hexagrammidae, and Liparidae, which were classified at the family level due to lack of confident discernment between two or more genera (Table 2). We eliminated seine sets for which set duration had not been recorded or for which there were other methodological issues, such as lengthy hang-ups on debris.

Table 2 Catch summary by taxonomic group and site, pooled across years. Numbers (N) are sample sizes of individuals caught for each taxon at each site, N tot. = total pooled among sites, % tot. is the percentage of all individuals comprised of that taxon, % sets is the percentage of seine sets in which at least one individual of that taxon was caught, and N sites is the number of sites in which at least one individual of that taxon was caught. All groups arranged taxonomically according to Mecklenburg et al. 2002. Site acronyms are defined in Table 1. The number of seine sets per site is reported in the Supplementary Information (Table S1)
Full size table

Catch-per-unit-effort (CPUE) was calculated as a measure of relative abundance for each species as the number of individuals caught per minute for each set individually, then averaged across sets for a given sampling event (site × date). The number of unique taxa was calculated for each seine set as the number of taxa caught per minute (standardized taxa richness). Summary statistics in Supplementary Information tables report catches per 5 min to approximate realistic seine set durations. We evaluated spatial and temporal patterns in taxonomic diversity and community composition using univariate and multivariate analyses, coupled with data visualization.

For five sites (Cowee Creek, Eagle River, Lemon Creek, Mendenhall River, Sheep Creek), we visualized temporal patterns in taxonomic composition and CPUE of dominant species to aid in interpretation of univariate and multivariate statistical results. First, we calculated the taxonomic composition of beach seine catches as the CPUE for each taxon (individuals/min, summed across sets) divided by the total CPUE of all fish across seine sets in a given month. Next, we plotted the relative abundance of dominant species (i.e., those caught in more than 30% of all sets), determined as the average CPUE per set in a given month.

Univariate analysis of taxa richness and relative abundance

To determine whether net type should be included as a covariate in statistical models, we used paired t-tests to compare standardized taxa richness (taxa/minute) and CPUE (individuals/minute) between white mesh and black mesh nets for paired sampling events in 2016–2017. There were significant differences between net types for both metrics (paired t-tests; CPUE: df = 61, t = 3.143, P = 0.003; taxa richness: df = 61, t = 6.076, P < 0.001), with the black net yielding higher CPUE and taxa richness. Therefore, we took net type into account for all further analyses.

We quantified seasonal trends in richness and relative abundance among the four most consistently sampled sites (Cowee Creek, Eagle River, Mendenhall River, Sheep Creek) using generalized additive mixed models (GAMMs) with a Gaussian distribution and identity link (“mgcv” package in R; Wood 2011, R Core Team 2019). Standardized taxa richness and CPUE for each set were separately modeled as a function of fixed effects day of year, site, and an interaction between day of year and site. Net type (i.e., small black mesh, large white mesh) and year were included as random effects to account for differences among gear types and interannual variation. Because of the expected nonlinear seasonal trends in taxa richness and CPUE, we included smoothing functions for day of year and day of year × site terms, with the degree of smoothing (k = 5) determined using generalized cross-validation (GCV; Wood 2017). Response variables were natural log-transformed to normalize the data. The full set of candidate models was comprised of all combinations of potential predictor variables using the “dredge” function from the “MuMIn” package in R (Bartoń 2020; R Core Team 2019) and model selection was performed using Akaike’s information criterion (AIC). Models within two units of the lowest AIC (∆AIC ≤ 2) were considered to perform equivalently (Burnham and Anderson 2002).

Multivariate analysis of community composition

Differences in community composition across sites, months, and years were assessed using multivariate ordination performed using data from the four most consistently sampled sites. Because of the net effect detected in univariate analyses, we performed multivariate analyses separately on data from the white mesh and black mesh nets. To reduce the effects of set-level variation, we used mean CPUE per sampling event (i.e., average of CPUE across sets for each sampling date × site combination) as the input for analyses. CPUE was square-root transformed to reduce the influence of rare species and used to generate an abundance-based Bray–Curtis similarity matrix (Clarke 1993). To test for differences in composition across sites, months, and years, we ran a 3-factor analysis of similarity (ANOSIM) on the Bray–Curtis matrix. After finding a non-significant effect of year (see “Results”), we then ran a 2-factor ANOSIM to test for differences across sites and months. Day of year was converted into month (April through September), as ANOSIM requires categorical factors.

ANOSIM produces both an R-value and a P-value. The R-values (0–1) indicate the degree of separation between groups; a value of 0 indicates groups of identical composition while a value of 1 indicates no overlap in the composition between the groups (Clarke et al. 2014). Following Creque and Czesny (2012), we interpreted R-values ≤ 0.25 as weak differences, R-values 0.26 to 0.50 as moderate differences, and R-values ≥ 0.51 as strong differences among the groups. The P-value produced by ANOSIM is permutation-based and affected by sample size (Clarke et al. 2014), and as such we did not rely on these values for inference. Following ANOSIM, we used similarity percentages (SIMPER) to identify which species contributed most to differences in composition among sites and months. Multivariate analyses were conducted in PRIMER V7.0.13 (Plymouth Marine Laboratories, UK, Clarke and Gorley 2015).

Results

We caught a total of 46,722 fish in 643 beach seines among all 137 sampling events (Table S1). A total of 30 different taxonomic groups was analyzed, comprising 27 genera and 3 families. Combining all seine sets, the most common taxa by total abundance were Pacific staghorn sculpin Leptocottus armatus (46%), starry flounder Platichthys stellatus (21%), juvenile salmon Oncorhynchus spp. (18%), and Dolly Varden Salvelinus malma (3%), while the remaining 26 taxa each accounted for less than 3% of the catch (Table 2). Less commonly encountered species included some that are commercially harvested as adults in Gulf of Alaska fisheries, including yellowfin sole (Limanda aspera; 1.2%), Pacific herring (Clupea pallasii; 1.1%), lingcod (Ophiodon elongatus; < 0.1%), and walleye pollock (Gadus chalcogrammus; < 0.1%, Table 2). Pacific staghorn sculpin was the most commonly occurring species, present in 90% of all beach seine sets, followed by starry flounder (88%), juvenile salmon (46%), rock sole Lepidopsetta spp. (34%), and Dolly Varden (32%); the remaining 25 taxa were found in fewer than 25% of seine sets. Three taxa were captured at all seven sites at least once over the duration of the study (Pacific staghorn sculpin, starry flounder, and juvenile salmon), while eight taxa were captured at only one or two sample sites (Table 2). Across all years and sites, we captured an average of 69 fish from 4.0 ± 2.2 (SD) taxonomic groups per 5-min seine set (Table S1, Table S2).

Patterns of taxa richness and relative abundance

For the five most frequently sampled sites, mean (± SD) taxa richness across sampling events ranged from 3.7 (± 2.0) taxa per 5-min set at Eagle River to 4.4 (± 2.3) taxa per 5-min set at Mendenhall River (Table S2). Seine catches at Cowee Creek, Eagle River, and Mendenhall River sites had a higher species richness across all sampling events combined (25, 25, and 24 unique taxa, respectively), followed by Sheep Creek (16 taxa; Table S2). Lemon Creek had the lowest overall species richness (6 taxa), with catches dominated by Pacific staghorn sculpin, juvenile salmon, and juvenile flatfishes (Table S2). Mean abundance (CPUE) per sampling event ranged from 42 fish per 5-min set at the Cowee Creek site to 108 fish per 5-min set at the Mendenhall River site (Table S1).

The best model explaining variation in standardized taxa richness (number of taxa caught per minute) included day of year, site, day of year × site interaction, year, and net type (Table 3, Table S3). For relative abundance (CPUE), the GAMM that included all predictors except day of year had the lowest AIC score but the full model performed equivalently based on ∆AIC < 2 (Table 3, Table S3). We visually inspected the estimated partial effects from the best fit GAMMs to assess changes in taxa richness and CPUE as a function of each predictor variable, when other sources of variation are held constant. Partial effects plots show increasing relative abundance and taxa richness through spring months to a peak in mid- to late-June before declining slightly and stabilizing at intermediate levels (Fig. 2a and 3a). The seasonal trend in taxa richness and timing of peak richness varied among sites, without a clear relationship to the watershed gradient (Table 3 and Fig. 4). Seasonal trends in relative abundance also varied among sites, exhibiting more dome-shaped curves with peaks in late June at the Sheep Creek and Cowee Creek sites (less glacierized) and approximately asymptotic growth to peak levels in late May at the Eagle River and Mendenhall River sites (more glacierized) (Table 3 and Fig. 5). Based on partial effects plots, relative abundance was lower and taxa richness was higher on average at the Sheep Creek and Cowee Creek sites compared to the Eagle River and Mendenhall River sites (Fig. 2c and 3c). Relative abundance and richness varied among years, with a peak in 2015 and lowest levels in 2013 and 2019, but the 95% confidence intervals of the estimated partial effects overlap among years (Fig. 2b and 3b). The GAMM results confirmed a significant net effect; specifically, relative abundance and taxa richness were higher on average in the black net compared to the white net (Fig. 2d and 3d).

Table 3 Generalized additive mixed model (GAMM) summaries for standardized taxa richness (taxa/min) and catch per unit effort (CPUE, individuals/min), which were natural log transformed and modeled separately. The full model estimated the fixed effects of day of year, site, day of year × site, and random effects of year and net type. Predictor variables included in the top three models for each response are indicated with an X. We also report the adjusted R2 (Adj. R2), percent deviance explained (% Dev.), degrees of freedom (df), log likelihood (logLik), Akaike’s information criterion corrected for small sample size (AICc), ∆AICc, Akaike weight (wi), and generalized cross validation score (GCV) for each model
Full size table
Fig. 2
figure 2

Partial effects of day of year (a), year (b), site (c), and net type (d) on standardized taxa richness, calculated as the number of unique taxa captured per minute of seining. Black lines show the estimated back-transformed predicted values from GAMMs fitted to natural log-transformed data. Gray bands are the 95% confidence intervals and tick marks show the distribution of individual observations. Sites are shown in separate panels, ordered from lowest percent glacier cover in the watershed (Sheep Creek; 0%) to highest percent glacier cover (Mendenhall River; 54%). Percent glacier cover is reported for all sites in Table 1

Full size image
Fig. 3
figure 3

Partial effects of day of year (a), year (b), site (c), and net type (d) on standardized relative abundance (CPUE), calculated as the number of individual fish captured per minute of seining. Black lines show the estimated back-transformed predicted values from GAMMs fitted to natural log-transformed data. Gray bands are the 95% confidence intervals and tick marks show the distribution of individual observations. Sites are shown in separate panels, ordered from lowest percent glacier cover in the watershed (Sheep Creek; 0%) to highest percent glacier cover (Mendenhall River; 54%). Percent glacier cover is reported for all sites in Table 1

Full size image
Fig. 4
figure 4

Seasonal trends in standardized taxa richness, calculated as the number of unique taxa captured per minute of seining. Black lines show the estimated back-transformed partial effects of day of year × site from GAMMs fitted to natural log-transformed data. Gray bands are the 95% confidence intervals and tick marks show the temporal distribution of individual observations. Sites are shown in separate panels, ordered from lowest percent glacier cover in the watershed (Sheep Creek; 0%) to highest percent glacier cover (Mendenhall River; 54%)

Full size image
Fig. 5
figure 5

Seasonal trends in standardized relative abundance (CPUE), calculated as the number of individual fish captured per minute of seining. Black lines show the estimated back-transformed partial effects of day of year × site from GAMMs fitted to natural log-transformed data. Gray bands are the 95% confidence intervals and tick marks show the temporal distribution of individual observations. Sites are shown in separate panels, ordered from lowest percent glacier cover in the watershed (Sheep Creek; 0%) to highest percent glacier cover (Mendenhall River; 54%)

Full size image

Variation in community composition

Results were qualitatively similar between net types, so we present only the results from the white mesh net here, as it was the net type used most consistently across sampling events. The 3-way ANOSIM indicated a weak, non-significant effect of year (R = 0.181, P = 0.058), a weak, significant effect of site (R = 0.215, P = 0.015), and a moderate, significant effect of month (R = 0.329, P = 0.003). A 2-way ANOSIM (with year removed) confirmed a weak, significant effect of site (R = 0.245, P = 0.001) and a moderate, significant effect of month (R = 0.353, P = 0.001).

Variation in community composition across months was visualized using a nonmetric multidimensional (nMDS) scaling plot (Fig. 6). The 3D stress value (0.15) indicated a reasonably good representation of the data in three dimensions, while the 2D stress value (0.21) indicated a moderate fit (Clarke and Gorley 2015). For ease of visual interpretation, we show the 2D ordination plot (Fig. 6), which supports the quantitative results indicating that the greatest differences were between early spring (April–May, light colors) and summer (August, darker colors) months. Across months, weak to moderate significant pairwise differences were found between Mendenhall River and Cowee Creek (R = 0.518), Eagle River and Sheep Creek (R = 0.385), and Eagle River and Cowee Creek (R = 0.240; Table 4). Across sites, moderate to strong, significant pairwise differences were found between April and June (R = 0.724), May and June (R = 0.561), and May and August (R = 0.499; Table 4). Based on SIMPER results, Pacific staghorn sculpin and starry flounder were the top two species contributing to differences across months (Table S4) and sites (Table S5). Pacific staghorn sculpin contributed a mean 24.7% dissimilarity among months (average of % dissimilarity for pairwise month comparisons) and 24.2% dissimilarity among sites, while starry flounder contributed a mean of 13.0% dissimilarity among months and 13.2% dissimilarity among sites (Table S4, Table S5). Juvenile salmon also contributed to the strong differences between spring and summer months (Table S4; mean of 15.9% dissimilarity among months).

Fig. 6
figure 6

Nonmetric multidimensional scaling (nMDS) ordination plot of species composition across sites (all sampling years included). Symbol color and shape represent sampling month. The ordination includes data from all years and sites, except for Berners Bay, Lemon Creek, and Sawmill Creek; only samples collected using the white mesh net were used in this analysis

Full size image
Table 4 Pairwise tests of differences between sites (bottom left, blue shading) and months (upper right, green shading) from the ANOSIM analysis, including R-statistics and p-values in parentheses. Significant differences are in bold print. Differences between years were all non-significant and are not shown. Site acronyms are defined in Table 1
Full size table

Visualization of taxon-specific seasonal patterns

We created visual summaries showing taxonomic composition of the catch across sites and months and relative abundance of five dominant species to aid interpretation of the statistical results. Pacific staghorn sculpin increased in proportional contribution to the catch (Fig. 7) and relative abundance (Fig. 8) from April to a peak around July at most sites. Sheep Creek was an exception, where Pacific staghorn sculpin continued to increase to a peak density in September. Pacific staghorn sculpin occurred in highest densities at Eagle and Mendenhall river sites. Starry flounder contributed a relatively consistent proportion of the catch throughout the summer and found in highest densities at Eagle and Mendenhall river sites (Fig. 7 and 8). A salmon pulse occurred in May, particularly at Sheep Creek, Mendenhall River, and Cowee Creek sites (Fig. 8). Dolly Varden peaked in June and their relative abundance was highest at Sheep and Cowee creek sites (Fig. 8). On average, a higher relative abundance of dominant species occurred at Eagle and Mendenhall river sites compared to Cowee and Sheep creek sites (Fig. 8).

Fig. 7
figure 7

Taxonomic composition of beach seine catches by site, calculated as the catch per unit effort (CPUE, N fish per minute seined) of each taxon divided by the total CPUE across all sets in a given month

Full size image
Fig. 8
figure 8

Relative abundance of dominant species, determined as the average catch per unit effort (N fish per minute seined) per set in a given month. Species caught in more than 30% of all sets are shown. Note the different y-axis scales for each taxon. Percent glacier cover is reported for all sites in Table 1

Full size image

Discussion

Consistent with estuaries in other regions, glacially influenced estuaries in Alaska experience large natural fluctuations in salinity, temperature, dissolved nutrients, and turbidity over tidal, seasonal, and interannual scales, which can affect trophic interactions (Whitney et al. 2017) and community dynamics (Allen 1982; Bechtol 1997; Robards et al. 1999; Abookire et al. 2000). The influence of freshwater, including glacier runoff, is massive within the inside waters of the Alexander Archipelago (Weingartner et al. 2009, Edwards et al. 2021). Differences in seasonal patterns of taxa richness and abundance between the most and least glacially influenced sites suggest that hydrological factors influence the structure of delta fish communities. However, strong mixing due to coastal circulation may have a homogenizing effect that dampens the influence of proximate watersheds on delta habitat and community structure, possibly explaining why differences in community structure among our sites were relatively small. In addition, glacier recession is only one of many landscape-scale drivers affecting watersheds and adjacent coastal ecosystems.

When examined alongside other nearshore surveys of Lynn Canal fish communities (e.g., Johnson et al. 2015), our findings suggest that delta fish communities may be unique from those in adjacent habitats with more complex structure (e.g., cobble, bedrock, vegetation). In particular, our sites were dominated by fewer, highly abundant species that are well adapted to low salinity environments. The relative importance of temporal and spatial factors in explaining variability in shallow-water marine fish communities in the Gulf of Alaska likely depends on scale. In this study, seasonal variation in fish communities was stronger than site effects; however, when Lynn Canal fish communities were compared with those in Kachemak Bay within the central Gulf of Alaska, regional differences were stronger than seasonal factors (Lundstrom et al. 2022). Our study provides an important baseline from which to compare shifts in fish communities over the coming decades as glacier recession is expected to accelerate, with consequences for downstream estuarine habitats.

Seasonal and interannual variation in fish communities

We found strong temporal variation in estuarine fish richness, abundance, and taxonomic composition, primarily due to seasonal effects. In general, both richness and abundance increased from April to a peak in May and June, followed by a decrease or plateau, depending on the site. Taxonomic composition also varied seasonally, with April and May samples exhibiting separation in multivariate space from those collected in June through September. This is explained, in part, by the influence of outmigrating, wild-born juvenile salmon and releases of hatchery-reared juvenile salmon that peak in abundance in the nearshore during May–June (Duncan and Beaudreau 2019; Lundstrom et al. 2022). With the exception of salmon, the early life history of fishes dominating our study sites is poorly understood, particularly the timing of recruitment for juvenile groundfishes. If the strong seasonality of primary production in high latitude estuaries (e.g., Carstensen et al. 2015) is a bottom-up driver of recruitment, we would expect to see higher densities of juvenile fish following the April–May plankton bloom (Ziemann et al. 1991). Catches of starry flounder generally increased gradually over the summer, while rock sole, Pacific staghorn sculpin, and Dolly Varden densities peaked in June–July; however, seasonal patterns in CPUE of individual species varied substantially across sites (Fig. 8). Seasonal shifts in estuarine fish assemblages and size structure are also a function of fish population dynamics (e.g., timing of recruitment) as well as movement of juvenile fishes to feed and/or seek optimal environmental conditions for growth (e.g., Nash 1988, Thiel et al. 1995). Even lower latitude estuaries with dampened temperature variability throughout the year exhibit seasonal changes in fish assemblages; for instance, only around half of the 70 fish species in a Georgia estuary were present year-round (Dalberg and Odum 1970) and abundances of four species fluctuated drastically among seasons in the Apalachicola Bay estuary in Florida (Livingston 1976). Given that our samples were collected during the warmest months, over the period of peak runoff from glacial systems, we would expect to see even stronger seasonality in fish assemblages sampled across the entire year. We did not sample during winter for the current study. However, a small number of seining efforts were conducted November 2017 through March 2018 at Eagle River, Sheep Creek, and Lemon Creek for a different project, and CPUE was nearly zero until early April (Bergstrom, unpublished data), suggesting that most fishes go into deeper water during the colder months.

While our strongest temporal signal was seasonal, there were also minor differences in fish abundance and richness across years, with lowest counts in 2013 and highest in 2015 (Fig. 2 and 3). However, this may have been partly a consequence of sampling error, as our lowest sampling efforts occurred in those years, both of which have relatively large 95% CIs in predicted partial effects from the GAMMs. Additionally, the most heavily sampled years of 2014, 2016, 2017, and 2019 exhibit only slight variability in relative abundance and taxa richness among them, and our multivariate analysis of community structure found an overall non-significant effect of year. Interannual variation in precipitation and temperature may influence estuary fish assemblages in the study area; for example, the summer of 2014 was relatively cool and wet compared to the summer of 2019, which was one of the warmest on record in Southeast Alaska (ACRC 2019; NOAA 2020). However, given only small interannual differences in community structure over the period of the study, we conclude that year-to-year variation at these sites is relatively low at present.

Differences among intertidal sites along the glacial to non-glacial gradient

Seasonal patterns in fish abundance and richness varied among sites, as indicated by the significant site*day of year effects in GAMMs. The more glacially influenced sites (Mendenhall River and Eagle River) reached a peak in relative abundance (CPUE) in May that persisted through September, while the less glacially influenced sites (Sheep Creek and Cowee Creek) peaked in May but decreased through late summer. Lundstrom et al. (2022) detected similar patterns in mean size of two dominant species, Pacific staghorn sculpin and starry flounder, at the same sites. If juvenile groundfish growth and distribution are strongly influenced by freshwater drivers, differences in the seasonal hydrograph of predominantly glacial and precipitation-fed streams could explain these site-specific seasonal trends. Estuaries at the outflow of glacial rivers experience persistent high discharge rates and cold temperatures through summer, while estuaries near less glacierized watersheds experience sharp decreases in flow rate and increases in temperature as rainfall decreases through the summer (Whitney et al. 2017; Lundstrom et al. 2022). The relative stability of summertime temperatures and salinities in glacial systems therefore provides extended periods of suitable habitat for the two most abundant species in this study, Pacific staghorn sculpin and starry flounder, which were more common in Mendenhall River and Eagle River sites. Subtle differences in seasonal patterns of taxa richness among sites are largely explained by differences in the proportions of dominant species across months. The distinct peak in taxa richness during June at Sheep Creek and Cowee Creek sites coincided with a higher proportion of salmonids, particularly Dolly Varden, during that month. Taxa richness at Eagle River and Mendenhall River sites showed more of a plateau during the summer, due to the higher dominance of Pacific staghorn sculpin at those sites from June onward.

Overall differences in taxa richness and abundance across sites were small compared to differences across months, but were significant and exhibited some patterns associated with degree of glacial influence. The more heavily glacially influenced sites, Eagle River and Mendenhall River (41% and 54% glacierized, respectively), had lower mean richness but higher mean abundance compared to Cowee Creek and Sheep Creek (10% and 0% glacierized, respectively), after controlling for other factors in the GAMMs. Increasing freshwater runoff and siltation from glacial meltwater may lead to habitat homogenization, resulting in a lower diversity of benthic species in delta habitats, as observed in Arctic coastal and fjord ecosystems (Wlodarska-Kowalczuk et al. 2007, Węsławski et al. 2011). Thus, there may be lower overall fish richness in estuaries more heavily influenced by glacial run-off (Fig. 2), but greater overall abundance (Fig. 3), driven by a small number of species that are tolerant to low salinities and high sedimentation rates. Fish community composition did not show as clear a separation between the most and least glacially influenced sites, although pairwise differences in community structure were significant between Cowee Creek and both the Eagle River and Mendenhall River sites (Figure S2). A previous study of stable isotopes in estuarine consumers at a subset of our study sites indicated that marine benthic invertebrates in the Cowee Creek estuary had a higher contribution of terrestrial-riverine derived carbon to their diets compared to those at the Mendenhall River site, which had a higher contribution of marine derived carbon (Whitney et al. 2018). Whitney et al. (2018) attributed this to higher retention of terrestrial organic matter in Cowee Creek delta habitat due to substantially lower flows and a higher percent forest cover in the Cowee Creek watershed, compared to the more glacierized Mendenhall River watershed. Another study found moderate differences in diet composition of Pacific staghorn sculpin between Cowee Creek and Mendenhall River estuaries during the period of peak glacial runoff in August (Whitney et al. 2017).

The relatively small number of sites included in this study warrants caution in our interpretation of the role that glacial meltwater plays in explaining ecological variation. A larger number of sites across the glacial to non-glacial watershed gradient would be necessary to confirm these fine-scale spatial patterns. In addition, while there is a gradient of glacial influence across the sites in this study, all of them are influenced by glacial run-off from adjacent systems to some degree, which could dampen variability in fish community structure among sites. Sheltered within the Alexander Archipelago, Lynn Canal has high freshwater retention throughout the spring and summer and relatively little exposure to the greater Gulf of Alaska (Weingartner et al. 2009). The spatial scale of this study was also relatively small compared to others in Alaska that have detected larger spatial differences in fish communities. For example, Lynn Canal estuaries in the eastern Gulf of Alaska are dominated by estuarine resident species like Pacific staghorn sculpin and starry flounder, while Kachemak Bay estuaries in the central Gulf of Alaska have greater proportions of marine species like Pacific sand lance (Ammodytes personatus), Pacific herring, and gadids (Abookire et al. 2000; Lundstrom et al. 2022). These regional differences are partly attributed to environmental conditions, particularly salinity, which is higher overall in Kachemak Bay estuaries that are more exposed to oceanic conditions compared to Lynn Canal estuaries (Lundstrom et al. 2022). An additional caveat is that differences in delta habitats and fish communities among sites is likely driven by variables other than just watershed composition (i.e., percent glacier cover). For example, estuary morphology, exposure, slope, substrate type, and circulation may all play a role in determining habitat features driving fishing community structure (e.g., Whitfield 1986, Carminatto et al. 2020). Sheep Creek and Eagle River are more exposed to prevailing north/south winds and tidal movement in Gastineau Channel and Lynn Canal, respectively, compared to the other sites. Additionally, Eagle River, Mendenhall River, and Sheep Creek sites are closer proximity to salmon hatchery release areas compared to Cowee Creek (Duncan and Beaudreau 2019).

Dominant species in glacially influenced sites

A small number of species had a large influence on our results. The multivariate analysis showed that the two most abundant and ubiquitous species in this study—Pacific staghorn sculpin and starry flounder—contributed the most to spatial and temporal differences in community structure. These dominant taxa exhibited considerable variability among months and among the most intensely sampled sites (i.e., Fig. 8), but occurred at higher relative abundances at Eagle and Mendenhall river sites overall compared to Cowee and Sheep creek sites. Pacific staghorn sculpin and starry flounder tolerate a wide range of temperature and salinity, including nearly freshwater conditions (Morris 1960; Mace 1983; Takeda and Tanaka 2007). These dominant species occur across the Gulf of Alaska, but seem to be particularly prevalent in Lynn Canal, especially in delta habitats (Johnson et al. 2015; Whitney et al. 2017; Duncan and Beaudreau 2019; Lundstrom et al. 2022). The fine sediments deposited at glacial river deltas may be ideal for starry flounder, which burrow into substrates to camouflage themselves (Moles and Norcross 1995). An analysis of starry flounder size structure at four of the estuary sites we sampled (Mendenhall River, Eagle River, Cowee Creek, Sheep Creek) indicated significant relationships between starry flounder mean length and temperature and turbidity, with length increasing with warmer temperatures and higher turbidity (Lundstrom et al. 2022).

The fish assemblage structure we observed at river deltas differs from findings of other intertidal fish community studies in our study region. For example, Johnson et al. (2015) found walleye pollock (Gadus chalcogrammus), Pacific sand lance, and Pacific herring to be the most abundant species in Southeast Alaska fish communities, based on beach seine surveys conducted from 1998 to 2004. However, their sampling was conducted in seagrass, bedrock, mud/sand, and cobble habitats using a larger seine net deployed by boat, which necessitated sampling deeper habitats than the shallow delta areas we sampled (Johnson et al. 2015). We did not fully characterize estuary fish communities at our study sites, as beach seines only sample a narrow band of intertidal habitat and, like all gear types, are not equally selective across species. Sampling adjacent subtidal areas using small trawls, for instance, would help to better resolve the community structure across a larger extent of estuarine and delta habitat. Despite these limitations, our study fills a gap in understanding of fish assemblages in delta habitats of Southeast Alaska, which may be unique and dominated by relatively few species with high physiological tolerance to fluctuations in salinity, temperature, and freshwater runoff.

Conclusions

Freshwater discharge into nearshore marine waters impacts physical habitat parameters (McLusky and Elliott 2004) and many biological and ecological characteristics of fish populations residing there (Sogard 1992; Nelson et al. 2015; Yetsko and Sancho 2015). In southeastern Alaska and the Gulf of Alaska, this annual discharge is massive (Hill et al. 2015) and comprised of both precipitation and glacial meltwater, the latter of which is predicted to increase drastically in the coming decades (Arendt et al. 2009; Hood and Berner 2009). Here, we characterized the structure of estuarine fish communities that are most likely to be affected by these hydrological changes across a spectrum of representative watershed types along the eastern Gulf of Alaska (Sergeant et al. 2020). Our findings improve knowledge about the current range of temporal and spatial variation in metrics of community structure (taxa richness, relative abundance) in glacially influenced estuaries of Alaska, thus providing a benchmark against which to track future ecological changes resulting from watershed transformation. As changing watersheds continue to alter stream and estuarine habitat, anadromous species may be particularly responsive to glacier retreat, which is expected to both create and eliminate habitat for Pacific salmon (Pitman et al. 2020). However, relative stability across a 7-year period with substantial interannual variation in precipitation and temperature suggests that estuarine fish communities in our study region are well adapted to variable conditions, which could confer resilience to environmental change, at least in the near term. Notably, we did not detect a response of fish communities at our sites to the longest persisting marine heat wave of the last decade (2014–2016), which resulted in abrupt ecological shifts across the Gulf of Alaska (Suryan et al. 2021). This lends credence to the expectation that increasing glacier runoff may serve as a cold-water buffer as stream and coastal waters warm (Pitman et al. 2020). Future studies would benefit from continuous measurements of water quality parameters (temperature, salinity, turbidity) to better resolve the mechanistic linkages between fish community structure and estuarine habitat conditions. Continued monitoring of these dynamic systems, especially those that capture intra-annual variability with multiple sampling efforts, will be important for tracking longer-term shifts in estuary habitats and communities affected by glacier retreat.