HACSim: an R package to estimate intraspecific sample sizes for genetic diversity assessment using haplotype accumulation curves

Algorithm functions

Our algorithm, HACSim (short for Haplotype Accumulation Curve Simulator), consisting of two user-defined R functions, HAC.sim() and HAC.simrep(), was created to run simulations of haplotype accumulation curves based on user-supplied parameters. The simulation treats species haplotypes as distinct character labels relative to the number of individuals possessing a given haplotype. The usual convention in this regard is that Haplotype 1 is the most frequent, Haplotype 2 is the next most frequent, etc. (Gwiazdowski et al., 2013). A haplotype network represents this scheme succinctly (Fig. 1).

Such an implementation closely mimics that seen in natural species populations, as each character label functions as a unique haplotype linked to a unique DNA barcode sequence. The algorithm then randomly samples species haplotype labels in an iterative fashion with replacement until all unique haplotypes have been observed. This process continues until all species haplotypes have been sampled. The idea is that levels of species haplotypic variation that are currently catalogued in BOLD can serve as proxies for total haplotype diversity that may exist for a given species. This is a reasonable assumption given that, while estimators of expected diversity are known (e.g., Chao1 abundance) (Chao, 1984), the frequencies of unseen haplotypes are not known a priori. Further, assuming a species is sampled across its entire geographic range, haplotypes not yet encountered are presumed to occur at low frequencies (otherwise they would likely have already been sampled).

Because R is an interpreted programming language (i.e., code is run line-by-line), it is slow compared to faster alternatives which use compilation to convert programs into machine-readable format; as such, to optimize performance of the present algorithm in terms of runtime, computationally-intensive parts of the simulation code were written in the C++ programming language and integrated with R via the packages Rcpp (Eddelbuettel & François, 2011) and RcppArmadillo (Eddelbuettel & Sanderson, 2014). This includes function code to carry out haplotype accumulation (via the function accumulate(), which is not directly called by the user). A further reason for turning to C++ is because some R code (e.g., nested ‘for’ loops) is not easily vectorized, nor can parallelization be employed for speed improvement due to loop dependence. The rationale for employing R for the present work is clear: R is free, open-source software that it is gaining widespread use within the DNA barcoding community due to its ease-of-use and well-established user-contributed package repository (Comprehensive R Archive Network (CRAN)). As such, the creation and disemination of HACSim as a R framework to assess levels of standing genetic variation within species is greatly facilitated.

A similar approach to the novel one proposed here to automatically generate haplotype accumulation curves from DNA sequence data is implemented in the R package spider (SPecies IDentity and Evolution in R; (Brown et al., 2012)) using the haploAccum() function. However, the approach, which formed the basis of earlier work carried out by Phillips et al. (2015), is quite restrictive in its functionality and, to our knowledge, is currently the only method available to generate haplotype accumulation curves in R because spider generates haplotype accumulation curves from DNA sequence alignments only and is not amenable to inclusion of numeric inputs for specimen and haplotype numbers. Thus, the method could not be easily extended to address our question. This was the primary reason for the proposal of a statistical model of sampling sufficiency by Phillips et al. (2015) and its extension described herein.

Algorithm parameters

At present, the algorithm (consisting of HAC.sim() and HAC.simrep()) takes 13 arguments as input (Table 1).

A user must first specify the number of observed specimens/DNA sequences (N) and the number of observed haplotypes (i.e., unique DNA sequences) (H*) for a given species. Both N and H* must be greater than one. Clearly, N must be greater than or equal to H*.

Next, the haplotype frequency distribution vector must be specified. The probs argument allows for the inclusion of both common and rare species haplotypes according to user interest (e.g., equally frequent haplotypes, or a single dominant haplotype). The resulting probs vector must have a length equal to H*. For example, if H* = 4, probs must contain four elements. The total probability of all unique haplotypes must sum to one.

