Size-spectra analysis in the estuary: assessing fish nursery function across a habitat mosaic

Estuaries act as nurseries for a wide variety of fish species, potentially providing vital foraging opportunities and refuge from predation for their juvenile residents. Yet, these dynamic environments are comprised of a mosaic of habitat types that span gradients of both salinity and physical habitat structure. Here, we present a novel use of size-spectra analysis to infer nursery habitat function across the estuary habitat mosaic. Interpreting slope and intercept values of abundance against body mass size-spectra regressions as indicators of predation risk and production, we constructed spectra for six distinct habitat types across the entire tidal influence of an unindustrialized estuary in coastal British Columbia. Based on catches of >200,000 individual fish representing 30 different species from April through September, the estuary rockweed mudflat habitat had the lowest size-spectra slope and highest intercept, consistent with lower predation risk and higher production. Size-spectra coefficients varied seasonally across the ecotone, indicating spatio-temporal variation in key nursery functions. Size-spectra can provide insight into key ecological processes of productivity and predation risk across dynamic aquatic habitats.


INTRODUCTION
Coastal wetlands, such as estuaries, salt marshes, mangroves, and seagrass meadows, are some of the most productive and ecologically significant ecosystems on the planet, yet they are also among the most threatened (van den Belt 2011). Estuarine environments are notably important to juvenile fishes, providing foraging opportunities, reduced predation risks, and a mosaic of habitats suitable for a variety of life stages (Dahlgren and Eggleston 2000, Peterson 2003. Juvenile fish should seek out habitats where there is the best potential to maximize growth and the lowest potential for mortality (Werner and Gilliam 1984, Dahlgren and Eggleston 2000, Cocheret de la Morini ere et al. 2002. As such, estuaries have been dubbed as nurseries by many scientists, conservation groups, managers, and the public at large (Beck et al. 2001). With nearly 60% of humanity residing within 100 km of the coastline (Vitousek et al. 1997) and associated rapid, largescale development in coastal habitats (Hughes et al. 2009), assessing fish nursery function within estuary environments can help prioritize protection and restoration (Beck et al. 2001, Peterson 2003.
With the rising popularity of ecosystem-based management, managers are asked to account for a multitude of interacting biological and physical factors while simultaneously uncovering the drivers and pressures which cause ecological change (Guiet et al. 2016). The increasingly complex ecological models being employed by ecosystem managers require rich datasets based on extensive field sampling for sufficient parameterization (Guiet et al. 2016). This is particularly challenging for estuaries where both the abiotic and biotic components of the ecosystem are in constant flux (Beck et al. 2001. Size distributions, socalled size-spectra, can be employed to infer ecosystem structure and function, while not being overwhelmed by complexity and associated model uncertainty (Giacomini et al. 2016, Guiet et al. 2016.
Size-spectrum theory is based upon the observation that in pelagic systems there are many more small individuals than large ones, but that total biomass is approximately equal across size classes (Sheldon et al. 1972). By plotting either biomass, abundance, or energy as a function of individuals' length, weight, or volume in logarithmic space, size-spectra summarize a complex suite of biophysical and trophic processes with simple linear regressions (Kerr andDickie 2001, Guiet et al. 2016). The coefficients of spectra regressions, slope (k) and intercept (Γ), can be used to infer different properties of community structure and ecosystem health (Kerr andDickie 2001, Guiet et al. 2016). Spectra slopes (k) measure the relative frequency of body sizes in a community and depend upon predator-prey dynamics and trophic transfer efficiencies (Platt and Denman 1978, Kerr and Dickie 2001, Giacomini et al. 2016, Perkins et al. 2018. When larger individuals are removed from the population, such as in the case of fisheries exploitation, spectra slopes (k) decrease and the abundance of smaller individuals increases, suggesting that the smaller size classes experienced a release from predation pressure (Gislason and Rice 1998, Bianchi et al. 2000, Blanchard et al. 2005. Thus, slope values can indicate the relative degree of predation risk, with steeper, more negative slopes being associated with lower predation risk due to the high number of smaller sizeclass fish relative to individuals in larger size classes. Intercepts (Γ) correspond to ecosystem production and are influenced by abiotic environmental factors such as nutrient availability and temperature (Boudreau and Dickie 1992, Maury et al. 2007, Guiet et al. 2016. Highly productive habitats, such as eutrophic lakes and upwelling zones, exhibit higher spectra intercept values than those with low overall productivity, such as oligotrophic regions (Bianchi et al. 2000, Guiet et al. 2016. That is, highly productive habitats would be expected to have higher fish abundances than less productive habitats which would be reflected in the spectra intercept. Size-spectra analysis has been used to provide insight into a variety of ecological processes, including predator-prey interactions (Thiebaux and Dickie 1992, Blumenshine et al. 2000, Kerr and Dickie 2001, effects of fisheries exploitation on community structure (Gislason and Rice 1998, Bianchi et al. 2000, Blanchard et al. 2005, and the effects of resource subsidies on food-web structure (Hocking et al. 2013, Perkins et al. 2018.
Here, we apply size-spectra analysis across an estuary to gain insight into the potential nursery functions of different habitats across this ecotone. We hypothesized that size-spectra along the estuary gradient from freshwater to saltwater would reflect spatial gradients of predation risk and productivity. First, we predicted that predation risk, as revealed by higher slopes (k), would increase closer to the ocean due to the relationship between habitat size and food chain length (McIntosh et al. 2018). Second, we predicted that production, as revealed by higher intercepts (Γ), would also increase in the seaward direction as temperate latitude oceans are more productive than their freshwaters (Gross et al. 1988). Thus, young fish would face trade-offs in predation risk and productivity across the estuary mosaic. However, it is also possible that particularly important nursery habitat areas break this tradeoff and have lower slopes and higher intercepts. We produced spectra based on abundance and body weight to describe the fish communities v www.esajournals.org across the full tidal influence of an unindustrialized estuary in the Great Bear Rainforest region of British Columbia, Canada (Fig. 1). Given that seasonal shifts in environmental conditions and biotic community structure can alter size-spectra (Pope et al. 1994, Maury et al. 2007, Guiet et al. 2016, we also examined how spectra coefficients fluctuated through time across the mosaic.

METHODS
The Koeye River estuary (51.7782°N, 127. 8737°W) is located on the Central Coast of British Columbia, Canada,~50 km south of the Heiltsuk village of Bella Bella (Fig. 1a). Sampling sites spanned the entirety of the tidal influence, ranging from the marine Fitz Hugh Sound to , the most frequently used beach seine sites (red dots), and the general salinity gradient across the estuary. The gray shaded area between reaches 2 and 3 is a narrow canyon with steep rock walls and swift currents that make beach seining nearly impossible and thus was not sampled. (b) A pictorial representation of the Koeye estuary size-spectra, showing approximate abundances of typical fish species caught in each size class. v www.esajournals.org 3 November 2020 v Volume 11(11) v Article e03291 ~5 km upstream of the river mouth. We divided the estuary into six reaches based on habitat characteristics ( Fig. 1a; Reshitnyk 2015): Reach 1 is typically fully saltwater with sandy beaches, Reach 2 is highly saline with dense eelgrass (Zostera marina) beds, Reach 3 has variable salinity with muddy substrate and large expanses of rockweed (Fucus spp.), Reach 4 has variable salinity and is the mainstem salt marsh, Reach 5 has variable salinity and is side-channel salt marsh, and Reach 6 is primarily freshwater with classic stream riffle-pool characteristics (Appendix S1: Table S1). Depths and bathymetry were variable across the estuary, with Reach 1 consisting of gentle slopes with a maximum seinable depth of~3 m, Reach 2 had steeper slopes with a maximum seinable depth of~4 m, Reach 3 had very low gradient slopes with a maximum seinable depth of~2 m, and Reaches 4, 5, and 6 were characterized by steep, almost vertical drops from the bank to a relatively flat river bottom with a maximum seinable depth of~1.5 m. Lower estuary sites were typically sampled during low tides to facilitate the effective seining of structural habitat elements such as eelgrass beds which grew in the midchannel and would not be able to be reached with our nets from shore during high tides. Middle and upper estuary sites could be extremely shallow and have very little wetted area during low tides, especially during July and August when river flows were at a minimum, so they were typically sampled during high tides when sufficient water depth was available to use our skiff and the larger 30-m seine rather than our small pole seine. Each reach had three to four sampling areas for a total of 19 sites across the estuary which were beach seined every two weeks from mid-April until the end of September 2018. In total, we performed 188 seine sets over 10 rounds of sampling in 2018.
During each seine set, we measured water chemistry parameters (temperature, salinity, dissolved oxygen, pH) using a YSI ProDSS multiparameter water quality meter (YSI Incorporated, Yellow Springs, Ohio, USA). The sensor was submerged to the middle of the water column at the deepest area of the seine and was calibrated prior to the start of each sampling round (Appendix S1: Table S1). Three different juvenile beach seines were employed based on site-specific bathymetry in order to ensure equal sampling coverage of the water column at all sites: a 22 9 3.1 m net was used primarily for deeper saltwater sites, a 30 9 1.8 m net was used for most other sites when water levels were high enough for boat operation, and a 13.7 9 1.2 m pole seine was used in shallow/confined areas. Seines were held on shore at one end, pulled out to their full extent, curled around until the other end was brought to shore forming a semi-circle, and then, the two ends were pulled in simultaneously, capturing all fish in the water column which was surrounded by the net. During each set, we identified and enumerated all fish species then measured the fork-length of up to ten randomly selected individuals of each species. We only measured a small subsample of each species during each seine set due to time and logistical constraints associated with the large variety and abundance of fish caught and our tight sampling schedule. With the exception of sand lance, which occurred infrequently but in hyper-abundance, we measured 22.7% of the individuals representing the other 29 species comprising the fish community (Appendix S1: Table S7). Most sampled fish were measured without the use of anesthetics, but larger individuals were briefly anesthetized in MS-222 (0.05 g/L) and allowed to recover in aerated water prior to release (University Animal Care Committee at Simon Fraser University protocol number 1270B-14). Direct mass measurements were not taken because of difficulties in obtaining accurate measures under variable field conditions.
To account for the wide variety of body plans exhibited by fish species in our study, we calculated our size-spectra based on individuals' body mass rather than length (Kerr andDickie 2001, Sprules andBarth 2016). To do this, we calculated the mass of sampled individuals of most species from their fork-lengths using speciesspecific Bayesian length À weight conversions found on FishBase (Froese et al. 2014, Froese andPauly 2018). Exceptions were coho salmon (Oncorhynchus kisutch), sockeye salmon (Oncorhynchus nerka), and Dolly Varden trout (Salvelinus malma) for which previously collected data on lengths and mass from Koeye were used and lingcod (Ophiodon elongatus) for which conversions were found in Forrester and Thomson (1969). Estimated masses were rounded up to the nearest 0.1 g. Returning adult salmon were v www.esajournals.org removed from the dataset at this point because they do not feed during this life stage and thus are essentially inert members of the estuary fish community.
To produce accurate size-spectra, we needed a sample set of body mass estimates equal to the total catch collected during each seining event.
To do this, we sampled with replacement from our measured individuals until we had a new working sample set of measures equal to the total number of fish caught. This was done for each individual species for each individual seining event. To account for our different net sizes, we corrected catch abundances from all sets to the 30-m net size using the ratio of the surface area seined which was estimated to be the area of a semi-circle with the outer circumference equal to the length of the net (net lengths = 30, 22, 13.7 m, respectively; surface areas seined = 143.24, 77.03, 29.87 m 2 , respectively; catch abundance correction ratio = 1:1.86:4.78). To produce our size distributions, body mass sample sets for all species were combined for each reach during each sampling round and individuals were sorted into log 10 equal body size classes, or mass bins (0.1-0.9, 1.0-9.9, 10.0-99.9, 100.0-999.9, and 1000.0-9999.9 g, respectively). Size-spectra were estimated for each reach-round combination with a simple linear regression relating the distribution of abundance across size classes (log 10 (n + 1)~log 10 (mass bin), with mass bins represented as 0.1, 1, 10, 100, and 1000; Fig. 2a). We then bootstrapped this process 1000 times to reduce any random sampling bias, retaining slope (k) and midpoint height (Γ H ) values for each of the regressions and using the mean values of the size-spectra coefficients found for each reach-round combination for the remainder of the analysis. Midpoint height (Γ H ) values were used instead of y-intercept (Γ) values to avoid correlation of the spectra coefficients which produces the statistical artifact of increasing intercept with decreasing slope (Daan et al. 2005, Guiet et al. 2016. Following best practices, we interpret these midpoint heights (Γ H ) in the same way as we would y-intercepts and so will refer to them simply as intercepts for purposes of clarity in the remainder of the paper.
To examine differences between habitat slope (k) and intercept (Γ H ) values across the six reaches, we performed the non-parametric Kruskal-Wallis test followed by the Conover-Iman test for stochastic dominance using the R package conover.test (Appendix S1: Tables S3, S4; Kruskal and Wallis 1952, Conover and Iman 1979, Conover 1999, Dinno 2017. We used these non-parametric tests because the relatively small sample size (10 spectra regressions per reach), did not allow consistent assessment of normality.
We also examined how size-spectra shifted seasonally across the ecotone. We first constructed whole estuary models for both slope (k; predation risk) and intercept (Γ H ; productivity) and compared six sequential polynomial models to determine the polynomial orders which best described the overall variation in coefficient values over the season (Appendix S1: Table S5). We then tested whether the different habitats showed synchronous, lm(coefficient~poly (week, order) + reach), or independent, lm(coef-ficient~poly(week, order) 9 reach), variation in spectra coefficients over time (Appendix S1: Table S6). Models were compared using Akaike's information criterion for small sample sizes (AIC c ).
All analyses and graphics were produced in the R statistical environment version 3.6.1 (R Development Core Team 2019) using the packages conover.test, AICcmodavg, dplyr, and ggplot2.

RESULTS
In total, we caught 216,691 fish of 30 different species, sampling 3986 individuals for forklength (Appendix S1: Tables S2, S7). Sand lance (Ammodytes hexapterus) made up the vast majority of our total catch (199,853 individuals caught or 92.2% of the total catch), most of which were captured during a single seine set (150,000 on 16 July 2018), but only appeared in 24 of the 188 seine sets (Appendix S1: Table S7, Fig. S1). With the exception of sand lance, we sampled 22.7% of the individuals representing the other 29 species comprising the fish community (Appendix S1: Table S7, Fig. S1). Juvenile coho salmon (Oncorhynchus kisutch) were the next most abundant (10,314 individuals caught or 4.8% of the total catch) and were also the most ubiquitous species across the estuary, appearing in 154 of the 188 seine sets (Appendix S1: Table S7, Fig. S1). Threespined stickleback (Gasterosteus aculeatus) and v www.esajournals.org shiner perch (Cymatogaster aggregata) were also major constituents of the catch (Appendix S1: Table S7, Fig. S1). The most abundant fish predators caught in the estuary were the staghorn sculpins (Leptocottus armatus), Dolly Varden (Salvelinus malma), starry flounder (Platchthys stellatus), and freshwater sculpins (Cottus asper/Cottus aleuticus), with a few cutthroat trout (Oncorhynchus clarkii) and resident rainbow trout (Oncorhynchus mykiss) in the upper estuary. Community structure and abundance varied greatly across the estuary, with overall richness and abundance increasing with salinity but diversity and evenness peaking in the middle estuary (Appendix S1: Table S2). For example, in the marine beach habitat we caught a total of 194,559 individual fish comprising 24 species, while in the tidal freshwater habitat, we caught a total of only 2040 individuals from 10 different species (Appendix S1: Table S2). The marine reach also had the greatest number of species, with five being found nowhere else in the estuary, including a juvenile salmon shark (Lamna ditropis; Appendix S1: Table S2).
Size-spectra varied across the estuary ecotone (Fig. 2). Contrary to the prediction that there would be contrasting spatial gradients across the ecotone, the intermediate estuary habitat (mudflat and rockweed habitat, Reach 3) had the lowest median slope (k; predation risk) and the highest median intercept (Γ H ; productivity) among the various habitat types (Fig. 2). The Conover-Iman test revealed that the mudflat reach was significantly different (P < 0.05) from all other habitats in both slope (k; predation risk) and intercept (Γ H ; productivity), with the exception of similar estimates of intercept (Γ H ) in the beach habitat (Reach 1; P = 0.5788; Appendix S1: Table S4) and a marginally non-significant difference in slope (k) with the mainstem marsh habitat (Reach 4; P = 0.0614; Appendix S1: Table S3).
Size-spectra shifted seasonally across the estuary (Fig. 3). AIC c comparisons of the whole estuary models indicated that seasonal patterns for slopes (k; predation risk) were best described by a 4th order polynomial. Intercepts (Γ H ; productivity) were best described by a 2nd order polynomial (Appendix S1: Table S5). Habitats varied synchronously for slopes (k; predation risk), decreasing from early spring to summer and then rising again into the fall ( Fig. 3a; Appendix S1: Table S6). In contrast, size-spectra intercepts (Γ H ; productivity) exhibited different seasonal patterns across the reaches, with the marine reaches having more variable intercepts representing more pulsed fish communities, whereas freshwater reaches were more consistent across the season ( Fig. 3b; Appendix S1: Table S6).

DISCUSSION
Here, we applied size-spectra to examine potential spatial gradients of predation risks and productivity and provide insight into estuary nursery function. Our analysis of size-spectra in the unindustrialized Koeye River estuary revealed that the middle of the estuary habitat, Reach 3, characterized by muddy substrate and large rockweed expanses, had size-spectra coefficients that indicated higher production (higher Γ H ) and lower predation risk (lower k) in the fish community (Fig. 2). The combination of high productivity and low predation risk has been suggested as a key reason why estuaries function as nurseries for so many species (Dahlgren and Eggleston 2000, Cocheret de la Morini ere et al. 2002. during each sampling round. Spectra are seen to vary between reaches, but also display temporal variation within each. Boxplots of the (b) slope (k) and (c) intercept (Γ H ) values found in each reach using mean values for each reach-round combination the mass sample extrapolation. Using the Conover-Iman test, the mudflats were found to be significantly different (P < 0.05) from all other reaches in slope (k) and intercept (Γ H ), except for marginally non-significant results for the main marsh slope (P = 0.0614; Appendix S1: Table S3) and with the beach intercept (P = 0.5788; Appendix S1: Table S4).
( Fig. 2. Continued) v www.esajournals.org Estuaries are comprised of a mosaic of habitat types that vary in structure, extent, and abiotic conditions (Peterson 2003, which likely contributes to the observed patterns in size-spectra and inferred estuary nursery function. For example, predator abundance, size, and trophic position tend to increase with habitat size and stability (Kushlan 1976, McIntosh et al. 2018. Indeed, we observed the overall largest fish (a juvenile salmon shark) as well as the highest abundance of upper size-class individuals in the more spatially expansive marine sites, while the fish communities of the more constrained, environmentally variable middle and upper estuary reaches were almost exclusively comprised of smaller size-class individuals ( Fig. 2a; Appendix S1: Table S1). Specifically, the reach that had the highest slope and intercept, the rockweed expanses of Reach 3, is extremely variable. This variability may serve to exclude Fig. 3. Model fits for the temporal variation associated with (a) slope (k), lm(k~poly(round, 4) + reach) (multiple R 2 = 0.6085, adjusted R 2 = 0.538, df = 50, P = 1.144 9 10 À7 ), and (b) intercept (Γ H ), lm(coefficient~poly (round, 2) 9 reach) (multiple R 2 = 0.7823, adjusted R 2 = 0.6941, df = 42, P = 4.981 9 10 À9 ). Habitats were found to vary synchronously in predation risk (slope [k]), but independently in productivity (intercept [Γ H ]). Model selection results can be found in Appendix S1: Tables S5, S6. larger predatory fishes which may not be as osmotically or thermally tolerant, might be at risk of becoming stranded in the shallows as tide waters recede. Intermediate reaches of estuaries also may be particularly productive due to a mixture of energy sources and local nutrient dynamics (Peterson 2003, Nelson et al. 2015, David et al. 2016; indeed, we anecdotally observed that this middle reach contained an abundance of benthic invertebrates, such as isopods and amphipods, and swarms of mysids providing a rich food source for rearing juvenile fish such as coho fry. Size-spectra are the emergent patterns that arise from this complicated and dynamic nature of fish communities and habitats in estuaries.
Spectra coefficients have been shown to reflect changes in environmental conditions (Maury et al. 2007, Guiet et al. 2016) and oscillate through time in response to shifting ratios in predator and prey size and abundance (Blumenshine et al. 2000, Law et al. 2009). In temporally dynamic systems, fish size cohorts are thought to track moving waves of productivity and predation risk, surfing the size-spectrum to maximize their foraging opportunities while minimizing their risk of predation (Pope et al. 1994). Seascape studies have shown that animals move among habitat patches daily, following tidal movements in search of food and shelter, but also show ontogenetic shifts in habitat use on longer time-scales as their resource needs change with growth (Dahlgren and Eggleston 2000, Cocheret de la Morini ere et al. 2002, Peterson 2003. For example, smaller-sized juvenile salmon may use shallow water habitats than larger-bodied individuals (Munsch et al. 2016). Here we present rare empirical evidence for seasonal shifts in size-spectra across multiple habitat types (Fig. 3), lending support to the surf-riding and mosaic seascape nursery hypotheses (Pope et al. 1994, Peterson 2003. Migration events, such as Pacific salmon smolt outmigration, bring large influxes of small size-class fishes to the estuary for a brief period, saturating the ecosystem with prey and presumably lowering the predation risk for any individual fish (May in Fig. 3a). Similarly, reproduction events within the estuary, such as those of shiner perch and sticklebacks, also create pulses in small fish abundance which affect individual predation risk (August in Fig. 3a). Such temporal dynamics in size-spectra are also likely linked to seasonal patterns of plankton blooms, plant growth, and macroinvertebrate abundance across the estuary habitat mosaic (Fig. 3b). Thus, estuary fish communities reveal spatially and temporally dynamic patterns of productivity and predation risk. Rather than consistent trade-offs between productivity and risk, the habitat with the least risk is not necessarily also the least productive (Fig. 3). In autumn, it appears that many fish move to overwintering habitats outside of the estuary and the community shifts back to its pregrowing season configuration, demonstrating the seasonality of nursery function in temperate estuaries (Fig. 3). Examining temporal shifts in spectra coefficients across the entire mosaic of estuary environments, allows for not only the identification of nursery function in discrete habitats (Beck et al. 2001), but also the study of more spatially and temporally dynamic processes shaping estuary fish communities (Peterson 2003.
By integrating abundance and relative body size frequencies, size-spectra analysis provides a relatively easy to sample and statistically simple method for studying the underlying ecological energetics and inferring potential underlying ecological processes of complex aquatic ecosystems (Kerr and Dickie 2001, Giacomini et al. 2016, Guiet et al. 2016, Sprules and Barth 2016. Thus, size-spectra analysis may prove to be an important new tool for understanding estuary fish ecology and habitat function. Here, we suggest that complex mosaics of habitats in estuaries present shifting trade-offs in predation risk and productivity that support different sizes of fish, a combination of intra-and inter-specific variation, across space, and time.

ACKNOWLEDGMENTS
This project was the product of the collaboration among SFU, the Hakai Institute, and the QQs Project Society from the Heiltsuik Nation. Special thanks to our Heiltsuk field technicians Richard Wilson-Hall, Jefferson Brown, and Cecil Brown. Kyle Wilson, Emma Hodgson, and Leithen M'Gonigle all provided key inputs during analysis. Karl Seitz received funding from the Tula-Mitacs Canada Grant IT09911, through internships No. FR23017 and FR23025. This project was supported by the Hakai Institute and Tula Foundation, the Liber Ero Chair of Coastal Science and Management, the Pacific Salmon Foundation v www.esajournals.org