Environmental DNA captures elasmobranch diversity in a temperate marine ecosystem

Abstract

Many sharks, skates, and rays (elasmobranchs) are highly threatened by the activities of commercial fisheries, and a clear understanding of their distributions, diversity, and abundance can guide protective measures. However, surveying and monitoring elasmobranch species can be highly invasive or resource‐intensive, and utilization of non‐invasive environmental DNA‐based methods may overcome these problems. Here, we studied spatial and seasonal variation in the elasmobranch community of the Western English Channel using environmental DNA (eDNA) collected from surface and bottom waters periodically over an annual cycle (2017–2018). In total we recovered 13 elasmobranch species within eDNA samples, and the number of transformed eDNA reads was positively associated with species (hourly) catch data resolved from 105‐year time series trawl data (1914–2018). These results demonstrate the ability of eDNA to detect and semi‐quantitatively reflect the prevalence of historically dominant and rare elasmobranch species in this region. Notably, eDNA recorded a greater number of species per sampling event than a conventional trawl survey in the same area over the same sampling years (2017–2018). Several threatened species were recovered within the eDNA, including undulate ray, porbeagle shark, and thresher shark. Using eDNA, we found differences in elasmobranch communities among sampling stations and between seasons, but not between sampling depths. Collectively, our results suggest that non‐invasive eDNA‐based methods can be used to study the spatial and seasonal changes in the diversity and abundance of whole elasmobranch communities within temperate shelf habitats. Given the threatened status of many elasmobranchs in human‐impacted marine environments, eDNA analysis is poised to provide key information on their diversity and distributions to inform conservation‐focused monitoring and management.

INTRODUCTION

Globally, elasmobranch diversity is threatened by human activities, and particularly by unsustainable fisheries. They are the most at-risk vertebrate class, after amphibians, with over one-third threatened with extinction (Dulvy et al., 2021; Pacoureau et al., 2021). Specifically, of the 1192 known elasmobranch species (Stein et al., 2018), 64 are considered by the IUCN to be at a high risk of extinction (Critically Endangered or Endangered), with a further 249 considered vulnerable or near threatened (IUCN, 2021). Given the increasing vulnerability of many elasmobranch species to unsustainable harvesting, there is a pressing need for management focused on the conservation of these species (Baum et al., 2003; Birkmanis et al., 2020; Pacoureau et al., 2021).

Effective management of elasmobranch species requires information on their distributions and abundance within ecosystems. Areas with high species richness or unique species compositions are often identified as conservation priority areas (Derrick et al., 2020), but the benefit of establishing marine protected areas (MPAs) for mobile or wide-ranging elasmobranch species can be unclear if their abundance cannot be reliably monitored. Non-invasive survey methods capable of capturing multispecies abundance data such as underwater visual censuses and baited remote underwater video (Juhel et al., 2018), as well as satellite imagery (Williamson et al., 2019), may overlook components of the full elasmobranch assemblage. Moreover, capture-based fisheries survey methods can have strong intrinsic vulnerability biases to survey gears (Young et al., 2019), and their use as survey methods are highly impactful and therefore unsuitable for protected areas. Meanwhile, passive acoustic monitoring (e.g., Rider et al., 2021) and satellite tracking of tagged individuals (e.g., Queiroz et al., 2019; Vedor et al., 2021) are capable of providing extremely detailed information on space use but, typically, they can only be used for a small number of individuals within a population. In the case of passive acoustic monitoring, it is only possible within the spatial context of a proximate hydrophone (acoustic receiver) array (Lea et al., 2016).

A potentially efficient and non-invasive method of capturing information on the whole community of elasmobranchs is to study their environmental DNA (eDNA). In the case of marine fishes, eDNA will be derived from multiple sources, including feces, urine, gametes, mucous or decomposing tissues. There are two main approaches employed to study the eDNA of fishes. The first approach is to design species-specific quantitative PCR (qPCR) assays for target DNA fragments, which enable the number of copies of that fragment in the eDNA sample to be determined. This can be an effective tool for studies focused on a small number of species, but it would be challenging to design complementary assays for a whole marine species assemblage. The second approach is to use metabarcoding, where “barcode” regions of multiple species are simultaneously PCR amplified from eDNA templates, sequenced using high-throughput technologies, and resultant sequences assigned to species using reference databases (Andruszkiewicz et al., 2017; Bista et al., 2017; Bohmann et al., 2014; Deiner et al., 2017). This approach has the advantage of being able to characterize the composition of whole communities, but the precise number of copies of target DNA of each individual species within the eDNA samples is not directly assayed, in part because of differences in amplification efficiencies of primers on different target DNA templates. Hence, metabarcoding methods are broadly considered to be only semi-quantitative approaches to assay target eDNA (Blabolil et al., 2021).

