Influence of Environmental Conditions on Mytilus trossulus Size Frequency Distributions in Two Glacially Influenced Estuaries

The Pacific blue mussel (Mytilus trossulus) is a foundation species in high-latitude intertidal and estuarine systems that creates complex habitats, provides sediment stability, is food for top predators, and links the water column and the benthos. M. trossulus also makes an ideal model species to assess biological responses to environmental variability; specifically, its size frequency distributions can be influenced by the environment in which it lives. Mussels that inhabit estuaries in high latitudes receive freshwater runoff from snow and glacial-fed rivers or can be under oceanic influence. These hydrographic conditions work together with local static environmental characteristics, such as substrate type, fetch, beach slope, distance to freshwater, and glacial discharge to influence mussel demographics. In 2019 and 2020, mussels were collected from two Gulf of Alaska ecoregions to determine whether mussel size frequencies change over spatial (local and ecoregional) and hydrographic scales and whether any static environmental characteristics correlate with this variability. This study demonstrated that mussel size frequencies were most comparable at sites with similar hydrographic conditions, according to the ecoregion and year they were collected. Hydrographic conditions explained approximately 43% of the variation in mussel size frequencies for both years, for the combined ecoregions. Mussel recruits (0–2 mm) were more abundant at sites with higher fetch, while large mussels (> 20 mm) were more abundant at more protected sites. Fetch and freshwater influence explained most of the variation in mussel size frequencies for both years and across both ecoregions, while substrate and slope were also important in 2019 and glacial influence in 2020. This study suggests that hydrographic and static environmental conditions may play an important role in structuring mussel sizes. Although differences in mussel size frequencies were found depending on environmental conditions, mussel sizes showed little difference across differing types of freshwater influence, and so they may be resilient to changes associated with melting glaciers.


Introduction
Warming climates are impacting hydrologic conditions in highlatitude estuaries. Glaciers are melting and receding, and precipitation is more commonly falling as rain than snow (Huss and Hock 2018;Jennings et al. 2018;Doumbia et al. 2020). These changing conditions may impact the biology and maintenance of ecosystem function of coastal marine organisms, such as mussels, which often play an integral and prominent role as ecosystem engineers and foundation species in intertidal systems throughout the world (Borthagaray and Carranza 2007;Buschbaum et al. 2009;Arribas et al. 2014). As ecosystem engineers, mussels create complex habitats by providing sediment stability through attachment to the substrate using byssal threads (Suchanek 1985). This network of byssal threads provides opportunities for settlement of other nearshore organisms (Borthagaray and Carranza 2007;Khalaman et al. 2021). As a foundation species, filter feeding mussels also connect the water column to the benthos by consuming phytoplankton and detritus from nearshore pelagic waters, in addition to other particles, and influence available nutrients through the release of undigested material through excretion known as pseudofeces (Prins and Smaal 1994;Bracken et al. 2012;Young et al. 1996). In addition, mussels serve as a food source to many higher trophic Communicated by Nathan Geraldi 1 3 level coastal species, such as whelks, sea stars, seabirds, sea otters, and humans (Carroll and Highsmith 1996;Sommer et al. 1999;Laidre and Jameson 2006;Coletti et al. 2016;Bodkin et al. 2018;Miller and Dowd 2019).
How mussels interact within their communities and the role they play in the ecosystem can depend on their size. Mussel size can have cascading ecosystem ramifications because some mussel predators are size selective. The ecosystem importance of mussel size results in size frequency studies that examine recruitment (Seed 1969), growth patterns (Wallace 1980), and mortality among cohorts (concomitant with growth and age models, see Kautsky 1982). For example, a size frequency distribution with a peak in the abundance of small mussels (0-2 mm) can be indicative of a recruitment event, and a steadily decreasing abundance of larger mussel sizes can indicate mortality as a result of predation, disease, or disturbance (Kautsky 1982;Witman 1987). Mussel size frequency distributions can also be controlled by environmental conditions, such as population changes in areas that are experiencing increased temperatures (Taylor et al. 2017), or by biological interactions, such as changes in nearshore trophic interactions .
Mussel demographics can be influenced by the environment at differing spatial and temporal scales, such as along environmental gradients and among seasons (Eckman 1996;Westerbom et al. 2002, Sanford andKelly 2011). These influences can be dynamic, where they fluctuate over various periods of time (such as daily or seasonal changes in temperature and salinity) or static, which generally do not change from year-to-year and are considered to be stable on ecological time scales (10 to 100 s of years). Static variables can also vary spatially (i.e., between ecoregions) and influence biological communities (Konar et al. 2016). Examples of static conditions can include variables like substrate type, fetch (a proxy for wave exposure), beach slope, distance to freshwater sources, and, in high-latitude systems, glacial influence (Warwick et al. 1991;Konar et al. 2016). Being mostly sessile, mussels may be exposed to a variety of environmental conditions, both dynamic and static (Seed 1969;Bodkin et al. 2018). Mussels are found in many different coastal water types, from those being primarily oceanic to those influenced by freshwater sources, and mussels living at high latitudes may experience strong spatial differences in these hydrographic conditions, such as exposure to glacial runoff (Neal et al. 2010). During times of glacial discharge, these estuaries experience fluctuations in temperature, salinity, nutrients, and particulate organic matter (Hood and Scott 2008;Hood and Berner 2009;Neal et al. 2010). As a result, glacial melt (mixed with other sources of freshwater runoff) carries terrestrial carbon matter (Hood and Scott 2008) into the intertidal zone, where marine sources of carbon are already present. This mixing of ocean and freshwater may affect the temperature as well as food quality and quantity available to mussels, which have been shown to select for marine-influenced carbon sources, as terrestrial sourced carbon may be more difficult for mussels to assimilate (Bracken et al. 2012;Mann 1988; but see Meerhoff et al. 2018;Schloemer et al. In Revision).
Substrate size is an important factor in substrate stability for sessile mussels, but can also have thermal effects on environmental conditions, with larger substrate (e.g., boulder) taking longer to heat (or cool) and providing more shade than smaller substrate (e.g., cobble) (Gedan et al. 2011). This concept was illustrated during the massive mussel die-off due to the June 2021 heat wave in Canada (Raymond et al. 2022), where blue mussels (Mytilidae) in areas with smaller substrates were more exposed to higher temperatures. The protection provided by larger substrates also aids in preventing dislodgment of organisms by wave action. Higher wave action often associated with sites of larger grain size can increase food availability and feeding time for suspension feeders, such as intertidal mussels (Ricciardi and Bourget 1999), as well as decrease emersion time (Harley and Helmuth 2003), resulting in higher recruitment and increased growth (Burrows et al. 2009). However, a shift in energy allocation from growth to attachment (as much as 10%) in mussels in areas with extremely high wave exposure can result in smaller mussels of the same age than mussels in areas with moderate wave exposure and similar food availability (Steffani and Branch 2003). For example, storm disturbance in exposed areas can cause size-related drag, which can supersede the strength of the byssal thread mussels' use for substrate attachment, causing preferential mortality in large mussels and resulting in a size limited mussel population (Witman 1987;Carrington et al. 2009). The impacts of storm disturbance and wave action can further differ depending on the slope of the beach. Steeper slopes typically reflect higher wave energy, whereas flat slopes dissipate wave energy, resulting in greater retention of food particles for suspension feeders, such as mussels (Ricciardi and Bourget 1999).
Since mussels are important ecological species in many different estuary types and ecoregions, and because their size can greatly influence their niche, it is important to understand if and how environmental conditions may affect their size frequencies. The purpose of this study was to gain a better understanding of how mussel size is influenced by certain static environmental variables that are considered to be stable on ecological time scales. Specifically, this study examined which static environmental variables correlate with mussel size frequency distributions from high-latitude estuaries in two ecoregions in Alaska, both with and without freshwater influence of varying degrees. The hypotheses that were tested include: (1) similar static characteristics, e.g., substrate grain size, fetch, beach slope, distance to freshwater sources, and percent glaciation in the watershed, will correlate with mussel size distributions regardless of ecoregion; (2) substrate size will have a positive correlation with mussel size; (3) fetch (potential for wave exposure) will have a negative correlation with mussel size; and (4) freshwater influence will have a negative correlation to mussel size.