The user can optionally input the fraction of observed haplotypes to capture p. By default, p = 0.95, mirroring the approach taken by both Zhang et al. (2010) and Bergsten et al. (2012) who computed intraspecific sample sizes needed to recover 95% of all haplotype variation for a species. At this level, the generated haplotype accumulation curve reaches a slope close to zero and further sampling effort is unlikely to uncover any new haplotypes. However, a user may wish to obtain sample sizes corresponding to different haplotype recovery levels, e.g.,  p = 0.99 (99% of all estimated haplotypes found). In the latter scenario, it can be argued that 100% of species haplotype variation is never actually achieved, since with greater sampling effort, additional haplotypes are almost surely to be found; thus, a true asymptote is never reached. In any case, simulation completion times will vary depending on inputted parameter values, such as probs, which controls the skewness of the observed haplotype frequency distribution.

The perms argument is in place to ensure that haplotype accumulation curves “smooth out” and tend to H* asymptotically as the number of permutations (replications) is increased. The effect of increasing the number of permutations is an increase in statistical accuracy and consequently, a reduction in variance. The proposed simulation algorithm outputs a mean haplotype accumulation curve that is the average of perms generated haplotype accumulation curves, where the order of individuals that are sampled is randomized. Each of these perms curves is a randomized step function (a sort of random walk), generated according to the number of haplotypes found. A permutation size of 1,000 was used by Phillips et al. (2015) because smaller permutation sizes yielded non-smooth (noisy) curves. Permutation sizes larger than 1,000 typically resulted in greater computation time, with no noticeable change in accumulation curve behaviour (Phillips et al., 2015). By default, perms = 10,000 (in contrast to Phillips et al. (2015)), which is comparable to the large number of replicates typically employed in statistical bootstrapping procedures needed to ensure accuracy of computed estimates (Efron, 1979). Sometimes it will be necessary for users to sacrifice accuracy for speed in the presence of time constraints. This can be accomplished through decreasing perms. Doing so however will result in only near-optimal solutions for specimen sample sizes. In some cases, it may be necessary to increase perms to further smooth out the curves (to ensure monotonicity), but this will increase algorithm runtime substantially.

Should a user wish to analyze their own intraspecific COI DNA barcode sequence data (or sequence data from any single locus for that matter), setting input.seqs = TRUE allows this (via the read.dna() function in ape). In such a case, a pop-up file window will prompt the user to select the formatted FASTA file of aligned/trimmed sequences to be read into R. When this occurs, arguments for N, H* and probs are set automatically by the algorithm via functions available in the R packages ape (Analysis of Phylogenetics and Evolution) (Paradis, Claude & Strimmer, 2004) and pegas (Population and Evolutionary Genetics Analysis System) (Paradis, 2010). Users must be aware however that the number of observed haplotypes treated by pegas (via the haplotype() function) may be overestimated if missing/ambiguous nucleotide data are present within the final alignment input file. Missing data are explicitly handled by the base.freq() function in the ape package. When this occurs, R will output a warning that such data are present within the alignment. Users should therefore consider removing sequences or sites comprising missing/ambiguous nucleotides. This step can be accomplished using external software such as MEGA (Molecular Evolutionary Genetics Analysis; (Kumar, Stecher & Tamura, 2016)). The BARCODE standard (Hanner, 2009) was developed to help identify high quality sequences and can be used as a quality filter if desired. Exclusion of low-quality sequences also has the advantage of speeding up compution time of the algorithm significantly.

Options for confidence interval (CI) estimation and graphical display of haplotype accumulation is also available via the argument conf.level, which allows the user to specify the desired level of statistical confidence. CIs are computed from the sample $\frac{\alpha }{2}100\text{\%}$ and $\left(1-\frac{\alpha }{2}\right)100\text{\%}$ quantiles of the haplotype accumulation curve distribution. The default is conf.level = 0.95, corresponding to a confidence level of 95%. High levels of statistical confidence (e.g., 99%) will result in wider confidence intervals; whereas low confidence leads to narrower interval estimates.

