Strong water stratification provides a refuge for rainbow smelt larvae Osmerus mordax in a sub-arctic estuary (Lake Melville, Labrador)

water stratification prevailing in upper Lake Melville provides a nursery for early life stages of fish, where they are protected from their predators. Ichthyoplankton aggregated just above and at the pycnocline, in the low salinity surface layer down to 25 m. Most adult pelagic fish occupied the bottom waters below the sharp pycnocline, although some ventured in the low salinity surface layer. Ten species of adult fish were captured in gill and fyke nets and 53 species were detected with eDNA. Larvae of rainbow smelt ( Osmerus mordax ) were ubiquitous in the surface layer in July and August and represented 100% of the ichthyoplankton assemblage sampled during these months. No fish larvae were detected in winter (February). We conclude that the low salinity surface layer provides a refuge for rainbow smelt larvae, a key forage species in the estuary. This study provides baseline information from which to assess future changes in biodiversity and distribution of fish in the Lake Melville estuary. It further supports the use of eDNA as a complementary tool for monitoring fish diversity in sub-arctic estuaries.


Introduction
The distribution of fish in estuarine environments primarily depends on the local hydrography, ontogeny, and their swimming ability.The latter varies between species, size, and developmental stages (Tzeng and Wang, 1993).Early life stages with limited swimming ability, such as eggs and larvae, often occupy the pelagic and surface waters and take advantage of the surface currents (Arevalo et al., 2023).Upon reaching a certain size, fish larvae and juvenile can position themselves vertically to avoid visual predation, find larger prey, or optimize their movement, and therefore their energy expenditure, using tidal stream transport (Miller et al., 1984;Arevalo et al., 2023).For example, herring (Clupea harengus; Fortier and Leggett, 1982) and gobies (Pomatoschistus spp.; Miró et al., 2022) descend near the bottom to take advantage of residual currents and avoid being advected outside the estuary (Arevalo et al., 2023).Sexually mature fish can also modify their distribution by choosing spawning areas and depths that will increase the fitness of their offspring (Iles and Sinclair, 1982;Grote et al., 2012;Sundby and Kristiansen, 2015).Despite the importance of ontogenetic vertical distribution and migrations, most studies on the use of estuaries as a habitat for fish have mainly considered horizontal habitats (e.g.Pihl et al., 2002).The vertical distribution of estuarine fish in relation to water masses remains less documented, especially across different life stages.
Buoyancy varies with temperature and salinity and strongly affects the vertical distribution and the dispersal of fish eggs (Sundby, 1990).Hydrography also drives the survival and dispersion of species entering the estuary at the larval or juvenile stages.Early life stages of most pelagic fish have evolved to adapt their specific gravity and vertical distribution to remain within the pycnocline, a layer where the water density increases rapidly with depth (Sundby, 1990;Sundby and Kristiansen, 2015).For instance, the larval dispersal of rainbow smelt larvae (Osmerus mordax), a forage species ubiquitous to temperate regions of the northern hemisphere, is closely linked to water stratification and seasonal changes in estuarine inflow and outflow (Bradbury et al., 2006).Positive buoyancy and passive swimming ability may explain surface aggregations of rainbow smelt larvae above the pycnocline.Previous studies have also revealed higher species richness and abundance of larvae in frontal habitats, such as haline tidal-mixing fronts (Sánchez-Velasco et al., 2012).Any changes in the hydrology or seasonal freshwater inflow in estuaries can thus have critical impacts on the development, survival, and diversity of fish larvae (Carassou et al. 2011;Pasquaud et al., 2012;Peterson et al., 2013;Rodrigues et al., 2022).In some estuaries, turbidity can also be an important driver of the distribution of early life stages of fish because it increases feeding success (Blaber and Blaber, 1980).
Lake Melville is the largest estuary of Labrador (northeastern Canada).Because of its relatively shallow and large embayment and its highly stratified water column, it is a typical sub-arctic fjard (Kamula et al., 2020).Freshwater, anadromous, and marine species cohabit in the Lake Melville estuary and migratory fish species, such as Arctic char (Salvelinus alpinus) and Atlantic salmon (Salmo salar), as well as forage species such as rainbow smelt, use riverine tributaries for spawning and for nursery habitats (Backus, 1957;McCarthy and Gosse, 2018).Lake Melville was identified as an Ecologically and Biologically Significant Area (DFO, 2013) but, despite its ecological importance, several knowledge gaps exist on the biodiversity and seasonal variation in the spatial distribution of fish within Lake Melville.Information on fish populations is limited to abundance and distribution of adult fish from net surveys conducted in summer (McCarthy and Gosse, 2018).Overall, the effects of changes in freshwater inflow across seasons and from anthropogenic stressors (Durkalec et al., 2016) on fish distribution remain poorly documented (Rytwinski et al., 2020).
In summers 2018 and 2019 and winters 2019 and 2020, we surveyed Lake Melville to assess the biodiversity and seasonal changes in horizontal and vertical distributions of pelagic fish in relation to hydrography.Here, we combine environmental and hydroacoustic data, ichthyoplankton and fish net samples, and water sampling for environmental DNA (eDNA) to test the hypothesis that the strong water stratification prevailing in this sub-arctic estuary provides a refuge for early life stages of fish.We further provide baseline information to assess the impacts of future changes in freshwater inflow and water stratification on the pelagic ecosystem of the region.

