Ecological indicators based on quantitative eDNA metabarcoding: the case of marine reserves

stands out. New genetic techniques such as environmental DNA (eDNA) metabarcoding have emerged and allow the detection of a wider range of species compared to conventional methods, but still fall short of providing reliable abundance estimations and subsequent ecological indicators. In this paper, we propose a combination of metabarcoding and quantitative polymerase chain reaction to obtain the quantity of eDNA molecules per species. This method was used inside and outside six no-take Mediterranean marine reserves to measure the effect of the protection on fish species and build a new indicator. Even if the total quantity of fish eDNA molecules was not different between the inside and outside of the reserves, we detected that cryptobenthic fish eDNA was significantly associated to the outside of reserves. Based on this observation, we propose a novel ecological indicator, the Demerso-pelagic to Benthic fish eDNA Ratio (DeBRa), taking advantage of the eDNA capacity to detect cryptobenthic reef fishes which are often missed by classical surveys. The DeBRa was significantly higher inside reserves, reflecting a higher relative quantity of eDNA molecules belonging to pelagic and demersal fishes under protection against fishing, therefore it appears to be a reliable eDNA-based indicator of human pressure. Furthermore, the DeBRa was not sensitive to habitat or environmental variations and does not require a complete reference database of eDNA sequences since it can rely on sequences assigned at the genus or family scale if possible and necessary.


