Relationship between brook charr ( Salvelinus fontinalis ) eDNA concentration and angling data in structured wildlife areas

Accurate management of exploited fish populations is essential to ensure their long-term sustainability. The use of eDNA as a tool for providing information on population relative abundance offers much potential although few examples in a fishery context have been documented. In this study, we collected 600 water samples from 30 lakes in Québec (Canada) to document the relationship between brook charr ( Salvelinus fon-tinalis ) angling data and their lake's eDNA concentration. Model selection with angling data and environmental parameters was used to find the best predictive model for eDNA concentration. We found a strong correlation between the average fish den - sity from current and previous years (fish harvested/ha, adj. R 2 = 0.76) with the mean eDNA concentration among lakes, supporting the growing trends in the literature. We observed very similar levels of correlation either when eDNA and angling data were from the same year or different years. We also found a pronounced inter-year difference in lakes' eDNA quantity measured in 2019 and 2020. We hypothesize that the main drivers for this difference were inter-seasonal variation including water temperature and associated variation in fish behavior. These results support the usefulness of eDNA as a quantitative tool for exploited fish populations.

While the magnitude and strength of the relationship between organism abundance (derived from traditional methods such as nets, hook and line, and electrofishing) and the concentration of eDNA in an environment exhibits substantial heterogeneity among studies (Biggs et al., 2015;Coulter et al., 2019;Horiuchi et al., 2019;Itakura et al., 2019;Pilliod et al., 2013;Schumer et al., 2019;Shelton et al., 2019;Yates et al., 2021), a generally positive trend between these two variables has been observed (Rourke et al., 2021;Yates et al., 2019).However, the extent to which the analysis of eDNA can be practically used to infer abundance in a fishery management context remains debated (Hinlo et al., 2018).Laboratory studies, where environmental variables can be controlled, have shown a stronger correlation between eDNA concentration and fish abundance than in natural ecosystems (Yates et al., 2019).There are still few empirical studies focused on exploited fish populations at a large scale in natural environments.Lacoursière-Roussel, Côté, et al. (2016) compared the CPUE (catch per unit effort) and BPUE (biomass per unit effort) from standardized gillnet surveys to the eDNA concentration of 12 natural lakes of exploited lake trout (Salvelinus namaycush).Relative indices of population abundance, expressed by CPUE and BPUE, were obtained using eight to 50 standardized gillnet sampling sites in each lake.They found that both CPUE (R 2 = 0.45) and BPUE (R 2 = 0.39) correlated with eDNA concentration.Yates et al. (2021) examined the relationship between eDNA concentration and absolute abundance data derived from mark-recapture studies on nine brook charr (Salvelinus fontinalis) populations experiencing little recreational fishing pressure.They also found a strong correlation with two indices of abundance: density (individuals/ha, adj.R 2 = 0.59) and biomass (kg/ha, adj.R 2 = 0.63) but an allometric correction substantially increased the level of correlation (R 2 = 0.78).Recently, Spear et al. (2021) demonstrated that walleye (Sander vitreus) density, estimated from the mark-recapture method from 24 exploited lakes, explained the majority of variance in eDNA observed across study lakes (R 2 = 0.62).The correlation was even stronger when incorporating the mean walleye size and accounting for variability among samples within a lake (R 2 = 0.81).Thus, eDNA could help inland fisheries management monitor fish abundance in a non-invasive way by lowering logistical and monetary costs and potentially complementing current conventional sampling methods (Evans et al., 2017;Garlapati et al., 2019;Harrison et al., 2019;Lim et al., 2016).However, further large-scale testing in a fishery context is needed to provide further confidence in the approach to managers and generalize its application.
This study focuses on documenting the relationship between brook charr eDNA concentration and angling data in structured wildlife areas in Québec, Canada.More precisely, angling data are collected annually in provincial wildlife reserves and national parks by the SÉPAQ (Société des établissements de plein air du Québec) on their exploited lakes in Québec, Canada (SÉPAQ, 2020).These statistics are primarily used to document the fishery status and temporal trends in fish populations, but also to evaluate the success of a given intervention or a new regulation (e.g., fish stocking or fishing quota), as well as to inform anglers (Arvisais, 2004;Plourde Lavoie, 2014).The leading indicators of abundance to classify lakes by the SÉPAQ are the fishing success or HPUE (Harvest per unit effort) expressed in the number of fish harvested per angler-days (fh/ad) and the fishing quality index (FQI) expressed in biomass per angler-days (kg/ad).Angler-days represent the total number of days spent by anglers on a lake.One person fishing on a lake will correspond to one angler-day for that specific day.There is no minimum time assigned to one angler-day.For example, one angler-day can correspond to either 1 hour or 4 hours of fishing on a lake per angler (Arvisais, 2004).Other indicators are also estimated, such as fish yield (total biomass harvested in a year per hectare [kg/ha]), fish density (total fish harvested in a year per hectare[fh/ha]), and fishing pressure (total anglers-days in a year per hectare [ad/ha]).Fishing success and fishing quality index are mainly used as a large-scale tool to indicate the fishery status and temporal trends of exploited populations while moderately correlating with abundance indices calculated from standardized surveys (e.g., gillnet surveys) (Jansen et al., 2013;Plourde Lavoie, 2014) However, in brook charr, this relationship is non-linear and tends to flatten out with higher fish abundance or biomass (Plourde Lavoie, 2014).Such angling data are expected to vary linearly with fish abundance only if the level of capturability remains the same among fish populations.Yet, the capturability of brook charr depends on their gregarious behavior and the angler's ability to search and find "hotspots" of greater capturability in a given lake.This may result in a non-linear relationship between fish abundance and fishing success (Lewin et al., 2006) whereby capturability increases as fish abundance decreases, which may in turn provide erroneous information to fishery managers and hide an accelerated declining population (Post, 2013).With over 11,000 fishing lakes exploited by the SÉPAQ over several years with different levels of fishing pressure (SÉPAQ, 2021), the fishing angler's database represents an excellent source of comparative information for testing environmental DNA as a method for estimating population abundance without depending on fish capturability.
The brook charr is a salmonid inhabiting fresh, well-oxygenated lakes and streams throughout North America (Bernatchez & Giroux, 2012).We chose this species for this study because of its socio-economic importance in Québec (with over 340 M$ CAD generated annually from recreational fisheries) and concerns about its overexploitation, with over 16 million captures with a catch release rate of 25% annually (MFFP, 2020).With over 3.5 million anglerdays, brook charr is the most sought-after fish in Québec by recreational anglers.In the eastern United States and eastern Canada, brook charr populations have been declining over the last century, partially due to overharvesting (Hudy et al., 2008).The integration of a new complementary (and cost-effective) method such as eDNA analysis could help reduce the declining rate of exploited brook charr populations throughout its native North American range.The eDNA analysis could potentially be used, for example, to determine whether populations are above critical fisheries conservation thresholds (Spear et al., 2021).Several eDNA studies have previously been conducted on this species, including several demonstrating positive relationships between eDNA and brook charr abundance (Baldigo et al., 2017;Lacoursière-Roussel, Rosabal, et al., 2016;Wilcox et al., 2016;Yates et al., 2021).More specifically, our goal was to document the relationship between angling data (as an abundance indicator) and the concentration of eDNA, and several environmental parameters for 30 exploited brook charr lakes to further refine the application of eDNA in a fisheries management context.
As sampled eDNA concentrations likely derive from recent eDNA sources, we primarily compared angling data from recreational fishing from the same year as the eDNA sampling.However, given the correlation that Lacoursière-Roussel, Côté, et al. (2016) observed previously between eDNA concentration and S. namaycush abundance collected in different years, we also explored the relationship of eDNA concentration with angling data from previous years to verify whether a single eDNA sampling event on exploited lakes tends to maintain the same level of explained variability over the years.