How does HACSim work?

Haplotype labels are first randomly placed on a two-dimensional spatial grid of size perms × N (read perms rows by N columns) according to their overall frequency of occurrence (Fig. 2).

The cumulative mean number of haplotypes is then computed along each column (i.e., for every specimen). If all H* haplotypes are not observed, then the grid is expanded to a size of perms × N* and the observed haplotypes enumerated. Estimation of specimen sample sizes proceeds iteratively, in which the current value of N* is used as a starting value to the next iteration (Fig. 2). An analogy here can be made to a game of golf: as one aims towards the hole and hits the ball, it gets closer and closer to the hole; however, one does not know the number of times to hit the ball before it lands in the hole. It is important to note that since sample sizes must be whole values, estimates of N* found at each iteration are rounded up to the next whole number. Even though this approach is quite conservative, it ensures that estimates are adequately reflective of the population from which they were drawn. HAC.sim(), which is called internally from HAC.simrep(), performs a single iteration of haplotype accumulation for a given species. In the case of real species, resulting output reflects current levels of sampling effort found within BOLD (or another similar sequence repository such as GenBank) for a given species. If the desired level of haplotype recovery is not reached, then HAC.simrep() is called to perform successive iterations until the observed fraction of haplotypes captured (R) is at least p. This stopping criterion is the termination condition necessary to halt the algorithm as soon as a “good enough” solution has been found. Such criteria are widely employed within numerical analysis. At each step of the algorithm, a dataset, in the form of a dataframe (called “d”) consisting of the mean number of haplotypes recovered (called means), along with the estimated standard deviation (sds) and the number of specimens sampled (specs) is generated. The estimated required sample size (N*) to recover a given proportion of observed species haplotypes corresponds to the endpoint of the accumulation curve. An indicator message is additionally outputted informing a user as to whether or not the desired level of haplotype recovery has been reached. The algorithm is depicted in Fig. 3.

In Fig. 3, all input parameters are known a priori except H_i, which is the number of haplotypes found at each iteration of the algorithm, and $R_i=\frac{H_i}{H^{\ast }}$, which is the observed fraction of haplotype recovery at iteration i. The equation to compute N* (1)${\displaystyle N_{i+1}^{\ast }=N_i+\frac{N_i}{H_i}\left(H^{\ast }-H_i\right)=\frac{N_i{H}^{\ast }}{H_i}=\frac{N_i}{R_i}}$

is quite intuitive since as H_i approaches H*, H^∗ − H_i approaches zero, $R_i=\frac{H_i}{H^{\ast }}$ approaches one, and consequently, N_i approaches N*. In the first part of the above equation, the quantity $\frac{N_i}{H_i}\left(H^{\ast }-H_i\right)$ is the amount by which the haplotype accumulation curve is extrapolated, which incorporates random error and uncertainty regarding the true value of θ in the search space being explored. Nonparametric estimates formed from the above iterative method produce a convergent monotonically-increasing sequence, which becomes closer and closer to N* as the number of iterations increase; that is, (2)${\displaystyle N_1^{\ast }\le N_2^{\ast }\le ...\le N_i^{\ast }\le N_{i+1}^{\ast }\to N^{\ast }}$

which is clearly a desirable property. Since haplotype accumulation curves are bounded below by one and bounded above by H*, then the above sequence has a lower bound equal to the initial guess for specimen sampling sufficiency (N) and an upper bound of N*.

Along with the iterated haplotype accumulation curves and haplotype frequency barplots, simulation output consists of the five initially proposed “measures of sampling closeness”, the estimate of θ (N*) based on Phillips et al. (2015)’s sampling model, in addition to the number of additional samples needed to recover all estimated total haplotype variation for a given species (N^∗ − N; Fig. 4) (Table 1).