Study area
Lake Melville is a 250 km long and 2100 km 2 sub-arctic estuarine fjard.In contrast to fjords, fjards are defined by irregular bathymetry and low catchment topographic relief (Kamula et al., 2020).The maximum depth of Lake Melville reaches 256 m with a mean depth of 84 m (Fig. 1).The Churchill River flows eastwards into Lake Melville via Goose Bay and the Goose Bay narrows (referred to as the narrows hereafter), a shallow (8 m) and narrow (2.5 km) sill.At its eastern end, Groswater Bay connects Lake Melville to nutrient rich waters from the Labrador shelf through another narrow (2.8 km) sill at the Rigolet narrows (not shown on the map).The 30 m deep Rigolet narrows restricts the inflow of seawater in Lake Melville (Durkalec et al., 2016).Lake Melville's hydrology is characterized by an estuarine circulation, where brackish surface waters from river runoff (salinity <10) flow seawards and the bottom saline waters of the Labrador Sea (salinity ~28) flow landwards.The estuary is seasonally ice-covered (December-May) and semi enclosed with a low salinity surface layer within the first 10-20 m (Bobbitt, 1982;Lu et al., 2013).Stratification of the fresher surface waters and saline bottom waters is more pronounced near the mouth of the Churchill River and gradually breaks down towards the Labrador shelf.The steep halocline persists throughout the year creating a consistent boundary layer between freshwater and saltwater (Durkalec et al., 2016).

Survey design
To optimize the sampling near the mouth of the Churchill River, the survey was concentrated in Goose Bay (GB) and the western part of Lake Melville (termed upper Lake Melville or ULM, hereafter), from 60 • 21.0′W to 59 • 36.6′W.We concentrated our sampling effort in summer Fig. 2. Section plots of conservative temperature (θ) and absolute salinity (S abs ) during summers 2018 and 2019 and winters 2019 and 2020 across Goose Bay (GB; west) and upper Lake Melville (ULM; east).For each panel, the density anomaly referenced to surface (σ 0 ; in kg m − 3 ) is plotted as solid light gray lines identified in the salinity panels.The brown shaded polygon represents bottom drawn using the depth at each station determined with hydroacoustic measurements.The gray dotted lines represent profile stations identified in the top of the salinity panels.White areas indicate the absence of data.The pycnocline is identified with a dashedred lines.The separation between Goose Bay (GB) and upper Lake Melville (ULM) occurs at the sill near station 5 and is identified in the salinity panel.Note that for winter 2019, a problem with the profiler prevented the presentation of data below ~20m.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) T. Small et al. and winter to contrast two seasons with different surface temperatures and freshwater inflow.In summer, the survey was conducted from the R/V Gecho II between June 30 to July 12 in 2018 and July 4 to July 14 in 2019.To ensure the acoustic detection of larval fish, summer surveys were conducted after the appearance of the swim bladder in juvenile rainbow smelt, occurring within the first 14 days after hatching (Belyanina, 1969).The hydroacoustic transects followed a stratified sampling with systematic samples nested within eastern and western strata (Parker-Stetter et al., 2009;Rudstam et al., 2009) (Fig. 1).Surveys were conducted during daylight only.Hydroacoustic transects were conducted at vessel speeds of 4-6 knots.Environmental data and net samples were collected at a maximum of 10 stations per season distributed from the head of Goose Bay towards the western end of upper Lake Melville.In winter, the hydroacoustic and environmental surveys were limited to the fixed stations because complete ice cover made continuous transects impossible.Not all stations were sampled in winter to avoid dangerous ice conditions (Supplementary Table 1).The study area was divided into two regions (GB and ULM) divided by the shallowest point of the narrows (Figs. 1 and 2).Permission to conduct this research was granted by the Nunatsiavut Government, Innu Nation, and the NunatuKavut Community Council prior to each sampling season.

Environmental data
Conductivity, temperature, and depth (CTD) profiles were collected at each station using a calibrated CastAway-CTD (summer), a vertical microstructure profiler (Model VMP-250) from Rockland Scientific International or an RBR Concerto 3 (winter 2019 and 2020, respectively).Turbidity and Chlorophyll-a (chl-a) data were also collected for all winter 2019 stations using the VMP's JFE Advantech Co. Fluorescence and Turbidity (JFE-FT) sensor.The manufacturer's calibration was used for this sensor and measurements are expressed in formazine turbidity units (FTU) for turbidity and part per billion (ppb) for chl-a (uranine reference).
Temperature, salinity, turbidity, and chl-a data were cleaned to omit the first meter of samples and to only include valid down casts.The conservative temperature (θ), the absolute salinity (S abs ) and the density anomaly referenced to the surface (σ 0 ) were derived using the TEOS-10 toolbox (McDougall and Barker, 2011).The data were then vertically binned (averaged every 25 cm).The depth of the pycnocline during each field season was determined by calculating the maximum difference in water density between two depths at each station and averaged by region (i.e., Goose Bay and upper Lake Melville).These depths were then used to distinguish the average temperature and salinities for the surface and bottom layers (i.e.above/below the pycnocline) for the two regions and for all sampling seasons.
Year-round environmental data were collected by a permanent mooring (53 • 22.2′N; 60 • 07.4′W; depth ~ 40 m) equipped with six HOBO TidbiT® v2 temperature loggers spread over the water column, in addition to a RBR Concerto 3 CTD and a Sea Bird Scientific ECO-FL Fluorometer targeting the near-surface layer.During the second year of the deployment, a Seabird CTD replaced the RBR CTD near the surface, and the latter was put near the bottom (see Supplementary Table 2 for instrument types and depths).The mooring was originally installed on July 2, 2019, and serviced 13-14 October, 2020 (the ECO-FL fluorometer stopped working on September 26, 2020).During the redeployment of the mooring, a piece of rope near the bottom-end of the mooring got tangled up with the acoustic release system bringing all instrument down by 3 m from their targeted depths (Supplementary Table 2).The second deployment acquired meaningful measurements until October 4, 2021.

