Groundtruthing of pelagic forage fish detected by hydroacoustics in a whale feeding area using environmental DNA

Conservation of whales, considered as umbrella species in marine environments, re - quires to be able to understand their relationships with ecosystem components such as prey species, including pelagic fish. However, studying such relationships in nature is a technical challenge. In this study, we used two noninvasive methods in combi - nation, namely hydroacoustics and environmental DNA (eDNA), to detect five pe - lagic or semipelagic fish species in the Saguenay–St. Lawrence Marine Park (Québec, Canada): the sandlance Ammodytes sp., the Atlantic herring Clupea harengus , the capelin Mallotus villosus , the rainbow smelt Osmerus mordax, and the redfish Sebastes sp. The Marine Park is a major summer


| INTRODUC TI ON
Conservation of marine biodiversity requires obtaining an accurate description of the trophic ecology between marine predators and their prey (Terborgh & Estes, 2013). Since the 70s, the global will to protect marine mammals, namely whale populations depleted by hunting, has specifically led to regulatory measures in whale catching for commercial purposes (Birnie, 1989). Although these measures favored recovery in some whale populations, others remain critically endangered (Attard, Beheregaray, & Möller, 2016;Clapham, Young, & Brownell, 1999). The increase in contemporary threats is thought to impede the recovery of depleted whale populations such as increased contamination in coastal areas, maritime traffic, overfishing of prey sought out by whale species, and the impacts of climate change. Conservation of whales, considered as umbrella species, has been a strong motivation for wider marine ecosystem conservation objectives. In this context, it has been proposed that regulatory measures must take into account the unique characteristics of each marine mammal population, including their relationships with ecosystem components such as prey species like other fish. However, tracking such relationships in nature remains a technical challenge (Birnie, 1989).
The Saguenay-St. Lawrence Marine Park (Québec, Canada) and surrounding waters are recognized as a major summer feeding ground for a wide diversity of marine wildlife species. Marine mammals converge to the head of the Laurentian Channel to feed from spring to fall, indicating that particular bathymetric features and oceanographic processes (e.g., tidal upwellings) have an influence on pelagic prey availability such as pelagic fish and euphausiids (Ménard, 1998;Simard, 2009). Schooling pelagic fish are major components of marine trophic webs that are consumed by higher trophic levels. Among the common wildlife species in the St. Lawrence Estuary (SLE), several species of whales, seals, and seabirds are pelagic fish predators. Changes in the availability of pelagic fish have been suggested as a possible explanation for the observed decline in fin whale abundance in their traditional feeding areas at the head of the Laurentian Channel observed since the early 2000s (Martins, Turgeon, Michaud, & Ménard, 2018). Since the same period, shifts in environmental regimes at the scale of the Gulf of the St. Lawrence Estuary, namely demersal and pelagic fish availability and composition, have also been proposed as a contributing factor in the decline of the St. Lawrence beluga whale population (Delphinapterus leucas) (Plourde et al., 2014). In this context, the need to understand the influence of pelagic prey availability in sustaining marine wildlife biodiversity in the SLE catalyzed the development of studies in the Saguenay-St. Lawrence Marine Park (Ménard, 2009).
Although significant work has been done in the SLE on zooplankton communities (Gagné et al., 2013), scarce research efforts have been dedicated to pelagic forage fish. The community of schooling pelagic forage fish is relatively low in diversity in the SLE, being essentially composed of five species: the sandlance Ammodytes sp., the Atlantic herring Clupea harengus, the capelin Mallotus villosus, the rainbow smelt Osmerus mordax, and the Atlantic mackerel Scomber scombrus.
These five fish species are considered as prey for cetaceans, although the cold upwelled waters at the head of the Laurentian Channel are considered inadequate habitat for both the rainbow smelt, found in surrounding warmer areas, and the Atlantic mackerel, found mainly in the lower SLE and Gulf of St. Lawrence. The redfish Sebastes sp. is a semipelagic species which has undergone a spectacular rebound of its population, with the presence of juvenile redfish abundance from 2011 to 2013 cohorts that demonstrate record biomass (Senay et al., 2019).
Redfish is not generally considered as an important prey species for cetaceans, although juveniles are pelagic and can form schools that could be consumed by whales.
In 2009, initiated largely for conservation purposes because of the observed decline in the number of fin whales and to better understand the ecological needs of the endangered St. Lawrence beluga whales, a monitoring of the spatial distribution and abundance of potential prey of marine wildlife predators, namely euphausiids and pelagic fish, was launched by Parks Canada in the Saguenay-St. Lawrence Marine Park.
This monitoring program relied on hydroacoustic technology and systematic simultaneous visual surveys of seabirds and marine mammals.
The objective was to better understand how prey spatiotemporal variability influenced marine wildlife predators in terms of abundance and diversity and to locate and characterize the environmental conditions of feeding areas for marine mammals.
Fisheries acoustic methods have been widely used for the monitoring of fish and zooplankton, and to estimate their abundance and distribution (Holliday, Pieper, & Kleppel, 1989;Masse, 1996). They offer the advantages of continuously sampling at high-resolution large volumes of water and being a noninvasive method. However, one of the main challenges and limitations in fisheries acoustics is species identification (Horne, 2000). It is now well documented that acoustic backscatter strength from fish and zooplankton is frequency-dependent. Differencing the mean volume backscattering strength (ΔMVBS) has been used as a multifrequency classification method to discriminate backscatter from fish to zooplankton (Madureira, Everson, & Murphy, 1993), between euphausiid species (McQuinn, Dion, & St. Pierre, 2013), and between fish species (Korneliussen & Ona, 2003). Since the swim bladder contributes to at least 90% of the sound scattered, fish species can be separated in two groups depending on the presence or absence of a swim bladder (Foote, 1980). With the exception of sandlance, a pelagic fish which is also characterized by its particular burrowing behavior (Robards, 1999), all the other species which form the pelagic fish communities at the head of the Laurentian Channel arbor an internal swim bladder. Previous studies documented that fish with swim bladder have a weakly decreasing volume backscattering coefficient (S v ) with increasing frequency, while fish without swim bladder have a S v essentially frequency independent at 18-120 kHz and peaked at 200 kHz (Korneliussen & Ona, 2002). However, species discrimination by multifrequency classification methods requires to be groundtruthed. This is traditionally performed by conventional capture methods such as trawling, which has the advantage of providing information on organism size but also the disadvantage of being invasive. Moreover, trawls are selective (Wileman, 1996), particularly for size (Bethke, Arrhenius, Cardinale, & Håkansson, 1999), and many fish can escape (Engås & Godø, 1989) or avoid nets altogether (Fonseca, Sanches Fernandes, Fontainhas-Fernandes, Monteiro, & Pacheco, 2016;Kaartvedt, Staby, & Aksnes, 2012). Due to local environmental conditions characterized by strong currents and a rugged seafloor, trawling has proven to be ineffective in the area studied here. Another method to obtain additional evidence for species identification, organism size and tilt angles, is imagery. Because of strong currents and dark and turbid waters, it might remain difficult or impossible to groundtruth acoustic signals using video imagery techniques. Thus, after numerous attempts to groundtruth acoustic signals by various fishing methods and imagery, there was a need to find a method that could effectively be used in combination with hydroacoustics to groundtruth the acoustic backscatter and improve the interpretation of acoustic data.
DNA is released by organisms in the environment from skin cells, mucus, metabolic waste, or gonads. In other words, eDNA is a signature of the organism's presence in the environment. Therefore, eDNA studies are used to detect DNA from organisms without the need to capture them by collecting samples from their environment (Lodge et al., 2012). eDNA allows to estimate fish diversity in aquatic environments (Li et al., 2018), as well as to predict the biomass of organisms (Coulter et al., 2019;Lacoursière-Roussel, Côté, Leclerc, & Bernatchez, 2016;Lacoursière-Roussel, Rosabal, & Bernatchez, 2016).
The objective of this study was to evaluate the efficiency of eDNA to identify the species composition of pelagic forage fish shoals detected by hydroacoustics in the Saguenay-St. Lawrence Marine Park. eDNA and hydroacoustics have been recently used in combination to assess the distribution and the abundance of the Japanese jack mackerel Trachurus japonicus (Fukaya et al., 2018;Yamamoto et al., 2016), and they have been described as an innovative approach to comprehensively detect aquatic species (Miksis-Olds & Watts, 2019). But to our knowledge, their ability to noninvasively detect several fish species simultaneously has never been tested so far.
Here, we used eDNA and hydroacoustics to detect five pelagic or semipelagic fish species (the sandlance Ammodytes sp., the Atlantic herring C. harengus, the capelin M. villosus, the rainbow smelt O. mordax, and the redfish Sebastes sp.). Our study, which was conducted over a 2-year period in different areas of the marine park, highlights the spatiotemporal variability in food resources that marine mammals encounter in the SLE.