These five quantities are given as follows: (1) Mean number of haplotypes sampled: H_i, (2) Mean number of haplotypes not sampled: H^∗ − H_i, (3) Proportion of haplotypes sampled: $\frac{H_i}{H^{\ast }}$, (4) Proportion of haplotypes not sampled: $\frac{H^{\ast }-H_i}{H^{\ast }}$, (5) Mean number of individuals not sampled: $N^{\ast }-N_i=\frac{N_i}{H_i}\left(H^{\ast }-H_i\right)$ and are analogous to absolute and relative approximation error metrics seen in numerical analysis. It should be noted that the mean number of haplotypes captured at each iteration, H_i, will not necessary be increasing, even though estimates of the cumulative mean value of N* are. It is easily seen above that H_i approaches H* with increasing number of iterations. Similarly, as the simulation progresses, H^∗ − H_i, $\frac{H^{\ast }-H_i}{H^{\ast }}$ and $N^{\ast }-N_i=\frac{N_i}{H_i}\left(H^{\ast }-H_i\right)$ all approach zero, while $\frac{H_i}{H^{\ast }}$ approaches one. The rate at which curves approach H* depends on inputs to both HAC.sim() and HAC.simrep(). Once the algorithm has converged to the desired level of haplotype recovery, a summary of findings is outputted consisting of (1) the initial guess (N) for sampling sufficiency; (2) the total number of iterations until convergence and simulation runtime (in seconds); (3) the final estimate (N*) of sampling sufficiency, along with an approximate (1 − α)100% confidence interval (see next paragraph); and, (4) the number of additional specimens required to be sampled (N^∗ − N) from the initial starting value. Iterations are automatically separated by a progress meter for easy visualization.

An approximate symmetric (1 − α)100% CI for θ is derived using the (first order) Delta Method (Casella & Berger, 2002). This approach relies on the asymptotic normality result of the CLT and employs a first-order Taylor series expansion around θ to arrive at an approximation of the variance (and corresponding standard error) of N*. Such an approach is convenient since the sampling distribution of N* would likely be difficult to compute exactly due to specimen sample sizes being highly taxon-dependent. An approximate (large sample) (1 − α)100% CI for θ is given by (3)${\displaystyle N^{\ast }\pm z_{1-\frac{\alpha }{2}}\left(\frac{{\hat{\sigma }}_H}{H}\sqrt{N^{\ast }}\right)}$

where $z_{1-\frac{\alpha }{2}}$ denotes the appropriate critical value from the standard Normal distribution and ${\hat{\sigma }}_H$ is the estimated standard deviation of the mean number of haplotypes recovered at N*. The interval produced by this approach is quite tight, shrinking as H_i tends to H*. By default, HACSim computes 95% confidence intervals for the abovementioned quantities.

It is important to consider how a confidence interval for θ should be interpreted. For instance, a 95% CI for θ of (L, U), where L and U are the lower and upper endpoints of the confidence interval respectively, does not mean that the true sampling sufficiency lies between (L, U) with 95% probability. Instead, resulting confidence intervals for θ are themselves random and should be interpreted in the following way: with repeated sampling, one can be (1 − α)100% confident that the true sampling sufficency for p% haplotype recovery for a given species lies in the range (L, U) (1 − α)100% of the time. That is, on average, (1 − α)100% of constructed confidence intervals will contain θ (1 − α)100% of the time. It should be noted however that as given computed confidence intervals are only approximate in the limit, desired nominal probability coverage may not be achieved. In other words, the proportion of times calculated (1 − α)100% intervals actually contain θ may not be met.

HACSim has been implemented as an object-oriented framework to improve modularity and overall user-friendliness. Scenarios of hypothetical and real species are contained within helper functions which comprise all information necessary to run simulations successfully without having to specify certain function arguments beforehand. To carry out simulations of sampling haplotypes from hypothetical species, the function HACHypothetical() must first be called. Similarly, haplotype sampling for real species is handled by the function HACReal(). In addition to all input parameters rquired by HAC.sim() and HAC.simrep() outlined in Table 1, both HACHypothetical() and HACReal() take further arguments. Both functions take the optional argument filename which is used to save results outputted to the R console to a CSV file. When either HACHypothetical() or HACReal() is invoked (i.e., assigned to a variable), an object herein called HACSObj is created containing the 13 arguments employed by HACSim in running simulations. Note the generated object can have any name the user desires. Further, all simulation variables are contained in an environment called ‘envr’ that is hidden from the user.