Sampling and processing
Scientific echosounders can measure the spatial distribution of larval and adult fish with very high vertical (cm to m) and spatial (seconds) resolutions, including in estuaries (Simmonds and MacLennan, 2005;Boswell et al., 2007).Here, summer hydroacoustic data were collected continuously (~8 h per day) using a multi-frequency echosounder (Simrad EK60; 38 kHz and 120 kHz) hull mounted on the port side 1 m below R/V Gecho II.The ping rate was set to 1 Hz and pulse durations and 3 dB beam angles for each acoustic transducer are detailed in Supplementary Table 3.
In winter, stationary hydroacoustic data were collected from holes drilled across the sea-ice.In winter 2019 we used a Wideband Autonomous Transceiver (WBAT, 38 kHz), and in winter 2020 hydroacoustic data were collected using a portable wideband echosounder system (Simrad EK80; 38 kHz Supplementary Table 3).In both years, transducers were fixed 1 m below the surface of the ice and recorded data for 45 min at each station.All echosounders were calibrated before or after each field season using the standard sphere method (Demer et al., 2015).
Hydroacoustic data were processed with Echoview (version 11.1; Echoview Software Pty Ltd.).The top 5 m below the surface were excluded to eliminate draft and the near-field transducer effects (Ryan et al., 2015) and to reduce the noise created from bubbles under the boat in summer.The sounder-detected bottom was corrected where necessary and a bottom exclusion line was created to eliminate data within 0.5 m from the bottom.The S v (Mean volume backscattering strength, dB re 1 m − 1 ) acoustic echograms were visually inspected and cleaned using Echoview's algorithms to remove background noise with a signal to noise ratio <10 dB, impulse noise, and attenuated noise signals (De Robertis and Higginbottom, 2007;Ryan et al., 2015).Clean hydroacoustic data were integrated into bins 0.50 nmi long and 1 m deep.As a proxy for abundance of pelagic organisms in both summer and winter, the Nautical Area Scattering Coefficient (NASC, s A in m 2 nmi − 2 ) was integrated at 38 kHz and 120 kHz and exported for post-processing.
To map the spatial distribution of larval and adult pelagic fish in summer and identify potential hotspots, we kriged the integrated NASC values using the "krig" function from the R package GSTAT (Pebesma, 2004).To achieve normality, the NASC was converted into Nautical area scattering strength (S A in dB re 1 m 2 nmi − 2 ).Variograms were chosen based on the best fit model (Spherical, Gaussian, Maternal, or Exponential) using the "fit.Variogram" function from the R GSTAT package (version 1.3.1056)(Pebesma, 2004).To investigate differences in vertical distribution, we also integrated and then kriged the S A values above and below 25 m, which corresponds to the maximum depth of the sound scattering layer observed on the echogram.
T-tests were used to check for significant differences between the mean S A above and below the pycnocline and between regions (upper Lake Melville and Goose Bay).Statistical tests were performed using R (version 1.3.1056).

Target Strength analyses
Target Strength (TS in dB re 1 m 2 ) analyses provide information on both the size and body composition of individual organisms.Here, Target Strength of single targets at 38 kHz were extracted using Echoview's single-echo detection algorithm for split beam echosounders (method 2) with a threshold of − 110 dB to detect both larval and adult fish (Supplementary Tables 4A and B).Echoview's fish tracking algorithm was then applied to the single target echograms to extract single fish tracks (Supplementary Table 5).Maximum number of pings between fish tracks were modified to 1 to account for higher densities seen on the echograms.We exported the mean compensated TS and the mean depth for each single target track.Percentages of total TS above and below 25 m, the maximum depth of the sound scattering layer, were calculated for pooled summers (2018-2019) and pooled winters (2019)(2020).Note that winter TS data come from narrowband in 2019 and from the nominal frequency from wideband signals in 2020, but that there was no significant difference when comparing both years (F = 1.06, num df = 2385, denom df = 8585, p-value = 0.06).

Fish sampling
While hydroacoustic surveys provide a high spatio-temporal resolution, their taxonomic resolution is poor.Net and trawl sampling is typically used to groundtruth the acoustic signal (e.g.Boswell et al., 2007;Herbig et al., 2023).In summer, a consistent sound scattering layer was present between ~10 and down to 25 m on the echograms.At each of the stations, a BONGO net, comprised of 2 adjacent 28-cm diameter ring nets with 335 μm mesh, was towed at 2 knots in the sound scattering layer for 10 min to capture ichthyoplankton.Samples were filtered through a 120 μm sieve.Larval fish were separated from zooplankton and were fixed in 90% ethanol that was replaced after the first 24 h.
Fish monitoring was completed using double bag fyke nets and experimental gillnets deployed in August 2018 and August-September 2019 to sample adult fish.All nets were set for 16 h to include both dusk and dawn periods when fish are active.Fyke nets sat on the bottom, ranged from 0.5 to 1.5 m in height, and were designed to capture small fish (between 40 and 50 mm in length) in coastal waters no more than 2 m deep.Experimental gillnets were designed to catch larger fish and consisted of four 50 ft panels (200 ft in total) with mesh size, 2″, 3″, 4″ and 5".Gillnets were deployed in bottom depth >2m.Abundance and biomass in Catch Per Unit Effort (CPUE) were calculated by dividing the number and weight of each individual species by the net soaking duration (16 h).Shannon-Weiner diversity index (Shannon, 1948) was calculated for each region (upper Lake Melville and Goose Bay) and year (2018-2019), except for Goose Bay in 2018 when no nets were deployed.