Several studies suggest that the quantity of the eDNA in marine environments, measured either through target DNA copy number (qPCR) or target read number (metabarcoding), can generally reflect the abundance of source individuals in the environment (Rourke et al., 2021; Salter et al., 2019; Sato et al., 2021; Stoeckle et al., 2021). This relationship is often weak, however, due to multiple contributing factors (Lamb et al., 2019). For example, there is evidence that fish breeding behavior can elevate the amount of eDNA in water (Bylemans et al., 2017), and that the direction and strength of currents will influence detectability at the site of production (Andruszkiewicz et al., 2019). Nevertheless, given the relatively rapid rate of eDNA degradation (Holman, Chng, et al., 2021, Holman, de Bruyn, et al., 2021), with half-life of eDNA in marine systems ranging from 18.2 to 71.1 h in seawater (Collins et al., 2018), there is confidence that the locations and times where eDNA is detected is likely to reflect the occurrence of species, and that eDNA quantities measured using qPCR or metabarcode read numbers can at least be partially indicative of fish abundance. Consequently, there is growing advocacy for the use of marine environmental DNA-based methods for routine monitoring, and for the derived data to inform marine management and policy decisions (Gilbey et al., 2021).

To date, research on marine elasmobranch communities using eDNA has primarily focused on subtropical or tropical environments (Bakker et al., 2017; Boussarie et al., 2018; Lafferty et al., 2018; Mariani et al., 2021; West et al., 2020). The ability of eDNA-based methods to detect and monitor abundance of elasmobranch communities in temperate waters has received less attention (see Weltz et al., 2017 for an example). This is notable because many threatened elasmobranch species are present in heavily fished temperate waters. In northeast Atlantic shelf seas these include the blue skate Dipturus batis, tope Galeorhinus galeus, angelshark Squatina, undulate ray Raja undulata, spiny dogfish Squalus acanthias, and porbeagle shark Lamna nasus (Heessen et al., 2015; Lawson et al., 2020). Therefore, eDNA-based methods could offer a much-needed enhancement of spatial and temporal monitoring practices for elasmobranch communities in temperate seas.

In this study, we used environmental DNA metabarcoding to describe the species composition of the elasmobranch community of the Western English Channel, near Plymouth, United Kingdom. The location was selected because the marine fish community of the region has been well characterized by over a century of survey trawls conducted by the Marine Biological Association of the UK, giving us knowledge of the rare species in the assemblage, and enabling us to test the ability of eDNA-based methods to detect them. By periodically sampling eDNA from surface and bottom waters at three stations over a full annual cycle, we were able to determine the diversity of species present, and the extent of spatial and temporal variation in the assemblage. The results are discussed with reference to the biology of focal species, and the potential for eDNA-based surveys to inform conservation-focused marine management.

MATERIALS AND METHODS Sampling

Environmental DNA samples were collected using Niskin bottles from three stations in the Western English Channel (L4, L5, E1; Figure 1), from both the surface and bottom, approximately once per month from February 2017 to April 2018. Surface water was collected at <1 m from the surface, to avoid detritus directly on the surface, while benthic water was collected close to the substrate without disturbing the sediment. Water was first strained through a 250 μm nylon mesh to remove large plankton and debris, before being transferred to Nalgene HDPE collection bottles pre-sterilized with a 10% bleach solution, and the samples were then placed on ice. During each sampling event, triplicate 2 L samples of seawater were collected at each depth. A complete list of samples collected is presented in Supporting Information Table S1.

FIGURE 1. Sampling stations in the English Channel. eDNA samples were taken from each station (L4, L5, and E1) in surface waters and bottom waters (L4, 50.25°N 4.22°W, 51–56 m depth; L5, 50.18°N 4.30°W, 60–65 m depth; E1, 50.03°N 4.37°W, 70–72 m depth). The yellow box indicates the distribution of the long-term survey trawls conducted by the Marine Biological Association of the UK (50.13–50.30°N and 3.92–4.65°W).

Within 5 h of collection, each 2 L sample was filtered through an 0.22 μm Sterivex-GP PES filter (SVGP01050; Merck Millipore) using a peristaltic pump. Across approximately 10% of events, 2 L of distilled water was taken into the field as a field negative control, and otherwise treated identically to the seawater samples.

DNA extraction, PCR amplification, and Illumina sequencing