Introduction
In marine and freshwater ecosystems, fish support a wide range of functions and contribute to nutrient cycles, carbon sequestration and the provision of vital nutrients to human populations (Holmlund and Hammer, 1999;Golden, 2016;Villéger et al., 2017;Mariani et al., 2020;Robinson et al., 2022). Through anthropogenic activities, fish and their habitats are under an ever-increasing number of threats such as overfishing or pollution at the local scale (Lloret et al., 2020;Cheminée et al., 2021) and climate warming at the global scale (Smale et al., 2019). These severe impacts jeopardize the nature's contributions to people provided by fish (Hicks et al., 2019;McLean et al., 2019;Robinson et al., 2019;Beger, 2021). To counteract these threats, Marine Protected Areas (MPAs) have become the main management tools on coastal ecosystems to maintain key habitats and viable fish populations (Grorud-Colvert et al., 2021). To reinforce conservation efforts in face of unprecedented biodiversity losses (Dulvy et al., 2021), international governing bodies have recently called for a new global target of 30% coverage by MPAs before 2030 (Jones et al., 2020).
Beyond coverage, there is a need to both comprehensively understand the effects of overfishing and to measure the ecological effectiveness of MPAs on fish communities. No-take or fully protected MPAs, also called marine reserves (Costello and Ballantine, 2015), have proven their ability to restore fish stocks and to produce more offspring from large individuals Ban et al., 2019). Yet, the effects of marine reserves on fish communities are not fully understood with some limited benefits  and even counter-intuitive results on biodiversity (Boulanger et al., 2021).
Whether they aim to measure management effectiveness or human impacts, various ecological indicators have been proposed (Meehan et al., 2020). However, most ecological indicators are based on surveys conducted with fishing gears, underwater visual censuses (UVCs), remote underwater videos (RUVs) and baited remote underwater videos (BRUVs), which (1) are selective and lead to biased estimates of biodiversity (Whitmarsh et al., 2017;Costello et al., 2017); (2) require highly advanced taxonomic expertise (Thomsen and Willerslev, 2015); or (3) are destructive or invasive (scientific trawling). The bias can be particularly striking when considering small species such as cryptobenthic fishes, hidden in the substrate, that cannot be correctly monitored with most fishing gears, UVCs or BRUVs (Smith-Vaniz et al., 2006;Alzate et al., 2014) while constituting the building blocks of reef ecosystems (Brandl et al., 2019). Large fish like sharks are also missed by most surveys since they are nocturnal, elusive, or rare (Boussarie et al., 2018). So new census methods with associated indicators are urgently needed to better assess the effects of protection, management and exploitation of coastal ecosystems.
Environmental genetic techniques such as environmental DNA (hereafter eDNA) metabarcoding have shown their ability to accurately detect the presence of a broad range of fishes in the marine environment, (Juhel et al., 2020) and elusive species (Boussarie et al., 2018), outperforming many classical techniques for fish inventories (Thomsen et al., 2012;Fediajevaite et al., 2021;Polanco Fernández et al., 2021). Yet, some studies highlight that eDNA provides poor estimates of benthic communities when considering a large variety of taxa and when sampling the water column (Antich et al., 2021;Laroche et al., 2020) but see Nichols et al. (2022). Moreover, detection ability depends on eDNA shedding and degradation rates that may vary within and among species (Thomsen et al., 2012;Sassoubre et al., 2016) but also with   (Harrison et al., 2019). eDNA detection currently relies on polymerase chain reaction (PCR), a process that may provide results biased in favor of certain species (Bellemain et al., 2010;Clarke et al., 2014;Kelly et al., 2019). Furthermore, species abundances, which inform various key indicators of ecosystem status (Stuart-Smith et al., 2013;Cinner et al., 2020), are still poorly estimated using eDNA in the marine environment; but see Stoeckle et al. (2021).
A recent review, based on 12 studies involving the eDNA metabarcoding of freshwater or marine fish species, reported a coefficient of determination (R 2 ), so a predictive power, ranging from 0.30 to 0.98 for fish counts and a R 2 of 0.75 for fish biomass in controlled environments (Rourke et al., 2021). In natural environments, the R 2 ranged from 0.24 to 0.68 for relative abundance, and from 0.26 to 0.45 for biomass. This difference was also reported in a meta-analysis on species-specific surveys, in which eDNA particle quantity explained 82% and 57% of the observed variation in fish abundance in laboratory experiments and natural environments, respectively (Yates et al., 2019). Indeed, many abiotic and biotic factors such as temperature, depth, activity, and metabolic rates render the estimation of species abundance from eDNA metabarcoding challenging (Rourke et al., 2021). Moreover, depending on possible primer mismatches, there is a chance to under-or overamplify the barcode sequence of some species or families (Kelly et al., 2019), which in turn would bias the relative abundance of eDNA molecules in a sample, after amplification. More research is needed to understand the extent of this bias and understand whether certain species, e.g. conservation-dependent species, could be underestimated by common markers. In general, this potential bias would be conserved throughout all samples and within families so it would not impact between-site comparisons.
In addition to the difficulty of linking the quantity of eDNA to species abundance, the distribution of eDNA molecules is often heterogeneously distributed in the environment (Turner et al., 2014;Turner et al., 2015). So, a large volume of filtered seawater (100 L) over long transects with many replicates may be required to obtain reliable estimates of species richness, but this sampling effort is often unrealistic in terms of deployment or cost (Stauffer et al., 2021). In light of these potential limitations, we need to test the ability of eDNA-based indicators to differentiate protected vs. unprotected marine areas using archetypal situations where fish communities are strongly affected by human pressure.
The Mediterranean Sea is amongst the most impacted marine areas (Cavan and Hill, 2021) with a strong decline of most commercial fish stocks (Colloca et al., 2013;Vasilakopoulos et al., 2014;Di Minin et al., 2019). MPAs of varying sizes and levels of protection have been implemented (Grorud-Colvert et al., 2021) and multiple studies have demonstrated the benefits of marine reserves on fish biomass (Giakoumi et al., 2017, Blowes et al., 2020. In this study, we conducted real-time quantitative PCR (qPCR) on eDNA samples filtered inside and outside six Mediterranean marine reserves to provide a quantity of eDNA molecules for each species in each sample. We make the hypothesis that the quantity of eDNA molecules is (i) higher inside protected areas owing to higher fish biomass and (ii) influenced by the position of fish in the water column with crypto-benthic species being less affected by fisheries (Boulanger et al. 2021) and less likely to release large eDNA quantities in the open sea than their demersal or pelagic counterparts. To test these hypotheses, we compared the quantity of eDNA molecules per species, and per category of species (i.e., cryptobenthic, benthic, demersal, and pelagic) inside and outside the reserves through two different approaches. First, we used statistical models to test the effect of protection and position in the water column on the quantity of fish eDNA molecules, while taking into account variability in fish body size. Second, we used a non-statistical approach based on Indicator Values to include both species eDNA quantities and occurrences within sites (Cáceres and Legendre, 2009). We then proposed the Demerso-pelagic to Benthic eDNA Ratio (DeBRa), a novel ecological indicator based on quantitative eDNA to measure the effect of protection or human pressure on fish communities, independently of habitat and environmental conditions.

Sampling protocol
eDNA samples were collected during the same season, in a twomonth period (May-June) with peristaltic pumps (1 L/min) one meter below the surface in six Mediterranean regions, each of them containing a fully protected marine reserve with a strict no-take policy known to provide higher fish biomass than fished areas (Sala and Giakoumi, 2018). In each region, two types of sites were sampled on the same day: inside and outside reserves (Fig. 1). We considered that contamination by dispersal between sampling sites, located 5 to 10 km apart, could be neglected given the strength of the currents in the area (1 km/day), eDNA persistence in the marine environment (24-48 h; Collins et al., 2018) and recent studies suggesting that eDNA signals display a high spatial fidelity (Murakami et al., 2019;West et al., 2020;Monuki et al., 2022;Jensen et al., 2022). Four replicates of 30L of seawater were sampled at each site along four 2 km-long transects using VigiDNA 0.2 µM filters. After filtration, the capsules were drained with air, filled with a CL1 buffer solution to conserve eDNA molecules, and stored at room temperature before eDNA extraction.

eDNA analyses
eDNA extraction, amplification, and sequencing were described in Boulanger et al. (2021). For the quantification of fish eDNA, the samples were amplified using the teleost-specific 12S mitochondrial rRNA primer (teleo) which has been extensively used for fish eDNA metabarcoding (Cantera et al., 2019;Cilleros et al., 2018;Coutant et al., 2021) and performs better (i.e. detects more species, with less bias and more specific amplification) than primers based on alternative loci (Collins et al., 2018;Weigand et al., 2019;Polanco Fernández et al., 2021). Then, the qPCR was performed in a final volume of 25 µL, which included 3 µL of DNA, 12.5 µL of SYBR® Green Master Mix (BioRad®), 7.3 µL of ddH2O, 0.5 µL of each "teleo" primer (10 mM), 1 µL of human blocking primer (100 µM; Valentini et al., 2016) and 0.2 g/μL of bovine serum albumin (BSA, Roche Diagnostic). Each sample was analyzed in 3 replicates. A dilution series of a synthetic gene of 10 -1 ng/µL, 10 -3 ng/µL, 10 -5 ng/µL, 10 -7 ng/µL (that correspond to 1,13E + 08, 1,13E + 06, 1,13E + 04 and 1,13E + 02 eDNA molecules respectively) was used as a standard for the qPCR analysis. The tubes containing the eDNA samples were sealed and the qPCR standards were added to the qPCR plate in a room separate from the eDNA extraction room. The qPCR cycle was as follows: 95 • C for 10 min, followed by 55 cycles of 95 • C for 30 s and 55 • C for 30 s. Melting curves were produced by plotting fluorescence intensity against temperature as the temperature was increased from 65 to 95 • C at a rate of 0.5 • C every 5 s. Samples were analyzed on a BIO-RAD® CFX96 Touch real-time PCR detection.

Incidence data
The OBITools package  was used to analyze sequence reads following the protocol described in Pont et al. (2018). The taxonomic assignment of reads was performed using the ecotag program with both an in-house local reference database and the sequences retrieved from the European Nucleotide Archive (ENA) database (Boulanger et al., 2021). The in-house database contains 320 fish species, corresponding to 41% of all teleost species in the Mediterranean Sea and to 75% of the regional species pool (Boulanger et al., 2021).
Reads showing 98% or more similarity with the reference database were kept. Taxa were preferentially assigned based on the local database, except if the similarity was higher for the ENA database. Considering the incorrect assignment of a few sequences to the sample due to tag-jumps (Schnell et al., 2015), we discarded all sequences with a frequency of occurrence < 0.001 per sequence and per library. The resulting dataset was manually checked to correct erroneous identifications and to remove foreign species and assignments that could not be assigned at the species level (Boulanger et al., 2021). The metabarcoding method provided a quantity of reads per species, i.e., the number of sequences that were detected for each species. Yet, this standardized quantity cannot be compared among samples because the initial eDNA quantity in each sample was yet unknown, but it allows to estimate species proportions in each sample.

eDNA quantity per species
We used a qPCR to obtain the total quantity of fish eDNA molecules in each sample, all species combined. This quantity was corrected in proportion to the quantity of reads that were excluded from the corresponding samples (foreign species, human contamination, unassigned taxa). For each sample, the corrected quantity of total molecules was then multiplied by species proportions in the corresponding sample obtained from the read metabarcoding outputs to calculate the quantity of eDNA molecules per species per sample. In more details, species proportions were obtained by dividing the number of metabarcoding reads of each species by the total number of reads before the bioinformatic analysis in each sample. Samples for which the qPCR did not provide a satisfying total quantity of eDNA (<10,000 molecules) were Table 1 Results of GLM that tested the effects of protection, species vertical distribution, length, trophic level, and the environment on the species detection probability and the quantity of DNA molecules. The reference for qualitative variables was the "cryptobenthic" category for vertical distribution, and "outside" for protection. We only show variables selected by the AIC criterion in the final models with statistical significance indicated by *: p < 0.05; **: p < 0.01; ***: p < 0.001. removed from the dataset.

Environmental variables
To describe the habitat and environmental variation among sites, we used the first four axes derived from a Principal Component analysis (PCA) (Boulanger et al., 2021). The PCA was computed using a set of habitat and environmental variables extracted for each site: habitat coverage, distance to land, sea surface temperatures (SST), mean bottom depths, and mean benthic and surface chlorophyll a. The PCA was applied to avoid multicollinearity between these variables and model overfitting (Fig. S1). The first four principal components explained 74% of the total variance in environmental conditions.

Modelling the quantity of eDNA
The quantity of DNA molecules per species and per transect was logtransformed (log10(x + 1)) before analyses. Our dataset at the species level contained a large proportion of zeros corresponding to either species absence or false negatives. We thus used a hurdle model, combining two subsequent models whose outputs were analyzed separately: a binomial generalized linear model (GLM) to predict species presence-absence and a gaussian GLM to predict the quantity of DNA molecules per species after excluding all zeros (Liu et al., 2019;Brown et al., 2021). We built these models to test the influence of protection and species categories on total quantity of eDNA molecules per transect (all species combined) and the total quantity of DNA molecules per transect per category of species: pelagic, demersal, benthic, and cryptobenthic species, as well as demerso-pelagic species (grouping demersal and pelagic species) and benthic species as a combination of cryptobenthic and benthic species. These categories were defined according to the vertical distribution of species in the water column, with pelagic species swimming in the water column, demersal species swimming near the bottom, benthic species living on the bottom and cryptobenthic species being smaller benthic species that mostly hide in the substrate (Brandl et al., 2018). Indeed, species with varying traits according to their morphological profiles (Friedman et al., 2020), and thus varying metabolic activities (Killen et al., 2010) may shed varying quantities of eDNA. Therefore, our models at the species level included the level of protection and fish vertical distribution but also species traits extracted from FishBase such as common length and trophic level to control for other sources of variability among species.
Because the age and size of the reserves ranged from 6 to 44 years of implementation, and from 65 to 1074 ha, we also included the age and size of the reserves. We normalized the age and size and the reserves such that the oldest or biggest value was 1, while fished areas were 0, following the path of Cinner et al. (2018). All models included the four first principal components of the previously mentioned PCA to control for habitat and environmental variability among sites and regions. To avoid collinearity between variables in the models, we used the variance inflation factor (VIF) to confirm that these variables could all be included in the model. We then performed a descending and ascending stepAIC (package MASS) to select the most parsimonious model to avoid overfitting and coefficient biases. Model fit and residuals were visually checked.

Indicator values and effect sizes
As a complementary approach to statistical models and to avoid the issue of zero-inflated distributions, we computed the Indicator Value index IndVal (Cáceres and Legendre, 2009) using the function multipatt (package indicspecies). IndVal combines two factors A and B, respectively corresponding to proxies of species specificity and fidelity for a given group of sites: A ij is defined as the mean abundance of the ith species within its jth protection level divided by the sum of mean abundances of this species over all protection levels, and B ij is defined as the number of sites where the species occurs inside its target protection level (Cáceres and Legendre, 2009). In other words, a species with 1 in specificity for protected areas will only be found in marine reserves, and a species with 1 in fidelity will be found in each reserve.
IndVal is, by construction, ranged between 0 and 100. The maximum value means that the species is only found in its target protection level and is always found in sites pertaining to this protection level. The minimum value means that the species is not found in sites of this protection level. Indicator values were computed for each species in protected (IndVal P ) and unprotected (IndVal U ) areas based on the quantity of eDNA. The following formula was then used to calculate the Effect Size of the reserves for the ith species: Because this method cannot account for environmental variability, we selected sites with similar habitat and environmental conditions by plotting the first three axes of the PCA and visually selecting a cluster containing both protected and unprotected sites with similar environments. The cluster represented sites displaying a value lower than − 0.5 on PC1 and comprised 8 protected sites and 28 unprotected sites. The effect sizes were computed on both the complete dataset and the cluster of environmentally similar sites. To test whether the Effect Size was significantly influenced by the species category while considering phylogenetic distances between species, we computed a Phylogenetic Generalized Least Squares model (pGLS) with the pgls function (package "caper"). The phylogeny used for this model was obtained via the Fish Tree of Life API (Rabosky et al., 2018).

The Demerso-pelagic to Benthic eDNA Ratio (DeBRa)
Considering that benthic species were recently found to be more present outside than inside Mediterranean reserves (Boulanger et al., 2021) and that the most threatened species, which are mainly demersal and pelagic, were found more present inside reserves , we computed the ratio between the log-transformed quantity of benthic fish eDNA molecules and the log-transformed quantity of demerso-pelagic fish DNA molecules to build a new indicator: The DeBRa was computed for each site, and the influence of protection was investigated using both a non-parametric rank test, and a gaussian GLM taking into account habitat and environmental variability.

Total quantity of eDNA per transect
The total quantity of fish eDNA molecules detected per transect ranged from 6,384 to 1,978,396 and followed a log-normal distribution Fig. 2a. The Wilcoxon test did not detect a significant difference between inside and outside marine reserves (W = 553, p-value = 0.343) (Fig. 3). The gaussian GLM revealed no effect of protection on the total quantity of eDNA molecules per transect but a strong and significant effect of environmental PC1 (Table 1).

Quantity of eDNA per species
The quantity of eDNA molecules per species per transect ranged from 0 to 1,965,527 with 94 different species detected throughout all samples ( Fig. S2, Table S1). Zero values represented 48% of the data, while positive values followed a log-normal distribution (Fig. 2b).
The binomial and gaussian GLM models revealed quite similar results, with significant effects of position in the water column and environmental PCA variables (Table 1). The probability to detect a species was similar inside than outside reserves (Fig. 4) while the quantity of eDNA molecules per detected species was markedly higher outside than inside reserves whatever the species category (Fig. 5). Regarding position in the water column, pelagic and cryptobenthic species had a lower detection probability (Fig. 4) and a lower quantity of eDNA molecules (Fig. 5) than benthic and demersal species. However, a significant interaction was found between protection and position in the water column in the binomial model. Pelagic species were more frequently detected inside marine reserves. The chance to be detected inside or outside reserve was equal for cryptobenthic species, and slightly higher inside reserves for benthic and demersal species, albeit non-significant. Furthermore, the gaussian model revealed that lower quantities of eDNA were detected for species with higher trophic level and larger body size. Reserve age had a positive effect on the quantity of eDNA but a negative effect on the detection probability.

Total quantity of eDNA per species category
The quantity of eDNA molecules per fish species category ranged from 0 to 1,967,876 and followed a log-normal distribution. Among all the models computed for each species category, none revealed a significant effect of protection on the quantity of eDNA (Table 1). After combining benthic with cryptobenthic species and demersal with pelagic species, no significant effect of protection was detected either. We only detected a strong effect of principal components corresponding to habitat and environmental variables in each of these models   (Table 1).

Indicator values and effect size
Indicator values ranged from 0 to 73.41 in reserves, and from 0 to 71.75 in fished areas, with uniform distributions for both. The Effect Size computed with these values ranged from − 0.50 to 0.64 and followed a gaussian distribution. The top 5% species with positive Effect Sizes are demersal or pelagic, while the top 5% species with negative values are cryptobenthic and benthic (Fig. S3). In the analysis restricted to a cluster of environmentally similar sites, the same pattern emerged: most positive values (higher inside reserves) correspond to demersal or pelagic species while most negative values (higher outside reserves) correspond to cryptobenthic or benthic species (Fig. 6). The pGLS successfully detected this pattern and indicated increasing Effect Sizes along the water column with significant or nearly significant effects of species category (Fig. 7).

Demerso-pelagic to Benthic eDNA Ratio (DeBRa)
Values for the DeBRa ranged from 1.0197 to 2.9410 and followed a normal distribution. The gaussian GLM revealed a significantly higher ratio of eDNA shed by demerso-pelagic fish than benthic fish in marine reserves compared to fished areas (p-value = 0.001). No significant effect PC2 (p-value = 0.1119) was detected, whilePC1, PC3 and PC4 were not retained in the model by the stepAIC (Fig. 8). Our model displayed a pseudo-R 2 of 0.16. This significant difference between both groups was also detected with a non-parametric rank test of Wilcoxon (p-value = 0.003), indicating a p-value more than one order of magnitude smaller than when taking only species richness into account (p-value = 0.05).

Discussion
Surprisingly, the GLM predicting the total quantity of eDNA molecules per transect did not reveal the expected effect of protection, but only a significant effect of environmental variables, while fully protected areas or marine reserves are notoriously known to promote higher fish abundance and biomass (Lester et al., 2009;Soykan and Lewison, 2015;Sala and Giakoumi, 2018). This result indicates that the total quantity of fish eDNA molecules cannot reliably be used as a proxy of fish biomass or as an indicator of management or fishing intensity.
At the species level, statistical models revealed a protection effect on the quantity of eDNA molecules. Even when accounting for environmental variability, which had a significant effect, there were more fish eDNA molecules outside than inside MPAs. Until recently, studies had only reported an increase in the abundance of pelagic over demersal species inside reserves (Libralato et al., 2010), but since classical surveys cannot offer such a species coverage for cryptobenthic species there are no previous records of increase in cryptobenthic abundance outside reserves. This finding could be explained by an hypoallometric relationship: a higher number of young and juvenile individuals may shed a higher eDNA quantity per gram of fish biomass than their adult counterparts (Maruyama et al., 2014;Takeuchi et al., 2019), and therefore provide a false signal of higher biomass, while decreasing the chances to amplify species with very low quantities of eDNA, which would explain higher quantities but lower detection rates outside marine reserves for some categories. Yet, our deep sequencing effort (1 million sequences) may partly mitigate this effect. Marine protected areas cover a relatively small area in the Mediterranean (6.81%) and considerable efforts have been made to create artificial nurseries far from predators (e.g., ports), which means that many coastal fish nurseries are outside MPAs (Med-PAN et al., 2016;Claudet et al., 2021). The large number of small individuals in these nurseries could explain higher quantities of eDNA molecules outside reserves. This result underlines the current complementarity of eDNA metabarcoding with conventional methods such as visual, acoustic or camera surveys that remain necessary to estimate fish body size.
The low quantities of eDNA molecules for pelagic species could be explained by their behavior since they often live in schools, are highly mobile and may have shed few and aggregated eDNA molecules in an area where they just passed. The differences in detection and in eDNA quantity between fishes with different positions in the water column could then partly be due to the variability in the shedding capacity between these categories, or simply to metabolic differences. This result is challenging to interpret and suggests that many hypotheses remain to be tested, notably in terms of biological trait influence on eDNA shedding. Finally, since (i) the substantial quantity of zeros demands to decompose our analysis in two subsequent models, each of which is missing a part of reality, (ii) our statistical models display a very low pseudo-R 2 and (iii) that no interaction could be detected in the gaussian model between protection level and position in the water column we suggest that working at the species level may not be an appropriate choice if the goal is to develop an ecological indicator based on the quantity of eDNA molecules. Yet, the protection effect was not significant either when computing models on each species category while the effect of environment was always significant as shown by Loiseau et al. (2021).
As an alternative, the Indicator Value approach, which is not based on statistical distributions and assumptions, has the advantage of using all the zeros in the dataset since it relies on both specificity and fidelity measures. In our case, after computing the Effect Sizes for each species, we observed a clear underlying pattern in the data: species with the highest effect size of protection were demersal or pelagic while those with the lowest effect size of protection were benthic or cryptobenthic. The pGLS further confirmed that the differences in effect size between species categories were significant while guaranteeing they were not due to phylogenetic similarities between species or other traits. Using Indicator Values we thus detected a marked pattern that could not be detected with species-based models.
Fishing seems to decrease the ratio between demerso-pelagic and benthic species eDNA quantity, meaning that fishing may favor benthic species to the detriment of demersal and pelagic species (Boulanger et al., 2021). In our study, the increase of the DeBRa indicator inside reserves indicates an increase in the number of demerso-pelagic fish eDNA molecules compared to benthic fish eDNA molecules, independently of habitat and environmental variables making this indicator transferable and comparable between ecosystems and seasons. Yet, this indicator does not provide information at the population level: we cannot determine if this increase is due to a simple increase in the density of adults inside reserves, or if it was caused by a massive increase in benthic juveniles outside during the spring season. Indeed, this hypothesis should be considered as many benthic species have pelagic early life stages, that could be less predated by demersal species, themselves under fishing pressure. Current research is exploring ways to use eDNA and environmental RNA (eRNA) to determine life-history stages from genetic material (Yates et al., 2021), so these methods could contribute to better understand whether unprotected areas may or not play an important and complementary role of nursery for some species.

Conclusion
Our findings suggest that eDNA metabarcoding coupled to quantitative real-time PCR can be used in field surveys to quantify eDNA molecules pertaining to each detected species. We show that Indicator Values could be used to detect some patterns that would not be detected by classical statistical models based on either species occurrences or abundances given the high proportion of zeros inherent to eDNA surveys. We also show that we cannot base our indicator on the total quantity of eDNA in a transect, and that the differences in quantities lie at the species level, with effects of many variables affecting the detected quantity of eDNA, such as body size, trophic level, environmental variables but also age of the reserve. More generally, our results emphasize the difficulty to evaluate benefits of MPAs on fish communities. We reveal antagonistic responses to protection between species categories, highlighting the need to base management strategies on multiple indicators and not only on the biomass of some targeted fishes. In the Mediterranean Sea, and more precisely in the Northwest regional pool, more than 75% of species are referenced in our genetic database, but in tropical regions which are species richness hotspots (Tittensor et al., 2010), a large proportion of species are missing in genetic reference databases . Since the DeBRa indicator only requires to determine if the species is cryptobenthic, benthic, demersal, or pelagic, it could be possible to assign unknown sequences at the genus or family level (for genera and families that have 95% of their taxa pertaining to the same category, for example), and thus confidently build this indicator without assignation at the species level. We expect that this new indicator can provide consistent results in other temperate and tropical regions to reveal human impacts or management benefits regardless of habitat, seasonal and environmental conditions.  Funding was provided by the Agence de l'Eau RMC, the Region Occitanie and the RESERVEBENEFIT project through the 2015-2016 BiodivERsA COFUND call with the ANR (France).

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Stephanie Manel reports financial support was provided by Rhone Mediterranee Corsica Region Water Agency. Stephanie Manel reports financial support was provided by French National Research Agency.
Stephanie Manel reports financial support was provided by Occitanie Region.