eDNA collection
We analysed eDNA in water samples to complement results from net catches.Over the last two decades, the use of eDNA has improved species detection and has increased our understanding of the dynamics and structure of communities, emerging as an essential tool for ecology and biodiversity studies (Yates et al., 2023).Environmental DNA is a cost-effective tool compared to other sampling methods.It can be used in environmental conditions where most methods are insufficient and, despite its own biases, eDNA systematically detects other species than just the common ones; species whose rarity, behavior, and size escape traditional methods (Cristescu and Hebert, 2018;Sard et al., 2019;Sigsgaard et al., 2015;Spear et al., 2015).As part of our study, water samples were collected for eDNA analyses of larval and adult fish during summer 2018.The depth of each sample varied depending on where the most backscatter was observed on the echosounder.Negative field controls (blanks) using distilled water were completed for every station.These controls were performed in the same field conditions as the samples and processed in the same way to estimate potential contamination.After collection, the samples were kept frozen and processed in laboratory at Université Laval.Unfortunately, no eDNA sampling was conducted in winter nor summer 2019.

eDNA extraction
DNA was extracted using a QIAshredder and DNeasy Blood and Tissue kit (Qiagen) following a modified manufactured protocol (Goldberg et al., 2011;Spens et al., 2017).Extraction was performed under a UV hood with all instruments either bleached and/or UV-treated.In all extraction batches, an extraction control was included to account for possible contamination.

Polymerase Chain Reaction (PCR) amplification and sequencing
eDNA amplification was conducted using MiFish primers (Miya et al., 2015; MiFish-U-F 5′-GTC GGT AAA ACT CGT GCC AGC-3′ and MiFish-U-R 5′-CAT AGT GGG GTA TCT AAT CCC AGT TTG-3′).These universal primers target a hypervariable region of the 12S rRNA gene (174 bp) designed specifically to amplify and disentangled at the species level all fishes.A single PCR reaction was conducted using dual index Illumina barcodes to reduce contamination.Each PCR reaction was composed of 12.5 μl of MasterMix (Qisgen), 2 μl of each primer (10 μM), 5.5 μl of diH 2 0, and 3 μl of DNA sample.PCR was run under the following conditions: 15 min at 95 • C and 35 cycles (30 s at 94 • C, 90 s at 65 • C, 60 s at 72 • C) followed by a final extension of 10 min at 72 • C. For each sample, including extraction control, PCR was done in five replicates that were pooled after amplification.In addition, a PCR negative control was included for each barcode combination.PCR products for each sample were pooled and verified on an 1% agarose gel.No amplification was observed in PCR negative controls.Pooled PCR products were then cleaned with AxyPrep Mag PCR Clean-Up kit (Corning) and quantify by fluorescence (Accuclear Ultra High Sensitivity dsDNA Quantification Kit (Biotium)).Samples were then pooled in equimolar proportion, cleaned and fragment size distribution checked on an Agilent 2100 Bioanalyser (High Sensitivity DNA kit).Sequencing was performed on a Illumina Miseq (Illumina, Nextera XT V3) at the genomic platform of the Institut de Biologie Intégrative et des Système (IBIS) at Université Laval.

Sequence cleaning and annotation
Raw sequencing reads were demultiplexed by the Miseq Control Software.Reads were filtered, merged, and annotated by the BARQUE pipeline (v1.5.4) (www.github.com/enormandeau/barque;accessed October 30, 2023).Species identification was performed using the Mitofish 12S database (Iwasaki et al., 2013), available from the latest BARQUE versions.A minimum of 97% similarity between the sequences of interest and the database sequences to assign taxonomic identification.
Non-fish sequences and out-of-range species (defined as species not documented in Eastern North America) were removed from the analysis.Individual species or sequences attributed to more than one species with numbers of sequences less than 50 were also removed from the analysis.Stations were divided into regions (upper Lake Melville and Goose Bay) and in layers above and below the pycnocline.