| Collection of acoustic data and water samples in fish shoals
A total of 20 shoals were sampled by acoustics and for environmental DNA, including 12 shoals in summer 2017 (Figure 1, orange circle; sampled between 18 July and 13 October 2017) and eight shoals in summer 2018 (Figure 1, red triangle; sampled between 17 July and 20 September 2018). The fish shoals were detected by hydroacoustics in areas known for fish aggregation and in different areas to insure variety in species composition. Information about each shoal (date of sampling, GPS coordinates, mean temperature (°C), mean salinity (PPM), and mean depth (m) measured at each location) is reported in Table S1.
Acoustic data were collected continuously and simultaneously to the water sampling from the Park's Canada vessel L'Alliance using a Simrad EK60 multifrequency acoustic system operating at 38, 120, and 200 kHz. Calibrations of the system were performed every year with the standard-target method (Foote, 1990;Simmonds & MacLennan, 2005), using a 38.1 mm tungsten-carbide (6%) sphere suspended from motorized tripoint lines. The transducers were installed at 1 m below the surface on the port side of the vessel. The acoustic system configurations were adjusted to sample equivalent volumes on all transducers (Korneliussen, Diner, Ona, Berger, & Fernandes, 2008). All frequencies were transmitted simultaneously once per second with a 1.0 ms pulse length. Raw volume backscattering strength (S v dB re 1 m −1 ) data were collected in Simrad format using the ER60 software. Additionally, a conductivity-temperaturedepth (CTD) profiler was deployed in all selected shoals to record environmental parameters ( Figure 2). Only the first CTD cast within each shoal was selected since it was the only one to respect the soaking time.
For environmental DNA, water samples were collected on the same boat with a sterilized Niskin bottle (1 L) deployed along the CTD cable at the desired depth of sampling in fish shoals ( Figure 2). No fish are being collected or processed on this boat or elsewhere by the members of the crew. Negative field controls (blanks) using distilled water were done for every shoal. These controls were performed on the boat in the field and treated in the same way as the real samples to estimate the possible contamination of the filtration material. Within a given shoal, water samples were taken at four locations (T1, T2, T3, and T4) at depths corresponding to the core of the shoal as detected by hydroacoustics, as well as at one location upstream of the shoal by taking into account the general direction of surface currents. The upstream water sample was taken outside of the shoal of interest (minimum distance of 323.71 m, mean distance of 1,275.05 m between the upstream location and the shoal of interest), in an area exhibiting highly similar physico-chemical properties (in terms of temperature, salinity, depth) as the core of the shoal (Table S1). Because physico-chemical properties were similar between the upstream location and the core of the shoal, we expected the eDNA concentrations not to be affected by these environmental factors. The