| Study species and angling data
Our study was performed in two wildlife territories managed by the SÉPAQ in Québec, Canada: Parc national des Grands-Jardins (47.67°N, 70.63°W) and the Réserve faunique des Laurentides (47.24°N, 71.23°W) (Figure 1).Thirty exploited populations of brook charr were selected in agreement with local SÉPAQ managers to cover a wide variety of lakes, both in terms of surface area, fishing success and fishing pressure (Table 1).All lakes are excellent brook charr habitats due to the cool, clear, and well-oxygenated waters of F I G U R E 1 Geographical locations of the 30 sampled lakes in the Parc National des Grands-Jardins and the Réserve Faunique des Laurentides, Québec, Canada.Red and black squares represent the lakes sampled in 2019 and 2020, respectively.Each acronym refers to the name of the lake associated in Table 1.
this region as well as the absence of competitor species and thus offer excellent opportunities for recreational fishing.Most of the lakes are allopatric for brook charr (meaning no other fish species are present), except for lake Boisvert (BOIS), which also supports Arctic charr (Salvelinus alpinus) in small numbers.Angling data of the 30 exploited brook charr populations in the last seven years were provided by the SÉPAQ (, 2020).

| Environmental DNA survey
The protocol used for the eDNA survey was inspired by the one developed and successfully applied in our laboratory for lake trout (S. namaycush) (Lacoursière-Roussel, Côté, et al., 2016)  Notes: Angling data include HPUE (fish harvested/angler-days), FQI (kg/angler-day), fish yield (kg/ ha), fish density (fish harvested/ha), mean mass of fish (g), and fishing pressure (angler-days/ha).Environmental parameters include temperature (°C), surface (ha), depth (m), dissolved oxygen (mg/L), pH, and mean eDNA concentration (copies/L) as the response variable for the 30 lakes studied in both 2019 and 2020.Dissolved oxygen and pH from 2019 could not be taken due to technical issue (see Methods).Lakes from 2019 and 2020 were sampled between the 10 and 19 of June and the 6 and 16 of July, respectively.Angling data information was obtained from SÉPAQ (2020).
TA B L E 1 Lake name (abbreviation), angling data, and environmental parameters for each lake included in the study collected, for a total of 600 samples.Water samples were collected between 0.5 and 2 m in depth and at 5 m to 10 m from the shoreline to collect in the areas most frequented by brook charr in early summer (Figure 2) (Bourke et al., 1997;Mucha & Mackereth, 2008).Cote et al. (2020) showed that brook charr can spend over 70% of their recorded time in less than 5 m depth.Similarly, Bourke et al. (1997) observed that only 18% of brook charr in a lake similar to those surveyed here spent their time away from the littoral zone in summer.
A systematic design was used to determine the location of the sampling points to obtain an average eDNA concentration as representative as possible of a given lake.The perimeter of each lake was divided by the number of samples (20) to determine sampling distance and to ensure equal spacing ranging from 28 to 330 m between samples.To reduce interannual and seasonal variation in eDNA concentrations among samples, the sampling period was scheduled for the 2nd week of June (June 10-19, 2019, andJune 8-17, 2020).However, due to complications associated with COVID-19 pandemic, sampling was delayed in 2020 from mid-June to July 6-16 of 2020.A ProDSS multiparameter digital water quality meter (YSI ProDSS) was used to measure the temperature, pH, and dissolved oxygen (DO) at the deeper location of each lake.Due to malfunctioning issues with the meter in summer 2019, pH and the DO measurements could not be taken.
To reduce the risk of contamination, the 2-L bottles used for sampling were previously decontaminated by autoclaving and rinsing with a 10% bleach solution, followed by a distilled water rinse.
To minimize eDNA degradation, water samples were kept chilled in coolers to be filtered on the same day with a vacuum water pump and a 1.2 μm glass microfiber filter (Whatman GF/C) for 2019 and 0.7 μm for 2020 (Barnes & Turner, 2016).Different filter pore sizes were used between years due to the COVID-19 shortage.In 2020, we used the closest available filter pore size smaller to the size used during the 2019 sampling to ensure we collected as much or more eDNA than the 2019 eDNA sampling (Barnes et al., 2020).However, we acknowledge that the 1.2 μm pore size filter could potentially collect more brook charr eDNA than the 0.7 μm filter as experimentally shown by Lacoursière-Roussel, Côté, et al. (2016).After water filtration, filters were stored in plastic bags dried with silica gel microbeads and frozen at −20°C to minimize eDNA degradation before laboratory extraction.To prevent cross-contamination between lakes, all sampling equipment (cooler, boat, life jacket, and boots) was decontaminated by spraying with a 10% bleach solution.In addition, all filtration equipment was also decontaminated with a 10% bleach solution and 30 minutes of UV light exposure before being individually sealed in a plastic bag to avoid contamination between samples.
To verify the effectiveness of decontamination procedures and the absence of contamination, two 2-liter bottles of distilled water per lake were used as a negative control during the field sampling and following the same filtration and extraction protocol.Two field crews were assigned to both filtration and field sampling to allow same-day filtration and to reduce filter contamination.Accommodation and logistical support were provided by the SÉPAQ, making sure to avoid any conflict with anglers on lakes.volume (800 μl) of AL buffer.The 5-ml tubes were then vortexed and incubated at 70°C in a water bath for 10 min.Afterward, 1 volume (800 μl) of 95% ETOH was added to the mix which was then vortexed.Then, 625 μl of the mix was loaded into a DNeasy spin column (Qiagen, Inc) and centrifuged for 30 s at 13,000 rpm, repeated until all the remaining lysate was transferred.The next step was to place each DNeasy spin column into a new collection tube and washed the column by adding 500 μl of AW1 buffer and centrifuged for 30 s at 13,000 rpm.Spin columns were then washed twice with 500 μl of AW2 buffer and centrifuged again for 30 s at 13,000 rpm both times.