Field Site Descriptions
This study took place in 9 watersheds within the Kachemak Bay and Lynn Canal ecoregions, which are located in the northwestern and southeastern Gulf of Alaska (GoA), respectively (Fig. 1). Kachemak Bay is separated into an outer and inner bay by the Homer Spit (Fig. 1), a land mass that juts southeastward into the bay for 7 km (Johnson 2021). The outer part of Kachemak Bay is exposed to oceanic influences with greater wave action, higher salinity, and colder water temperatures (Spurkland and Iken 2011). The inner bay is more protected, and the nearby Harding Icefield and multiple rivers provide variable freshwater input (Spurkland and Iken 2011). The general circulation pattern in Kachemak Bay consists of several gyres forming at the mouth, which circulate oceanic water in those areas. A second gyre in the middle of the inner bay brings freshwater out from the inner bay and along the northern coast of the bay via currents before heading out and northward into Cook Inlet (Johnson 2021). Study sites were located in a variety of settings throughout Kachemak Bay, from the outer coast, to islands, to smaller inlets in the inner bay. The oceanic sites located in the southern portion of the outer bay were Port Graham, Outside Beach, Elephant Island, and Cohen Island (Fig. 1), and are influenced by oceanic water flowing into the bay. Two sites, Bishop's Beach and Bluff Point, were also in the outer bay ( Fig. 1), but because of their northern shore locations, they are influenced by the combined outflow from Kachemak Bay, including freshwater influences. Although Jakolof Bay is an outer bay site (Fig. 1), it is adjacent to a non-glacial fed river. The inner bay sites were in estuaries with varying amounts of glacial cover, given here as percent of glacial cover in the watershed (Jenckes et al. In Revision). These sites were at the heads of Tutka Bay (8%), Halibut Cove (16%), off the Wosnesenski River (within Neptune Bay; 27%), and off the Grewingk River (60%) (Fig. 1).
Lynn Canal consists of a system of large fjords located about 130 km inland from the GoA (Bruce et al. 1977) (Fig. 1). Compared to Kachemak Bay, this region has less water currents and oceanic influence and higher precipitation (Weingartner et al. 2009). As a result, freshwater runoff plays a prominent role in Lynn Canal. The mountains strongly influence the regional hydrological cycle, because they support glaciers, such as the Eagle and Mendenhall glaciers, and their narrow and steep estuaries respond rapidly to the heavy precipitation in the region (Weingartner et al. 2009). Like Kachemak Bay, the study sites in Lynn Canal were off estuaries with varying amounts of glacial cover: Cowee Creek (10%), Eagle River (41%), and Mendenhall River (54%). Similar to Jakolof Bay, Sheep Creek was downstream of a watershed with no glacial influence and located next to a non-glacial fed river (Fig. 1). There were no sites characterized as oceanic due to the distance of Lynn Canal to the open ocean.