Environmental parameters
The Goose Bay/upper Lake Melville ecosystem was characterized by two water density layers separated by a strong haline stratification across all seasons (Fig. 2).The depth of this interface between the two layers, the pycnocline, varied between seasons and regions (Supplementary Table 6).In Goose Bay, the pycnocline was deeper in the summer (10.4 and 12.9 m for both years) compared to the winter (6 and 7 m).In upper Lake Melville, the depth of the pycnocline underwent larger variations, with averaging at 7 and 12 m in the summer and and 11 m in the winter for both years, respectively.
The upper Lake Melville was consistently colder and more saline than Goose Bay for both layers (Fig. 2 and Supplementary Table 6).Above the pycnocline a low salinity surface layer (LSSL) was present all year long, with summer salinities as low as 0.8 g kg − 1 and 2.4-2.7 g kg − in Goose Bay and upper Lake Melville, respectively.During the winter, the salinity in the LSSL decreased from 2.1 g kg − 1 to 0.5 g kg − 1 between winter 2019 and 2020 in Goose Bay, and from 10 g kg − 1 to 4 g kg − 1 in upper Lake Melville, although the latter presented relatively large variations (Supplementary Table 6).In the bottom layer, salinity was much higher, ranging from 16 g kg − 1 to 24.4 g kg − 1 across both regions and seasons.The bottom salinity decreased in Goose Bay and increased in upper Lake Melville between winters 2019 and 2020.
In the summer, the LSSL was consistently warmer (up to 12.9 • C in Goose Bay and 13.1 • C in upper Lake Melville on average) than the bottom layer.The latter was relatively stable in temperature throughout the seasons (between 1.6 and 2.2 • C in Goose Bay and 0.5 and 1.1 • C in upper Lake Melville).In winter, the temperature of the LSSL was near the freezing point everywhere.The sill separating Goose Bay from upper Lake Melville limits exchanges of bottom water between both regions, as suggested by the contrasted properties of the bottom layer described above (see also Fig. 2).
Turbidity and chl-a concentration measurements were available for the top ~20m of the water column in winter 2019 (Fig. 3).The turbidity was higher in the LSSL and lower in the bottom waters.It was also higher in Goose Bay than in upper Lake Melville, which supports existing literature that has suggested a high turbidity region near the mouth of the Churchill River (Durkalec et al., 2016).Similarly, the concentration of chl-a in the surface layer was higher in Goose Bay than in upper Lake Melville.
The vertical distribution of temperature throughout the water column, as well as the salinity and chl-a fluorescence at discrete depths, were also measured nearly continuously with a mooring near the mouth of the Churchill River between July 2019 and October 2021 (Fig. 4).A warmer surface layer developed annually at depth <15 m between June and November, with weekly averaged temperature peaking at about 17 • C, before decreasing to freezing point in winter (Fig. 4A).In the bottom layer, water temperature remained relatively constant below 4 • C throughout the year.Note that the temperature in the bottom layer was colder (<2 • C; darker shades of purple) during the winter compared to the winter 2021.This can be explained by the fact that 2021 was a relatively mild winter in the region (Cyr et al., 2022).
The salinity in the surface layer (measured at 8 m and 11 m during both deployments, respectively) remained between 10 and 17 g kg − (2019-2020) and 17-19 g kg − 1 (2020-2021) from the late summer to the following spring when it rapidly dropped near 0 g kg − 1 at the onset of the spring freshet and the melting season in late May (Fig. 4B).The greater values and weaker variability of the salinity during the second year are partially explained by the deeper location of the CTD during the 2020-2021 deployment.
Chlorophyll-a concentrations show clear spring blooms occurring in late May in both years (Fig. 4C).These blooms corresponded to the rapid freshening of the surface layer discussed above (Fig. 4B) and likely associated with sea ice melt and further stratification of the water column.Later in the year, Chlorophyll-a concentrations slowly decreased from the summer values until the next spring, although small increase in the fall (i.e.fall blooms in August-September) are visible.Again, differences in concentration between the two years may be partially explained by differences in the depth of the sensors.

Spatial distribution of larval and adult pelagic fish
In summer, a strong sound scattering layer was consistently observed above and within a depth range close to the pycnocline, between ~10 m and a maximum depth of 25 m (Fig. 5A).Although the pycnocline was present year-round, the absence of the sound scattering layer in winter (Fig. 5B) supports the assumption that it was mainly of biological origin rather than purely physical (Fig. 5).Moreover, the pycnocline was consistently present across the study area while the distribution of backscatter was patchy with a varying maximum depth, as expected for biological targets.Net sampling within the sound scattering layer confirmed the occurrence of rainbow smelt larvae at all stations in summer, which composed 100% of the ichthyoplankton community.Between 13 and 241 rainbow smelt larvae were caught at each station, for a total of 770 and 1129 larvae in summer 2018 and 2019, respectively.No ichthyoplankton was detected on the echosounder nor Fig. 4. A) Temperature contours at the mooring location from two deployments (2019-07-15 to 2020-10-14 and 2020-10-13 to 2021-10-04), from 6 TidbiTs temperature sensors (position identified with dashed-gray lines), one RBR CTD located at 8 m during the first deployment (dashed-blue line) and 38.4 m during the second deployment (dashed-red line) and one Seabird CTD located at 11 m during the second deployment (dashed-green line).Note that during the second deployment, the mooring release got tangled bringing all instruments down by about 3 m from their targeted depth.B) Corresponding absolute salinity from the CTDs mentioned above.The gray lines correspond to the hourly averages while the thicker colored lines are the weekly averages.C) Chlorophyll-a concentration (using the calibration from the manufacturer) from a Sea Bird Eco-FL fluorometer.The gray lines correspond to the hourly average while the thicker colored lines are the weekly averages.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)sampled in winter.During both summer and winter, the echosounder also detected larger fish, probably older juveniles or adults, as scattered targets in the bottom water layer (Fig. 5A and B and 6).
During both summer and winter, weaker targets from larval fish were distributed in the top 25 m and stronger targets, most likely from larger fish, were concentrated in the bottom layer at depths >25m (Fig. 6).In summer, single targets detected at 38 kHz within the top 25 m had a dominant Target Strength (TS) mode at − 71 dB (re 1 m 2 ) and a second mode at − 52 dB (Fig. 6A).Target strengths below the pycnocline in summer were also bimodal with the highest rate of occurrence centered at − 55 dB and a second, lower, mode at − 71 dB.In winter, few targets were detected in the top 25 m, where the mode was centered at − 56 dB (Fig. 6B).During that season, most targets were distributed below 25 m and also had a TS distribution centered around − 56 dB.

Biodiversity of pelagic fish
A total of 275 adult fish were sampled in the gill and fyke nets between August 11-15, 2018, and 211 were caught between August 19-September 27, 2019 (Table 1).Ten species co-occurred: Salvelinus fontinalis (brook trout); Catostomus Catostomus (longnose sucker); Osmerus mordax (rainbow smelt); Microgadus tomcod (tomcod); Catostomus commersoni (white sucker); Couesius plumbeus (lake chub); Gasterosteus acusleatus (threespine stickleback); Pseudopleuronectes americanus (winter flounder); and Cottidae (mottled and spiny sculpin).All taxa were sampled both years except for winter flounder which was present exclusively in 2018 and sculpin which was present in 2019 only.The most abundant adult fish in the nets deployed in upper Lake Melville in 2018 were white sucker and longnose sucker (n = 77 and 59, respectively).In 2019, longnose sucker was the most abundant fish in the nets in upper Lake Melville (n = 68) (Table 1).Fyke and gill nets were not deployed in Goose Bay in 2018, but Lake chub, rainbow smelt, threespine stickleback, and tomcod were the only species that were found in both Goose Bay and upper Lake Melville in 2019.We observed net selectivity, with fyke nets capturing more longnose sucker, tomcod, lake chub, and threespine stickleback, while gill nets captured more white sucker, brook trout, and rainbow smelt.Net selectivity was consistent from one year to the other (Table 1).Fish diversity measured from net samples was higher in upper Lake Melville than in Goose Bay (Table 2).

