Genomics versus mtDNA for resolving stock structure in the silky shark (Carcharhinus falciformis)

Introduction

Many elasmobranchs around the globe have experienced devastating population declines due to overfishing in both target and non-target fisheries (Musick et al., 2000; Clarke et al., 2006; Ferretti et al., 2010; Heupel et al., 2014; Dulvy et al., 2014; Oliver et al., 2015; Dulvy & Trebilco, 2018). These species are especially vulnerable to overfishing due to life history traits such as late maturity, slow growth, low fecundity, and high juvenile mortality, which collectively result in low intrinsic rate of population increase (Baum et al., 2003; Dulvy et al., 2008). Elasmobranch populations take decades to recover from overfishing, and only if fishing pressure is relieved for an extended period (Stevens et al., 2000). Furthermore, many threatened and endangered elasmobranchs have little to no population genetic data that would assist in the resolution of management units (reviewed in Domingues, Hilsdorf & Gadig (2018)).

Genetically distinct populations are isolated management units known as stocks; however, stocks can be defined on a smaller scale than genetic populations through other criteria, such as an exclusive economic zone boundry (Carvalho & Hauser, 1994; Ovenden et al., 2015). Reduced gene flow indicates that if a population is overfished it will not be replenished by immigrants from surrounding populations. This is why managing on a genetic stock-by-stock basis is essential for successful maintenance of exploited species and is sorely needed for over-harvested elasmobranchs (Dizon et al., 1993; Heist, 2004; Tallmon et al., 2010).

For the past two decades the standard for examining population structure in elasmobranchs has been a section of the mitochondrial genome, usually the control region (mtCR) (Duncan et al., 2006; Hoelzel et al., 2006; Keeney & Heist, 2006; Castro et al., 2007; Whitney et al., 2012; Clarke et al., 2015; reviewed in Domingues, Hilsdorf & Gadig (2018)). Though recent studies are moving towards multi-marker approaches (Momigliano et al., 2017; Pazmiño et al., 2018: Green et al., 2019), there is still a large body of literature focusing on mtCR. The mitochondrial genome has a higher rate of mutation than most of the nuclear genome (Brown, George & Wilson, 1979; Charlesworth & Wright, 2001; Neiman & Taylor, 2009) and this rate of mutation is a key advantage in vertebrates with slowly-evolving genomes (Avise et al., 1992; Martin, Naylor & Palumbi, 1992). Elasmobranch mtDNA studies to date have been successful in elucidating population partitions and evolutionary divergences, but the maternal inheritance of mtDNA can limit conclusions about gene flow in cases of sex-biased (usually male) dispersal. Both mtDNA and nuclear markers often have concordant results in sedentary species (Lavery, Moritz & Fielder, 1996; Avise, 2004; Zink & Barrowclough, 2008; DiBattista et al., 2015) but, when examined alone, may miss key components of population structure, particularly in migratory fauna (Pardini et al., 2001; Bowen et al., 2005; Toews & Brelsford, 2012). When highly mobile elasmobranchs are examined with both mtDNA and nuclear markers (usually microsatellites), a different picture often emerges in which females are more resident and males are dispersive (Pardini et al., 2001; Schultz et al., 2008; Portnoy et al., 2010; Karl et al., 2011; Daly-Engel et al., 2012; Portnoy et al., 2015: Bernard et al., 2017; Domingues et al., 2018). Identifying outlier SNPs in the nuclear genome can highlight genes possibly under selection, or show functional responses to environmental changes that have important management consequences (Jones et al., 2012; Fischer et al., 2013; Barrio et al., 2016; Guo, Li & Merilä, 2016). Therefore, the combination of mitochondrial and nuclear markers can yield fundamental ecological and evolutionary insights.

High-throughput sequencing is a powerful tool for revealing fine-scale population structure that may be missed by single locus studies (Andrews et al., 2016; Hohenlohe et al., 2018). However, this method can be costly, especially when examining many individuals as is typical of population genetic or phylogeography studies, and the perceived cost may prevent some from considering a high-throughput sequencing approach. For population genetics approaches based on differences in allele frequencies among populations, equimolar pooling of samples before sequencing is an affordable and accurate strategy for large-scale genetic analysis (Schlötterer et al., 2014). Several studies have successfully resolved population structure using a pooled site-associated DNA approach known as pool-seq, including some in commercially valuable marine species (Gautier et al., 2013; Mimee et al., 2015). Pool-seq provides estimates of allele frequencies for thousands of loci distributed across the genome simultaneously, which in some cases gives greater statistical power that can actually exceed the accuracy of allele frequency estimates based on individual sequencing (Futschik & Schlötterer, 2010, but also see Anderson, Skaug & Barshis, 2014). Therefore, a comparison of results between the standard mtCR analysis and high-throughput pool-seq is informative in evaluating the relative power and cost of the two approaches for examining population structure.