Equal haplotype frequencies

Figure 5 shows sample graphical output of the proposed haplotype accumulation curve simulation algorithm for a hypothetical species with N = 100 and H* = 10. All haplotypes are assumed to occur with equal frequency (i.e., probs = 0.10). Algorithm output is shown below.

Algorithm output shows that R = 100% of the H* = 10 haplotypes are recovered from the random sampling of N = 100 individuals, with lower and upper 95% confidence limits of 100–100. No additional specimens need to be collected (N^∗ − N = 0). Simulation results, consisting of the six “measures of sampling closeness” computed at each iteration, can be optionally saved in a comma-separated value (CSV) file called ‘output.csv’ (or another filename of the user’s choosing). Figure 5 shows that when haplotypes are equally frequent in species populations, corresponding haplotype accumulation curves reach an asymptote very quickly. As sampling effort is increased, the confidence interval becomes narrower, thereby reflecting one’s increased confidence in having likely sampled the majority of haplotype variation existing for a given species. Expected counts of the number of specimens possessing a given haplotype can be found from running max(envr$d$specs) * envr$probs in the R console once a simulation has converged. However, real data suggest that haplotype frequencies are not equal.

Unequal Haplotype Frequencies

Figures 6 and 7 show sample graphical output of the proposed haplotype accumulation curve simulation algorithm for a hypothetical species with N = 100 and H* = 10. All haplotypes occur with unequal frequency. Haplotypes 1–3 each have a frequency of 30%, while the remaining seven haplotypes each occur with a frequency of c. 1.4%.

Note that not all iterations are displayed above for the sake of brevity; only the first and last two iterations are given. With an initial guess of N = 100, only R = 83.3% of all H* = 10 observed haplotypes are recovered. The value of N* = 121 in the first iteration above serves as an improved initial guess of the true sampling sufficiency, which is an unknown quantity that is being estimated. This value is then fed back into the algorithm and the process is repeated until convergence is reached.

Using Eq. (1), the improved sample size was calculated as $N^{\ast }=100+\frac{100}{8.3291}\left(10-8.3291\right)=120.061$. After one iteration, the curve has been extrapolated by an additional N^∗ − N_i = 20.06099 individuals. Upon convergence, R = 95.4% of all observed haplotypes are captured with a sample size of N* = 180 specimens, with a 95% CI of 178.278–181.722. Given that N = 100 individuals have already been sampled, the number of additional specimens required is N^∗ − N = 80 individuals. The user can verify that sample sizes close to that found by HACSim are needed to capture 95% of existing haplotype variation. Simply set N = N* = 180 and rerun the algorithm. The last iteration serves as a check to verify that the desired level of haplotype recovery has been achieved. The value of N* = 188.7623 that is outputted at this step can be used as a good starting guess to extrapolate the curve to higher levels of haplotype recovery to save on the number of iterations required to reach convergence. To do this, one simply runs HACHypothetical() with N = 189.

Lake Whitefish (Coregonus clupeaformis)

An interesting case study on which to focus is that of Lake whitefish (Coregonus clupeaformis). Lake whitefish are a commercially, culturally, ecologically and economically important group of salmonid fishes found throughout the Laurentian Great Lakes in Canada and the United States, particularly to the Saugeen Ojibway First Nation (SON) of Bruce Peninsula in Ontario, Canada, as well as non-indigenous fisheries (Ryan & Crawford, 2014).