Environmental DNA
A total of 8.38 million of the 9.40 million sequences obtained from the Illumina Miseq platform were annotated to a species by the BARQUE pipeline.A total of 53 unique fish species were identified by eDNA, and all larval and adult fish species collected in the nets were also detected by eDNA.Overall, 30 marine species, 13 freshwater species, and 10 anadromous species were detected (Table 3).Forty-nine species were detected in the LSSL (above the pycnocline), and 44 in bottom waters below the pycnocline.Thirty-nine species were detected in Goose Bay and 48 species were detected in upper Lake Melville.Freshwater species made up the top 10% of sequences of each species found in upper Lake Melville above the pycnocline, while below the pycnocline the highest number of sequences came from Gadidae with the second most being freshwater species.In Goose Bay, the greatest number of sequences above and below the pycnocline came from Gadidae, followed by Pleuronectidae and freshwater species (Table 3).

The low salinity surface layer and pycnocline of Lake Melville as a habitat for rainbow smelt larvae
In July and August, rainbow smelt larvae formed a ubiquitous sound scattering layer above and within a depth range related to the sharp pycnocline separating the low salinity surface layer (LSSL) from the colder and saltier bottom waters.Larvae of other species might develop in freshwater (e.g.trout and pike), in the outer estuary (e.g.capelin), or closer to the bottom (e.g.plaice, flounder, sculpins), which would explain why no other ichthyoplankton species was collected.The sound scattering layer was absent in winter, despite the stratification prevailing year-round in Goose Bay and upper Lake Melville.Very few small targets, such as fish larvae, were acoustically detected in the bottom layer (Fig. 6C and D), and we deduced that most rainbow smelt larvae remained above and within the pycnocline.The TS mode at − 71 dB prevailing in the LSSL and pycnocline in summer would correspond to a 7.4 mm rainbow smelt, based on an equation developed for rainbow smelt at 70 and 120 kHz from Rudstam et al. (2003) (derived from Love, 1977).This is within the size range of the rainbow smelt we sampled (4.9-13.1 mm with an average of 7.7 mm; Sutton, 2022) and supports the assumption that the sound scatter layer originates from larval rainbow smelt.The TS analysis suggests that larger juvenile and adult fish mostly remained in the bottom layer (Fig. 6).This is supported by eDNA observations from which 57% of the species detected were purely marine and likely prefer to remain in the bottom waters while in the estuary.Based on these results and those from a companion paper (Sutton, 2022), we conclude that, in Lake Melville, rainbow smelt larvae hatch in late June -early July at river mouths and move downstream into the estuary where they remain in the low salinity surface layer until they become juveniles at around 90 days post hatch.This behaviour is similar to that from other rainbow smelt populations, for instance on the south east coast of Newfoundland (Bradbury et al., 2006).
While the pycnocline can represent a physical barrier preventing zooplankton and ichthyoplankton from crossing water mass boundaries (Röpke et al., 1993), several ecological advantages could also explain the aggregation of rainbow smelt larvae at the pycnocline.First, the size segregation suggests that predator avoidance partly drives the distribution of smaller fish in the top layers.Similar distribution and behavior has been observed near river mouths where negative geotaxis by rainbow smelt larvae was linked to predator avoidance in natal river systems (Bradbury et al., 2004).Yet, the occurrence of some targets with a stronger TS mode of − 52 dB in the LSSL in summer (Fig. 6A) suggests that some of the larger juvenile and adult fish ventured in the top layer, presumably to feed on larval fish and zooplankton.This would also explain the occurrence of adult fish in the shallow coastal waters, as seen from the fyke and gill net samples.The design of the survey did not allow sampling at night, but we can speculate that rainbow smelt larvae remained in the top layer and that some larger fish performed diel vertical migrations and ventured in the top layer at night.
Secondly, an abundance of palatable zooplankton at the pycnocline could increase prey encounters for first feeding larvae and cause the larvae to aggregate at the surface.After hatching in fresh or brackish areas, rainbow smelt larvae are passively transported towards the surface and start feeding on Mysis, Daphnids and copepod nauplii just after the resorption of their yolk sack at 7 days old (Sirois et al., 1998;Stritzel Thomson et al., 2011).These zooplankton were abundant in the LSSL and pycnocline in summers 2018 and 2019 (Maxime Geoffroy, Memorial University, unpublished data).Zooplankton generally aggregate at the pycnocline because of the consistent density differences of the water masses (e.g.Meerhoff et al., 2013;Geoffroy et al., 2017;Pérez-Santos et al., 2018).Additionally, zooplankton are generally abundant in estuarine turbidity maximum (ETM) (Blaber and Blaber, 1980;Fuji et al., 2010) and our results demonstrate higher turbidity above the pycnocline.By foraging at the pycnocline rainbow smelt larvae thus increase their feeding success and, ultimately, larval survival and recruitment (Leggett and Deblois, 1994;Pepin et al., 2015).These observations emphasize the importance of the freshwater surface layer on the distribution of rainbow smelt larvae.If the seasonal freshwater inflow changes in the Lake Melville estuary, it may impact the distribution and survival of the zooplankton prey, for example by modifying the estuarine circulation.In turn, it could potentially result in a mismatch between first feeding larvae and the occurrence and distribution of their main prey (Cushing, 1990).
In addition to predator-prey interactions, the freshwater refuge hypothesis could partly explain the aggregation of rainbow smelt above and within the pycnocline.Developed for Arctic cod, this hypothesis stipulates that for small pelagic fish hatching in cold waters, such as rainbow smelt in Lake Melville, areas with high freshwater discharge provide a warm thermal refuge for survival of early life stages (Bouchard and Fortier, 2011).A similar situation could be at play in Lake Melville where the larvae would have better growth and survival rates in the warmer waters prevailing near the surface from June to October relative to the colder bottom waters below the pycnocline (Fig. 2).These larvae would reach a larger size at the onset of the winter, when ice starts forming, which would result in higher recruitment and selectivity for these individuals (Bouchard and Fortier, 2008).

