Shark and ray diversity, abundance and temporal variation around an Indian Ocean Island, inferred by eDNA metabarcoding

Sharks embody several major aspects of modern marine management: they are traditionally antagonized, exploited or by‐caught by humans, typically vulnerable to extirpation, pursued as luxury food, yet valued as wildlife and essential as top‐down regulators in marine food webs. Due to their generally large size, elusiveness, high mobility, and potentially dangerous nature, elasmobranchs pose substantial technical challenges to biodiversity monitoring, which prompted recent efforts to harness the power of environmental DNA (eDNA) as a noninvasive survey method for these taxa. Here we deployed an elasmobranch‐specific metabarcoding assay to characterize shark and ray diversity around Reunion Island, during the austral summer, detecting at least 14 species and a strong overall correlation between frequency of species detection and read abundance. Over 90% of sequence reads belonged to three large predators: scalloped hammerhead (Sphyrna lewini), tiger shark (Galeocerdo cuvier) and bull shark (Carcharhinus leucas). While the importance of tiger and bull sharks is well established in Reunion Island, and a major focus of the local shark control program, the prevalence and abundance of scalloped hammerhead has so far been grossly neglected. We also confirm the absence of typical “tropical reef sharks” around the island and reveal an important temporal fluctuation in tiger shark during the study period. Collectively, results show how eDNA can help circumvent barriers, bias and drawbacks associated with monitoring shark populations using visual and capture‐based techniques, and generate spatial and temporal biodiversity data on these species for rapid consideration by marine environmental managers.