The colonization of refugia during Pleistocene glaciation is thought to have resulted in high levels of cryptic species diversity in North American freshwater fishes (Hubert et al., 2008; April et al., 2011; April et al., 2013a; April et al., 2013b). Overdyk et al. (2015) wished to investigate this hypothesis in larval Lake Huron lake whitefish. Despite limited levels of gene flow and likely formation of novel divergent haplotypes in this species, surprisingly, no evidence of deep evolutionary lineages was observed across the three major basins of Lake Huron despite marked differences in larval phenotype and adult fish spawning behaviour (Overdyk et al., 2015). This may be the result of limited sampling of intraspecific genetic variation, in addition to presumed panmixia (Overdyk et al., 2015). While lake whitefish represent one of the most well-studied fishes within BOLD, sampling effort for this species has nevertheless remained relatively static over the past few years. Thus, lake whitefish represent an ideal species for further exploration using HACSim.

In applying the developed algorithm to real species, sequence data preparation methodology followed that which is outlined in Phillips et al. (2015). Curation included the exclusion of specimens linked to GenBank entries, since those records without the BARCODE keyword designation lack appropriate metadata central to reference sequence library construction and management (Hanner, 2009). Our approach here was solely to assess comprehensiveness of single genomic sequence databases rather than incorporating sequence data from multiple repositories; thus, all DNA barcode sequences either originating from, or submitted to GenBank were not considered further. As well, the presence of base ambiguities and gaps/indels within sequence alignments can lead to bias in estimate haplotype diversity for a given species.

Currently (as of November 28, 2018), BOLD contains public (both barcode and non-barcode) records for 262 C. clupeaformis specimens collected from Lake Huron in northern parts of Ontario, Canada and Michigan, USA. Of the barcode sequences, N = 235 are of high quality (full-length (652 bp) and comprise no missing and/or ambiguous nucleotide bases). Haplotype analysis reveals that this species currently comprises H* = 15 unique COI haplotypes. Further, this species shows a highly-skewed haplotype frequency distribution, with a single dominant haplotype accounting for c. 91.5% (215/235) of all individuals (Fig. 8). The output of HACSim is displayed below.

From the above output, it is clear that current specimen sample sizes found within BOLD for C. clupeaformis are probably not sufficient to capture the majority of within-species COI haplotype variation. An initial sample size of N = 235 specimens corresponds to recovering only 73.8% of all H* = 15 unique haplotypes for this species (Fig. 9).

A sample size of N* = 604 individuals (95% CI [601.504–606.496]) would likely be needed to observe 95% of all existing genetic diversity for lake whitefish (Fig. 10).

Since N = 235 individuals have been sampled previously, only N^∗ − N = 369 specimens remain to be collected.

Deer tick (Ixodes scapularis)

Ticks, particularly the hard-bodied ticks (Arachnida: Acari: Ixodida: Ixodidae), are well-known as vectors of various zoonotic diseases including Lyme disease (Ondrejicka et al., 2014). Apart from this defining characteristic, the morphological identification of ticks at any lifestage, by even expert taxonomists, is notoriously difficult or sometimes even impossible (Ondrejicka, Morey & Hanner, 2017). Further, the presence of likely high cryptic species diversity in this group means that turning to molecular techniques such as DNA barcoding is often the only feasible option for reliable species diagnosis. Lyme-competent specimens can be accurately detected through employing a sensitive quantitative PCR (qPCR) procedure (Ondrejicka, Morey & Hanner, 2017). However, for such a workflow to be successful, wide coverage of within-species haplotype variation from across broad geographic ranges is paramount to better aid design of primer and probe sets for rapid species discrimination. Furthermore, the availability of large specimen sample sizes for tick species of medical and epidemiological relevance is necessary for accurately assessing the presence of the barcode gap.