DNA was extracted from filters using the DNeasy PowerWater Kit (Qiagen). Using PCR we amplified a ~182 bp fragment of the mitochondrial 12S gene using the Elas02 elasmobranch primer pair (Taberlet et al., 2018), which are modified versions of the MiFish-U/E primers (Miya et al., 2015). A list of the PCR efficiencies of these primers against target sequences of species in the regional elasmobranch community, as determined by the decipher v2.22.0 package (Wright, 2016) is provided (Supporting Information Table S2). These PCR primers were adapted with unique 8-mer sample-identifying barcode tags identical on both the forward and reverse primer and incorporating 2–4 random 5` bases to increase sequencing heterogeneity. A total of eight PCRs were performed on each extracted eDNA template. Each PCR was conducted in a 20 μl volume comprising: 10 μl AmpliTaq Gold 360 Master Mix (4398876; Applied Biosystems); 0.16 μl Bovine Serum Albumin (B14; ThermoFisher); 1 μl forward primer (5 μM); 1 μl reverse primer (5 μM); 5.84 μl molecular grade water; and 2 μl eDNA template. Thermocycling parameters comprised: polymerase activation at 95°C for 10 min; 40 cycles of 95°C for 30 s, 59°C for 30 s, 72°C for 60 s; and a final extension of 72°C for 7 min. Alongside the extracted 209 samples (Supporting Information Table S1) we included seven filtration negative controls, eight extraction negative controls, 18 negative no-template PCR controls, and three positive PCR controls using 0.04–9.1 ng genomic DNA extracted from two non-UK species (fin tissues, spinner shark Carcharhinus brevipinna, and rough ray Raja radula) (Supporting Information Table S3). The eDNA extractions, pre-PCR preparations, and post-PCR procedures were carried out in separate rooms.

PCR products were checked by gel, and then pooled and purified using the MinElute PCR Purification Kit (28004; Qiagen) following the manufacturer's protocol. Illumina sequencing adapters were attached to the amplicons using the NEXTflex PCR-Free kit (5142-01; PerkinElmer) following the manufacturer's protocol. A total of 3 libraries using unique indexes were created. Libraries were then quantified using a NEBNext (E7630S; New England Biolabs) qPCR assay and sequenced on an Illumina NextSeq using v2 (2 × 150 bp paired-end) chemistry and 10% phiX spike-in.

Bioinformatic analyses

Raw sequencing reads were processed using the meta-fish-pipe v1.0 bioinformatics module (https://doi.org/10.5281/zenodo.5083336), following Collins et al. (2019). In brief, the following steps were carried out: (i) reorientation and demultiplexing of reads using cutadapt v3.4 (Martin, 2011); (ii) denoising, merging, removal of chimeric reads, and dereplication using dada2 v1.20 (Callahan et al., 2016); (iii) homology filtering using hmmer v3.1 (Eddy, 1998); (iv) first pass taxonomic assignment using sintax (Edgar, 2016) and NCBI RefSeq v206 reference library obtained using refseq-reflib v1.0 (https://doi.org/10.5281/zenodo.5083346); and (v) exhaustive taxonomic assignment using the meta-fish-lib v243 custom UK fish reference library (Collins et al., 2021), blastn v2.11.0 (https://blast.ncbi.nlm.nih.gov), and EPA-ng v0.3.8 (Barbera et al., 2019). Resulting amplicon sequence variants (ASVs) were cross-referenced with those from concurrent lab projects to control for laboratory contamination, and all non-elasmobranch species were removed for downstream statistical analyses.

Trawl surveys

Our eDNA collections took place in a marine region where the demersal fish assemblage has been intensively surveyed by the Marine Biological Association of the UK (MBA) since 1911. Hence, we used these quantitative survey data to describe the known elasmobranch community of the region, enabling us to determine if the prevalence of species in the eDNA sample we collected was broadly reflective of the known prevalence of species in the historical trawl survey data. Because of the vastly differing timescales of the eDNA and trawl surveys, we do not use the trawl survey data to translate eDNA read abundance into metrics of actual elasmobranch abundance in the region.

The MBA trawl survey took place between 50.13 to 50.30°N and 3.92 to 4.65°W during the years 1911, 1913–1914, 1919–1922, 1950–1958, 1967–1979, 1983–1994, 2001–2010, and 2016–2018. For much of this period the focal location has been the L4 sampling site. In total eight vessels have been used for sampling (1911–1919 SS Oithona, 1920–1922 RV Salpa, 1950–1952 RV Sabella, 1952–1973 RV Sula, 1974–2003 RV Squilla, 1979 RV Sarsia, 2004–2015 RV Plymouth Quest, 2015–2018 RV Sepia). Records suggest that the survey trawls have been broadly comparable throughout the series, being conducted at the same speed (ca. 4 knots), with gear of similar dimensions (headline length range, 16.2–19.8 m; groundrope length range, 19.8–27.4 m; main net stretched mesh diameter, 75–270 mm, and all vessels used a fine-mesh cod end or a cod-end cover). We conducted our analyses on data from the survey years where elasmobranch individuals were reliably identified to species level, namely 1914, 1919–1922, 1953–1958, 1976–1979, 1983–1994, 2001–2018. Within this subset of data (1037 trawls), the average number of hauls during sampling years has been 23 (range 1–45), and the average duration of each haul has been 49 min (range 14–180 min). Records of smooth-hound Mustelus mustelus in the trawl data were considered to be starry smooth-hound Mustelus asterias, given genetic analyses suggesting all Mustelus in this region are starry smooth-hound (Farrell et al., 2009). Average catch per unit effort (CPUE) of each species was calculated as the average catch in numbers of each species per hour trawling. Using these same data, the frequency of occurrence of species across all 1037 hauls was also calculated (% of hauls in which at least one individual of the species was recorded as present).

Analyses of eDNA data

Community-level analyses of eDNA samples were conducted in R v3.6.2 (R Core Team, 2019). We removed all samples where elasmobranchs were absent, leaving 174 of the 209 samples, and two dataframes were generated using data transformations previously proposed as appropriate for eDNA-derived metabarcode data. First, we generated a matrix comprising Wisconsin transformed data (following Kelly et al., 2019), using the “wisconsin” function in vegan v2.5.7 (Oksanen et al., 2020). Second, we generated a matrix comprising Hellinger transformed data (following Laporte et al., 2021), using the “hellinger” function in vegan.

To test for differences among stations, between sampling depth and between sampling months, we used PERMANOVA with the “adonis2” function in vegan, with 10,000 permutations. To test for post-hoc differences between sample groups, we used the “pairwise.adonis” function, with 100,000 permutations (Martinez Arbizu, 2020). To identify species associated with statistically significant differences among sampling stations and depths we used the multilevel patten analysis “multipatt” function in indicspecies v1.7.9 (De Cáceres & Legendre, 2009) on the transformed data. To ordinate differences among samples, we used Principal Coordinates Analysis implemented with the “pcoa” function in ape v5.0 (Paradis & Schliep, 2019) in each of the two matrixes, using the resulting primary axes of variation to visualize differences among sampling stations and between sampling months.

To explore seasonal changes in abundance at the species level, we calculated the Hellinger standardized number of reads per sampling month across stations, and generated a heatmap of abundance. Then, again using Hellinger standardized data, we quantified seasonal variation by fitting generalized additive models (GAMs) to the data for the six species that were most abundant in the eDNA metabarcode reads (small-spotted catshark Scyliorhinus canicula, starry smoothhound M. asterias, small-eyed ray Raja microocellata, thornback ray Raja clavata, spotted ray Raja montagui and blonde ray Raja brachyura). Models were generated using mgcv v1.8.33 (Wood, 2011), using the following predictors: the smooth factor sampling month (k = 5), and fixed factors of sampling year, sampling depth, and sampling station. Response data were assumed to have a negative-binomial distribution, and models were fitted using the REML smoothing parameter estimation method.

Comparison of eDNA data to trawl data

To compare the total eDNA read abundance across the 209 samples, and the composition of the elasmobranch communities using the trawl surveys (average CPUE across the 1037 hauls, and frequency of occurrence across the 1037 hauls), we used linear regression employing the lm function in base R (R Core Team, 2019). For these analyses we 4^th root transformed the total eDNA read abundance and 4^th root transformed the average CPUE data. To compare the numbers of species encountered as a function of sampling effort (number of sampling events, either trawl or eDNA) we calculated sample-level species accumulation curves using the “specaccum” function in vegan. We undertook this analysis for all survey trawls across all time periods where individuals in the trawl were identified to species level, and for the 22 survey trawls that took place during the eDNA sampling period (February 2017 to April 2018). Although these 22 survey trawls temporally overlapped with eDNA survey period, they were not conducted on the same days, and since only three elasmobranch species were caught across 22 survey trawls, they were not suitable for use to undertake analyses attempting to calibrate eDNA metabarcode read number against abundance measurements. Finally, we compared the species richness resolved through eDNA and the survey trawls to the diversity to all species of elasmobranchs encountered in proximity to Plymouth (Start Point in Devon to Looe in Cornwall, southward to the outer Channel grounds) using the Plymouth Marine Fauna (Marine Biological Association, 1957).

RESULTS eDNA sequencing statistics

In total, 209 samples from the three locations (L4: 66 samples, E1: 71 samples, L5: 72 samples) were collected during 2017–2018. A total of 161,183,652 raw sequencing reads were generated for the three libraries across the samples and controls. After quality-control filtering and taxonomic assignment, 58,684,923 reads were remaining (Supporting Information Table S4). Then, following removal of non-elasmobranch and control species from samples, a total of 38,615,907 reads were assigned to native elasmobranch species within our 209 samples (Supporting Information Table S5).

Community composition

In total 13 species were recovered in the eDNA samples from 2017–2018, and included seven sharks and six skate species (Figure 2a). Shark species were thresher shark (Alopias vulpinus), tope, small-spotted catshark, nursehound (Scyliorhinus stellaris), starry smooth-hound, spiny dogfish, and porbeagle shark. Skate species recovered were spotted ray, thornback ray, small-eyed ray, blonde ray, cuckoo ray (Leucoraja naevus), and undulate ray (Figure 2a). By contrast only three elasmobranch species were encountered in survey trawls that took place during the 2017–2018 sampling eDNA sampling period, specifically small-spotted catshark, spotted ray, and thornback ray.

FIGURE 2. (a) Total read abundance of the 13 species of elasmobranch identified across all 209 field samples. (b) Association between the number of eDNA reads in samples (y-axis, 4th root transformed) and the catch per unit effort of 16 taxa recovered in survey hauls (1911–2018). (c) Association between the number of eDNA reads in samples (y-axis, 4th root transformed) and the frequency of occurrence of 16 taxa recovered in survey hauls (1911–2018)

In total, in the 1037 survey trawls of the region between 1914 and 2018 that we analyzed, 14 taxa have been recorded, including 12 of the 13 species recovered in the eDNA. The only species present in the eDNA but absent in the trawl survey was thresher shark. Meanwhile, the taxa absent from eDNA, but present in the long-term trawl time series were shagreen ray (Leucoraja fullonica) and angelshark. These two taxa were extremely uncommon in the trawl survey: shagreen ray (captured in three survey hauls in 1921) and angelshark (captured in two survey trawls, one in 1921 and one in 1957). Across species, the catch per unit effort in trawl surveys of the region between 1914 and 2018 was significantly positively associated with the numbers of eDNA reads of species recovered in the 2017–2018 sampling across the 209 eDNA samples (linear model, both CPUE and total read abundance 4^th root transformed, n = 15, F_1,13 = 7.258, p = 0.035, r² = 0.298, Figure 2b). Frequency of occurrence of species within 1914 and 2018 trawl surveys was not significantly associated with the numbers of eDNA reads assigned to those species (n = 15, total read abundance 4^th root transformed, F_1,13 = 4.493, p = 0.054 r² = 0.257, Figure 2c).

Spatial and temporal patterns of community structure

Overall, there were significant differences in elasmobranch community structure among sampling stations, and between sampling depths (Figure 3; Table 1). In post-hoc comparisons, the spatial differences were most striking between the inshore site L4 and the two sites further offshore (L5 and E1), while there was no evidence of any spatial differentiation between L5 and E1 (Figure 3; Table 2). From indicator analyses, we found that the only significant differences among sites corresponded only to a greater abundance of small-eyed ray and nursehound at L4 relative to L5 and E1. (Supporting Information Table S6). Indicator analyses provided no evidence of significant abundance differences between the sampling depths (Supporting Information Table S6).

FIGURE 3. Differences in elasmobranch community structure between the three sampling stations as resolved from environmental DNA reads. (a) Wisconsin-transformed data, (b) Hellinger-transformed data. Each point is one sample

TABLE 1 Statistical significance of differences in elasmobranch community structure from eDNA reads among sampling stations, between sampling depths, and among sampling months, as resolved using PERMANOVA

Data	Factor	Df	SS	r2	F	p
Wisconsin transformed	Station	2	1.308	0.0193	2.264	0.006
	Depth	1	0.622	0.009	2.152	0.035
	Month	9	12.439	0.183	4.784	<0.001
	Station * Depth	2	0.403	0.006	0.698	0.809
	Station * Month	15	10.465	0.154	2.415	<0.001
	Depth * Month	9	2.853	0.042	1.097	0.278
	Station * Depth * Month	14	4.746	0.070	1.173	0.115
	Residual	121	34.961	0.512
	Total	173	67.797	1
Hellinger transformed	Station	2	1.002	0.017	1.993	0.024
	Depth	1	0.541	0.009	2.154	0.042
	Month	9	11.817	0.197	5.226	<0.001
	Station * Depth	2	0.344	0.006	0.684	0.799
	Station * Month	15	8.944	0.149	2.373	<0.001
	Depth * Month	9	2.614	0.044	1.156	0.196
	Station * Depth * Month	14	4.199	0.070	1.194	0.107
	Residual	121	30.403	0.508
	Total	173	59.864	1

Note:

Bold indicates p < 0.05.

TABLE 2 Statistical significance of differences in elasmobranch community structure from eDNA reads between pairs of sampling stations, using pairwise PERMANOVA

Dataset	Site	F	r2	p
Wisconsin	L4 vs. E1	2.496	0.021	0.013
	L4 vs. L5	2.120	0.018	0.033
	E1 vs. L5	0.309	0.003	0.973
Hellinger	L4 vs. E1	2.201	0.019	0.036
	L4 vs. L5	1.938	0.016	0.064
	E1 vs. L5	0.152	0.001	0.998

Note:

Bold indicates p < 0.05.

We found highly significant differences in community structure among sampling months (Table 1; Figure 5a). Focusing on the primary axes of variation in the Principal Coordinates Analysis (PCoA-1 or PCoA-2), a seasonal cyclic pattern was present irrespective of the data transformation method (Figure 4). Exploring temporal changes in eDNA abundance using GAM models indicated strong differences among species in seasonal read abundance, and evidence of significant seasonal (between month) variation was present in four of the six species we considered (Supporting Information Table S7). Specifically, small-spotted catshark showed read abundance peaks in March-May, thornback ray peaked in September-November, and small-eyed ray and starry smooth-hound peaked in November-April (Figure 5b–g). We found no significant seasonal (between month) variation in blonde ray or spotted ray. Temporal differences were also apparent in reads of rarer species, with porbeagle shark, thresher shark, and tope shark present in July and November, while spiny dogfish was present between October and April (Figure 5a).

FIGURE 4. Seasonal shifts in elasmobranch community structure, summarised using Principal Coordinates Analysis scores. Each point is one sample arranged by sampling month, for 2017 (dark blue) and 2018 (light blue). Illustrated are the scores from the primary axes of variation PCoA-1(a–b) and PCoA-2 (c–d), for each data transformation Wisconsin standardized (a,c), Hellinger standardized (b,d). Superimposed are 3rd order polynomial curves with the shaded area illustrating one standard error. Variance captured by each PCoA axis is reported in Figure 3

FIGURE 5. (a) Reads per sample per sampled month (Hellinger transformed) between February 2017 (2017_02) and April 2018 (2018_04). (b–g) Generalized additive model plots illustrating associations of Hellinger transformed species-level read abundance (standardized, y-axis) in relation to sampling month (x-axis)

DISCUSSION

This study provides strong evidence of the ability of eDNA-based methods to generate information on the structure of both pelagic and demersal components of a temperature marine elasmobranch assemblage. In total, eDNA samples collected monthly over one year detected nearly all species recovered in demersal trawl surveys over a century-long scale. The exceptions were species that are typically rare in the region, and have only been caught by trawl surveys on a small number of occasions. We did not aim to calibrate eDNA as a tool for the fully quantitative assessment of elasmobranch communities. Such validation would require a very substantial amount of trawl effort, as the 22 trawls undertaken over the same timescale as the eDNA sampling caught individuals of only three species. Nevertheless, it was notable that the total number of reads we recovered of species was significantly positively associated with their average catch per unit effort in the region derived from the full trawl survey data collected over the previous century. Thus, our results support the concept that eDNA metabarcoding can provide semi-quantitative information pertaining to dominant and rare species that may help to map distributions and primary habitat of elasmobranchs across marine regions. This is of importance, as current knowledge of elasmobranch assemblages is often based on visual or capture-based survey methods that can be strongly biased toward the species that are more abundant, more easily captured, and/or less cryptic (Boussarie et al., 2018).

Spatial structure

Spatial differences in the abundance of elasmobranchs can be explained by differences in core ecological niches of the species (Humphries et al., 2016), as well as vulnerability to local fisheries (Brander, 1981). We found evidence of significant differences in the spatial distribution of elasmobranch species, over a spatial scale of 30 km. Contrasts were most apparent between the inshore site (L4), relative to offshore sites (L5, E1), with the inshore sites being characterized by a greater abundance of small-eyed ray and nursehound in particular. This is supportive of eDNA being capable of resolving some of the fine-scale differences in depth and substrate preferences of the UK elasmobranch fauna, that have previously been resolved through trawl surveys (Kaiser et al., 2004) and tracking of tagged individuals (e.g., Kaiser et al., 2004; Humphries et al., 2016; Simpson et al., 2021).

In marine systems, eDNA has been able to resolve spatial differences in fish community structure over small (<5 km; Jeunen et al., 2019; West et al., 2020), moderate (20-100 km; Lafferty et al., 2018; Mariani et al., 2021), and broader spatial scales (>100 km; Sigsgaard et al., 2020; Fraija-Fernández et al., 2020; Holman, Chng, et al., 2021, Holman, de Bruyn, et al., 2021; Valdivia-Carrillo et al., 2021; West et al., 2021). However, the extent of resolution is likely to depend partly on the degree of eDNA transport, which in turn depends on both the rates of eDNA persistence and the amount of horizontal advection of the eDNA (Andruszkiewicz et al., 2019). Therefore, the observed heterogeneity in spatial patterns of eDNA abundance may be related to oceanographic differences among the three sampling locations, as well as the relative rates of eDNA persistence at locations. Notably, all locations have similar oceanographic properties, including seasonal stratification and seasonal nutrient profiles (Smyth et al., 2010). Nevertheless, there are some differences linked to the proximity of L4 to the coast, most notably seasonal surface freshening, linked to freshwater input (Smyth et al., 2010). Thus, L4 may be more heavily influenced by eDNA from proximate shallow water inshore sites than the further offshore L5 and E1, potentially explaining some differences in eDNA composition. Experimental work around Plymouth has indicated eDNA to be detectable for around 48 h (Collins et al., 2018); however, the rate of decay was 1.6× faster at inshore waters (Sutton Harbour in Plymouth Sound) than the offshore waters (E1). Therefore, offshore sites may be more homogenous, perhaps due to longer eDNA persistence providing more opportunity for mixing.

Differences in marine communities resolved using eDNA-based methods over depth gradients have been reported, for example, Jeunen et al. (2020) who studied a depth-temperature-salinity gradient in a New Zealand fjord, and Canals et al. (2021) who studied a 2000 m open ocean depth gradient, encompassing both the epipelagic and mesopelagic zones. By contrast we found no clear-cut differences between elasmobranch communities resolved using eDNA from the surface and bottom waters over a distance of 50 m, which may be related to the mobility of the focal species. Several of the shark species are pelagic and therefore have the capability to move rapidly between surface and bottom waters and indeed have been tracked doing so in the southwest UK region (e.g., porbeagle; Pade et al., 2009). Moreover, benthic species, including the thornback ray, blonde ray, spotted ray, and small-eyed ray, undertake diel vertical migrations from deep benthic habitat during the day to shallow benthic habitat during the night (Humphries et al., 2017). However, absence of any clear differences between the surface and bottom waters may also simply reflect intrinsic mixing of waters in the sampling region. Irrespective of the causes, our results are notable as they suggest reliable information on shelf sea pelagic, benthopelagic, and benthic faunas may be sourced from eDNA collected from surface waters of temperate shelf seas, even when those species do not have pelagic larval dispersal phases.

Temporal structure

Seasonal shifts in fish community structure have been reported using marine and estuarine eDNA (Djurhuus et al., 2020; Stoeckle et al., 2017, 2021). Similarly, our data were characterized by temporal shifts in read abundance, likely linked to seasonal differences in habitat use. Seasonal occurrence of pelagic shark species in our data can be explained by seasonal migration. Porbeagle shark were present in our data in October, and Biais et al. (2017) report movements of satellite-tagged porbeagle northwards into UK waters during warmer summer-autumn months, before moving further north and west in autumn-winter and returning to southerly waters of the Iberian Peninsula during the coldest winter periods. Thresher shark were recorded in the eDNA in October, and these are typically recorded in UK waters in summer months (Stevens, 1976). Populations on the western Atlantic undertake north-south migrations, being further north in summer-autumn, and returning south in winter-spring (Kneebone et al., 2020). Spiny dogfish were recorded in eDNA in November and February. This is a species known to spend winter in the Western English Channel and move northwards into the more northerly European shelf waters during summer (Vince, 1991). We also found the starry smooth-hound to be most abundant in the winter months. Tagging studies have shown this is a highly migratory species, and consistent with our results, it has a general pattern of overwintering in the English Channel and Bay of Biscay, where pupping takes place, before spending summer months in the North Sea (Brevé et al., 2016, 2020; Griffiths et al., 2020).

Benthic shark and skate species tended to be more consistently present in the data throughout the year, but with some seasonal peaks in abundance. For these species, it is less clear if migrations can explain the variation observed, as although seasonal migrations in benthic elasmobranchs are known (e.g., thornback ray; Hunter et al., 2006), most tagged individuals exhibit strong philopatry, being recaptured within 50 km of the immediate release site (Bird et al., 2020; Rodrıguez-Cabello et al., 2004; Simpson et al., 2020). Thus, if seasonal migrations do explain the variation, then they must be reflecting more modest within-region shifts in habitat use, perhaps related to reproduction. This explanation is supported by the movement patterns of species and within-species sexes recorded by long-term acoustic tracking of skates in the Western English Channel off Plymouth (Simpson et al., 2021). Equally, it is possible that eDNA abundance is linked to activity levels (de Souza et al., 2016; Thalinger et al., 2021), perhaps associated with reproduction. For example, we found a peak in abundance of the small-spotted catshark during late spring (March-May), corresponding to peak breeding season in Plymouth in April (Sumpter & Dodd, 1979). The spotted ray, small-eyed ray, and thornback ray all showed peak eDNA abundances in late winter and spring, during which time these species have been recorded to start egg laying in northern European waters (Clark, 1922; Holden, 1975; Koop, 2005).

Detection of threatened species

Of the 13 species recovered within eDNA reads, several are rare species of particular conservation concern. Tope shark is listed as Critically Endangered by the IUCN, and was recovered in one eDNA sample (total 4404 reads). This is a benthopelagic species that is widespread in the Eastern Atlantic, undertakes long migrations throughout the region (Holden & Horrod, 1979), but is rarely caught in trawl surveys (Heessen et al., 2015). In MBA survey trawls they have been encountered in seven sampling events between 1921 and 2005, consistent with occasional presence in the region. The undulate ray was also present in multiple eDNA samples yet was represented by a relatively low number of reads (total 34,254 reads) in comparison to other skate species. The undulate ray is IUCN listed as Endangered, and in northern European waters is abundant in parts of the English Channel away from Plymouth (e.g., Jersey), at depths less than 100m (Heessen et al., 2015). In Plymouth, the species was historically “not uncommon” at depths of “20m or more S. of Eddystone” (Marine Biological Association, 1957), but has only rarely been caught in MBA survey hauls, consistent with the species having been sporadically present, but also suggesting our sampled region is not core habitat for the species.

Notably, the eDNA analyses failed to record two species that historically have been encountered in MBA survey trawls, but not during recent decades. The first species is the shagreen ray, which is most abundant in deep water trawls >70 m (Heessen et al., 2015), so may be expected to be absent from our sampling locations. The second species is the angelshark, listed as Critically Endangered by the IUCN. This is a species that was formerly common in Plymouth waters (Marine Biological Association, 1957), but last sampled in an MBA trawl survey in 1957. The absence of angelshark in our eDNA samples is compatible with records suggesting the species is now extirpated from the English Channel region, with the nearest extant population in inshore waters at Cardigan Bay, West Wales (Ellis et al., 2020; Hiddink et al., 2019).

Collectively, our results show that eDNA metabarcoding can reliably capture the diversity of the proximate elasmobranch assemblage. By contrast trawl surveys undertaken over the same timescale were only able to capture a small number of species, most likely due to the intrinsic rarity of most elasmobranch species in the survey area. Trawl surveys are the most commonly used – and destructive – method of surveying fish assemblages in temperate seas. We have shown that eDNA-based methods have potential to reveal part of the assemblage that would otherwise be unrepresented within contemporary surveys (Figure 6). However, there are additional species that have historically been recorded within Plymouth waters that were not encountered in eDNA or survey trawls, and these species remain as “dark diversity” (following Boussarie et al., 2018), In the Western English Channel dark diversity would include pelagic species such as basking shark Cetorhinus maximus and blue shark Prionace glauca, as well as benthic feeding species that have been occasionally recorded, such as common eagle ray Myliobatis aquila, common stingray Dasyatis pastinaca, and blue skate Dipturus batis (Figure 6). Although most plausibly the absence of these species from our results is that they are either absent or intrinsically rare in the year we sampled (e.g. basking shark were common off Plymouth between 1995 and 2006; Sims, 2008), there remains the additional possibility that their ability to be detected may in part be related to volumes of water sampled, number of PCR replicates, depth of sequencing, or technical aspects of the assay such as primer efficiency and bioinformatic filtering (Díaz-Ferguson and Moyer, 2014, Pilliod et al., 2014).

FIGURE 6. The elasmobranch community of the Western English Channel near Plymouth. Shown are a rarefied sampling curves for species found in the 2017–2018 eDNA survey (orange, total 13 species, 209 sampling events), for species caught in the trawl survey from 1914–2018 (green, total 14 species, 1037 sampling events), and for species caught in the trawl survey only during 2017–2018 (blue, total 3 species, 22 sampling events). Orange circles indicate presence in eDNA, green circles indicate presence in survey trawls. Also shown is the “dark diversity” (sensu Boussarie et al., 2018), which are species that have been recorded in the Western English Channel near Plymouth (Marine Biological Association, 1957), but were not recovered in either trawl records analyzed or eDNA (14 species). Images from FAO, and are not to scale

To conclude, this study has shown the ability of eDNA to illuminate the species richness of a temperate elasmobranch community, and its spatial and temporal structure, which fits with expectations from habitat features and species life histories. The results suggest that eDNA could be used for mapping and routine monitoring of elasmobranch assemblages, enabling semi-quantitative assessments of the effectiveness of marine management objectives. Further refinement of methodological aspects, especially pertaining to eDNA transport and the associations between eDNA metabarcoding read number and organismal abundance, will play a major part in facilitating the transition of eDNA monitoring from emerging tool to established practice in marine science.

ACKNOWLEDGEMENTS

We thank the captain and crew of the RV MBA Sepia at the Marine Biological Association for assistance collecting eDNA samples, and the numerous researchers and crew who have collected Marine Biological Association long-term survey data.

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.

AUTHOR CONTRIBUTIONS

Conceptualization. MJG, SM, DWS, ZL. Acquisition of field samples and field data. SR, RB, AMG. Acquisition of laboratory data. ZL, RAC, CB. Analysis of data. ZL, RAC, MJG. Writing of manuscript. ZL, MJG, RAC.

DATA AVAILABILITY STATEMENT

All code to reproduce the analyses in this study can be obtained from https://doi.org/10.5281/zenodo.5497562. The raw sequencing reads generated during this study are available from the Sequence Read Archive (BioProject PRJNA808852).