| Environmental DNA extraction
Spin columns were centrifuged for 1 min 30 s at 13,000 rpm, and the remaining flow-through was discarded (column drying).Finally, 80 μl of warm nuclease-free water was added and columns were incubated at room temperature for 10 min.Then, each column was centrifuged for 2 min at 13,000 rpm and split into two storage tubes of 40 μl each before moving them into a freezer dedicated to eDNA projects at −20°C.Laboratory controls not containing eDNA were analyzed with the same steps explained here without a filter to act as a negative control during extraction manipulation.

| eDNA amplification
We used qPCR to estimate eDNA concentrations in MicroAmp fast were then compared with the standard curve known concentration to estimate their eDNA concentration as copies per liter.The limit of detection (LOD) of the assay was calculated by diluting synthetic DNA (gBlocks™) into nine concentrations (1000, 500, 250, 50, 10, 4, 2, 1, and 0.5 molecules/μl) with 10 replicates each (Forootan et al., 2017).A 95% threshold approach was applied to analyze the LOD as recommended by Klymus et al. (2020).

| Model construction
All statistical analyses were performed using the R 4.0.2software (R Core Team, 2022 was constructed to establish the threshold where models become less plausible to correctly predict eDNA quantity.The mixed models were built using the package glmer.nbwith the variable "lake" as a random effect accounting for the variability between lakes.Angling data from two lakes (WOOD and BOIS) in 2020 were not available for their sampling year.Angling data of a given lake tend to stay consistent over the years (Figure 5), thus we replaced the missing data of WOOD and BOIS with their respective 2019 angling data.
Negative binomial mixed models were used because it has been shown that such distribution reproduces well-observed eDNA concentrations (Chambert et al., 2018;Furlan et al., 2016).
Indeed, it is expected to observe more extreme values than a normal distribution due to the stochastic capture of eDNA clusters, stretching the observed eDNA distribution to the right (Lacoursière-Roussel, Côté, et al., 2016;Lacoursière-Roussel, Rosabal, et al., 2016;Wilcox et al., 2016).Also, outliers may be due to clusters of highly concentrated eDNA from a recent secretion event that did not have time to disperse to the rest of the lake, causing a high concentration of eDNA locally (Klymus et al., 2015).
The collinearity of the explanatory variables was also tested by performing a Pearson coefficient test between each pair of them.
In addition, variance inflation factors (VIFs) were evaluated to ensure the absence of multicollinearity between the variables.A VIF >5 is recognized as a strict value for the presence of multicollinearity (Legendre & Legendre, 1998).To respect the assumptions, explanatory variables with a Pearson correlation > 0.5 or a VIF > 5 were used in separate models.As a result, both fish density and fish yield were extremely collinear and showed almost identical linear regression with eDNA.The same observation was made with HPUE and FQI (Gaudet-Boulay et al., 2022).Thus, we only provide results of the linear regression with fish density and HPUE to avoid redundancy but kept fish yield and FQI in the model selection analysis as they are still used as indicators of fish abundance in these brook charr populations.Model assumptions were checked on the residuals of the global model for each year, which includes all model parameters.

| Model selection
Three selection criteria were used to determine the best model.First, the models within a Bayesian information criterion (BIC) value <2 of the best model were retained as the most plausible.Secondly, the conditional R 2 (variance explained by both fixed and random effects) and marginal R 2 (variance explained by fixed effects only) were also examined during model selection using the variable lake as a random effect, as recommended by Nakagawa and Schielzeth (2013) for mixed models using the r2 function from the package performance.
Thirdly, we examined the coefficient of variation (R 2 adj ) and the significance level (p-value) from the linear regression between the mean lake eDNA concentration and the best variables retained from the most plausible models from 2019 and 2020.We compared eDNA concentration to angling data from their respective sampling season but also for each year between 2014 and 2020 and the mean value across those years.In this way, we compared angling data to a more representative value of eDNA concentration per lake.Assuming similar correlation across years (and it was the case, see Results), this approach would also be more practical from a manager's point of view to monitor fish abundance with only one mean eDNA value to compare between each lake.Finally, we tested the impact between the fishing pressure (angler-days/ha) and the relationship between eDNA concentration and the fishing success (fish harvested/anglerday). Higher fishing pressure, involving a more extensive sampling of the whole lake's area, should give a more representative fishing success and could thus enhance the relationship with eDNA concentration.

| Brook charr eDNA detection and concentration
Brook charr eDNA was detected in all 30 lakes surveyed (Figure 3).
In 2019, the average concentration of eDNA collected ranged from 313 copies/L in Lake Chavaudray to 185,033 copies/L in Lake Sirois.
In 2020, the average eDNA concentration collected per sample was much lower than in 2019, ranging from 60 copies/L in Lake Craine to 8114.1 copies/L in Lake Malbaie (Table 1).No amplification from both laboratory and field controls was observed except for one of the six replicates of the Lake NOEL field control (B1) (521.1 copies/L).Some qPCR sample replicates did not detect brook charr eDNA (13 from LAND, 18 from NOEL, 1 from LILY and 6 from CRAI) and were discarded from further analysis.Furthermore, LILY eDNA concentration appeared to have two distinct groups of eDNA concentration with samples from 1 to 10 showing the lowest eDNA concentration of 2019 lakes (mean = 647 copies/L) while samples from 11 to 20 showed a higher concentration of eDNA (mean = 8781 copies/L).
Synthetic DNA was successfully detected in all positive controls.
The R 2 value of the qPCR standard curve ranged from 0.993 to 1, and the efficiency ranged from 94.3% to 106.2%.Inhibition of the polymerase during qPCR can lead to more cycles to obtain twice as many molecules as in the previous cycle, which can potentially affect efficiency (BioSistemika, 2021).It is accepted that the accuracy can vary between 90% and 110%, which is what we observed (Thermofisher, 2021).The limit of detection (LOD) for qPCR assays was at 4 DNA copies per reaction.

| Environmental DNA variance among lakes
We observed a pronounced difference in eDNA concentrations between 2019 and 2020 and a high level of variability between lakes within the same year.Thus, the average eDNA concentration in 2019 was 33,977 ± 37,530 copies/L (mean ± SE), while in 2020 the average was 1049 ± 1127 copies/L (Table 2).High standard deviations were observed within lakes due to the right-skewed distribution of eDNA.
The range of coefficients of variation (CV) was similar in both years even if in 2019 two lakes had a CV greater than 100% (ESPE, LAND), ranging from 14% (SIRO) to 191% (ESPE) in 2019 and from 26% (MARA) to 97% (MALB) in 2020 (Table 2).Lake temperature in 2019 ranged from 9.5 to 16.8°C with an average of 13.2 ± 2.0°C, which was significantly lower (t-test, p < 0.001) than the 2020 lake temperatures, which ranged from 18.8 to 22°C with a mean of 20.5 ± 1.0°C (Table 1).The pH and DO measurements in 2020 (not available in 2019) varied little among lakes with an average of 6.7 ± 0.6 and 7.9 ± 0.6 mg/L, respectively (Table 1).

| Mixed-model performance
For lakes sampled in 2019, both fishing success (HPUE) and fishing quality index (FQI) were not significantly correlated with observed eDNA concentrations (Table 3).The mixed-model analysis found that the only variable that had a significant effect on observed eDNA concentration in 2019 was "Temperature," which was negatively correlated with lake eDNA concentration and explained 15% of the observed variation among lakes (Table 3).Similarly, fish density and fish yield derived from angling data was not correlated with eDNA concentration in 2019.The best mixed model included "Temperature" (BIC = 6390, Marg.R 2 = 0.15), but most of the variation in eDNA was explained in combination with the random effect "lake" (Cond.R 2 = 0.81) (Table 3).The predictive power of the best-selected model was not strong since it was no more than two BIC units apart from the "Null" model (Table 3).  2 In 2020, the mixed-model selection revealed again that HPUE and FQI were not correlated with the observed eDNA concentration (Table 4).Unlike 2019, the parameter Temperature did not explain the observed eDNA variation (Table 4).Among the mixed models from 2020 lakes, the best models retained were those with fish density (BIC = 4465, Cond.R 2 = 0.73, Marg.R 2 = 0.44), and fish yield (BIC = 4468, Cond.R 2 = 0.73, Marg.R 2 = 0.38) (Table 4).The fish density model was the best predictive in 2020 being above 7.8 BIC units from the random model and 2.9 BIC units above from the second-best model (fish yield).

| Linear regressions against mean lake eDNA concentrations
The average fish density from 2014 to 2020 strongly correlated (p < 0.001) with the mean eDNA concentration of the 2020 sampled lakes, explaining 76% of the observed variance (Figure 4a).The average fishing success (HPUE) from 2014 to 2020 correlated weakly but significantly (p = 0.03) with the mean eDNA concentration of the 2020 sampled lakes (adj.R 2 = 0.25) (Figure 4b).Both fish density and HPUE from 2014 to 2020 were weakly correlated with the observed eDNA concentration in 2019 (Figure 5a).In contrast, linear regression with fish density from 2014 to 2020 explained variation in eDNA concentration among lakes quite well, with 70% of the variance being explained in 2020 (Figure 5b [G]).Linear regressions between eDNA concentration from 2020 and angling data from 2014 to 2020 showed a consistent correlation over the years (Figure 5b).
In 2020, both HPUE and FQI showed a significant interaction with fishing pressure (Table 4).To examine this interaction, lakes were assigned to the closest value of the following three categories of fishing pressure: low (1 ad/ha), medium (2.5 ad/ha), and high (5 ad/ ha) using the Visreg function.Value of each category was selected to avoid a disproportionate number of lakes in one category and to ensure a normal distribution between them; 4 lakes in low (MARA, CHEN, MICH, and CINT) 7 in medium (BOIS, CHAR, ANNE, CRAI, CARR, STYM, and BOUI) and 4 in high (WOOD, FAGU, BLAN, and MALB).For lakes in the high-pressure category, more pronounced fishing pressure increased both the correlation between HPUE and FQI with the average eDNA concentration of 2020 lakes, while low fishing pressure lakes showed little or no correlation between eDNA concentration and the same predictors (Figure 6).

| DISCUSS ION
The main objective of this study was to document the relationship between angling data and eDNA concentration among several exploited brook charr populations.Whereas this relationship is typically explored using standardized surveys (Yates et al., 2019), no previous study has documented relationships between eDNA concentration and angling data.Our results revealed a significant correlation between mean eDNA concentration with fish density (adj.R 2 = 0.70) and fish yield from lakes sampled in 2020, while no correlations were observed with the two other angling data used as abundance indicators, namely HPUE and FQI (Gaudet-Boulay et al., 2022).The same tendencies were found with the average value of these angling data from 2014 to 2020, revealing consistency over the years using a single eDNA sampling event.Poor relationships were found in 2019 for all angling variables, while Temperature explained some of the observed variation of eDNA concentration.Our most plausible models did not incorporate environmental variables although they likely affect eDNA concentrations (Boivin-Delisle et al., 2021;Caza-Allard et al., 2021;Harrison et al., 2019;Saito & Doi, 2021).
The modest intra-year variation in environmental variables among lakes could have diminished their predictive power and could explain the absence of temperature influence on eDNA in 2020.However, in both years the most plausible mixed models included high and constant conditional R 2 (cond.R 2 = 0.81 in 2019 and cond.R 2 = 0.73 in 2020) while having lower marginal R 2 (marg.R 2 = 0.15 in 2019 and marg.R 2 = 0.44 in 2020) (Tables 3 and 4).Thus, the random factor (lake) explained a greater proportion of the variance observed among lakes, with unique combinations per lake of both abiotic and biotic characteristics (e.g., temperature (Caza-Allard et al., 2021), chlorophyll, pH, dissolved oxygen (DO) (Boivin-Delisle et al., 2021), fish spatial distribution (Thalinger et al., 2021), and foraging behavior (Klymus et al., 2015)) that could explain the relatively high and similar Cond.R 2 within all models.Our results also revealed a pronounced temporal variation in eDNA concentration between 2019 and 2020 likely due to the different sampling periods and filter pore size.Although one of the initial objectives of this study was to develop a predictive model with all 30 lakes, the constraint of having to use different filter pore sizes between the two years as well as the inter-year difference in eDNA concentration (possibly due to environmental and fish behavioral differences) prevented us to incorporate all lakes into a single model.Thus, we treated the two years independently to document their relationship with angling data and investigated the variables explaining the large inter-year difference.

| Relationship of eDNA with angling data as abundance indicators
Fish density (fh/ha) and fish yield (kg/ha) correlated strongly with eDNA concentration measured in 2020 but weakly in 2019.The results from 2020 agree with previous studies that investigated the relationship between these two parameters and eDNA concentration (Spear et al., 2021;Yates et al., 2021).However, these studies did not compare eDNA to angling data, instead employing standardized surveys.Notably, in a gillnet survey, Boivin-Delisle et al. ( 2021) highlighted the utility of eDNA as a monitoring tool by demonstrating a significant positive relationship (R 2 = 0.71) between eDNA concentration and walleye density.The few studies on brook charr populations to date showed moderate-to-high correlation with eDNA.Yates et al. (2021) revealed a strong correlation between eDNA concentration with fish density (fish/ha, R 2 = 0.59), biomass 26374943, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/edn3.341 by Cochrane Canada Provision, Wiley Online Library on [15/11/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License (kg/ha, R 2 = 0.63), and an even stronger relationship when adjusting with brook charr allometric scaling (ASM/ha, R 2 = 0.78).Baldigo et al. (2017) also showed a moderate level of correlation with stream brook charr eDNA concentration with population density (R 2 = 0.44) but lower correlation with population biomass (R 2 = 0.25) using electrofishing surveys to estimate fish abundance.Wilcox et al. (2016) found a similar correlation between brook charr eDNA and fish density (R 2 = 0.59) based on electrofishing surveys in small streams.Our results are most directly comparable to those of Yates et al. (2021), as we examined eDNA concentrations in lentic systems where horizontal dispersion is limited (Goldberg et al., 2018;Takahara et al., 2013).
Our results also provide evidence, at least in some situations such as the lakes we surveyed here in 2020, that a snapshot survey of brook charr eDNA correlates with usual angling data (fish density and fishing yield) collected over multiple years.This finding could be useful for managers as one eDNA sampling event could inform on the temporal stability of differences in fish relative abundance among different lakes being exploited.The relatively strong correlation we observed may also be due to the relatively large scale of the study which covers two sets of 15 lakes with 20 per lake.Indeed, to the of our knowledge, this is one of the largest studies to document the relationship of eDNA concentration with abundance indicators in a fishery management context.
In contrast to fish density and fish yield, the fishing success (HPUE) and the fishing quality index (FQI) correlated poorly with the mean eDNA concentration of exploited brook charr populations in both years as well as through 2014 to 2020.To our knowledge, there is no available study that explored the relationship between HPUE/  correlation (Spear et al., 2021;Yates et al., 2021).Interestingly, all of these studies have shown strong correlation of eDNA with population density and biomass.
Since brook charr is a gregarious species usually found in its localized preferred habitats (Bourke et al., 1997;Cote et al., 2020), anglers will tend to actively seek out those areas with the highest catch rate.If fish harvested in those preferred habitats are rapidly replaced by others, then fishing success will remain high while the abundance of the lake population decline (Dassow et al., 2020;Post, 2013).behavior can also affect the catchability of brook charr, where less experienced anglers (catching less fish per unit of effort) will tend to avoid lower-density lakes (Post, 2013).In contrast, experienced anglers often tend to exploit lower-density lakes where there are larger fish sizes and more restrictive regulations (Hutt & Bettoli, 2007), leading to higher than expected catch success in such lakes.These factors may result in a state of socalled hyperstability, where populations are thought to be abundant due to stable or high fishing success and no management action is taken despite declining and at risk of collapsing (Erisman et al., 2011;Post et al., 2002).Therefore, catchability may not be density-dependent (Zyll et al., 2002;Curry et al., 2003;)  explain the weaker relationship between HPUE and FQI with standardized surveys at higher fish abundance as discussed in the introduction (Plourde Lavoie, 2014).The poor correlations between HPUE and FQI with eDNA concentration could also be explained by the fact that eDNA is a density-dependent parameter as many studies have shown (Baldigo et al., 2017;Lacoursière-Roussel, Côté, et al., 2016;Shelton et al., Takahara et al., 2012;Yates et al., 2021).Yet, fish density estimated from angling data also has its limitation because it depends on angling fishing pressure.Thus, similar and sufficient fishing pressure is needed to obtain reliable fish density angling data (Arvisais, 2004).When this assumption is respected, fish density derived from angling data represents a good abundance indicator of exploited brook charr populations, and it could therefore mean that eDNA is also a good abundance indicator.However, we showed that even if is variation in fishing pressure over the years, fish density regressions against  eDNA concentration remain constant (Figure 5).With approximately 50% of exploited brook charr populations in Québec being overexploited (MFFP, 2020), one cannot rule out the possibility that this state of overexploitation could be due to misleading information by the HPUE and FQI, which overestimates fish abundance, hence delaying management actions.The interaction observed between HPUE/FQI and fishing pressure (angler-days/ ha) also supports this hypothesis (Figure 6).With an increase in fishing pressure (ad/ha), the chances of sampling the whole lake area are increased in contrast to lakes with low fishing pressure where areas with high catch rates are prioritized potentially inflating the estimated brook trout abundance.Fishing success and fishing quality index are currently more commonly used to inform management decisions, so it would be important to take into consideration the biases of these parameters discussed above.

| Heterogeneity in eDNA concentration between 2019 and 2020
Abiotic and biotic factors may act in concert influencing fish distribution in lakes, determining, in that way, the spatial distribution of their eDNA.Thus, an increase in temperature generally reduces the amount of eDNA present in the environment by increasing the rate of eDNA degradation (Caza-Allard et al., 2021;Jo et al., 2019;Stewart, 2019).Our results may reflect this phenomenon, with lower eDNA concentration found in the warmer 2020 lakes compared with the cooler 2019 lakes; temperature could be the main driver explaining the difference in eDNA concentration in both years (Table 1).thermal refuges at greater depth that corresponds to their metabolic preferences (Cote et al., 2020).The horizontal transport of eDNA in lentic systems can be limited (Goldberg et al., 2018).Consequently, the heterogeneous spatial distribution of brook charr between seasons could affect the distribution of lake eDNA across seasons, given that we observed lower surface water eDNA concentrations in July 2020 compared with June 2019 (Thalinger et al., 2021).Additionally, eDNA sampling in 2020 occurred after the prime time for brook charr angling in those lakes (mid-May to end of June), which would remove a considerable number of fish and the associated potential source of eDNA just before sampling.Changes in the methodology could have affected the observed difference in eDNA concentration between 2019 and 2020.However, in a study on brook charr, Lacoursière et al. (Lacoursière-Roussel, Rosabal, et al., 2016) compared the filtration efficiency of multiple sizes of filter between two temperatures and demonstrate that more eDNA was captured with the 1.2 μm than the 0.7 μm at 14°C but not at 7°C.We used the same type of 1.2 μm filter size at a similar temperature (~13°C), whereas the 0.7 μm was used at a higher temperature (~20°C).Also, Zhao et al. (2021) demonstrated that most eDNA particles collected were larger than 1.2 μm with very few eDNA copies being collected at smaller filter pore sizes, suggesting that the difference in filter pore sizes cannot easily explain the difference in eDNA concentration we observed between 2019 and 2020.Indeed, the summer 2020 eDNA concentrations we detected are more similar to those reported in other brook charr lentic (Yates et al., 2021) and lotic studies (Baldigo et al., 2017;Wilcox et al., 2016) than values we observed in 2019.For example, Yates et al. (2021) found that mean eDNA concentration ranged from 592 copies/L to 7805 copies/L in brook charr lakes sampled between June 30 and July 13 2018, closely matching our 2020 eDNA concentrations and sampling period (Table 1).Overall, these observations tend to support our interpretation that environmental conditions, more than the change in filter pore size, were the primary cause of variation in eDNA concentration between years.

| CON CLUS ION
In summary, our results provide further support to the possibility of obtaining a strong relationship between eDNA concentration and fish abundance indicators by using angling data as a large-scale example in a fishery management context.Our results also illustrate the importance of sampling during the same period of the year and using the same filter pore size to standardize surveys among lakes and years.Despite eDNA method limitations, we show that it represents a valuable and reliable tool to assist abundance estimation of fish exploited populations on temperate lakes.Additionally, the possible impact of both biotic and abiotic variables needs to be taken into consideration while doing eDNA study design and analyzing eDNA data (Evans & Lamberti, 2018;Pilliod et al., 2014).To better understand eDNA dynamics in nature and improve its use as a survey tool, it would be useful to document interannual eDNA variations in lakes.This information will contribute to improve sample design and will inform how eDNA concentration evolves with the community and environmental changes.Overall, our results show that the quantitative eDNA approach is a promising tool that could provide information on exploited fish population abundance at a lower cost and in a non-invasive way.

CO N FLI C T O F I NTE R E S T
There is no conflict of interest to declare.

Environmental
DNA was extracted by following the DNeasy Blood and Tissue kit extraction protocol from Qiagen Inc(Goldberg et al., 2011) adapted for Whatman GF/C filters.Extractions were performed under a UV cabinet in a dedicated laboratory for eDNA projects.All equipment (forceps, scissors, and rack) was sterilized with a 10% bleach solution and 30 minutes under UV.Pipettes were also dedicated to each project and sterilized under UV lamps for 30 min each day before manipulation.Gloves were changed between each sample to avoid contamination.First, each filter was cut in half and both of them were placed in a 5 ml tube with 720 μl of ATL buffer and 80 μl of Proteinase K (Qiagen, Inc.).The solution was then mixed using a vortex mixer and incubated at 56°C overnight.The next day, the cut filter and lysate were moved into a QIAShredder spin column (Qiagen, Inc.) and centrifuged for 30 s at 13,000 rpm.The remaining flow-through was transferred in a new 5 ml tube with 1 F I G U R E 2 Conceptual representation of the spatial lake design used for the eDNA sampling protocol and hypothesized distribution of brook charr.Bird view illustration (a) of twenty eDNA samples (red circle) that were equally spaced (perimeter/20) with one environmental profile sample taken at the deepest location of the lake (green circle) for each lake.Lateral illustration of the same lake (b) showing that each eDNA sample was collected between 0.5 and 2 m deep (red lines) while the environment profile sample was collected on the entire water column (green line).Blue squares (b) represent hypothesized eDNA distribution following brook charr habitat preference near the surface in the littoral zone in the spring.
optical 96-well plates with the 7500 Fast Real-Time PCR system thermal cycler (Applied Biosystems).The qPCR assay specific to brook charr (amplicon = 140 bp) used for eDNA quantification was designed byWilcox et al. (2013) with the following primer sequences; forward primer (BRK2_F) 5'CCACAGTGCTTCACCTTCTATTTCTA 3′, reverse primer (BRK2_R) 5' GCCAAGTAATATAGCTACAAAACCTAATAGATC 3′, probe (SAFO_Probe3) 5' ACTCCGACGCTGACAA 3′.This assay was designed to maximize the base pair difference in the primerbinding regions to enhance target specificity between closely related species to brook charr.Each eDNA sample and extraction control was mixed with the Environmental Master Mix 2.0 TaqMan solution (Applied Biosystems) as well as the SPUD assay to reduce and detect inhibition(Nolan et al., 2006).In each well, 2 μl of the sampled eDNA was amplified in a total reaction volume of 20 μl including 10 μl of Environmental Master Mix, 1.8 μl of each primer, 0.5 μl of probe, 1.2 μl of both SPUD primer, 0.5 μl of SPUD probe and 1 μl of target DNA following these stages: 2 min at 50°C, 10 min at 95°C and 50 cycles of 15 s at 95°C and 60 s at 60°C.A standard curve was used to quantify the number of DNA molecules amplified with 3 technical replicates for each known number of copies of synthetic DNA using gBlocks™ Gene Fragments (100,000, 10,000, 1000, 100, and 10 copies per reaction) designed from the COI sequence.The standard curve concentration was prepared identically as the other eDNA samples, except the SPUD was replaced by 3.9 μl of sterile water.Six replicates per eDNA sample were made as well as six negative control wells containing everything except the eDNA and one positive control well with synthetic DNA.The highest threshold line detected among all qPCR analyses of each year was set as the new respective threshold line to ensure a higher precision in the quantification of eDNA.Cycle threshold (Ct) values of each replicate

26374943, 0 ,
Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/edn3.341 by Cochrane Canada Provision, Wiley Online Library on [15/11/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License F I G U R E 3 Boxplot illustrating the distribution of eDNA concentration (copies/L) for (a) 2019 lakes and (b) 2020 lakes.Red diamonds represent the mean, and the dash lines in boxes represent the median.The complete names of each lake can be found in Table

FQI
from angling data and eDNA concentration.However, several studies examined the relationship between CPUE and BPUE with standardized surveys and eDNA.These returned mixed results, with some showing moderate(Lacoursière-Roussel, Côté, et al., 2016) to high correlation(Boivin-Delisle et al., 2021) and others showing low

F
I G U R E 4 Log-transformed mean eDNA concentration from 2020 lakes as a function of (a) the average fish density (mean of fish harvested /ha, adj.R 2 = 0.76) and (b) HPUE (mean of fish harvested/angler-days, adj.R 2 = 0.25) from 2014 to 2020.Regression lines are represented in blue with a 95% confidence intervals (CI) gray band and the p-value from the F-test.Confidence intervals (95%) of the mean eDNA concentration and both angling data are represented by vertical and horizontal error bars.Vertical gray dots represent each individual eDNA sample of a given lake associated with an angling data value.F I G U R E 5Log-transformed mean eDNA concentration (log(copies/L)) of (a) 2019 and (b) 2020 lakes as a function of each fish density (fish harvested/ha) and HPUE (fish harvested/angler-days) through 2014 and 2020.Regression lines are represented in blue with a 95% confidence intervals (CI) gray band, their respective adjusted R 2 value and p-value from the F-test.

Furthermore, the study
by Lawson Handley et al. (2019) showed strong seasonal differences in the distribution of eDNA in lentic systems, which closely matched the ecology of several lentic species (e.g., S. alpinus, Salmo trutta, and Perca fluviatilis).The study by Littlefair et al. (2021) also revealed a strong seasonal stratification of eDNA during summer and mixing during autumn, which reflected the thermal preferences of fishes in temperate lakes.Water layers become stratified and warmer in summer, and brook charr tend to prefer cooler temperatures.As a result, during warm periods they may seek

F
Interaction between the 2020 (a) HPUE (fish harvested/angler-days) and fishing pressure (angler-days/ha) (p = 0.02) and the 2020 (b) FQI (biomass/angler-days) and fishing pressure (ad/ha) (p < .01)from the 2020 lakes.The bands correspond to the 95% confidence intervals, and the triangle represents the value of each lake.Blue represents lakes with a high fishing pressure (~ 5ad/ha), green with a medium fishing pressure (~ 2.5ad/ha), and red with a low fishing pressure (~1ad/ha) (see Results for details).
major contributions to the design of the study; the acquisition, analysis, and interpretation of the data; and writing of the manuscript.EGM has made major contributions in the conception and design of the study; the acquisition and interpretation of the data; and writing of the manuscript.ML has made major contribution in the analysis, interpretation and writing of the manuscript.MY has made major contributions in the interpretation and writing of the manuscript.BB and CH have made major contributions in the acquisition of the data.GC made major contributions in the acquisition and the interpretation of the data, and the writing of the manuscript.AG made major contribution in the logistic and the acquisition of the data.LB has made major contributions to the conception and design of the study; the interpretation of the data; and writing of the manuscript.ACK N OWLED G M ENTSWe would like to thank all the field technicians of MFFP for their help in the eDNA and water sampling.A special thank you for all the lab members that help in every stage of this project: A. Perreault-Payette, C. Babin, and everyone else.Thank you Stéphanie Gagné of MFFP for her comments and suggestions in the making of the manuscript.We are also grateful to Associate Editor C. Jerde and 2 anonymous reviewers for their very valuable and constructive comments on an earlier version of this manuscript.This project was supported by the Ministère de la Faune, des Forêts et des Parcs (MFFP), the Société des Établissements de Plein Air du Québec (SÉPAQ), and a Strategic Partnership Grants for Projects from the Natural Sciences and Engineering Research Council of Canada (NSERC) to LB.This study is part of the research program of Ressources Aquatiques Québec (RAQ).

Lakes (abbreviation) Mean eDNA (copies/L) Median eDNA (copies/L) Standard deviation CV (%)
Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/edn3.341 by Cochrane Canada Provision, Wiley Online Library on [15/11/2022].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License and could Variables are defined in the Methods section.p-value represents the significance of the parameter on the response variable (significance threshold: p-value < 0.05), Cond.R 2 represents the conditional R 2 and the Marg.R 2 is the marginal R 2 (see methods for details).∆BIC corresponds to the difference in BIC with respect to the best-performing model (i.e., the one with lowest BIC) and models within 2 units from it are the most plausible.BICWt is the BIC weight of each model.The "Random" model represents a control model without predictors, which is only the random effect.If other models have a ∆BIC within 2 units or more of this model, then they are considered to have poor predictive power.
Variables are defined in the Methods section.p-value represents the significance of the parameter on the response variable (significance threshold: p-value < 0.05), Cond.R 2 represents the conditional R 2 and the Marg.R 2 is the marginal R 2 (see methods for details).∆BIC corresponds to the difference in BIC with respect to the best-performing model (i.e., the one with lowest BIC) and models within 2 units from it are the most plausible.BICWt is the BIC weight of each model.The "Random" model represents a control model without predictors, which is only the random effect.If other models have a ∆BIC within 2 units or more of this model, then they are considered to have poor predictive power.