Fish diversity and distribution
Although we cannot determine relative abundance based on the number of DNA sequences, by combining the eDNA results and net catches we can make inferences on presence-absence and diversity of all life stages of fish (Lodge et al., 2012;Thomsen et al., 2012;Evans et al., 2017).As expected, the highest number of sequences for freshwater fish, such as longnose sucker (Catostomus catostomus), white sucker (Catostomus commersoni), and longnose dace (Rhinichthys cataractae), were found above the pycnocline and the highest number of sequences for marine species, such as American plaice (Hippoglossoides platessoides), Atlantic cod (Gadus morhua), and Greenland cod (Gadus ogac), were detected below the pycnocline (Table 3).Yet, the DNA of both freshwater and marine species were detected throughout the whole water column, likely due to fish migrations and to mixing and drifting of biological tissues via currents (Deiner and Altermatt, 2014).
Estuaries are complicated water interfaces, and one can expect that eDNA vaguely distinguishes differences in communities.However, García-Machado et al. (2022) found clear changes in community composition between fluvial, brackish, and marine sectors of the St. Lawrence estuary, and Hallam et al. (2021) similarly distinguished between fluvial and estuarine communities of the Thames River.In addition, Berger et al. (2020) demonstrated a strong bi-dimensional pattern of distribution of fish communities associated with different water mass characteristics and Littlefair et al. (2021) demonstrated that fish behavior and thermal stratification determine the eDNA distribution in the water column.Our methodology to sample and process eDNA samples, including field and laboratory controls, multiple PCR, and bioinformatic filtering, furthermore diminishes the impacts of contamination from exogenous sources.Yet, under certain circumstances, eDNA transported by currents can confound the detection of community differences and the distribution of species.To our knowledge, no method of detection is immune to catchability biases, and combining several approaches decreases the risk of false negatives and positives.
The eDNA results confirmed the presence of 48 species that have previously been documented within Lake Melville and surrounding tributaries (Backus, 1957;Wells et al., 2017;McCarthy and Gosse, 2018).Of these 48 species, 11 have rarely been observed in the estuary: American sand eel; burbot; herring; fish doctor; snakeblenny; rock gunnel; Atlantic halibut; Allegheny pearl dace; Atlantic redfish spp.; fourline snakeblenny; and Pacific sandlance.Common dab were regularly detected in our eDNA samples, but existing literature reports that the northern limit of the species was assumed to be in the Strait of Belle Isle (Backus, 1957).Other species detected with eDNA that were reported for the first time in Lake Melville comprised blackspotted  Table 3 Number of sequences of environmental DNA of all species above and below the pycnocline in upper Lake Melville and Goose Bay divided by habitat (blue = Marine, gray = Freshwater, green = Anadromous).
stickleback, haddock, prickly sculpin, and mountain whitefish.Although the use of eDNA is not without constraints, our results support previous studies which concluded that eDNA is a useful complementary tool for monitoring fish diversity and species richness in northern estuaries (Berger et al., 2020;Lacoursière-Roussel et al., 2018;Gibson et al., 2023).The net and eDNA data included freshwater, anadromous and marine species, and indicated that upper Lake Melville was more biologically diverse than Goose Bay, at least for adult fish.The hydroacoustic surveys also detected higher backscatter from larger fish in the bottom layer of upper Lake Melville than Goose Bay in summer.Despite potential exchange within the top layer across both regions, the larger inflow of Atlantic water from the Labrador Sea in upper Lake Melville than in Goose Bay, due to the sill at the narrows (Durkalec et al., 2016), likely explains higher abundances of marine fish species.The Atlantic water inflow could also result in higher abundance of lipid rich copepods, which in turn would attract fish (Head et al., 2003;Pepin, 2013).This hypothesis is supported by the higher acoustic backscatter at 120 kHz, a frequency at which meso-and macrozooplankton dominate the acoustic signal, in upper Lake Melville than in Goose Bay.As for the freshwater and anadromous species present in upper Lake Melville, they may be related to the freshwater inputs from the Northwest and Kenamu rivers.For example, both rivers are known by residents as prime fishing grounds for salmonids (Doug Blake, Northwest River, personal communication).
The higher abundance of large pelagic fish in Goose Bay than upper Lake Melville in winter (Fig. 8) is more puzzling.A recreational ice fishery for adult rainbow smelt occur in Goose Bay during winter and spring and it is possible that mature rainbow smelt migrate to Goose Bay to spawn at the mouths of the Churchill and Goose rivers at this time of the year.In winter, because the salinity of Goose Bay is lower the freezing point is higher, which results in higher temperatures (Fig. 2 and Supplementary Table 6).Coupled with higher turbidity and primary production (Fig. 4), higher winter temperatures in Goose Bay could provide a better habitat for zooplankton and fish (North and Houde, 2003).Further studies should investigate the importance of Goose Bay as a winter habitat for estuarine fish.