| INTRODUCTION
Elasmobranchs (i.e., "sharks and rays") remain the least thoroughly assessed vertebrates on Earth (Dulvy et al., 2014), despite their established ecological significance (Ferretti, Worm, Britten, Heithaus, & Lotze, 2010) and their diverse, conflicting value to human societies (Dulvy et al., 2017). Sharks and rays, and in particularly the largest, most charismatic species, epitomize the current relationship between humans and the oceans, in that different interests-as disparate as seeing the same animals as both a lucrative food commodity (Cardeñosa et al., 2018) and an awe-inspiring tourist attraction (Mustika, Ichsan, & Booth, 2020)-converge on the same resource with seemingly irreconcilable positions (Booth, Squires, & Milner-Gulland, 2019).
One emblematic issue with elasmobranch conservation resides within the boundaries of tourism and recreational activities: sometimes, in a few particular regions, a handful of shark species pose a risk to humans (Chapman & McPhee, 2016). Unfortunately, this situation may lead to human-shark conflicts, for which shark culling is often identified as a solution (Meeuwig & Ferreira, 2014), despite recent global analyses showing this to be the main contributor to reef shark population collapse (MacNeil et al., 2020). While alternative solutions to culling exist to protect humans (O'Connell, Andreotti, Rutzen, Meÿer, & Matthee, 2017), reconciling shark conservation with human safety remains challenging due to the complex interaction of factors influencing risk perception (Le Busque, Dorrian, & Litchfield, 2021), especially in communities where shark bite incidents have recently occurred.
One of the world's shark bite hotspots (McPhee, 2014) is the small oceanic island of Reunion, situated over 600 km off the eastern coast of Madagascar, where shark bite incidence rate has risen sharply over the last decade (Chapman & McPhee, 2016), in part as a result of the rapid development of tourism and the increased popularity of surfing (Lagabrielle et al., 2018). Concerns over the potential negative impacts of sharks on watersports and other coastal leisure activities have simultaneously prompted the activation of a shark control programme, and an unprecedented level of research effort on the status of the elasmobranch fauna around the island (Blaison et al., 2015;Guyomard et al., 2019), mostly focused on the two species responsible for the unprovoked bites, bull shark (Carcharhinus leucas) and tiger shark (Galeocerdo cuvier) (Ballas, Saetta, Peuchot, Elkienbaum, & Poinsot, 2017;Guyomard et al., 2020). Such a scenario presents a major conservation challenge, because management strategies are being introduced reactively, in an environmental context with large gaps in ecosystem knowledge, and before scientists are in a position to form a comprehensive understanding of the broader ecological consequences of these management actions.
Recent studies have shown how newly developed environmental DNA (eDNA) metabarcoding techniques have the ability to rapidly generate large species inventories, which can readily portray differences in marine habitats and community types (Sigsgaard et al., 2019;West et al., 2020). Of particular relevance to sharks, these approaches have been used to compare regions under varying levels of anthropogenic impacts (Bakker et al., 2017), to detect previously unrecognized biodiversity features (Boussarie et al., 2018), and to follow seasonal migrations (Postaire, Bakker, Gardiner, Wiley, & Chapman, 2020). Here we set out to characterize the shark and ray community around the coasts of Reunion Island through the retrieval of elasmobranch DNA from temporally replicated water collections. We specifically aimed to: (a) identify the dominant elasmobranch species using a non-invasive method; (b) evaluate the effect of sampling depth on community composition; (c) examine spatial variation in the community; (d) explore possible temporal patterns over a 13-week sampling period covering the austral summer. Results amass several lines of evidence that can be of immediate relevance for the management of the island's coastal marine environment.

| Sample collection
We identified 16 sites as representative of the shelf reef habitats (Pinault et al., 2014) in Reunion Island, with a special emphasis on coral reefs, which are more prevalent on the west side of the island (Figure 1), and are the most vulnerable ecosystems, where increasing anthropogenic pressures have been observed in the last decade (Lagabrielle et al., 2018;Magnan & Duvat, 2018). Sites were located between 50 m and 1 km from the coast and were accessed by boat between December 2018 and March 2019. Each site was sampled between two and five times (with at least one week elapsing between events), taking two samples on each date, at both 10 m and 40 m depth, using 5 L Niskin bottles operated by a crane system from the boat (details on collections are shown in Table S1). We also collected three samples from the Reunion Aquarium in Saint Gilles, as positive controls, as well as to verify the potential influence of the Aquarium recirculation system on the samples from the corresponding coastal area of Roches Noires.
At the end of each collection day, water samples were transferred to a sterilized environment and filtered using a vacuum pump. For each sample, 4 L of seawater were passed through a cellulose filter (Millipore, 47 mm diameter and 0.45 μm pore size), after which each filter was folded, covered in silica beads, and stored at −20 C in 2 mL screw-cap tubes. On five occasions, 4 L of commercial mineral water were also filtered and processed in the same way, as "field blanks". In total, we obtained 114 filtered samples for later laboratory analysis (Table S1).

| Laboratory procedures
All samples were transferred to the eDNA facilities at the University of Salford, where different rooms are used for eDNA extraction, pre-PCR preparations and post-PCR procedures, in order to minimize the risk for cross-contamination. Full suits, masks and overshoes were used in the DNA extraction room, which was fully sterilized using UV lights and bleach at the end of each session. DNA was extracted from the filters using the Qiagen PowerSoil DNA Isolation Kit. During each session, a DNA extraction "blank" was also created using only extraction buffer, in order to measure the extent of contamination associated with the procedure. A total of 11 extraction blanks were processed. DNA concentration of extracts was checked using a Qubit fluorometer (ThermoFisher).
DNA metabarcoding was carried out by targeting a 171-bp 12S ribosomal fragment of the mitochondrial genome (Miya et al., 2015), using the slightly modified Elas02 primer sequences (Taberlet, Bonin, Zinger, & Coissac, 2018), which were predicted to have the greatest universality for elasmobranchs than any other primer set (Collins et al., 2019), while circumventing the nonspecific nature of COI markers (Bakker et al., 2017). For the purpose of sample identification following multiplexed Illumina sequencing, samples were amplified using primer pairs specific to each sample, containing F I G U R E 1 (a) Map of Reunion Island, showing the steep bathymetric profile, the positioning of all sampling sites and the prevalence of reefs to the West. The color-coding of sampling sites (blue: North, green: West, red: South, yellow: East) illustrates the groupings used to explore community composition differences around the island. b) Histogram representing the overall elasmobranch taxonomic composition (proportions of square-root-transformed read abundance) recovered at each site over the study period; sites are coded and labelled according to the four-group structure mentioned above 8-base oligo-tags attached to the forward and reverse primers (which differed from other samples' tags in at least three positions), and preceded by two-to-four fully degenerate positions, to inflate sequence diversity and generate well-balanced libraries. The full sequenced products were therefore >190 bp long. Three PCR "blanks," loaded with UV-irradiated nucleotide-free H 2 O in place of DNA template, were also produced as further negative controls.
Each of the 128 samples (106 from natural seawater, 3 from the aquarium tanks, 5 field blanks, 11 extraction blanks, and 3 PCR blanks) was amplified in 20 μL total volume, composed of 10 μL AmpliTaq Gold Master Mix (ThermoFisher), 0.16 μL BSA, 1 μL of 5 μM forward primer, 1 μL of 5 μM reverse primer, 10 ng in 2 μL of eDNA template, and 5.84 μL of molecular biology grade water. Thermocycling conditions included an initial denaturing step of 95 C for 15 min, followed by 40 cycles of 94 C for 1 min, 54 C for 1 min, 72 C for 1 min, and a final extension step of 72 C for 5 min. To reduce PCR stochasticity, three technical amplification replicates per sample were carried out and subsequently pooled. Amplification success was checked by electrophoresis, loading 5 μL of PCR product through a 1.5% agarose gel.
Two dual-indexed Illumina libraries were prepared using the KAPA HyperPrep PCR-free library preparation kit, which were quantified by means of qPCR using the KAPA library quantification kit (Roche). Libraries were then pooled in equimolar concentrations along with 10% PhiX (Illumina) serving as a positive sequencing quality control, and sequenced on an Illumina MiSeq platform using v2, 2×150 paired-end chemistry.

| Bioinformatics
The bioinformatic pipeline was essentially based on the OBITools metabarcoding software suite (Boyer et al., 2016), and can be summarized as follows: first, the two libraries were demultiplexed with "bcl2fastq v. 2.20" (Illumina), then FastQC was used to assess read quality using base scores > Q30, and "illuminapairedend" used to retain paired-end alignments with a quality score > 40. Samples were demultiplexed using "ngsfilter" and sequences length-filtered via "obigrep," to include fragments between 140 and 200 bp. Finally, "obiuniq" was used to dereplicate reads and "uchime" to remove chimaeras. Molecular Operational Taxonomic Units (MOTUs) were clustered using "swarm" (Mahé et al., 2015), with a clustering threshold d = 3 (corresponding to >98% sequence identity), and then assigned using "ecotag" (Boyer et al., 2016), against a custom-made database retrieved from GenBank, generated through in silico PCR using the 12S primers via ecoPCR. The final database was composed of 77,416 12S sequences, which included also other vertebrates.

| Data analysis and mining
General description of taxon diversity and read abundance was provided after using blanks to remove potential contamination noise. Extraction and PCR controls never contained more than 2 reads from any one taxon (Table S2) and therefore showed no impact of contamination on laboratory procedures. Field blanks showed presence of Sphyrna lewini (3 out of 5 blanks) and G. cuvier reads (1 in 5 times), and in those instances, those reads were subtracted from any other sample collected in those sampling dates.
Environmental DNA detections were compared with existing information available for Reunion Island, which included faunal lists in public databases (i.e., www. fishbase.org), scientific surveys (Fricke et al., 2009;Letourneur et al., 2004), and two information sources on fishing activity available for the study period (December 2018 to March 2019): (a) catch records from the local shark control program (Guyomard et al., 2019) and (b) a social media survey of photographs of catches regularly posted by anglers (e.g., https://www.facebook.com/ groups/241858615965673).
We explored αand β-diversity using both presence/ absence and quantitative approaches, in line with several recent empirical studies that suggest a significant quantitative value of eDNA data (Djurhuus et al., 2020;Shelton et al., 2019), especially when the focus is on taxonomically similar and dimensionally comparable species, as it is the case with aquatic vertebrates screened at 12S ribosomal markers (Russo et al., 2020). Site diversity across the sampling period was described using species richness, S, and Shannon's H index. Overall community composition was explored using Principal Component Analysis (PCA) through a covariance matrix of square-roottransformed read abundance data, and non-metric Multi-Dimensional Scaling (nMDS) based on Jaccard distance. The effect of sampling depth and geographic location (considering both a four-group "North/West/South/East" design, see Figure 1, and a simpler "East/West" one, where all North/West/South samples were pooled into one group) was tested using PERMANOVA using both Jaccard and Bray-Curtis dissimilarity measures.
To further investigate the relationship between the two dominant species in the data set, we used a residual randomized permutation procedure (RRPP) with "tiger shark" as a response variable (continuous, fourth-rooted read number) and "scalloped hammerhead" (continuous, fourth-rooted read number) and "time" (categorical, sampling weeks) as predictor variables.

| RESULTS
After bioinformatic filtering and the removal of taxa with sequence identity below 0.7, we obtained 9,159,609 reads; of these, 4,022,414 belonged to the Class Elasmobranchii and 1,838,347 to bony fishes. The remainder of the reads (>3 × 10 6 reads) almost entirely belonged to mammals, 95% of which was human DNA (Table S2). Over 140 species of bony fishes could be identified with this assay, despite it being designed for elasmobranchs ( Figure S5).
A total of 19 elasmobranch taxa were identified, of which 10 were Selachii ("sharks") and nine Batoids ("rays"). Five of these species were almost exclusively detected in the aquarium samples: Carcharhinus melanopterus (65,133 reads), Triaenodon obesus (44,527 reads), Leucoraja erinacea (13,810 reads), Rhinoptera sp. (6,171 reads) and Mustelus sp. (210 reads), which were therefore removed from further analysis. Carcharhinus melanopterus was also detected in one shallow sample at Roches Noires, but that was considered most likely derived from the aquarium. On the other hand, Rhyncobatus sp. (found in two natural samples) and Pteroplatytrygon violacea (found in one sample at Boucan Canot) were considered potentially true detection, as the two species are frequently observed in Reunion Island coastal waters. Naturalistic survey records accumulated for over a century (Fricke et al., 2009;Letourneur et al., 2004) report a total of 64 elasmobranchs in Reunion Island waters, 41 of which are reported in FishBase. Local catches during the study period included 12 species, only two of which (the Arabian smoothhound, Mustelus mosis, and the short-nose spurdog, Squalus megalops) were not detected by eDNA (Table S3).
Overall, elasmobranch sequences were recovered in 71 out of 106 natural samples (67%), with over 91% of sequence reads assigned to three major large predatory sharks: scalloped hammerhead (Sphyrna lewini, 45% of reads), tiger shark (Galeocerdo cuvier, 36% of reads) and bull shark (Carcharhinus leucas, 7% of reads). These were also the most frequently encountered species, being detected in respectively 63, 30 and 19 of the 71 positive samples. Interestingly, we found a strong positive correlation between frequency of occurrence and read abundance of species (Spearman R = 0.91, p = .000013). Scalloped hammerhead, tiger and bull shark reads were also detected in the aquarium samples, despite no individuals of these species being kept in the facilities, suggesting that the DNA of abundant species may also be pumped into the aquarium recirculation system from the outside.
The most diverse sites over the sampling period were Boucan Canot and Roches Noires in the North (species richness, S = 8) and Saint Pierre in the South (S = 7). Most sites appeared dominated by S. lewini, but where this was not the case, G. cuvier was the most abundant species (Figure 1). Indeed, β-diversity across all positive samples was essentially shaped by the relative proportions of scalloped hammerhead and tiger shark (Figure 2a,b), but no major differences were apparent between areas around the island when exploring community structure using presence/absence data (Figure 2c). PERMANOVA analyses carried out to investigate the influence of eDNA sampling depth and location showed that depth had no influence on elasmobranch taxon composition, whether binary or quantitative data were used (Table 1). Sampling area effect was also not significant using binary data but became significant when quantitative data were considered (two-area design: F = 0.9062, p = .028; four-area design: F = 1.527, p = .032).
When we assessed the relationship between G. cuvier and S. lewini in function of time, we found that tiger shark abundance was inversely correlated to both scalloped hammerhead abundance (F = 5.7, p = .02) and time (F = 31.7, p = .001), but most of the variance was explained by time (R 2 = 0.33), rather than the abundance of S. lewini (R 2 = 0.06).

| DISCUSSION
Gathering diversity, abundance and ecological data on elasmobranchs, rapidly and extensively, is broadly viewed as a fundamental step towards their conservation. Our study strengthens the position of eDNA analysis as an agile, accurate and efficient way to monitor diversity patterns in the oceans, particularly in large, elusive and potentially dangerous organisms such as sharks and rays. We characterized the elasmobranch fauna of a tropical oceanic island using a set of primers that offer optimal return for this type of surveys. We detected at least 14 elasmobranch species around the island, while never visually recording any shark or ray during sampling operations for 13 weeks. The number of species detected was greater than what recorded by fishermen during the sampling period, and almost the same as recorded by Sigsgaard et al. (2019) (N = 15) in the Arabian gulf, over a larger area, a greater diversity of habitats, and using one order of magnitude greater sequencing depth.
Elasmobranch diversity in Reunion Island appears comparable with what a previous eDNA metabarcoding study found across the Caribbean (Bakker et al., 2017), albeit using a less efficient COI marker and a considerably lower sequencing depth. Collins et al. (2019) recently indicated that non-specific amplification using a range of COI markers would severely impair eDNA metabarcoding analysis of fish communities, while cautioning about the lower taxonomic resolution of the 12S region. Here we show that the 12S "Miya" region offers, for sharks and rays specifically, the best combination of sequencing depth (millions of reads, and accumulation curve reaching saturation, with 90% of species detected with <70 samples; Figure S4) and assignment success: of the 19 elasmobranchs identified in total (five of which deemed to reflect outflow from the aquarium), 12 were assigned to species and 7 to genus. Elasmobranchs are less speciose than teleosts and are generally characterized by greater inter-specific genetic divergence, providing the grounds for immediate applicability, even when regionspecific, curated sequence libraries are not available (such as in the present study).
Another valuable property of the assay employed in the present study is its remarkable "bonus" of non-target teleost sequences obtained. These allowed the identification of 143 taxa (Table S2), reflecting the assemblage of the studied habitats, with four surgeonfishes (Fam. Acanthuridae), two triggerfishes (Fam. Balistidae), one wrasse (Fam. Labridae) and a couple of small pelagic species among the dominant top 10 species ( Figure S5). Sigsgaard et al. (2019) recovered a remarkably similar level of teleost diversity in the Arabian Gulf: 148 species from a much heftier multi-marker, deep-sequencing effort, which suggests that a single-marker approach and a moderate sequencing effort may be sufficient to capture the main biodiversity features of fish assemblages, even in tropical seas. This also raises the possibility that elasmobranch eDNA surveys across the globe would intrinsically yield valuable additional data for reconstructing and monitoring trophic interactions, especially in diverse tropical reef habitats (Rizzari, Bergseth, & Frisch, 2015;Roff et al., 2016).
Environmental DNA data portray Reunion Island, during the austral summer, as an oceanic hotspot of elasmobranch diversity, whose waters appear dominated, as anticipated, by tiger (G. cuvieri) and bull (Carcharhinus leucas) sharks, but chiefly, and less expectedly, by scalloped hammerhead (S. lewinii), which was present and abundant in all of the 16 locations, and in 89% of the positive samples. These three most frequently recorded species were also the top three in terms of eDNA read abundance, and more generally, a very strong association was observed between frequency of detection and read abundance (R 2 > 0.82) across all species. The extent to which eDNA datasets can be used to extract measures of taxon abundance remains a major discussion point F I G U R E 2 Representations of β-diversity around Reunion Island, in function of sampling depth and location. (a) Principal component analysis (PCA) joint plot on square-root-transformed read abundance data of 71 samples that contained elasmobranch DNA. Principal components PC1 and PC2 together explain respectively 76.5% of the variance in the data matrix. Orange dots are samples collected at <10 m, while turquoise dots refer to samples collected at 40 m. Vectors for only the most important species variables are visualized. The scatter along the vectors of Sphyrna lewini (scalloped hammerhead) and Galeocerdo cuvier (tiger shark) reflect the abundance of these species in the samples. (b) Correlation coefficients associated with PC1 for all the species detected, scalloped hammerhead and tiger shark are coloured in purple to highlight their importance. (c) nMDS plot (stress: 0.21) based on Jaccard similarity index, with data points and convex hulls coloured according to geographic location (following coding in Figure 1). One sample that only recorded the single detection of Loxodon macrorhinus was removed from the visualization. The total number of data points is <71 because some samples contained the same few species and therefore appear completely overlapped in the plot (Deiner et al., 2017), especially in metabarcoding studies. Several authors conservatively opt to use frequency of detection as a proxy for the abundance of a given species across time and space (Boussarie et al., 2018;Postaire et al., 2020); here, we show that occurrence frequency and read abundance are strongly correlated, suggesting that, at least when PCR-amplification is efficient (Kelly, Shelton, & Gallego, 2019), DNA copy abundances can be used to examine relative proportions and temporal variation of target taxa that are taxonomically, dimensionally and physiologically comparable. This may prove momentous in shark ecology and conservation, given the substantial biases associated with the selectivity of all other monitoring methods. Blaison et al. (2015) showed that, although both tiger and bull sharks are present in Reunion Island waters all year round, longline catches and acoustic telemetry detections provide conflicting results regarding their seasonal peak abundances; eDNA might just offer a less burdensome and least-biased approach for the study of shark temporal variation. We found a strong temporal signal in tiger shark, for example ( Figure 3), with the peak of abundance in December mirroring the longline catch maximum in Blaison et al. (2015); on the other hand, the comparatively lower incidence of bull shark during the study period (December to March) appears more consistent with T A B L E 1 Full result table of Permutational Multivariate Analysis of Variance (PERMANOVA) carried out to explore community changes in relation to sampling Depth and Location, as well as their Interaction. Results are reported for both binary (Jaccard's coefficient) and quantitative (Bray-Curtis distance) data and considering two different area designs (see Methods and Figure 1 for details about the spatial groupings)

Jaccard
Bray-Curtis acoustic telemetry data, which indicate peak detections in this species during the austral winter (Blaison et al., 2015). It is hoped that future studies, covering a longer time interval, may provide a more robust portrayal of elasmobranch seasonal variation. The high frequency and abundance of scalloped hammerhead around Reunion Island is remarkable. Until about a decade ago the presence of this species around the isle was still somewhat unconfirmed (Fricke et al., 2009), and S. lewini is still not included in several high profile faunal inventories, such as FishBase. Recently, Guyomard et al. (2019) noted that this species accounted for nearly 9% of the captures by SMART ("Shark Management Alert in Real-Time") drumlines deployed in the West and South-West of the island, but eDNA evidence suggests that this species may be significantly more important in Reunion Island than currently assumed, at least in the period investigated. Sphyrna lewini has recently been assessed as critically endangered (Rigby et al., 2019), with several populations under severe threat from the finning industry (Fields et al., 2020), so the identification of regional hotspots that may play a significant role during the life cycle of this species is of paramount conservation urgency. In this sense, it is worth noting how scalloped hammerhead is demonstrably the species most vulnerable to the deterrents currently trialed in Reunion Island to limit the impact of shark bites: while tiger and bull shark survival rate after capture on SMART drumlines is 94% and 96% respectively, the rate drops to 46% in scalloped hammerhead (Guyomard et al., 2019).
Despite the notable detection of zebra shark (Stegostoma fasciatum)-at Roches Noires once, and on two separate events, in December and March, at Saint Leu-our eDNA analysis highlights a general lack of typical "reef" sharks (e.g., Triaenodon obesus, Carcharhinus melanopterus, C. amblyrhynchos, C. galapagensis), especially in association with the coral reefs on the Western side of the island. This result is consistent with recent drumline surveys (Guyomard et al., 2019), which show that typical reef sharks accounted for just 7% of recorded captures, and among them the majority were tawny nurse sharks (Nebrius ferrugineus). The amplification of DNA from C. melanopterus in this study, and of other reef sharks in other ongoing projects using these primers, indicate that the non-detection of these species is unlikely to stem from failed amplification. Some authors have hypothesized that this may reflect historical overfishing of these species, which created niche availability for the expansion of the more adaptable and generalist tiger and bull sharks (Guyomard et al., 2019), and possibly also contributing to the spate of shark bites observed in the last decade (Lagabrielle et al., 2018). In reality, it is hard to provide a robust reconstruction of events, as elasmobranch populations in Reunion Island have mostly been investigated in recent times, specifically in response to the rapid increase of shark bites, and with a focus on the two main dangerous species (Blaison et al., 2015). Fortunately, we can now profit from eDNA metabarcoding approaches, which can readily generate biodiversity baselines for the purpose of short-/long-term monitoring of community shifts in response to local and/or global environmental processes.
While the value of eDNA metabarcoding approaches continues to be explored in a myriad of environmental management contexts, with strides towards numerous practical applications being made at varying speeds (Djurhuus et al., 2020;Holman et al., 2019;Jacobs-Palmer et al., 2020), the case for upscaling the use of these methods in elasmobranch ecology and conservation is becoming particularly compelling. In the face of major global threats (MacNeil et al., 2020;Queiroz et al., 2019) and the inherent cost, difficulties and drawbacks associated with monitoring shark populations using visual and capture-based techniques, the DNA fragments available in the oceans offer a universal source of qualitative and quantitative spatio-temporal information that can jumpstart the next phase of marine megafauna conservation.

ACKNOWLEDGMENTS
This study was part of projects Direc-ADNe, funded by the DEAL-Réunion (BOP 113), and SeaDNA (grant NE/N005759/1), funded by the UK Natural Environmental Research Council. J. Gadenne, E. Hoarau from Centre Sécurité Requin and T. Poirout from Université de La Réunion contributed to the sampling. P. Bocchiardo and T. Joubert from the Aquarium de La Réunion allowed the collection of water samples in the aquarium tank.

CONFLICT OF INTEREST
The authors have no competing interests.

AUTHOR CONTRIBUTIONS
Sebastien Jaquemet and Stefano Mariani conceived the study, Chloe Fernandez collected samples with Helene Magalon and carried out laboratory analyses with Charles Baillie. Charles Baillie led bioinformatic and ecological analyses with help from Chloe Fernandez and Stefano Mariani. All authors discussed and interpreted the results. Stefano Mariani drafted the manuscript, with contributions from all other authors. Sampling summary table (Table S1), filtered sequence data (Table S2), full summary of species records (Table S3), taxon accumulation curve ( Figure S4) and teleost read abundance histogram ( Figure S5) are available online. The authors are solely responsible for the content and functionality of these materials. Queries (other than absence of the material) should be directed to the corresponding author.

ETHICS STATEMENT
Research protocols followed approved standards of the University of Reunion. No ethical review was required as the data were obtained from water samples.