| Acoustic data processing
Acoustic data were processed using Echoview Software version 6.1 (Myriax Software) and standard cleaning techniques (Simmonds & MacLennan, 2005). The seafloor was automatically detected, and an offset of 0.5 m was applied. Then, the bottom line was manually edited to exclude all false bottom detection. Data shallower than 5 m, which is twice the transducer near-field, were also removed from analysis to eliminate the transmit pulse and minimize backscatter from surface bubbles (Simmonds & MacLennan, 2005). In order to minimize the effects of background noise, we estimated the signalto-noise ratio (SNR) at each frequency using a previously developed method (De Robertis & Higginbottom, 2007) and excluded all data with SNR <10 dB for either frequency. In addition, a background threshold of −80 dB was applied. Finally, the echograms were manually edited to remove logging artifacts (e.g., surface bubbles, side lobes). For each fish shoal sampled, a postcalibration was applied to the acoustic data in order to modify the sound speed according to the CTD profile. One shoal (ID: 2017_02) was completely excluded from the acoustics analysis since the most of the backscatter was included in the side lobes. isms. Finally, density function was applied using the library ks to define 70% volume contour with kernel density estimate in order to compare fish shoals acoustic multifrequency response in the same two-dimensional space (Duong, 2014). The Kernel density estimate is a common statistical method for density estimation (Silverman, 1986). The ΔS v,120-38 threshold previously described (McQuinn et al., 2013) of 5 dB for krill and fish discrimination was used for data visualization.

