Marine protected areas (MPAs) have been increasingly developed to protect fish communities and restore their ecological services. Shallow water (above -30 m) fish populations may also find a natural refuge in the depths of the mesophotic reefs to escape overfishing and climate change. While a lot of knowledge has been acquired in last decades about the effect of MPAs and more recently about mesophotic ecosystems, thanks to advances in monitoring methods, little is known about the effect of the interaction between protection and depth. Using environmental DNA (eDNA) metabarcoding, we sampled rocky fish assemblages of the Cap Roux no-take marine reserve and the surrounding fished areas, covering several depth strata ranging from the surface to mesophotic depths. We found that protection, depth and their interaction had a significant effect on fish diversity and assemblage composition. The vast majority of the 66 fish taxa identified in this study were detected inside the MPA. On the other hand, depth had a negative impact with mesophotic reefs harbouring the lowest diversity. A reserve effect on species richness was observed for shallow reefs but not for the mesophotic zone. Shallow protected locations were composed of rich fish assemblages including targeted and threatened species that were rarely detected or even undetected outside the MPA and/or at mesophotic depths. In this study, depth does not seem to act as an efficient refuge for shallow communities. We suggest that the only effective protection comes from the MPA, highlighting the importance of these management tools for the long-term conservation of fish populations.

Study area