Potential impacts of anthropogenic stressors
This study complements other baseline studies of Lake Melville (Lu et al., 2013;McCarthy and Gosse, 2015;2018;Durkalec et al., 2016).Unlike previous studies, our survey coincided with the completion of the second phase of the Churchill River hydroelectric development (i.e., Muskrat Falls) on the main tributary of the estuary, and the first winter freshwater release in 2020.In this context, the observed 4-fold decrease in Goose Bay surface salinity between winters 2019 and 2020 is worth noting (Supplementary Table 6).For the same period, a little more than a two-fold salinity decrease was also observed in upper Lake Melville surface waters, although with relatively large variations.In the bottom layer, an increase in salinity was observed in upper Lake Melville in the winter 2020 compared to 2019, which is consistent with increased estuarine circulation and inland pumping of oceanic waters at depth associated with increased river discharge (Saucier et al., 2009).While it can be hazardous to draw any conclusion from such a short time period, these baseline observations emphasize the importance of closely monitoring hyrographic changes in the ecosystem, especially in the context of planned alteration of the natural debit of the Churchill River.In other regions, hydroelectric dams have been attributed to habitat fragmentation and decline in biodiversity (Wu et al., 2003;Liu et al., 2019aLiu et al., , 2019b)).In the light of the findings made here, it was decided to leave the mooring at the mouth of the Churchill River to start such monitoring, which hopefully can bring answers to some questions raised in this study.

Conclusions
This study demonstrates that 1) rainbow smelt is one of the most abundant forage fish in Lake Melville; 2) the distribution of rainbow smelt larvae is closely tied to the physical oceanographic conditions; and 3) eDNA can be used as a complementary tool to monitor fish diversity in sub-arctic estuaries.Rainbow smelt were consistently sampled in the ichthyoplankton nets, gill and fyke nets, and detected in the eDNA samples.The species is known to funnel energy between zooplankton and top predators such as trout, ringed seals, and seabirds (McCarthy andGosse, 2015, 2018).It is also an important cultural species for local communities and a yearly recreational ice fishery targets rainbow smelt.The recent completion of the Muskrat Falls hydroelectric project on the lower Churchill River, the main tributary of Lake Melville, could modify the seasonal inflow of freshwater.This process could already be underway, with an observed decrease in salinity in the winter of 2020 compared to 2019.Our results suggest that, as a counterpart, if the freshwater inflow is reduced in summer it could impact the distribution and, potentially, survival of rainbow smelt larvae.Given the importance of the species for higher trophic levels, any change in the abundance and ecology of rainbow smelt will have cascading impacts on the whole ecosystem of Lake Melville.We thus suggest monitoring the abundance, distribution, and condition of rainbow smelt as a sentinel species for the ecosystem of Lake Melville.Monitoring fish biodiversity, for instance using eDNA, would also indicate potential changes in the ecosystem of the estuary.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Maxime Geoffroy reports financial support was provided by The Marine Environment Observation, Prediction and Response Network (MEO-PAR).Frederic Cyr reports financial support was provided by Fisheries and Oceans Canada.Maxime Geoffroy reports financial support was provided by Natural Sciences and Engineering Research Council of Canada.James McCarthy reports financial support was provided by WSP E&I Canada Ltd. James McCarthy reports a relationship with Nalcor Energy Ltd that includes: consulting or advisory.

Fig. 1 .
Fig. 1.Bathymetry of Goose Bay and upper Lake Melville within Newfoundland and Labrador (insert).The continuous gray line indicates the hydroacoustic transects followed in summers 2018 and 2019.Environmental stations (dots) were samples in summer and winter.The locations of fyke and gill net sampling are denoted by triangles (2018) and crosses (2019).The star indicates the position of the long-term oceanographic mooring.

Fig. 3 .
Fig. 3. Turbidity (in Formazine Turbidity Unit; FTU) and chl-a (in part per billion; ppb) section plot in the top 20 m for Goose Bay (GB) and upper Lake Melville (ULM) during winter 2019.For each panel, the density anomaly referenced to surface (σ 0 ; in kg m − 3 ) is plotted as solid light gray lines identified in the chl-a panel.The vertical dashed lines indicate the location of the stations.The sill between GB and ULM is identified in brown.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5 .
Fig. 5. Examples of raw S V echograms at 38 kHz at station 2, where the bottom depth was 32 m, during A) summer 2018 as measured with a hull mounted Simrad EK60; and B) winter 2019 as measured with a Simrad Wideband Autonomous Transceiver (WBAT).The dashed red lines indicate the sound scattering layer from larval fish, which in this case was located between 9 and 15 m.Corresponding temperature (red) and salinity (black) profiles for station 2 in C) summer 2018 and D) winter 2019.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 6 .
Fig. 6.Percentage of total Target Strengths (TS) at 38 kHz for A) pooled summer data (2018-2019; n = 122,924); and B) pooled winter data (2019-2020; n = 18,175) at depth ≤25 m (orange) and >25 m (gray).The dashed lines indicate the modes.Average depth of each target and their corresponding TS for targets detected in C) summer and D) winter.The maximum range of the WBAT and portable echosounder in winter was set to only 100 m. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 8 .
Fig. 8. Nautical area scattering strength (S A in dB re 1 m 2 nmi − 2 ) at 38 kHz above and below 25 m.The dashed line delineates Goose Bay (GB) and upper Lake Melville (ULM).

Table 1
Mean abundance (CPUE) in fyke and gill nets with corresponding mean biomass (CPUE) and average length for each species.

Table 2
Shannon-Weiner diversity Index (H) by region and year calculated from catches in gill and fyke nets deployed in shallow water (<2m).