| eDNA extraction
eDNA extraction was performed using a QIAshredder and DNeasy Blood and Tissue Kit (Qiagen) according to a previously developed protocol (Goldberg, Pilliod, Arkle, & Waits, 2011). The extracted eDNA was stored at −20°C until amplification by quantitative realtime PCR (qPCR) method. For each extraction batch, an extraction negative laboratory control was added to account for possible contamination.

| Design of primers and probes specific to the five fish species and optimization
In order to identify the five fish species of interest in our study (Ammodytes sp., C. harengus, M. villosus, O. mordax, and Sebastes sp.), we used DNA sequences from the mitochondrial gene cytochrome c oxidase subunit 1 (COI). COI is a gene that has been very successfully used to discriminate between closely related fish species based on sequence divergence (e.g., April, Mayden, Hanner, & Bernatchez, 2011;McCusker, Denti, Van Guelpen, Kenchington, & Bentzen, 2013).
COI sequences from the five species were obtained from the bold system (Barcode of Life Database http://www.bolds ystems. org/index.php/). The COI sequences from related species present in F I G U R E 2 Measuring hydroacoustics and eDNA in the Saguenay-St. Lawrence Marine Park. Left: The conductivity-temperature-depth (CTD) profiler used in the study, with a sterilized Niskin bottle (1 L) deployed along the CTD cable. Right: 38 kHz echogram showing the water sampling into a fish shoal, the line to the left indicating the cable descending into the shoal and the right line while it is ascending. The color scale illustrates the variation in volume backscattering coefficient (S v ) that ultimately allows to discriminate shoals, their composition, and their density the same environment were also used to design the specific primers (Table S2). Each primer pair and probe was developed to maximize the number of mismatches between the targeted species and the related species. The sequences were aligned using the software Geneious (https ://www.genei ous.com/). From aligned COI gene sequences, specific primers and probes were successfully designed also using Geneious and verified using the software Primer Express 3.0 (Life Technologies). In addition, the Primer-BLAST tool (https :// www.ncbi.nlm.nih.gov/tools/ primer-blast/ ) was used to verify the specificity of the amplification on other species that could be present in the same environment. The sequences of the designed primers and probes are presented in Table 1. Note that for the primers designed for O. mordax, there was a polymorphism on the 3′ end of the reverse primer for the subspecies O. m. dentex, which would likely prevent the use of these primers to target this specific subspecies.

| Test of primer and probe specificity by qPCR method
Specific primers and probes were tested on DNA tissues of targeted and related species (list of tested species are available in Table S2).
Tissues of one sandlance, one herring, and one redfish were obtained from Parks Canada. Tissues of capelin and rainbow smelt and other related species were available in our laboratory at the Institut de Biologie Intégrative et de Systèmes (IBIS), Université Laval, Québec. No tissue for related species of sandlance was available. Genomic DNA was extracted using a modified version of a salt extraction protocol (Aljanabi & Martinez, 1997). In order to identify the species of fish of interest, the COI was amplified by PCR using the FishCOIF primers (forward, F and reverse, R, Table S3) (Lara et al., 2010)  were aligned and compared to the BOLD fish database using the software Geneious and the BOLD system (http://www.bolds ystems.org).
The qPCR method was used to test specificity and efficiency of each primer pair on DNA of targeted species and related species (Table S2) Specificity and efficiency of all primer pairs were first tested with FAST SYBR Green method, and once amplification of targeted species only was achieved, the TaqMan assay was performed with the species-specific probes for further specificity enhancement.

| Assay sensitivity
A standard curve experiment was performed following the reaction and condition protocol using TaqMan Environmental Master Mix 2.0 as described above. A synthetic DNA template of 500 base pairs including the target amplicon sequence was designed from the COI gene sequence and prepared by Integrated DNA Technologies Inc.
From the stock at 1.00E + 10 copies/µl, an eight-level dilution series (2,000, 1,000, 500, 100, 20, 8, 4, 2, and 1 copies per reaction) was prepared in sterile yeast tRNA at 10 µg/µl. Ten replicates of each dilution series were run to determine the amplification efficiency and the limit of detection defined as the lowest copies per reaction with >95% amplification success (Bustin et al., 2009). The threshold of one molecule was detected at 40 C t for each primer pair. Therefore, all signals exceeding 40 C t were considered as potential artifact of PCR and were eliminated from further analyses.

| eDNA samples analysis by qPCR method
To analyze the collected eDNA samples, the TaqMan qPCR method was used with the addition of the SPUD to the reaction, as well as a standard curve. SPUD is a universal system for rapid quality control of nucleic acid templates before qPCR (Nolan, Hands, Ogunkolade, & Bustin, 2006). DNA presence of each targeted species was tested on six technical replicates for each eDNA sample TA B L E 1 Species-specific qPCR primers used for the amplification of different regions of the mitochondrial COI gene designed to identify targeted shoal species The amplification was performed in a final volume reaction of 20 μl including 1.8 μl of each primer (10 μM), 0.5 μl of probe (10 μM), 10 μl of Environmental Master Mix 2.0 (Life Technologies), 3.9 μl of SPUD, and 2 μl of DNA following these conditions: 2 min at 50°C, 10 min at 95°C, 50 cycles of 15 s at 95°C, and 60 s at 60°C.
All qPCR results were quantified using a standard curve of known DNA quantities which allowed us to quantify positive PCR amplification in number of molecules. The average of the number of molecules measured in the six replicates of a given sample was calculated. For each shoal, except for the five shoals for which no upstream location could be sampled, the average number of molecules at one location sampled within the shoal was corrected by subtracting the value obtained at the upstream location. Thus, the average number of molecules obtained at the upstream location was always scaled back to zero. If the corrected average of number of molecules measured for a species was zero or below zero in a given sample, then we concluded that the DNA of this species was not present in that sample. Alternatively, detection of a species was inferred when the corrected average of number of molecules measured for a species at a specific location was above zero (Table S1).

| Statistical analyses for environmental DNA
Statistical analyses were performed using the R software version 3.5.2 (R Core Team, 2017). First, we used the package lme4 (Bates, Mächler, Bolker, & Walker, 2015) to determine whether the location within a shoal (upstream, T1, T2, T3, or T4) or the location of a shoal compared to other shoals impacted eDNA concentration.
Distinct analyses were performed for the shoals with or without the upstream location sampled. We used a linear mixed model including the depth, the temperature, and the salinity of the shoals as fixed effects and the IDs of the shoals as a random effect. We checked for homogeneity of variances using a Levene test and also tested normality of residuals using graphical assessment (Q-Q plot and normal distribution histogram) and a Shapiro-Wilk test. To obtain p-values, we performed pairwise comparisons using locations within/between shoals with the package lmerTest (Kuznetsova, Brockhoff, & Christensen, 2017). Then, we used a nonparametric Mann-Whitney test to compare eDNA concentrations of a specific species between shoals and years.

| Comparison of hydroacoustics and eDNA data
The comparison between both datasets was examined using prin-

| Hydroacoustic approach
The density analysis of the acoustic data (EIUs, 25 × 2 m) in a two-dimensional space (ΔS v,120-38 , ΔS v,200-120 ) revealed significant differences in the acoustic responses of the different shoals ( Figure 3; Figure S1; Table S1) with different CTD profiles ( Figure   S2).The 70% volume contours for each sample illustrate two main  Results of the specific amplification by qPCR are presented in Figure S3. Specific primers amplified the desired species, and no amplification was detected for the related species. For the primers specific to M. villosus, we observed an amplification with a DNA mix containing three species, but this amplification had a very low efficiency and was under the threshold of detection.

| Specificity of primers and probes
Based on the standard curve experiment, the assays for the targeted species had an amplification efficiency between 95.5% and 114.5% and a limit of detection between 20 mtDNA copies/rxn (O. mordax) and 4 mtDNA copies/rxn (A. hexapterus) depending on the targeted species (Table 2).

Detection of targeted species in shoals
An environmental DNA-qPCR-based approach was applied to infer the presence of five target fish species in 20 shoals sampled in the Saguenay-St. Lawrence Marine Park. Typically, two or more targeted species were identified in the eDNA samples from all shoals (Table 3).
Only one shoal (ID: 2017_04) had a positive amplification for one species only, namely the capelin M. villosus (Table 3 panel A).
Furthermore, two species were detected in seven shoals comprising different assemblages of the following three species: capelin, redfish, and herring (Table 3 panel B). In total, 10 shoals showed a positive amplification for at least three targeted species, including sandlance, herring, capelin, or redfish (Table 3 panels C and D). These four species were all identified in three shoals sampled in 2017 (Table 3

Capelin eDNA detection
The presence of capelin was detected in all shoals but three: only three shoals did not have a positive amplification for the capelin. replicates. Therefore, we decided to tentatively consider these shoals as positive for the presence of capelin (despite negative numbers of eDNA molecules after correction with the upstream location) (Table 3 panel B).

Redfish eDNA detection
In 2017, nine out of the 12 shoals that were sampled included  (Table 3 panels B and C).

Sandlance eDNA detection
As reported for the redfish, the identification of the sandlance  (Table 3 panels B-D).

Absence of the rainbow smelt eDNA detection
The presence of the rainbow smelt was not found in any of the shoals of 2018, the only year this species was tested for eDNA.

Variations in the level of eDNA detection
DNA quantity and the number of positive detections tended to be variable depending on shoals (Table 3 and Table S4), on the position within a shoal (upstream, T1, T2, T3, or T4) (Figure 4), and on the targeted species (Table 3, Table S4 and Figure 4). However, using a linear mixed model, we did not find any significant effect of the location within a shoal or of the location of a specific shoal on eDNA concentration. Among the 60 samples coming from 15 shoals for which the presence of all species but rainbow smelt was estimated, 38 of them showed the lowest DNA quantity at the upstream location compared to within the shoal (locations T1, T2, T3, or T4).

Negative controls and positive amplification sequencing
Positive amplifications were detected in thirteen-field negative controls (in six controls from shoals sampled in 2017 and in seven controls from shoals sampled in 2018). However, quantities were much lower than actual shoal samples; therefore, the results can be taken with confidence. All extraction negative controls showed no positive amplification indicating the absence of contamination during sample extraction.

| Combining hydroacoustic and environmental DNA
Principal component analysis and hierarchical clustering on principal component (HCPC) analysis allowed comparing clusters obtained with both datasets (Figure 5 and Figure S4). The variables used were the three frequency differences (ΔS v,200-120 , ΔS v,200-38, and ΔS v,120-38 ) and Cluster #1 (yellow) appeared to be more associated with fish with a swim bladder as negative values were computed for 200 and 120 kHz frequencies with acoustics dataset, and a positive v-test value was TA B L E 2 Percentage of amplification efficiency (eff%), limit of detection, intercept (y-inter), and r 2 results for each standard curve developed with a synthetic DNA template detected for redfish variable with eDNA data. In the opposite, cluster #2 (green) was more likely associated with fish without swim bladder as samples with a positive response to ΔS v,200-120 and a negative response to redfish were detected in the group. In cluster #3 (gray), all samples overlapped with the 5 dB ΔS v,120-38 threshold ( Figure 3) and almost all of them were located near the mouth of the Saguenay Fjord ( Figure 1). Clusters #4 (red) and #5 (magenta) computed with eDNA dataset contained only one sample (ID: 2018_07 and 2017_11, respectively). Cluster #4 was attributable to a negative response to capelin after applying the correction for the mean eDNA concentrations.
Cluster #5 was associated with a strong positive response for both sandlance and herring. The similarity between clusters (#1, #2, and #3) for both methods ranged from 50% to 75% of the samples. tect Ammodytes sp. in the water column (Sigsgaard et al., 2017), it is important to keep in mind that the sandlance alternates between swimming pelagically (for feeding) and lying buried in the sediments (Behrens & Steffensen, 2007;Robards, 1999). This particular behavior could have restricted the detection of eDNA molecules for the sandlance as our study was based on the analysis of samples from the water column. But perhaps more important is the possibility that more than one species of the genus Ammodytes coexist in the system and that the primers we developed based on the A. hexapterus specimen that was collected in the studied system cannot efficiently amplify DNA from A. americanus, which could be the predominant Ammodytes species in this system. Clarifying this issue will have to await the collection of more specimens of both species in order to thoroughly compare their sequences and possible redesign new primers. Finally, while herring was detected in 2017 and 2018, eDNA from rainbow smelt was not found in any of the samples. It may be explained by the oceanographic processes (e.g., tidal upwellings) near the head of the Laurentian Channel that prevent the transport between the rainbow smelt populations of the middle and lower estuary (Bernatchez & Martin, 1996;Ingram & El-Sabh, 1990). In order F I G U R E 5 Comparison between (a) differences in mean volume backscattering strength (200, 120, and 38 kHz) and (b) eDNA concentrations combined with environmental parameters (depth, salinity, and temperature) of the sampling sites. Hierarchical clustering on principal component (HCPC) was performed for each dataset, allowed in grouping samples into distinct clusters (3 and 5). V-test was computed to determine the significant variables (p < .05) for each cluster, and the v-test values were associated with the overall mean of the variable. The hierarchical trees were cut at suggested level (dashed black line). The symbol (¤) was added to samples when eDNA concentrations were not corrected by the eDNA value obtained at the upstream location. One shoal (ID: 2017_02) was completely excluded from the analysis to obtain a better estimate of the rainbow smelt distribution, future studies may concentrate sampling efforts on the north and south shores of the SLE, where more suitable habitat for this species is found.

| D ISCUSS I ON
The use of eDNA methods allowed to groundtruth acoustic data, and their combined integration led to specifically define five clusters in the Saguenay-St. Lawrence Marine Park. The first cluster was associated with fish with a swim bladder and was primarily located on the north shore of the marine park. The second cluster was more associated with swim bladderless fish and was generally found in the center of the SLE. The third cluster composition was mixed and was detected near the mouth of the Saguenay. This zone is known for being highly dynamic with localized upwelling currents and tidal mixing (Marchand, Simard, & Gratton, 1999 shoal was the southernmost shoal located in the middle estuary. Interestingly, the highest levels of eDNA molecules for sandlance and herring were measured in this shoal (respectively, 20/24 and 24/24 eDNA replicates were positive for these species, with a corrected number of molecules detected of 0.6 and 16.6). The middle estuary therefore appears as a promising place to perform sampling and try to increase the detection levels of the fish species, especially for the sandlance.
Combining hydroacoustics and eDNA helped to obtain a picture of the dynamics of pelagic forage fish species in the Saguenay-St.
Lawrence Marine Park over a 2-year period. However, the limitations of these methods, as well as the complexity of this marine ecosystem, have to be taken into account when interpreting the collected data. First, hydroacoustics allowed to discriminate between shoals including fish with swim bladder or fish without swim bladder, mixed shoals, and shoals including euphausiids. Nonetheless, this nonselective method remains inadequate for species identification and for collecting samples near the surface and the bottom of the water column (Scalabrin, Marfia, & Boucher, 2009;Thorne, 1983). Second, the level of eDNA molecules detected by qPCR was used to infer and compare the relative abundance of each species in the shoals of the marine park, as it was previously done in other aquatic systems (Pilliod, Goldberg, Arkle, & Waits, 2013 In these cases, eDNA detection represents probably traces of the species carried at the surface by strong currents rather than direct presence in the shoal.
To get a better estimate of fish abundance, we measured for each shoal, an upstream location for which no or very low levels of species detection were expected, and which could therefore be used as an external reference. Indeed, it has previously been reported that spatial location is an important factor when interpreting eDNA data (Evans et al., 2017). However, for five shoals, we were not able to collect upstream samples to correct the eDNA values measured. Correction with the upstream location primarily allowed obtaining more stringent results: In more than 80% of the samples, the correction with the upstream location simply reduced the levels of eDNA molecules measured for the species detected in a specific shoal, without influencing their presence/absence. In this context, the eDNA values measured for the five uncorrected shoals are less stringent, but may nevertheless represent a good indicator of species presence. Furthermore, for one-third of the samples, the eDNA values measured upstream were higher than in the actual shoal samples, suggesting that the upstream locations were not sampled far enough from the shoals to prevent species detection. In particular, for two shoals sampled in 2018 (IDs: 2018_02 and 2018_07), the concentrations obtained upstream were very high suggesting that fish were still present at high density in this external location. This especially led to underestimate the presence of the capelin inside these two shoals. Interpreting eDNA data in terms of fish presence and abundance therefore requires to choose appropriate external, negative controls, as well as appropriate locations for the sampling sites. Future studies in the SLE could, for example, select control sites more distant to the shoals of interest (at least 1 km from the core of the shoals) and sample water at different depths within the shoals. We chose to sample 1 L of water at each position within each shoal. The choice of this volume was appropriate in regard to other studies that generally use 1 L to maximize the recovery of eDNA (Hinlo, Gleeson, Lintermans, & Furlan, 2017). However, sampling a more important volume could be an efficient way to increase species detection.
For example, it was suggested in river systems that five 1 L water samples may be required to enable the detection of multiple species (Shaw et al., 2016). Also, hydroacoustic data could be used to estimate total fish biomass (Boswell, Wilson, & Wilson, 2007) in relation to eDNA concentrations. Finally, metabarcoding performed on eDNA could offer promising opportunities to simultaneously identify all species present in the shoals of interest . This could help refining the estimate of species proportions with or without swim bladder in each shoal, as well as determining whether other pelagic fish species than the ones targeted in this study could contribute to the diversity in the SLE.
To our knowledge, our study is the first to combine hydroacoustics and eDNA to document the occurrence of different whale forage species and that suggested that prey composition may vary annually in this highly dynamic marine environment that represents a major feeding ground for marine mammals. Our study was conducted during a 2-year period, and sampling efforts should be extended over years in order to get a complete picture of the situation in the area. Climatic changes, which already appear to contribute to habitat degradation of St. Lawrence belugas (Plourde et al., 2014), are most likely to exacerbate these fluctuations in fish resources, as it was already described in other systems of the Atlantic (Godo, 2003).
Furthermore, our study highlights the necessity to sample appropriate external controls in eDNA studies to get a precise estimation of fish abundance. We hope that our study will ultimately facilitate the implementation of conservation and fisheries management policies that will be appropriate to the ecological variability reported in the SLE. Hopefully also, it will stimulate the interest in performing similar combined studies between eDNA and hydroacoustics in other systems.

ACK N OWLED G M ENTS
The authors would like to highlight the contribution of the members of Parks Canada conservation team at the Saguenay-St.
Lawrence Marine Park for their participation in data collection.
In particular, Francis Lapointe (boat operator) and Laurence

CO N FLI C T O F I NTE R E S T
None declared.