The silky shark (Carcharhinus falciformis (Müller & Henle, 1839)) is the second most commonly harvested shark on Earth (Rice & Harley, 2013; Oliver et al., 2015). They are one of the top contributors to the shark fin trade and the most common elasmobranch bycatch species in tuna purse-seine fisheries around the world (Clarke et al., 2006; Oliver et al., 2015; Cardeñosa et al., 2018). This pelagic shark, formerly abundant in all tropical oceans, has declined by an estimated 85% in the last 20 years, and is now listed as vulnerable and declining by the International Union for the Conservation of Nature (Rice & Harley, 2013; IUCN, 2017). Currently silky shark population assessments are conducted at the scale of regional fishery management organization, and conservation management measures are implemented at this scale in the absence of genetic or movement data to define population boundaries. Clarke et al. (2015) surveyed silky sharks across these regional management regions and found the western Atlantic was strongly differentiated from the Indian Ocean, but the North Atlantic, Gulf of Mexico, and Brazil could not be differentiated and appeared to comprise a single population. In contrast, using the same mtCR marker, Domingues et al. (2017) examined five regions across the Western Atlantic and found the North Western Atlantic was distinct from the South Western Atlantic. The difference between the two studies results from additional sampling in the South West Atlantic from further south than Clarke et al. (2015).

In an era where wildlife management needs far exceed the financial resources to address them, many seek to find the most accessible, robust, and economical means to define management units. In this study, we provide a direct comparison of population genetic analysis methods between Sanger sequencing of the mtCR region and high-throughput sequencing of regional pools of individuals. The same individuals from Clarke et al. (2015) were re-sequenced using pool-seq approaches. Regions re-sequenced included Gulf of Mexico, North West Atlantic, and Brazil, as well as one geographically distant location in the Red Sea (Fig. 1). We focused this analysis on SNPs from the mitochondrial DNA as well as nuclear DNA. We did not analyze any microsatellite loci because they were not a part of Clarke et al. (2015). We then evaluate the economics of conducting pool-seq relative to conventional Sanger sequencing of these same individuals. Ecological and management implications will be addressed in a subsequent companion paper.

Materials and Methods

Results

A total of 30.8 million reads were generated for the four geographic regions, which averaged 7.7 ± 3.0 million reads per pooled library. Results from the MutliQC assessment showed fairly homogenous output between libraries in regard to sequence quality scores, GC and per base sequence content, sequence length distributions, duplication levels, overrepresented sequences, and adapter content. Once assembled, aligned, and mapped, 5,792 SNPs were resolved across the mitochondrial and nuclear genomes combined. There were 4,103 biallelic SNPs, 168 multialleleic SNPs and 48 insertions and deletions (INDELs). INDELs and multiallelic SNPs remain a challenge for quantification software, so we restricted our analysis to biallelic loci (Fracassetti, Griffin & Willi, 2015). AssessPool creates visualizations of F_ST values and allows for visual outlier inspection. No visual outliers were present and given these SNPs are distributed haphazardly across the genome, therefore they are assumed to be putatively neutral.

Mitochondrial genome

Analysis of the complete mitochondrial genome (17,774 bp) revealed 804 variable sites: 681 biallelic and 17 multiallelic SNPs. Because coverage in this dataset was fairly low on average, most of these SNPs did not meet the filter threshold. After further filtering for the highest quality markers, 30 SNPs were selected to calculate allele frequencies. Pairwise F_st values were all significant (Fig. 2; Table S2). The Red Sea had much higher F_st values (ranging from 0.367 to 0.745) than any inter-Atlantic comparison (ranging from 0.029 to 0.135). However, all comparisons within the Atlantic still showed significant F_ST values, the highest being between the North Atlantic and Brazil, and the lowest between Brazil and Gulf of Mexico (Fig. 2; Table S2).

Cost analysis