Notably, the deer tick (Ixodes scapularis), native to Canada and the United States, is the primary carrier of Borrelia burgdorferi, the bacterium responsible for causing Lyme disease in humans in these regions. Because of this, I. scapularis has been the subject of intensive taxonomic study in recent years. For instance, in a recent DNA barcoding study of medically-relevant Canadian ticks, Ondrejicka, Morey & Hanner (2017) found that out of eight specimens assessed for the presence of B. burgdorferi, 50% tested positive. However, as only exoskeletons and a single leg were examined for systemic infection, the reported infection rate may be a lower bound due to the fact that examined specimens may still harbour B. burgdorferi in their gut. As such, this species is well-represented within BOLD and thus warrants further examination within the present study.

As of August 27, 2019, 531 5′-COI DNA barcode sequences are accessble from BOLD’s Public Data Portal for this species. Of these, N = 349 met criteria for high quality outlined in Phillips et al. (2015). A 658 bp MUSCLE alignment comprised H* = 83 unique haplotypes. Haplotype analysis revealed that Haplotypes 1–8 were represented by more than 10 specimens (range: 12–46; Fig. 11).

Simulation output of HACSim is depicted below.

The above results clearly demonstrate the need for increased specimen sample sizes in deer ticks. With an initial sample size of N = 349 individuals, only 78.7% of all observed haplotypes are recovered for this species (Fig. 12). N* = 803 specimens (95% CI: 801.551-804.449) is necessary to capture at least 95% of standing haplotype variation for I. scapularis (Fig. 13) . Thus, a further N* − N = 454 specimens are required to be collected.

Scalloped hammerhead (Sphyrna lewini)

Sharks (Chondrichthyes: Elasmobranchii: Selachimorpha) represent one of the most ancient extant lineages of fishes. Despite this, many shark species face immediate extinction as a result of overexploitation, together with a unique life history (e.g., K-selected, predominant viviparity, long gestation period, lengthy time to maturation) and migration behaviour (Hanner, Naaum & Shivji, 2016). A large part of the problem stems from the increasing consumer demand for, and illegal trade of, shark fins, meat and bycatch on the Asian market. The widespread, albeit lucrative practice of “finning”, whereby live sharks are definned and immediately released, has led to the rapid decline of once stable populations (Steinke et al., 2017). As a result, numerous shark species are currently listed by the International Union for the Conservation of Nature (IUCN) and the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES). Interest in the molecular identification of sharks through DNA barcoding is multifold. The COI reference sequence library for this group remains largely incomplete. Further, many shark species exhibit high intraspecific distances within their barcodes, suggesting the possibility of cryptic species diversity. Instances of hybridization between sympatric species has also been documented. As establishing species-level matches to partial specimens through morphology alone is difficult, and such a task becomes impossible once fins are processed and sold for consumption or use in traditional medicine, DNA barcoding has paved a clear path forward for unequivocal diagnosis in most cases.

The endangered hammerheads (Family: Sphyrnidae) represent one of the most well-sampled groups of sharks within BOLD to date. Fins of the scalloped hammerhead (Sphyrna lewini) are especially highly prized within IUU (Illegal, Unregulated, Unreported) fisheries due to their inclusion as the main ingredient in shark fin soup.

As of August 27, 2019, 327 S. lewini specimens (sequenced at both barcode and non-barcode markers), collected from several Food and Agriculture Organization (FAO) regions, including the United States, are available through BOLD’s Public Data Portal. Of these, all high-quality records (N = 171) were selected for alignment in MEGA7 and assessment via HACSim. The final alignment was found to comprise H* = 12 unique haplotypes, of which three were represented by 20 or more specimens (range: 28–70; Fig. 14).

HACSim results are displayed below.

Simulation output suggests that only 82.6% of all unique haplotypes for the scalloped hammerhead have likely been recovered (Fig. 15) with a sample size of N = 171.

Further, HACSim predicts that N* = 414 individuals (95% CI [411.937–416.063]) probably need to be randomly sampled to capture the majority of intraspecific genetic diversity within 5′-COI (Fig. 16). Since 171 specimens have already been collected, this leaves an additional N* − N = 243 individuals which await sampling.