Mytilus trossulus Size Frequency
Mytilus trossulus size frequency distributions were constructed from samples taken in May 2019 and June 2020 from the study sites in both ecoregions described above. At each site, mussels were scraped from the substrate from ten randomly placed 0.25 m × 0.25 m quadrats along a permanent 50-m transect. These transects ran through the middle of the mussel zone within the high rocky intertidal. Mussels were placed into labeled bags and returned to the lab. Shell length (± 0.1 mm) from the umbo to the posterior tip of the shell was then measured using electronic calipers. These lengths were then rounded to the nearest mm, following the general rounding convention (≥ 0.5 rounded up), for the purpose of binning mussel lengths into 2 mm size classes to assess size frequency distributions by site and year.

Static Environmental Correlates
Substrate grain size, fetch, beach slope, distance to freshwater sources, and percent glaciation in the watershed were determined at each site. Substrate was determined by visually estimating percent cover of the various substrate types (using a modified Wentworth scale, Blair and McPherson 1999) in ten 25 cm × 25 cm or 50 cm × 50 cm quadrats along the permanent 50-m transect and then averaging all quadrats for each site. The modified Wentworth scale included mud (< 0.25 mm), sand (0.26-2 mm), gravel (2.1-64 mm), cobble (64.1-256 mm), boulder (256.1 mm-1 m), and bedrock (> 1.1 m). Wave action, as a measure of physical disturbance, can be characterized by the exposure of the location to varying wind speeds and can be represented by using fetch as a proxy (Burrows et al. 2008;Burrows 2012;Konar et al. 2016). Fetch was calculated using the "spoke pattern" method, where vertices were created every 10° for 360° around the center of a site, with a vertex length of up to 200 km (modified from Konar et al. 2016). These vertices were clipped whenever they met with land, using a combination of available data layers (shapefiles) and spatial analyses (waver package) in R Studio. The sum of the clipped vertices was used to estimate the total fetch at each site. Slope was calculated by measuring the rise (vertical) and run (horizontal) at five points along each 50-m transect for every 1 m elevation from mean low water to 4 m tidal elevation (or the supratidal margin). These individual slope measurements were added together and then divided by the total number of measurements to obtain the mean slope for each site. The distance from each study site to the nearest freshwater source was used as a proxy for salinity and was measured as a straight-line distance from the mouth of the nearest river or stream to each study site in ArcGIS Pro version 2.7.1 (ESRI, Redlands, CA). Sources of freshwater were downloaded from the Alaska Hydrography Database (https:// www. arcgis. com/ home/ item. html? id= 9143a 3fc2c 9b41d ca6f9 50bde 6f02a 69) and were added as shapefiles along with glaciers and shoreline into the software. Using the geoprocessing "Generate Near Table" tool, distances from sites to the nearest freshwater source were calculated, only across water bodies and excluding land masses (Konar et al. 2016). Percent glaciation was calculated by taking the amount of glaciated land within the watershed boundary at each site and dividing it by the total area within each watershed (Jenckes et al. In Revision).

Data Analyses
For all analyses, study sites were categorically grouped a priori based on their water type. Water type is a general description of hydrographic conditions, based on the influence of fresh or oceanic water for each site. The different water types are described in this study as "glacial" for those sites adjacent to a glacial fed river; "riverine" for those sites adjacent to non-glacial fed rivers; and "oceanic" for the Kachemak Bay sites under oceanic influence (Fig. 1). The oceanic sites were further grouped into northern and southern oceanic, based on the currents in Kachemak Bay (Johnson 2021). A PERMANOVA analysis in Primer-e version 7 (v6, Plymouth Marine Laboratories, Anderson et al. 2008; hereafter Primer) tested for variation in size frequencies among mussel populations according to region, water type, site, and year using a crossed and nested design (region: fixed; water type: fixed; site: random nested in region and water type; year: fixed). Interactions among factors that were considered high (p ≥ 0.25), and whose contributions were close to zero or negative, were combined or pooled, one at a time following the conventions in Underwood (1997). Mussel size frequency distribution histograms were constructed for each site using 2-mm bin length classes, and the counts of mussels for each length class were plotted using the "breaks" code in R Studio version 1.3.1093 (Becker et al. 1988; R Core Team 2020). The modes of these histograms were calculated in R Studio and used to describe and compare histograms by region, water type, site, and year. Mussel size bins were further expressed as cumulative percentages, and these size class percentages were used to construct a dissimilarity matrix (Manhattan distance, to preserve the distance relationships among the size classes, Clarke et al. 2014). Size class cumulative percentages were further used to construct a non-metric multi-dimensional scaling plot (nMDS) using Primer to visualize size frequency distributions by region, water type, site, and year.
Static environmental variables (i.e., substrate, fetch, slope, percent glaciation, and distance to freshwater) were analyzed for autocorrelation using Spearman rank correlations, resulting in draftsman plots (Appendix Table 4). Static environmental variables were normalized in Primer by subtracting the mean and dividing by the standard deviation for each individual variable for all sites (as in Konar et al. 2016).

3
A resemblance matrix was created based on the similarities in static variables among sites using Euclidean distances on the normalized data, and a principal component analysis (PCA) was used to visualize these similarities by displaying the first two axes (PC1 and PC2) of variation in the static variables explained among the sites.
A distance-based linear model (DISTLM in Primer) was used to determine the strength of the correlations of environmental variables to mussel size frequencies. The DIS-TLM routine is used to analyze and model the relationship between a multivariate biological dataset and one or more predictor variables. This process works as a regression on multiple variables in multivariate space and allows predictor variables to be fitted individually or grouped together (e.g., substrate type, fetch, slope, distance to freshwater, and percent glaciation) (Anderson et al. 2008). The Akaike information criterion (AIC) was used to determine the best model fit for the data. A distance-based redundancy analysis (dbRDA) performs an ordination of fitted values from the DISTLM model and was used to visualize the results of the DISTLM. An eigenanalysis of the fitted data was performed on the resemblance matrix (Manhattan distances on the mussel size classes, see above) and was constrained to find linear combinations of the predictor variables that explained the greatest variation in the data cloud (Anderson et al. 2008). The output plot represents the first two axes that explained the greatest amount of variation; each dbRDA axis represents the relationships between the static environmental variables and the mussel size frequencies.

Mussel Size Frequency
With all sites and years included in the PERMANOVA analysis, mussel size frequency distributions differed significantly among water types (oceanic north and south, riverine, and glacially influenced), with water type explaining approximately 43% of the variation in mussel size frequencies (Table 1). Although region and year collectively explained ~22% of model variation, water type explained ~3× more variation than region and ~5× more variation than year, respectively ( Table 1). As such, the influences of region and year may not necessarily be biologically meaningful. Additionally, while water type was significantly different within both years (p = 0. 0001), mussel size frequencies were not significantly correlated with the variable glacial coverage in 2019 (referred to as percent glaciation, DISTLM, p = 0.09), but they were in 2020 (DISTLM, p = 0.04, Table 2). However, the mussel size frequency distribution means did vary by hydrographic condition (Fig. 2). Glaciated and riverine sites had the highest mean mussel sizes for both years (Fig. 2), while oceanic sites (northern and southern) had the smallest (Fig. 3). There was a high amount of variation in the mussel size ranges among the glacial and riverine sites for 2019 compared to 2020, while the oceanic sites were more consistent in sizes for both years (Figs. 2 and 3). Overall, there was no apparent trend in mussel size that corresponded with the amount of glacial cover in a watershed (Fig. 2).
While there was visually no clear pattern in the size frequency distributions among the glacial and riverine sites for both years, there was a difference between ecoregions; Kachemak Bay glacial and riverine sites were mostly unimodal; and Lynn Canal glacial and riverine sites were mostly bimodal (Fig. 2). The oceanic sites, which were all in Kachemak Bay, were mostly unimodal, with smaller sizes compared to the glacial and riverine sites in that ecoregion (Fig. 3). Within the Table 1 PERMANOVA results testing differences in the mussel size frequency data by ecoregion (Kachemak Bay and Lynn Canal), water type (oceanic north and south, riverine, and glacial, nested in region), year (2019 and 2020), and site (nested in region and water type, not included, because the site is the level of replication; therefore, no test is performed) Pooled (1) refers to the pooled interaction terms site (nested in region × water type) + region × water type, and pooled (2)   1 3 oceanic sites, there was a large number of small mussels at the northern compared with the southern sites, as seen at Bishop's Beach and Bluff Point in 2019, which together had over 30,000 mussels for all analyzed quadrats combined (Fig. 3). In contrast, abundance of recruitment-sized mussels (≤ 2 mm) at the southern oceanic sites was at least an order of magnitude lower than the northern oceanic sites (Fig. 3).
In 2019, all the Kachemak Bay oceanic sites and the Lynn Canal riverine site (Sheep Creek) had an initial peak in recruitment-sized mussels followed by a steady decline in frequency with size ( Figs. 2 and 3). This recruitment event was not seen in 2020 at the Lynn Canal riverine site, but was observed at several of the oceanic sites. When visually comparing size frequency distributions by water type, it was apparent that the glacial sites had a higher frequency of larger mussels (> 10 mm) than the other water types (Figs. 2 and 3). In addition, an nMDS based on the mussel size frequencies across all sites and years further supported the organization of mussel size frequencies by glacial and oceanic water types (Fig. 4). A stress value of 0.02 indicates a very good representation, with little chance of misinterpretation (Clarke et al. 2014). Glacial sites aggregated in spite of the inherent differences in size frequencies along the glacial gradient and oceanic sites grouped based on north or south affiliation. The riverine site from Kachemak Bay (Jakolof Bay) grouped with the glacial sites, and the Lynn Canal riverine site (Sheep Creek) clustered with the southern oceanic sites (Fig. 4).

Static Environmental Correlates
Spatial patterns in static environmental variables were only analyzed for 2019, because we assumed that they do not substantially vary over the short time period of this study. None of the correlation coefficients among environmental variables was greater than 0.95; so by convention (see Clarke and Ainsworth 1993), no variables were eliminated during further analyses. A PERMANOVA showed that region and water type were both significant factors in the static environmental data that were examined, which included substrate grain size, fetch, beach slope, distance to freshwater sources, and percent glaciation (Table 2). A principal component analysis (PCA) Fig. 2 Size frequency histograms of mussels from glacial and riverine sites in Kachemak Bay and Lynn Canal for the years A 2019 and B 2020. Percent glaciation is given in parentheses next to the site name, and the mode for mussel size is represented by a dashed line, while the mean is indicated as a solid line on the normalized environmental data created components of variation based on the differences among static variables and displayed these components as different axes (Fig. 5). The PCA plot displayed the first two axes (PC1 and PC2) of variation explained among the sites. PC1 explained approximately 41% of variation, and PC2 explained approximately 20%, totaling 61%. Based on their static environmental variables, sites grouped together according to water type, with the oceanic sites clustering together on the negative side of PC1, and the glacial and riverine sites aggregating on the positive side (Fig. 5). PC2 represented a gradient of substrate and slope, with coarser substrate (i.e., cobble and boulder) and greater slope being associated with positive axis values (Fig. 5). The northern oceanic sites were characterized by flat beaches and no glacial influence (see Appendix Table 5) and correlated with smaller substrate (Fig. 5). The southern oceanic sites had a high correlation with large (bedrock and boulder) substrate and increased distance to freshwater (Fig. 5) and were characterized by steeper beaches and no glacial influence (Appendix Table 5). The riverine and glacial sites grouped on the positive side of the PCA ordination, but differed in their environmental correlates. For example, the Kachemak Bay riverine site (Jakolof Bay) correlated with cobble substrate and steep slope, while the Lynn Canal riverine site (Sheep Creek) was more influenced by muddy substrate, percent glaciation, and gradual slope (Fig. 5). Both riverine sites were characterized by low fetch and were in close proximity to a freshwater source (Appendix Table 5). Sites with the highest percent of glacial influence correlated with gravel substrate (Fig. 5), and these sites were best described by a gradual slope (Appendix Table 5). The rest of the glacial sites (with the exception of Cowee Creek, which correlated with sand) correlated with cobble substrate and steep slope (Fig. 5).

Environmental Influence on Mussel Size Frequencies
For both years and regions combined, the significant environmental variables that contributed to variation in mussel size frequencies included fetch (p = 0.003) and distance to freshwater (p = 0.011, Table 3). Substrate and slope contributed to the distance-based linear model in 2019 (p = 0.043 and 0.141, respectively), though slope was not a significant variable. Distance to freshwater and percent glaciation were significant in 2020 (p = 0.024 and 0.043, respectively, Table 3). Mussel size frequencies clustered according to water type in both years, with the Kachemak Bay oceanic (north and south) sites and the Lynn Canal riverine site correlating to higher negative values for the first axis (dbRDA1), and all glacial sites and the Kachemak Bay riverine site on the positive side (Fig. 6). Overall, the first axis (dbRDA1) explained 93.2% of the fitted variation (86% of the total variation) in 2019 and 98% of the fitted variation (89% of the total variation) in 2020 (Fig. 6). The second axis (dbRDA2) explained relatively little, only 5.3% of the fitted variation (4.8% of the total variation) in 2019, and 2.4% of the fitted variation (2.2% of the total variation) in 2020 (Fig. 6). Sites with the highest frequency of recruitment-sized mussels (0-2 mm) in 2019 and 2020 were the Kachemak Bay northern oceanic sites, which were mostly influenced by high Fig. 4 A non-metric multidimensional scaling (nMDS) plot showing the rank order relationships in mussel size data for each site by water type, labeled by regional association (KB, Kachemak Bay, and LC, Lynn Canal). Both years (2019 and 2020) are shown for each site The circle visualizes the correlations between the original dataset features and the principal components, which are displayed as coordinates. Principal component 1 (PC1) represents approximately 41% of variation in static variables among study sites, and principal component 2 (PC2) represents approximately 20% of variation, for a total of 61% variation explained. Sites are labeled according to water type (legend symbols) and region (KB, Kachemak Bay; LC, Lynn Canal) fetch (p = 0.003), with the largest negative values for dbRDA1. One northern oceanic site (Bluff Point) correlated with large substrate (bedrock) for both years, while the other (Bishop's Beach) correlated with smaller substrate; substrate was also a significant environmental variable in 2019 (p = 0.043). The Lynn Canal riverine site (Sheep Creek) correlated with negative values on the dbRDA1 axis and, as such, clustered with the oceanic-influenced sites (Fig. 6). Glaciated sites had the highest frequencies of larger-sized mussels (≥ 10 mm) in both years and correlated with higher positive dbRDA values representing higher percent glaciation, middle-sized substrate types (cobble and gravel), and gradual slope. The positions of the oceanic (north and south) and the glacial sites at opposing sides of the dbRDA axis indicate that these mussel size distributions were most dissimilar from each other; mussels growing at glacially influenced sites, or in areas with low salinities, tended to be larger than at sites with oceanic influence, or in areas with high salinities.

Discussion
Mytilus trossulus is exposed to a multitude of environmental conditions in nearshore ecosystems around the world. These conditions are likely influencing their recruitment, growth, mortality, and size distribution, which may then influence the role they play in the overall community. In high latitudes, for example, the GoA, glacial presence may also shape M. trossulus size frequency distributions. This study demonstrated that M. trossulus size frequencies were generally comparable at sites with similar water types and sampling years. Most of the same static variables were significant and had similar effects on size frequencies regardless of ecoregion (Kachemak Bay and Lynn Canal), in agreement with the first hypothesis that similar static characteristics would correlate with M. trossulus size distributions across ecoregions. This result agrees with other studies that have shown that the influence of environmental variables, as well as species demographics, can be similar across regional scales (Konar et al. 2016;Bodkin et al. 2018).
We also hypothesized that M. trossulus size frequency would be positively correlated with substrate size. In 2019, our hypothesis was supported, but not in 2020, although it was still included in the best-fitting model based on the AIC criterion, indicating it may have some influence on M. trossulus sizes. This is in agreement with a study from New England, showing that larger rock substrates stayed cooler and facilitated greater survival of another prominent nearshore invertebrate, the barnacle Semibalanus balanoides (Gedan et al. 2011). Additionally, M. trossulus beds located in areas with smaller substrate (e.g., mud or cobble) could also be exposed to higher temperatures due to the lack of thermal mass and, consequently, ideal microhabitat conditions in small substrate (Gedan et al. 2011).

Increases and decreases in temperatures can result in smaller
Mytilus sizes, as mussels exposed to temperatures outside their tolerated range of ~5-32 °C (Seed and Schuanek 1992;Braby and Somero 2006) will allocate energy away from growth and not grow as large as mussels living within their tolerated range (Helmuth 1998). For example, M. trossulus in Alaska grow larger in water temperatures of ~9 °C (max. length = ~60 mm, sea surface temperature|National Marine Ecosystem Status (noaa.gov)) than M. trossulus in California, where water temperatures reach ~32 °C (max. length = ~50 mm, Braby and Somero 2006). Exposure to very high temperatures can also result in mortality events, such as in British Columbia during the June 2021 heat wave in Canada (Raymond et al. 2022). Further, mussel populations in Greenland occur in protected microhabitats that are not exposed to sub-freezing temperatures outside of their tolerance range, and thus, mussels can grow between 40 and 73 mm (Thyrring et al. 2017). Temperature was not included in the present study; therefore, it is inconclusive if substrate played a role in temperature modulation, and if the mussels in this study were smaller at sites with higher temperatures. However, it is difficult to infer growth from size, as mussels in an area with high disturbance can be smaller than mussels of the same age in a more protected area (Seed 1969).
In addition, substrate size is often related to wave exposure; areas with higher wave exposure generally have larger substrate sizes compared to protected areas, as smaller substrates get swept away or are frequently overturned by wave action (Skilbeck et al. 2017). In this study, M. trossulus at two sites with the smallest-and largest-sized substrate and higher fetch were typically smaller in size when compared to M. trossulus from sites with medium-sized substrate and lower fetch. These results illustrate how substrate size and fetch can concurrently affect M. trossulus size frequencies, although it is difficult to conclusively say that this study supports hypothesis 2: substrate size did correlate positively with mussel size. In this study, fetch was negatively correlated to mussel size in both sampling years and appeared to play a larger role than the other static variables in structuring M. trossulus size frequency distributions, supporting our hypothesis 3: fetch correlates negatively to mussel size. Smaller M. trossulus found on more exposed beaches are probably partially due to byssal thread strength of larger individuals not being able to withstand the extreme wave action due to size-related drag (Carrington et al. 2009). Studies have found that mussels can outgrow the strength of their byssal threads (Witman 1987;Carrington et al. 2009) and this may be what is limiting mussel size at exposed beaches in this study.
Glaciation also contributed to the mussel size variation in both sampling years. Contrary to our hypothesis 4: freshwater influence will have a negative correlation to mussel size, mussels growing at freshwater-influenced sites (glacial and riverine) tended to be larger than at sites with no freshwater influx. The salinity tolerance range for mussels is typically 10-30 PSU, with the preferred salinity of ∼ 25 PSU (Riisgård et al. 2012). Because larger M. trossulus individuals were found in the freshwaterinfluenced sites, this may indicate that the salinity values at these sites are ideal for this species. A lack of predators at these sites Fig. 6 Distance-based redundancy analysis (dbRDA) for the environmental correlates and the mussel size frequency distributions for A 2019 and B 2020 at all sites (labeled by region: KB, Kachemak Bay, and LC, Lynn Canal), based on the BEST selection procedure in DISTLM. Sites are designated by region and water type. Axes represent dbRDA coordinate scores; dbRDA1 in 2019 accounted for 93.2% of the explained variation of the fitted model, and 85.7% of the explained variation of the total variation; dbRDA1 in 2020 accounted for 98% of the explained variation of the fitted model and 88.8% of the explained variation of the total variation. Vectors for each variable are overlaid, their length representing both the strength and direction of the effect of environmental variables on the dbRDA axes scores. Correlation circles are pictured in both plots; sites closer to the center of the circle and farther away from the circle are not well represented by the data; sites close to the circumference are best represented by the data could also alter size frequencies; however, sea stars, sea otters, and predatory gastropods are common (pers. obs., McCabe and Konar 2021). Instead, the link between glacial melt and freshwater discharge and the output of nutrients could explain some of the positive effect of percent glaciation on mussel size (Hood and Scott 2008;Arimitsu et al. 2017). For example, a study in glacially influenced estuaries in Chile reported that an increase in glacial melt resulted in more available nutrients for benthic organisms (Meerhoff et al. 2018). Similarly, studies done in glacially influenced Alaskan waters found an increase in nitrogen and carbon available to nearshore organisms (Hood and Berner 2009;Arimitsu et al. 2017). Another Alaskan study found that ∼ 7% of dissolved organic carbon and ∼ 38% of dissolved organic nitrogen from terrestrial dissolved organic matter were bioavailable to marine microbial communities on short, 4-6-day time scales (Sipler et al. 2017). If these freshwater-delivered terrestrial materials, e.g., in form of microbial particles, are being used by the filter feeders (i.e., mussels), it could indicate high-latitude sources of freshwater are, in fact, delivering more labile nutrients to nearshore organisms in high-latitude estuaries, as was seen in a Washington State study (Howe et al. 2017), although this is contrary to other studies (Bracken et al. 2012, Thyrring et al. 2017. Another consideration is that static variables can include dynamic components and can impact dynamic variables (Konar et al. 2016). For example, the percent of glaciated area in an estuary is static, but the discharge off the glacier is variable depending on season as well as watershed and glacier size. In another example, the distance from a study site to a freshwater source is static, but the amount of freshwater input can be temporally variable depending on the timing of sampling and the river/stream size. As seen in this study, the influence that an environmental variable has may also vary from year to year. For example, slope was important in 2019, but not in 2020. In 2020, sampling occurred a month earlier than in 2019, and so we may have missed the recruitment window that was seen in 2019. This suggests that the influence of slope is higher when mussel recruitment is elevated, as mussels tend to settle in places with increased food availability, such as in areas with gentle slopes (Ricciardi and Bourget 1999). This is an environment where algae also typically settle (Adami et al. 2004), as mussels often settle on filamentous and other types of algae (Seed 1969;Seed and Suchanek 1992). Slope and fetch can impact where mussels and algae settle (i.e., areas with a steep slope and high fetch will have low settlement), and filamentous algae typically do not survive in areas with extreme wave action (Sousa 1979), thereby affecting the ability for mussels to settle at those locations (Adami et al. 2004).
Mussel recruitment events are marked by an abundance of small mussels (< 2 mm, Kautsky 1982;Westerbom et al. 2002), and can occur multiple times throughout the year, mainly during spring and late summer (Suchanek 1985). Annual variability in mussel recruitment has been shown in other studies (Kautsky 1982;Seed and Suchanek 1992;Hunt and Scheibling 1996) due to an increase in fecundity (reproductive output), food availability, a decrease in mussel predation (Kroeker et al. 2016), or a shift in the coastal currents that assist in juvenile mussel dispersal and settlement (Connolly et al. 2001, except, see Bodkin et al. 2018, which showed consistent recruitment at some sites across the Gulf of Alaska but did not include Kachemak Bay). In this current study, M. trossulus had high recruitment in 2019 at Kachemak Bay's Bishop's Beach and Bluff Point sites and in Lynn Canal at Sheep Creek. While it is unknown why recruitment was higher in 2019 than 2020, this study has shown that recruitment events can be spatially and temporally variable in high-latitude estuaries. In other areas with little environmental fluctuation, mussel recruitment has occurred continuously (Seed 1969;Seed and Suchanek 1992;Hunt and Scheibling 1996). Though temporal variation in mussel recruitment was not addressed in this study, spatial variation was observed between ecoregions. An example of this may be due to the more consistent hydrographic conditions at the southern oceanic sites in Kachemak Bay compared to more variable glacial conditions in both regions. In this study, size frequency histograms for the southern oceanic sites demonstrated little variation, though there was evidence of some recruitment. Most of the populations also had right-skewed distributions, with mussel frequencies decreasing with increasing mussel size, suggesting growth of cohorts in these populations (Kautsky 1982;Connolly et al. 2001). Alternatively, the observed distribution pattern could also suggest predation by size-selective predators (see below).
Growth of cohorts can be seen in the size frequency histograms for glacial and riverine water types, indicated by peaks in multi-modal distributions (Kautsky 1982;Connolly et al. 2001). The glacial sites in Lynn Canal in 2019 had bimodal distributions, suggesting growth as cohorts (Seed 1969;Connolly et al. 2001). The glacial sites in Kachemak Bay in 2019 had more unimodal distributions, with a lower occurrence of recruitment-sized mussels and a higher frequency of large mussels (> 20 mm). This may be due to a decrease in intraspecific competition. For example, mussel populations dominated by recruitment-sized mussels can experience high intraspecific competition, as large numbers of new mussels compete for food (Kautsky 1982).In the absence of many small mussels, larger mussels in protected areas may capitalize on the lack of competition and the added resources, resulting in faster growth. This concept follows the widely known Kleiber's law, where larger organisms have lower metabolic turnover and energy demand per unit mass than smaller-bodied organisms that are growing and have a higher energy demand (Sukhotin et al. 2020). Therefore, mussels in glacially influenced estuaries may be growing quickly because of higher food availability (Kroeker et al. 2016), both high-quality marine and terrestrial material (Hood and Berner 2009;Neal et al. 2010). Combined with the positive effect of being protected from wave action, mussels in glacially influenced estuaries may experience ideal conditions for growth, resulting in more mussels that reach the size preferred by some mussel predators (Sommer et al. 1999;Doroff et al. 2012;Robinson et al. 2019).
Differences in mussel size availability can have ecosystem repercussions. For example, sea otters and black oystercatchers select for mussels > 20 mm as prey (Laidre and Jameson 2006;Miller and Dowd 2019). Sea otters in Kachemak Bay have seasonal patterns when selecting foraging habitats, preferring protected bays during winter storms (Gill et al. 2009). Sea otter spraint (feces) analysis in Kachemak Bay revealed that M. trossulus is a dominant prey source (Doroff et al. 2012), indicating that sea otters frequent areas where mussels grow larger. This coincides predominantly with the protected glacial sites in this study, though sea otters in Kachemak Bay are also frequent in the non-glaciated areas (Gill et al. 2009). Black oystercatchers are typically found in protected bays where they have demonstrated prey size fidelity (30-45 mm, Meire and Ervynck 1986) when feeding chicks such that rather than selecting smaller mussels or different prey species for their chicks, adults will travel farther looking for larger mussels (Robinson et al. 2019). Mussel selectivity by black oystercatchers may be contributing to a top-down control on the mussel populations in protected bays in this study, resulting in a lower frequency of large mussels compared to the frequency of small mussels at the oceanic sites. A general relationship between mussel size and sea star predation has also been demonstrated, with sea stars selecting mussels in relation to their own body size (Gooding and Harley 2015). However, sea star wasting decimated sea star populations along the Pacific coast, including Alaska (Konar et al. 2019), although they are starting to repopulate (pers. obs.). As mussels are prominent competitors for space, the loss of this important mussel predator may have contributed to an increase in mussel cover (Weitzman et al. 2021). These predator-prey relationships indicate that if a shift occurs in mussel size frequencies, it could have effects on nearshore predators. For example, a higher frequency of large mussels can be beneficial for nearshore predator species.

Conclusion
There are many environmental factors that may explain variation in mussel size frequencies. This study has demonstrated that hydrographic conditions in addition to some static environmental variables (specifically substrate type, fetch, percent glaciation, and to a lesser degree slope and distance to freshwater) may influence mussel size frequencies in high-latitude estuaries. As air temperatures increase, glaciers in the GoA will recede at an increasing rate (Dyurgerov and Meier 2000;Hood and Berner 2009;Kendall et al. 2010). Seasonal variation in precipitation adds another significant freshwater source to the coastal environments in the GoA (Bieniek et al. 2014). This increase in glacial melt can lead to changes in water temperature, salinity, and available nutrients, resulting in changes in the nearshore environment (Hood and Berner 2009). In future studies, examining the roles of these more dynamic environmental drivers, such as seasonal discharge or water temperatures, would be valuable. Additionally, examining these patterns over a longer time series would be useful to determine if the patterns in mussel size frequency found in this study are consistent over time. In this study, M. trossulus at glacial and riverine sites grew larger, suggesting that these mussels may do well in future high-latitude estuaries influenced by freshwater sources. Because larger mussels are preferred by some top predators, more larger mussels could support more top predators. This research suggests that hydrographic conditions and other static attributes are important in determining mussel sizes but that the lack of mussel size response with increasing glacial melt indicates that mussels may be resilient to differing freshwater influences associated with melting glaciers. Table 4 Spearman Rank correlation for all static environmental variables in this study, showing no correlation greater than 0.95 among the variables, which was used as the cut-off value for high correlation (Clarke and Ainsworth 1993) Substrate sizes and units are in parentheses next to the substrate types. "Slope 0.5" indicates that slope was measured at each site at the 0.5 m tide level. "Fetch (200 m)" was measured up to 200 m out from each site, every ten degrees, around 360°. "Dist. to. fresh(m)" is the distance to a freshwater source from each site, measured in meters. %Glaciation is the percent of glaciated area inside each watershed, where each sampling site is located. All data were normalized by subtracting the mean and dividing by the standard deviation for each individual variable for all sites (as in Konar et al. 2016 Substrate sizes and units are in parentheses next to the substrate types. "Slope 0.5" indicates that slope was measured at each site at the 0.5 m tide level. "Fetch (200 m)" was measured up to 200 m out from each site, every ten degrees, around 360°. "Dist. to. fresh(m)" is the distance to a freshwater source from each site, measured in meters. %Glaciation is the percent of glaciated area inside each watershed, where each sampling site is located Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.