The findings for cost analysis indicate that pool-seq reaches a threshold at approximately 300 individuals, after which this approach offers cheaper results than individual Sanger sequences. Furthermore, the cost is only twice as expensive at just over 100 individuals (Fig. 3A). The pool-seq approach provides a far higher ratio of information for the cost, yielding greater population resolution. This cost assessment does not include analytical time, labor, or effort associated with pool-seq analyses such as access to computer resources and expertise with bioinformatic pipelines. However, these costs are likely to decrease in the near future as bioinformatic pipelines are improved and become more widely available, for example as applications deployed via cloud based platforms such as Galaxy (https://usegalaxy.org/) or CyVerse (https://cyverse.org/). It is also important to note that the choice of pool-seq methodology has many caveats, which are discussed in greater detail in the ‘considerations on pool-seq’ section of the discussion below.

Discussion

Elasmobranchs are being harvested at unsustainable levels in several commercial fishing industries around the world. A fundamental step in successful management of any species is resolving population boundaries so they can be managed on a genetic stock-by-stock basis. As genetic sequencing technologies advance, there is greater opportunity to detect even small-scale genetic differences between populations. When these differences amount to statically significant allele frequencies at the population level, this indicates limited exchange among distinct stocks.

Here, we validate the utility of pool-seq using the same individuals as a previous study (Clarke et al., 2015) and show that pool-seq recovers additional population structure relative to Sanger sequencing of the mtDNA control region. Pool-seq was able to detect isolated populations between the Gulf of Mexico, Western Atlantic, and along the Brazilian coast, where Clarke et al. (2015) found no population structure. As expected, the Red Sea population was highly isolated from Atlantic conspecifics using both approaches.

One advantage of this pool-seq approach is that we recover SNPs throughout the entire mitochondrial genome along with thousands of additional nuclear loci that together provide greater statistical power to detect finer scale population structure (Ryman & Palm, 2006; Larsson et al., 2009; Kurland et al., 2019). The pool-seq approach yielded significant genetic structure among inter-Atlantic regions in both mtDNA and nuclear loci, whereas Sanger sequencing of the mtCR lacked power to resolve significant differences among the same populations. The congruence between the mitochondrial genome and nuclear loci reinforces the conclusion of population structure among all regions sampled in this study.

In this case, pool-seq lived up to the promise of increased power to detect fine-scale structure, but does it live up to the promise (Ferretti, Ramos-Onsins & Pérez-Enciso, 2013; Schlötterer et al., 2014) of being cost-effective? Individual extraction costs remain fixed across both approaches and Sanger sequencing generally has a flat rate per individual, including PCR primers and reagents, and sequencing per individual per locus. In contrast, pool-seq has a flat sequencing cost determined by the number of reads generated from the high-throughput sequencing platform, plus a small additional cost per pool for the exact quantification of DNA for equimolar pooling and the library preparation for high-throughput sequencing. Comparing costs at our institution between a single Sanger sequencing marker and pool-seq on the Illumina MiSeq platform indicates pool-seq becomes less expensive when sample size of the study rises above 300 individuals. Although the cost per pool is essentially fixed, when higher numbers of individuals are included per pool, the price per individual analyzed is further reduced (Fig. 3B). Our comparison here is limited to 12 pools due to the maximum number of reads per lane produced on the MiSeq platform. Therefore, analyzing more than 12 pools would require additional sequencing runs and result in a step increase in the cost per individual/pool, although this would differ among other Illumina machines (such as the HiSeq, NextSeq or NovaSeq) or other high-throughput sequencing platforms (such as the PacBio Sequel II). Larger numbers of pools could be run on some of these machines, but with differing individual read lengths and sequencing depths, which also introduce other trade-offs. Likewise, samples can also be run with individual barcodes, therefore gaining the individual information lost by pooling specimens, but with increased initial setup and sequencing costs. There are so many options by which to apply these methods that we cannot possibly consider them all here, and the availability, cost, and trade-offs associated with each should be ideally considered by individuals when designing high-throughput sequencing projects. In our case, we considered only the options currently available to us through our campus sequencing core, and all these pool-seq price comparisons are to a single Sanger-sequenced marker. Thus, when considering the information acquired from pool-seq compared to the cost from traditional single mitochondrial marker the price per individual advantage is massively amplified.

Conclusions

The finding of population structure on the scale of North Atlantic/Gulf of Mexico/Brazil is nearly unprecedented for a pelagic shark. Population structure in globally distributed sharks is typically detected on a scale of ocean basins (Atlantic versus Indo-Pacific, Castro et al., 2007; Graves & McDowell, 2015) and a few pelagic fishes have no population structure on a global scale (e.g., Basking shark, Cetorhinus maximus, Hoelzel et al., 2006; Blue shark Prionace glauca, Veríssimo et al., 2017; Wahoo, Acanthocybium solandri, Theisen et al., 2008). The resolution of isolated populations on the scale of North Atlantic Ocean is more typical of coastal species than pelagic species. The silky shark seems to be a pelagic species with a somewhat coastal population structure. This has strong implications for international management because smaller stocks imply smaller populations which are more readily depleted. At a minimum, these data require rethinking a single population management approach for the Atlantic, and this pattern needs to be investigated for this species across the Indo-Pacific as well.

Overall this study demonstrates pool-seq is a powerful and cost-effective tool for analyzing large portions of the genome which the methods traditionally used for elasmobranchs could not supply. Sharks and rays are an imperiled group of species that could benefit from advanced genomic studies to outline appropriate management units. Finally, although the technology is becoming cheaper and easier to apply, it is a common pitfall to assume everyone in the field can afford, or must use, these approaches to produce defensible science. Bowen et al. (2014) advocate judicious rather than wholesale application of genomic approaches as the most robust course of study, particularly when considering the global inequities in available research budgets. Sanger sequencing is still more cost effective for small numbers of individuals, but as the number of individuals included in a study rise, the cost per individual reaches the point where high throughput sequencing studies can be cheaper than sequencing a single mitochondrial marker from each individual. We provide an example of just such a case here, and highlight the potential advantage of cost savings together with increased power for resolution of fine scale population structure. Though there is still additional cost of using cluster computer servers and bioinformatics programs, these cost are dropping as technology advances. When study organism and sampling strategies are assessed and implemented into the study design, pool-seq has great promise for augmenting the scientific foundations for management of marine recourses.

Supplemental Information