The Cap Roux MPA is a no-take marine reserve located in the North-Western part of the Mediterranean Sea between Cannes and Saint-Raphael, France (Fig. 1). Established in 2003, this MPA was initiated by the professional fishermen from the “Prud’homie de pêche” of Saint-Raphel, to sustain local small-scale fisheries. Located along the Estérel coast, a natural area with little urbanisation, the Cap Roux MPA covers an area of 450 ha spanning from the shore to the -80 m isobath (Bottin et al., 2020, www.medamp.org). This reserve features diversified fish communities and typical Mediterranean marine habitats such as Posidonia oceanica meadows, macro-algae dominated rocky reefs and corraligenous reefs. Outside the MPA border, fish populations are harvested by recreational (e.g., angling and spearfishing) and professional fishing practices (e.g., gillnet), targeting both shallow and mesophotic areas (information provided by the “Comité des pêches maritimes et des élevages marins des Alpes Maritimes”, the Alpes-Maritimes Fishing and Marine Farming Committee, https://www.cdpmem06.fr). Although surveillance is lacking, previous studies aiming to evaluate Cap Roux MPA effectiveness found a significant reserve effect (Seytre et al., 2013; Seytre & Francour, 2009). This part of the Mediterranean Sea is also heavily threatened by climate change with an increasing frequency and intensity of marine heat waves causing mass mortality events of benthic communities down to -30 m (Bramanti et al., 2023). Considering all the aforementioned characteristics, this reserve constitutes an ideal testing ground to address our questions.

Metabarcoding: water sampling 

Field sampling was carried out between September 26 and 28, 2023. The eDNA samples were collected in three depth ranges: one meter below the surface (hereafter: surface), between -20 and -30 m (hereafter: mid-depth) and between -50 and -60 m (hereafter: mesophotic). Hereafter, the term “shallow” will be used to refer to depths located above -30 m, so it includes surface and mid-depth sampling stations. For each depth, two sites were sampled inside the MPA and two sites outside. Outside sites were located 2 km from the MPA borders (Fig. 1). For each site, two replicates were simultaneously sampled for a total of 24 eDNA samples collected in this study (Table S1).

The eDNA samples were collected using two onboard diaphragm pumps (Argaly; flow: 1.0 L/min) for surface samples and an underwater double-head pump (Subspace; flow: 1.0 L/min) for mid-depth and mesophotic samples. The pumps were connected to a filtration capsule (eDNA water filter, Waterra; 600 cm2; Polyethersulfone; 0.45 µm), allowing on-site live filtration. Each sample consisted of 30 L of filtered seawater. The filtration was conducted from fixed point directly above rocky reefs (< 10 m from the substrate), either macroalgae-dominated or corraligenous reefs depending on the depth. The underwater pump was manually deployed from the boat using a rope, with a delayed start so the filtration process only began when the target depth was reached. A diving computer was placed on the pump to register depth and temperature, as well as a camera to verify afterward that the filtration occurred above rocky habitats (Fig. 2).

Immediately after sampling, 50 mL of Longmire buffer solution (Longmire et al., 1997) was injected into the eDNA capsules to allow the long-term conservation of eDNA before the laboratory procedure. eDNA samples were always manipulated with gloves to avoid contamination. Back at the lab, capsules were shaken vigorously and the eDNA extract was stored at room temperature in the dark until extraction (Roblet et al., 2024b).

During the field campaign, several negative field controls were performed to check for contamination that could have occurred during capsule handling on the boat. These controls consisted of 1 L of ultrapure water filtered with a capsule connected to the surface pump. They were treated the same way as true field samples.

Metabarcoding: lab processing 

DNA extraction and PCR amplification were performed by Argaly (Sainte-Hélène-du-Lac, France), in laboratories dedicated for handling eDNA samples. The 50 mL falcon tubes were centrifuged for 15 min at 15,000 g. DNA from the pellets was then extracted, using the NucleoSpin eDNA water kit protocol (Macherey Nagel), following the manufacturer’s instructions.

For PCR amplifications, three primer sets selected from our recent pilot study were used (Roblet et al., 2024b): Fish16S, targeting the 16S rRNA locus and specific to Actinopterygii species; AcMDB07, also Actinopterygii specific and targeting the 12S rRNA gene; and Vert16S, which targets the 16S locus but is Vertebrate specific (Table S2). These primer sets are highly effective to detect Teleostean and Chondrichthyan fishes, and are complementary, which justifies combining them (Roblet et al., 2024b). 

For each of these three primer sets, the extracted DNA from each sample was amplified in 12 replicates. Each PCR replicate was uniquely identified by a combination of two eight-base tags appended to the PCR primer at the 5’ end. These tags were used during bioinformatics analysis to assign sequences to the corresponding replicate. Following amplification, all samples were purified with the MinElute purification kit (Qiagen). Library constructions and sequencing were then performed by Fasteris (Geneva, Switzerland). The libraries were prepared according to the Metafast protocol, designed to minimize sequencing artefacts. The libraries were then sequenced in several Illumina MiSeq runs with paired-end reads of 2 x 250 bp (for AcMDB07 and Vert16S) or 2 x 150 bp (for Fish16S).

Various quality controls were conducted at each step of the protocol to identify potential contamination, ensuring an accurate interpretation of the results. For each primer set, the following controls were performed: three negative extraction controls, two negative PCR controls, two positive controls and 96 bioinformatic controls. The success of the amplifications and purifications was confirmed on a 2 % agarose gel (E-Gel Power Snap, Invitrogen).

Metabarcoding: Bioinformatics

Argaly conducted the bioinformatic steps, using the following procedure: the raw sequence data for each primer were analysed using the suite of OBITools programs (https://pythonhosted.org/OBITools/welcome.html; Boyer et al., 2016), which is designed specifically for metabarcoding data analysis. For each primer set, paired sequences were assembled (“illuminapairedend” command). Then, sequences with an alignment score >= 40 (i.e., corresponding to an overlap of 10 bases) were assigned to the corresponding amplification replicate, thanks to the tags inserted in the 5' end of the primers (“ngsfilter” command). The resulting dataset was dereplicated (“obiuniq” command), then filtered (“obigrep” command) to remove low-quality sequences (i.e., containing at least one N), sequences whose length did not belong to the length range observed in silico for the target group, and singletons (i.e., sequences observed only once in the dataset). Sequences sharing 97 % identity were grouped into clusters using SumaClust (Mercier et al., 2013). The abundances of sequences within each cluster were summed for each PCR replicate. The cluster head, representing the most abundant sequence in the cluster, was chosen as the representative sequence. Clusters appearing less than 10 times in a PCR replicate were deleted. A taxonomic assignment of the cluster heads was performed using the “ecotag” command to obtain a list of MOTUs (Molecular Operational Taxonomic Units). The reference sequences used for this taxonomic assignment were obtained by performing an in-silico PCR on the public sequence database GenBank (v.249) using the ecoPCR program (Ficetola et al., 2010). This in silico PCR was conducted using the PCR primers associated with each marker, allowing a maximum of three mismatches per primer and retaining only sequences assigned at least at the family level. The R package “metabaR” (Zinger et al., 2021) was then used to remove artefactual sequences that are present in low abundance in the metabarcoding data, but which may influence the ecological conclusions that can be drawn from them (Calderón‐Sanou et al., 2020). This process included removing i) MOTUs with sequence similarity to any sequence in the reference database below 0.95, as they are potential chimeras; ii) MOTUs whose frequency over the entire dataset is maximum in at least one negative control ("max" method of the “contaslayer” function), because they are potential contaminants; and iii) MOTUs with a relative frequency < 0.03% within a PCR replicate (“tagjumpslayer” function), because they are potentially artefacts generated during library construction for sequencing (i.e., "tag jumps"; Schnell et al., 2015). PCR replicates with a sequencing coverage < 100 sequences were also removed and then the remaining PCR replicates were aggregated by sample using the “aggregate_pcrs” function. Finally, MOTUs observed less than 10 times in a sample were recoded as absent in that sample.

A manual taxonomic cleaning was performed for each primer set dataset to remove non-fish MOTUs, freshwater fish MOTUs and MOTUs assigned to marine fish that do not occur in the study area. The assignation of the remaining MOTUs was carefully checked and modified, if necessary, based on biogeographic data or in case of incorrect taxonomy. As we were not interested in intraspecific diversity, MOTUs assigned to the same taxa were merged. The datasets of the three primer sets were then transformed into presence-absence data and pooled into a single common dataset.

Statistic analyses 

All statistical analyses were conducted in R v4.4.0 (http://www.R-project.org; R Core Team, 2024). The effect on fish assemblages of two main factors: protection status (2 levels: protected / fished) and depth (3 levels: surface / mid-depth / mesophotic), as well as their interaction, were assessed in this study. Statistical analyses were conducted on the full list of fish taxa detected, and for some analyses on a subset including only taxa targeted by professional or recreational fishermen in the region (hereafter: target taxa) (Table S3). Fish taxa were also pooled into trophic categories based on their trophic level following the classification of Aglieri et al. (2021) (Table S4).

First, two Venn diagrams (“VennDiagram” package; Chen & Boutros, 2011), a stacked bar plot and a heatmap were constructed to visualize differences in fish richness between depths and protection statuses. Then, the effects of these two main factors and their interaction on the average species richness and target species richness per sample were assessed using an analysis of variance (ANOVA). This ANOVA was performed after fitting a generalized linear model (glm) with a Poisson distribution (‘glmmTMB’ function of “glmmTMB” package; Brooks et al., 2017). Issues of overdispersion and zero-inflation were evaluated using the ‘testDispersion’ and ‘testZeroInflation’ functions (“DHARMa” package; Hartig, 2020). A residual check was conducted with the function ‘simulateResiduals’ from the same package. Three variables were used to take into account environmental variability between samples: seawater temperature, distance to shore and habitat diversity (i.e., number of habitats in buffer of 100 m radius around each sampling station, according to the “Donia expert” classification, https://plateforme.medtrix.fr). Seawater temperature was highly correlated to sampling depth and was therefore removed from the model (Spearman’s rank correlation coefficient: rho = -0.92, p < 0.001). The most parsimonious model, selected using the Akaike Information Criterion (AIC), only retained the two main predictors and their interaction but no environmental factor.

For beta diversity, a Principal Coordinate Analysis (PCoA) based on the Jaccard distance matrix was first computed to visualize the ordination of the data grouped by protection status and depth (“vegan” package; Dixon, 2003). A permutational multivariate analysis of variance (PERMANOVA) was then run on this distance matrix to evaluate the effects of these two main factors and their interaction on community composition (‘adonis2’ function of “vegan” package), followed by post hoc pairwise comparisons (‘pairwise.adonis2’ function of “pairwiseAdonis” package; Martinez Arbizu, 2020). The homogeneity of group dispersions between protection statuses and between depths was also analysed (‘betadisper’ function of “vegan” package). An Indicator Species Analysis (‘multipatt’ function of “indicspecies” package; Cáceres & Legendre, 2009) was performed to identify species significantly associated with a given depth or protection level.

To compare the compositional change of species across protection status and depth and to understand the underlying ecological processes driving it, we further calculated the Jaccard dissimilarity index (βjac), and its turnover (βjt) and nestedness (βjn) components (“betapart” package; Baselga & Orme, 2012). Turnover is a measure of taxon replacement among samples while nestedness quantifies the degree to which sample compositions are subsets of one another (Baselga & Orme, 2012; Legendre & De Cáceres, 2013). The βjac index ranges from 0, indicating the samples have identical species composition, to 1, indicating their species compositions are completely different. Finally, we tested whether βjac differed between paired depths using a Kruskal–Wallis rank sum test.