ORIGINAL RESEARCH article
Sec. Evolutionary and Population Genetics
Volume 10 - 2019 | https://doi.org/10.3389/fgene.2019.00011
Cryptic Diversity, but to What Extent? Discordance Between Single-Locus Species Delimitation Methods Within Mainland Anoles (Squamata: Dactyloidae) of Northern Central America
- 1Department of Biology, Indiana University of Pennsylvania, Indiana, PA, United States
- 2Department of Biological Sciences, Clemson University, Clemson, SC, United States
- 3Department of Biology, Central Michigan University, Mount Pleasant, MI, United States
- 4Senckenberg Forschungsinstitut und Naturmuseum, Frankfurt, Germany
- 5Federación Hondureña de Deportes de Montaña Y Escalada, Departamento de Francisco Morazán, Tegucigalpa, Honduras
- 6Centro Zamorano de Biodiversidad, Escuela Agrícola Panamericana Zamorano, Departamento de Francisco Morazán, Tegucigalpa, Honduras
Single-locus molecular barcoding is a useful method for identifying overlooked and undescribed biodiversity, providing the groundwork for further systematic study and taxonomic investigation. A variety of methods for delimiting species from barcoding libraries have been developed and applied, allowing for rapid estimates of species diversity in a broad range of taxa. However, tree-based and distance-based analyses can infer different group assignments, potentially over- or underestimating the number of putative species groups. Here, we explore diversity of mainland species of anole lizards from the Chortís Block biogeographical province of northern Central America using a DNA barcoding approach, generating and analyzing cytochrome oxidase subunit I (COI) sequences for over 400 samples assignable to 33 of 38 (86.8%) native and one introduced mainland species. We subsequently tested the effects different models of nucleotide substitution, different species-delimitation algorithms, and reducing our dataset had on species delimitation estimates. We performed of two distance-based (ABGD, RESL) and three tree-based (bPTP, mPTP, GMYC) analyses on both the full dataset and a dataset consisting only of unique halotypes. From 34 nominal taxa, analyses of the full dataset recovered between 34 and 64 operational taxonomic units (OTUs), while analyses of the reduced dataset inferred between 36 and 59. Reassigning individuals to either mPTP-inferred or ABGD clustered (7.2% threshold) groups improved the detection of a barcoding gap across three different models of nucleotide substitution, removing overlap between intra- and interspecific distances. Our results highlight the underestimated diversity of mainland Chortís Block anoles, but the lack of congruence between analyses demonstrates the importance of considering multiple analytical methods when dealing with single-locus datasets. We recommend future studies consider the effects of different models of nucleotide substitution on proposed barcoding gaps, as well as the effect reducing a dataset to unique haplotypes may have on proposed diversity estimates.
The majority of fields in the biological sciences rely on the proper identification of species (Wheeler, 2004; de Queiroz, 2007); thus, the efficient and reliable identification and delimitation of taxa is of pivotal importance. A valuable tool for characterizing the taxonomic diversity of large clades of organisms with phenotypically conserved taxa both quickly and efficiently is DNA barcoding (Meyer and Paulay, 2005). The use of an approximately 650 basepair (bp) sequence of mitochondrial DNA (mtDNA) from the cytochrome c oxidase subunit I (COI) gene has been recognized as an effective “barcode” for animals (Hebert et al., 2003), though other fragments can be effective for species identification also (Hajibabaei et al., 2006). This target region can be sequenced easily from an unknown specimen and referenced against a library of known taxa to discern whether the specimen falls out within a recognized lineage or not (Meyer and Paulay, 2005).
Accelerating the pace of species discovery is a secondary goal of DNA barcoding (Stoeckle, 2003; Hebert et al., 2004), and the discovery of cryptic lineages using single-locus delimitation analyses can have profound impacts on a range of fields (Bickford et al., 2007). This approach can highlight areas where more careful taxonomic study, integrating morphological examination with multi-locus molecular analyses, is needed (Hajibabaei et al., 2007; Padial et al., 2009, 2010), especially in areas where genetic diversity is masked by conserved morphology (Hebert et al., 2003). Single-locus barcoding analyses, however, do not provide enough information to formally describe new species (DeSalle, 2006). Discordance between the boundaries of putative species inferred by different methods of single-locus delimitation can lead to uncertainty in diversity inferences, due to either over- or under-estimating the true number of lineages present in a sample (Blair and Bryson, 2017).
DNA barcoding has been applied to numerous vertebrate (e.g., Hebert et al., 2004; Álvarez-Castañeda et al., 2011; Mendoza et al., 2016; Ramirez et al., 2017; Barman et al., 2018) and invertebrate (e.g., Barrett and Hebert, 2005; Moura et al., 2011; Barco et al., 2016; Jin et al., 2018) taxa, but the barcoding of herpetofauna largely has lagged behind. The initiation of the “ColdCode” campaign (Murphy et al., 2013) has led to a more concentrated effort to provide reference libraries for reptiles and amphibians, including those focused on specific taxa (e.g., Dang et al., 2015; Liu et al., 2015) and broad-scale regional assessments (e.g., Nagy et al., 2012; Guarnizo et al., 2015; Hawlitscheck et al., 2016; Vasconcelos et al., 2016; Deichmann et al., 2017).
In Northern Central America, the Chortís Block biogeographic province is increasingly being recognized as an important and distinctive area for herpetofaunal diversity (Figure 1; Townsend, 2014). This geologically and ecologically-delimited region, also referred to as “Eastern Nuclear Central America” (Campbell, 1999; Wilson and Townsend, 2007; Townsend, 2009), is comprised of eastern Guatemala, all of mainland Honduras and El Salvador, and northern Nicaragua, including Isla del Tigre, the Honduran Islas de la Bahía, and the Nicaraguan Cayos Miskitos (Townsend, 2014). Over 400 reptile and amphibian species are found here, and importantly, over 37% of those recognized species are endemic to the region, including more than half of the area’s amphibians and over 30% percent of its squamates (Townsend, 2014).
Figure 1. Map of the Chortís Block Biogeographical Province; boundaries are denoted by a dashed line.
Anoles (Squamata: Dactyloidae) are particularly strong representatives of the Chortís Block’s squamate diversity, the region’s high levels of endemism, and the need for increased study and conservation work. As currently considered, thirty-eight species of anole lizards are native to the mainland Chortís Block, 20 (52.6%) of which are considered endemic (Table 1; Townsend, 2014; McCranie and Köhler, 2015; Köhler et al., 2016; Hofmann and Townsend, 2017, 2018). Two additional species, Anolis allisoni and Norops sagrei, have been introduced to the northern coast of mainland Honduras (McCranie and Köhler, 2015). They inhabit a variety of physiographic regions, habitats, and ecological roles throughout the Chortís Block (McCranie and Köhler, 2015); however, their often-conserved morphology combined with the paucity of published data on these species has led to a multitude of taxonomic complications.
Table 1. Summary of the mainland anoles of the Chortís Block: Species, number of samples (N) and localities used in this paper, and conservation scores.
Here, we highlight the discordance between putative species limits inferred by several popular single-locus delimitation analyses. To this end, we test two distance-based and three tree-based methods on a DNA barcode library for anoles of the mainland Chortís Block of Northern Central America—consisting of sequence data for 33 of the 38 (86.8%) native species and one introduced to the mainland (otherwise native to Cuba and to the Islas de la Bahía, as well as numerous other islands and keys in the Caribbean). Our goals are to provide an efficient molecular reference for species identification and draw attention to cryptic lineages in need of further taxonomic investigation, while simultaneously stressing the importance of incorporating multiple tests before drawing conclusions of cryptic species-level diversity.
Materials and Methods
Justification of Nomenclature
We follow Nicholson et al. (2018) in their application of a rank-based, multi-genera taxonomy for anoles (valid under the International Code of Zoological Nomenclature; ICZN) based on the clade names proposed by Poe et al. (2017). As such, we refer to the monophyletic grouping of beta anoles as Norops (sensu Nicholson et al., 2018), while recognizing the criticisms of the multiple-genera taxonomy (e.g., Poe, 2013). All native, mainland Norops in the Chortís Block are beta anoles and members of clade Draconura [Poe et al., 2017 (auratus group of Nicholson et al., 2012)], one of the three clades within Norops, all of which share the synapomorphy of anterolaterally directed transverse processes on their caudal vertebrae (Etheridge, 1959; Nicholson, 2002) as well as numerous molecular characters (Nicholson et al., 2012; Poe et al., 2017).
Between June 2006 and January 2016, we sampled 96 localities in the Chortís Block, including localities within 17 of 18 departments in Honduras and four of 15 departments in northern Nicaragua (Table 1 and Supplementary Table 1). Tissue samples were collected from vouchers and stored in SED buffer (250 mM EDTA/20% DMSO/saturated NaCl; Seutin et al., 1991; Williams, 2007). Voucher specimens were preserved in 10% formalin solution, and later transferred to 70% ethanol for permanent storage. Vouchers were deposited in the Florida Museum of Natural History (FLMNH), University of Florida (UF); the Museum of Vertebrate Zoology, Berkeley (MVZ); the National Museum of Natural History, Smithsonian Institution (USNM); the Senckenberg Forschungsinstitut und Naturmuseum (SMF); and the natural history collection of the Universidad Nacional Autónoma de Nicaraga-León (UNAN). For one taxon where samples from within the Chortís Block were not available (Norops uniformis) and another with uncertain species boundaries (N. cupreus), samples from outside the Chortís Block were included (Table 1). Three additional samples representing populations of three species from El Salvador were received from the Herpetology collections of the University of Kansas Biodiversity Institute (KU). All specimens were identified a priori based on external morphology, using the keys of Köhler (2008) or McCranie and Köhler (2015). We follow McCranie and Köhler (2015) in considering Norops dariense a synonym of Norops cupreus, until further taxonomic study is completed.
COI Amplification and Sequencing
DNA extraction, amplification, and sequencing prior to August 2015 was carried out at the Smithsonian Institution Laboratory of Analytical Biology (Suitland, MD, United States) following standardized DNA Barcode of Life (BOLD) protocols, while all other sequences were generated in the Townsend Lab at the Indiana University of Pennsylvania (Indiana, PA, United States). Whole-genome DNA was extracted from tissue using PureLink Genomic DNA Kits (Life Technologies). The mitochondrial gene cytochrome oxidase subunit I (COI), the standard vertebrate barcoding gene, was amplified by polymerase chain reaction (PCR) using the primers dgLCO-1490 and dgHCO-2198 (Meyer, 2003). PCR products were visualized via 1.5% agarose gel electrophoresis. Unincorporated nucleotides were removed from the PCR product using 2 μL of ExoSAP-IT® per sample, and sequencing (forward and reverse) was carried out via Eurofins SimpleSeq DNA Sequencing service (Louisville, KY, United States) following manufacturer’s protocols. Chromatograms were checked manually and sequences were assembled using Geneious v.7.1.7 (Kearse et al., 2012). Sequence alignment was carried out in MEGA7 (Kumar et al., 2016) using the ClustalW algorithm (Thompson et al., 1994) and sequences were manually checked to ensure there were no indels or stop codons. Two datasets were formed: a “full” dataset consisting of all sequences, and a “reduced” dataset consisting only of unique haplotypes. The numbers of conserved, variable, and parsimony-informative sites (not including the outgroup Anolis allisoni) were calculated in MEGA7.
Sequence divergences were estimated using uncorrected p-distances and under the K2P model, a model extensively used in barcoding studies (e.g., Liu et al., 2015; Hawlitscheck et al., 2016), with 1000 bootstrap estimates in MEGA7. Recent publications have discussed the widespread misuse of the K2P model in barcoding analyses and its tendency to underestimate diversity (e.g., Barley and Thomson, 2016; Zinger and Philippe, 2016). Therefore, we selected the best-fit model of nucleotide substitution for the dataset using jModelTest 2.0 (Darriba et al., 2012) based on the Bayesian Information Criterion. TrN+G was selected (Tamura and Nei, 1993; Δ = 4.65) as the best-fit model (TIM1+I+G) was not implementable in MEGA7. We then estimated sequence divergences a third time using this model. Between group and within group averages were taken from initial a priori identifications and by categorizing sequences to lineages as assigned by ABGD at a 7.2% threshold and the mPTP analysis (see below). These three groupings were chosen as conservative exemplars of a priori, distance-based, and tree-based delimitation methods, together with the fact that these produced congruent results between the “full” and “reduced” datasets (see below). We generated histograms of pairwise distances for all three models in R v3.2.3 (R Core Team, 2015) using the “checkDNAbcd” function in the package ad hoc (Sonet et al., 2013). We tested for possible substitution saturation and plotted transitions and transversions against K2P distances (Kimura, 1980) using DAMBE v.6.3.3 (Xia and Lemey, 2009; Xia, 2013). The first and second codon positions, the third codon position, and all codon positions were tested separately.
Evolutionary model-based hypotheses of phylogenetic relationships were estimated for these data. The goal of these analyses were not to infer evolutionary histories of the taxa using only single-locus datasets, but to provide a framework in which to test species delimitation methods. Both the reduced and full dataset were subjected to maximum likelihood (ML) phylogenetic analyses using the program RAxML v7.2.8 (Stamatakis, 2014), performed with 1,000 bootstrap replicates under the GTR+GAMMA substitution model with data partitioned by codon position. Samples of Anolis allisoni were included as the outgroup. To create an ultrametric gene tree for GMYC analyses (see below), we analyzed both datasets in BEAST v.2.4.7 (Bouckaert et al., 2014), using a strict clock, yule tree prior, a GTR+GAMMA model of substitution, and all other priors left at default values. Analyses were performed for 100 million generations, sampling every 5000. Tracer v.1.6 (Rambaut et al., 2014) was used to assess convergence and adequate posterior sampling (ESS > 200), and a maximum clade credibility tree was created using TreeAnnotator v.2.4.7 (Bouckaert et al., 2014) using mean heights for annotation.
Species Delimitation Analyses
We tested three tree-based methods of species delimitation on both the full and reduced datasets: bPTP (Zhang et al., 2013), multi-rate PTP (mPTP; Kapli et al., 2017), and the single threshold general mixed Yule coalescent model (GYMC; Pons et al., 2006; Fujisawa and Barraclough, 2013). bPTP analyses were performed on the online server1 using the ML trees from RAxML. We ran 500,000 generations with a thinning of 500 and a burn-in of 0.1, then assessed convergence visually using the MCMC iteration v. log-likelihood plots generated automatically. Next, we applied the recently introduced mPTP method, which improves upon the Poisson Tree Processes (PTP; Zhang et al., 2013) for single-locus species delimitation, to our datasets. Instead of all species sharing the same rate of evolution (λ) as in PTP, each species branch has its own λ in the mPTP model. This method determines which number of species fits best to the given data by utilizing the Akaike Information Criterion (rather than a p-value test as in PTP) because of the different number of parameters. mPTP has been shown to be consistent and very effective for species delimitation in datasets with uneven sampling (Blair and Bryson, 2017). Using the ML trees from RAxML, we performed four simultaneous Markov Chain Monte Carlo (MCMC) runs of 100,000,000 steps, sampling every 10,000 steps, to assess Average Support Values (ASV), the confidence of the ML delimitation for each species. Convergence of the runs was assessed visually using the outputted likelihood plot of the combined runs (created using the “--mcmc_log” command) and the Average Standard Deviation of Delimitation Support Values (ASDDSV), which approaches zero as the multiple MCMC runs converge on the same delimitation distribution. Finally, we incorporated the single threshold GMYC model for our full and reduced datasets using the summarized BEAST ultrametric trees. These analyses were performed using the R package ‘splits’ (Ezard et al., 2009).
To compare tree-based and distance-based methods of species delimitation, as well as statistically detect the barcode gap in our data, we then performed two distance-based methods: the Refined Single Linkage algorithm (RESL; Ratnasingham and Hebert, 2013) and the Automated Barcode Gap Discovery method (ABGD; Puillandre et al., 2012). RESL was implemented directly in the Barcode of Life Datasystem (Ratnasingham and Hebert, 2007, 2013)2, and used to assign sequences to OTUs. ABGD infers a model-based confidence limit for intraspecific divergence based on prior intraspecific divergences, clustering similar haplotypes together as “species.” The barcode gap is the first significant gap beyond the intraspecific divergence limit, and therefore two samples taken from distinct clusters will have a distance between them larger than the barcoding gap. This inference and gap detection is then continuously applied to the previous clusters until a final partition is reached (Guarnizo et al., 2015). Our alignments were processed in ABGD web3 using the Kimura two-parameter substitution model, prior for maximum value of intraspecific divergence between 0.001 and 0.1, 15 recursive steps, and a gap width (X) of 1.0.
In this formula, Nmatch refers to the number of delimited species that exactly match a taxonomically defined morphospecies (not including taxa split as cryptic lineages), Ndelimited refers to the total number of lineages delimited by an analysis, and Nmorph refers to the total number of morphologically defined species (morphospecies). The use of a match ratio provides a better comparable value for different analyses than simply reporting the Nmatch, where splitting and lumping species cancels out match values.
We successfully amplified 410 samples representing 33 native species of Norops and two samples of Anolis allisoni for COI, generating a full dataset of 412 sequences. Twenty-nine of the 34 total nominal species were represented by two or more samples, with five species represented by a single sample. Sequences ranged in length from 521 to 654 bp, with 395 sequences (95.9%) longer than 600 bp. Three hundred fifty-two sites (53.8%) were conserved, 302 (46.2%) were variable, and 281 (43.0%) were parsimony-informative. After removing redundant identical haplotypes, the reduced dataset consisted of 290 unique sequences (288 ingroup, 2 outgroup).
For tests of substitution saturation on all codon positions and codon position 1 and 2, index of substitution saturation (Iss) values were less than Issc.Sym (critical index of substitution saturation assuming a symmetrical topology) or Issc.Asym (critical index of substitution saturation assuming an asymmetrical topology) for all numbers of species simulated (NumOTU), suggesting little substitution saturation has occurred (Supplementary Table 2 and Supplementary Figure 1). Tests on only codon position 3, however, resulted in Iss values higher than Issc.Asym values for NumOTU values of 16 and 32. This suggests there potentially is saturation of the 3rd codon position for the dataset. DAMBE is limited to NumOTU ≤ 32, but tests random subsets of 4, 8, 16, and 32 OTUs multiple times in order to circumvent this limitation (Xia and Lemey, 2009).
Summaries of uncorrected pairwise distances, K2P-corrected distances, and TrN+G corrected distances are shown in Table 2 and Figure 2, with the full data available in Supplementary Data Sheet 1. Generating pairwise distances using the better-fit model (TrN+G) resulted in larger values than both uncorrected pairwise and K2P models, increasing the mean distance between samples from 0.171 (uncorrected p) to 0.211. Assigning sequences to lineages recovered by ABGD at a 7.2% threshold tightened the average intra- and interspecific distances, greatly improving the delineation of a barcoding gap threshold regardless of the model implemented. Assignment of samples by mPTP similarly tightened the range of average intraspecific distances, but led to a wider range of interspecific distances.
Table 2. Summaries of uncorrected pairwise distances, K2P-corrected distances, and TrN+G-corrected distances for the dataset, including averages based on a priori identification (morphospecies), ABGD clustering (at 7.2% threshold), and mPTP delimitations.
Figure 2. Histograms of pairwise distances under different models of nucleotide substitution: uncorrected (top), K2P (middle), and TrN+G (bottom).
Tree- and distance-based methods of species delimitation did not produce congruent results, often inferring different numbers of species with the same method applied to the full and reduced datasets (Table 3 and Supplementary Tables 3, 4). Methods applied to the reduced dataset (290 sequences) inferred between 36 (ABGD threshold of 7.2%) and 59 (bPTP) species, while the same methods applied to the full dataset (412 sequences) inferred between 34 (ABGD threshold of 10%) and 64 (bPTP) species (Table 3 and Figures 3, 4). Match ratios of analyses performed on the reduced dataset were higher than their complement performed on the full dataset in all cases except the ABGD threshold of 7.2% (higher in full) and RESL (equal). Across both datasets, ABGD match ratios were the highest of all analyses, while mPTP match ratios were highest among tree-based analyses. Results of ABGD at 7.2% (36 inferred species), RESL (51), and mPTP (46; confidence interval 44–49) were congruent across both datasets. GMYC inferred 50 (full dataset; confidence interval 50–55; match ratio 0.548) and 51 (reduced dataset; confidence interval 50–54; match ratio 0.565) species. bPTP recovered an unreasonably high number of species with wide confidence intervals in both the full (64: 53–80) and reduced (59: 53–78) datasets, and had the lowest match ratios of any analysis (full: 0.429; reduced: 0.452). Several previous studies have recommended the use of 10% threshold when interpreting ABGD results (Kekkonen and Hebert, 2014; Blair and Bryson, 2017). This threshold inferred 34 nominal species and the highest match ratio (0.824) when applied to our full dataset, failing to delimit N. unilobatus from N. wellborane or N. limifrons from N. zeus, but delimiting two OTUs within both N. rubribarbaris and N. yoroensis (Figures 3, 4).
Table 3. Number of putative species (OTUs) inferred by each delimitation method from input datasets (“Full”: 412 sequences; “Reduced”: 290 unique sequences) of 34 nominal taxa, with exact matches (Nmatch) and match ratios.
Figure 3. Maximum likelihood phylogeny of the full COI sequence dataset, with lineage assignments from the three tree-based (mPTP, bPTP, GMYC) and two distance-based (RESL, ABGD at three thresholds) methods. ML bootstrap support and Bayesian posterior probabilities shown when ≥ 50 and 0.50, respectively. Gray bars span all samples assigned a priori to named taxa (morphospecies). Note (∗): JMS71 was assigned to its own lineage by GMYC, apart from all other samples assigned to N. cupreus (i.e., there were two inferred lineages within N. cupreus, not three).
Figure 4. Maximum likelihood phylogeny of the full COI sequence dataset, with lineage assignments from the three tree-based (mPTP, bPTP, GMYC) and two distance-based (RESL, ABGD at three thresholds) methods (continued from Figure 3). ML bootstrap support and Bayesian posterior probabilities shown when ≥ 50 and 0.50, respectively. Gray bars span all samples assigned a priori to named taxa (morphospecies).
Species Delimitation Method Performance
Numerous studies have compared and contrasted the single-locus delimitation methods tested here across empirical and simulated datasets (e.g., Talavera et al., 2013; Dellicour and Flot, 2015, 2018; Ahrens et al., 2016; Blair and Bryson, 2017; Luo et al., 2018). Our results are consistent with many empirical studies showing that different methods often produce different delimitation scenarios when using single-locus data. Widely used methods of single-locus species delimitation tested here are each subject to the potential biases and differing conditions inherent in empirical datasets.
Despite recovering largely incongruent numbers of OTUs, the methods were consistent in recovering more OTUs than the number of species originally considered (Figures 3, 4 and Table 3). Given the rapid evolution of mainland anoles and the lack of clarity regarding the relationships of all populations analyzed, it is not surprising that more lineages were inferred than are currently recognized. While some of these OTUs might correspond to undescribed cryptic species, it is also likely that some structure is a result of genetic drift or isolated populations currently undergoing speciation. In some cases, however, it is clear that the analyses over-split taxa. As in Blair and Bryson (2017)’s analyses of Phrynosoma lizards, our bPTP analyses produced unreasonable delimitations with wide confidence intervals; some of these clusters do not reflect relationships as understood with better molecular sampling (Hofmann and Townsend, 2017) and others were separated into numerous lineages despite little-to-no divergence between them (e.g., N. cupreus; Figure 3). bPTP is known to be sensitive to different mutation rates, but unlike in the simulated data of Dellicour and Flot (2018), here it produced the least-accurate delimitations.
Several factors might have influenced the incongruent results among the other algorithms. ABGD is known to overlump, performing poorly on more speciose datasets with faster speciation rates compared to smaller, more slowly speciating populations (Dellicour and Flot, 2015, 2018). ABGD’s conservative tendencies have also been shown across a variety of loci and taxa, including amphibians (Guarnizo et al., 2015), other lizards (Blair and Bryson, 2017), brittle stars (Boissin et al., 2017), and insects (e.g., Song et al., 2018). Here, ABGD recovered the fewest inferred species and had the highest match ratios at all levels, but in no instance was the delimitation completely congruent with current taxonomy.
In contrast to ABGD’s overlumping, single-threshold GMYC is known to oversplit species (Talavera et al., 2013), as higher substitution rates, uneven sampling (including singletons or the inclusion of identical sequences), variation in population size among species, ongoing gene flow, or unresolved nodes could bias results (Esselstyn et al., 2012; Tang et al., 2014; Ahrens et al., 2016; Blair and Bryson, 2017; Luo et al., 2018). Five of the species sampled herein were represented by only a single sequence; in combination with variable effective sample sizes of these anoles, these samples might have contributed to the more liberal GMYC delimitations (Ahrens et al., 2016). This method recovered more clusters than distance-based analyses across our data, as well as in studies of Socotran reptiles (Vasconcelos et al., 2016), Madascincus lizards (Miralles and Vences, 2013), and Ophiomorus geckos (Korniolios et al., 2018). Interestingly, Blair and Bryson (2017) found it to be the most conservative method applied to their data.
The use of multiple models of nucleotide substitution (uncorrected pairwise, K2P, and TrN+G) resulted in a wide range of observed distances within the samples, with TrN+G estimating the widest range of distances. A difference of greater than 9% at the upper limit of interspecific differences is substantial enough to suggest the incorporation of multiple models as a standard practice in barcoding investigations, including those traditionally used models (uncorrected p and K2P), as well as the best fit model, as suggested by several authors (Srivathsan and Meier, 2012; Barley and Thomson, 2016; Zinger and Philippe, 2016). Regardless of the model of nucleotide substitution used, a clear barcoding gap was recovered by assigning samples based on delimitation methods rather than a priori (Table 2). The use of either ABGD or mPTP delimitations improved distance analyses, tightening intraspecific average distances across the board. A proposed COI barcoding gap between 3.7 and 7.1% could be considered for these anoles, though we caution against a strict cutoff as it is clear further taxonomic investigation is necessary in order to more clearly delineate a barcoding gap in these taxa.
The incongruous results of different methods in many empirical studies emphasizes the importance of using multiple species delimitation methods for single-locus data; using one line of evidence might serve only to increase the confusion regarding species boundaries. While one method might be advantageous due to performance or speed, we agree with Dellicour and Flot (2018) that one should compare the results of several approaches when attempting to draw conclusions from single-locus data. Similarly, researchers should carefully consider the models of nucleotide substitution used, as well as the effect of reducing their datasets down to unique haplotypes has on their analyses and inferences. Here, we demonstrate that the same method applied to a full dataset can return different results compared to one consisting only of non-redundant haplotypes, in some cases overlumping or oversplitting lineages.
Analyses of a single mitochondrial locus should not be the sole line of evidence used to draw taxonomic conclusions (DeSalle, 2006), just as using any single character to dictate systematic changes is ill-advised (Will et al., 2005). Analyzing a single gene barcode as a first step in identifying previously overlooked lineages, however, complements taxonomic investigation, providing a roadmap for groups of taxa in need of more complete study. Furthermore, barcoding analyses can be useful in identifying areas of population structure and potential areas where gene flow is not occurring between populations (Tavares and Baker, 2008). In poorly represented taxa, a single gene can test the efficacy of groupings based solely on morphology.
Across all analyses, no method inferred species limits completely congruent with current taxonomy, even under the most conservative considerations. Based on the performance of the multiple methods on both datasets, we consider ABGD at a 7.2% threshold to be the most reasonable distance-based estimate and mPTP to be the most reasonable tree-based methods of delimitation for this dataset; these methods inferred 36 and 46 lineages, respectively, and were congruent between full and reduced datasets. These estimates suggest that species diversity in Chortís Block anoles could be underestimated by as much as 26%. The recently described N. caceresae (Hofmann and Townsend, 2018) was supported by every analysis. Several recognized species were inferred by these methods to represent two or more OTUs, each of which appear to represent reciprocally monophyletic populations often separated biogeographically. These methods congruently inferred multiple geographically isolated species level lineages within Norops mccraniei, N. morazani, N. rubribarbaris, and N. yoroensis. mPTP additionally inferred multiple lineages within N. biporcatus, N. heteropholidotus, N. laeviventris, N. lemurinus, N. rodriguezii, and N. uniformis. Conversely, ABGD at a 7.2% threshold inferred N. limifrons and N. zeus as the same lineage, but as distinct lineages at a 5.2% threshold; N. unilobatus and N. wellbornae were similarly considered a single lineage by this method. These results suggest that further taxonomic investigation into these lineages, including more thorough molecular sampling, is warranted.
In this study, we utilized a COI-barcode dataset of mainland anoles to test inferences of several commonly utilized single-locus delimitation methods. Our results provide further support for the necessity of utilizing multiple methods in barcoding studies, as results differed between tree- and distance-based analyses, as well as between the same analyses applied to a full dataset and one with the redundant haplotypes removed. Additionally, we highlight previously unrecognized diversity in need of further taxonomic and systematic investigation. Finally, a goal of this study was to provide a foundation for further studies of Chortís Block anoles by providing a reference library of barcodes to allow for the efficient identification of difficult-to-assign specimens, particularly females and juveniles. While a tremendous amount of work remains to be done if our understanding of mainland anoles will ever rival that of Caribbean anoles, the stage is set for a generation of researchers to undertake that effort.
Data Availability Statement
The sequence data generated for this study has been deposited in GenBank. See Supplementary Table 1 for accession numbers.
This study was carried out in accordance with the recommendations of Indiana University of Pennsylvania Institutional Animal Use and Care Committee (IUP-IACUC) and the University of Florida IFAS Animal Research Committee (UF-ARC). Protocols were approved by the IUP-IACUC and UF-ARC.
EH and JT designed the study. All authors collected and contributed samples. EH, IL-M, MM-F, and JT generated the molecular data. EH analyzed the data. EH and JT wrote the first draft of the manuscript. All authors reviewed, edited, and contributed to the final manuscript.
Results included in this paper were obtained with support from the National Science Foundation (DEB-0949359, to KN), Critical Ecosystem Partnership Fund (to JT), Indiana University of Pennsylvania (IUP) Department of Biology (to EH and JT), IUP College of Natural Sciences and Mathematics (to JT), IUP School of Graduate Studies and Research (to EH), IUP Faculty Senate (to JT), Commonwealth of Pennsylvania University Biologists (to EH), and a Pennsylvania State System of Higher Education Faculty Professional Development Grant (to JT).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Fieldwork in Honduras was carried out under a series of research permits issued by the Instituto Nacional de Conservación y Desarrollo Forestal, Áreas Protegidas y Vida Silvestre [ICF] (Resolución GG-MP-055-2006 and Dictamen DAPVS 0091-2006; Resolución DE-MP-086-2010 and Dictamen DVS-ICF-045-2010; Resolución DE-MP-095-2014 and Dictamen ICF-DVS-112-2014), and in Nicaragua under research permits issued in 2007, 2008, 2013, and 2016 by the Ministerio del Ambiente y los Recursos Naturales (MARENA), Managua, Nicaragua. Logistical support in Honduras was provided in 2008, 2009, and 2015 by the Centro Zamorano de Biodiversidad, and we thank Jorge Iván Restrepo and José Mora (2008–2009) and Oliver Komar and Karla Lara (2015) for facilitating our fieldwork. We thank Carol L. Spencer, Jimmy A. McGuire, and David B. Wake (MVZ) and Rafe M. Brown, Rich E. Glor, and Luke J. Welton (KU) for tissue loans from their respective institutions, and are grateful to them, as well as Steve Rogers and Kaylin Martin (Carnegie Museum of Natural History [CM]) and Roy W. McDiarmid and Steve Gotte (National Museum of Natural History [USNM]) for facilitating timely deposition and accession of specimens associated with this study. Amy Driskell and Dan Mulcahy (Smithsonian Institution Laboratory of Analytical Biology) generated some COI sequence data, in collaboration with Roy McDiarmid (USNM), as part of the “Barcoding the Herpetofauna of Eastern Nuclear Central America” project. We are grateful for field assistance provided by Jason M. Butler, Thomas J. Firneno, Luis Herrera, Alexander Hess, Michael Itgen, Mariah Kenney, Catherine Krygeris, Lesster Munguía, Jorge Luis Murillo, Lenin Obando, Fatima Pereira-Pereira, Sandy Pereira-Pereira, Javier Sunyer, Scott Travers, Josué Vasquez, Hermes Vega, and Kayla D. Weinfurther. We also thank Waldir M. Berbel-Filho and Celia A. May whose comments and suggestions significantly improved this manuscript.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00011/full#supplementary-material
Ahrens, D., Fujisawa, T., Krammer, H. J., Eberle, J., Fabrizi, S., and Vogler, A. P. (2016). Rarity and incomplete sampling in DNA-based species delimitation. Syst. Biol. 65, 478–494. doi: 10.1093/sysbio/syw002
Barco, A., Raupach, M. J., Laakmaan, S., Neumann, H., and Knebelsberger, T. (2016). Identification of North Sea molluscs with DNA barcoding. Mol. Ecol. Resour. 16, 288–297. doi: 10.1111/1755-0998.12440
Barman, A. S., Singh, M., Singh, K. S., Saha, H., Singh, Y. J., Laishram, M., et al. (2018). DNA barcoding of freshwater fishes of Indo-Myanmar biodiversity hotspot. Sci. Rep. 8:8579. doi: 10.1038/s41598-018-26976-3
Bickford, D., Lohman, D. J., Sodhi, N. S., Ng, P. K. L., Meier, R., Winker, K., et al. (2007). Cryptic species as a window on diversity and conservation. Trends Ecol. Evol. 22, 148–155. doi: 10.1016/j.tree.2006.11.004
Blair, C., and Bryson, R. W. Jr. (2017). Cryptic diversity and discordance in single-locus species delimitation methods within horned lizards (Phrynosomatidae: Phrynosoma). Mol. Ecol. Resour. 17, 1168–1182. doi: 10.1111/1755-0998.12658
Boissin, E., Hoareau, T. B., Paulay, G., and Bruggemann, J. H. (2017). DNA barcoding of reef brittle stars (Ophiuroidea, Echinodermata) from the southwestern Indian Ocean evolutionary hot spot of biodiversity. Ecol. Evol. 7, 11197–11203. doi: 10.1002/ece3.3554
Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C-H., Xie, D., et al. (2014). BEAST 2: a software platform for bayesian evolutionary analysis. PLoS Comput. Biol. 10:e1003537. doi: 10.1371/journal.pcbi.1003537
Campbell, J. A. (1999). “Distribution patterns of amphibians in Middle America,” in Patterns of Distribution of Amphibians: A Global Perspective, ed. W. D. Duellman (Baltimore, MD: The Johns Hopkins University Press), 111–120.
Dang, N. X., Sun, F. H., Lv, Y. Y., Zhao, B. H., Wang, J. C., Murphy, R. W., et al. (2015). DNA barcoding and the identification of tree frogs (Amphibia: Anura: Rhacophoridae). Mitochondr. DNA Pt. A. 27, 2574–2584. doi: 10.3109/19401736.2015.1041113
Deichmann, J. L., Mulcahy, D. G., Vanthomme, H., Tobi, E., Wynn, A. H., Zimkus, B. M., et al. (2017). How many species and under what names? Using DNA barcoding and GenBank data for west Central African amphibian conservation. PLoS One 12:e0187283. doi: 10.1371/journal.pone.0187283
Dellicour, S., and Flot, J.-F. (2015). Delimiting species-poor data sets using single molecular markers: a study of barcode gaps, haplowebs, and GMYC. Syst. Biol. 64, 900–908. doi: 10.1093/sysbio/syu130
Esselstyn, J. A., Evans, B. J., Sedlock, J. L., Khan, F. A. A., and Heaney, L. R. (2012). Single-locus species delimitation: a test of the mixed Yule-coalescent model, with an empirical application to Philippine round-leaf bats. Proc. R. Soc. B-Biol. Sci. 279, 3678–3686. doi: 10.1098/rspb.2012.0705
Ezard, T., Fujisawa, T., and Barraclough, T. G. (2009). SPLITS: SPecies’ Limits by Threshold Statistics. R Package Version, 1. Available at: http://r-forge.r-project.org/projects/splits/
Fujisawa, T., and Barraclough, T. G. (2013). Delimiting species using single-locus data and the generalized mixed yule coalescent approach: a revised method and evaluation on simulated data sets. Syst. Biol. 62, 702–724. doi: 10.1093/sysbio/syt033
Greenbaum, E., and Komar, O. (2010). “A conservation assessment of Salvadoran protected areas: priorities and recommendations based on amphibian and reptile distributions,” in Conservation of Mesoamerican Amphibians and Reptiles, eds L. D. Wilson, J. H. Townsend, and J. D. Johnson (Eagle Mountain, UT: Eagle Mountain Publishing LC), 437–459.
Guarnizo, C. E., Paz, A., Muñoz-Ortiz, A., Flechas, S. V., Méndez-Narváez, J., and Crawford, A. J. (2015). DNA barcoding survey of anurans across the eastern Cordillera of Colombia and the impact of the Andes on cryptic diversity. PLoS One 10:e0127312. doi: 10.1371/journal.pone.0127312
Hajibabaei, M., Singer, G. A., Hebert, P. D., and Hickey, D. A. (2007). DNA barcoding: how it complements taxonomy, molecular phylogenetics and population genetics. Trends Genet. 23, 167–172. doi: 10.1016/j.tig.2007.02.001
Hajibabaei, M., Smith, M. A., Janzen, D. H., Rodriguez, J. J., Whitfield, J. B., and Hebert, P. N. (2006). A minimalist barcode can identify a specimen whose DNA is degraded. Mol. Ecol. Notes 6, 959–964. doi: 10.1111/j.1471-8286.2006.01470.x
Hawlitscheck, O., Morinière, J., Dunz, A., Franzen, M., Rödder, D., Glaw, F., et al. (2016). Comprehensive DNA barcoding of the herpetofauna of Germany. Mol. Ecol. Resour. 16, 242–253. doi: 10.1111/1755-0998.12416
Hofmann, E. P., and Townsend, J. H. (2017). Origins and biogeography of the Anolis crassulus subgroup (Squamata: Dactyloidae) in the highlands of Nuclear Central America. BMC Evol. Biol. 17:267. doi: 10.1186/s12862-017-1115-8
Hofmann, E. P., and Townsend, J. H. (2018). A cryptic new species of anole (Squamata: Dactyloidae) from the Lenca Highlands of Honduras, previously referred to as Norops crassulus (Cope, 1864). Ann. Carnegie Mus. 85, 91–111.
IUCN (2017). The IUCN Red List of Threatened Species. Version 2017-02. Available at: http://www.iucnredlist.org
Jin, Q., Hu, X. M., Han, H. L., Chen, F., Cai, W. J., Ruan, Q. Q., et al. (2018). A two-step DNA barcoding approach for delimiting moth species: moths of Dongling Mountain (Beijing, China) as a case study. Sci. Rep. 8:14256. doi: 10.1038/s41598-018-32123-9
Kapli, P., Lutteropp, S., Zhang, J., Kobert, K., Pavlidis, P., Stamatakis, A., et al. (2017). Multi-rate poisson tree processes for single-locus species delimitation under maximum likelihood and markov chain monte carlo. Bioinformatics 33, 1630–1638.
Kearse, M., Moir, R., Wilson, A., Stones-Havas, S., Cheung, M., Sturrock, S., et al. (2012). Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649. doi: 10.1093/bioinformatics/bts199
Köhler, G., Townsend, J. H., and Peterson, C. B. P. (2016). A taxonomic revision of the Norops tropidonotus complex (Squamata, Dactyloidae), with the resurrection of N. spilorhipis (Álvarez del Toro and Smith, 1956) and the description of two new species. Mesoam. Herpetol. 3, 7–41.
Korniolios, P., Kumlataş, Y., Lymberakis, P., and Ilgaz, C. (2018). Cryptic diversity and molecular systemtatics of the Aegean Ophiomorus skinks (Reptilia: Squamata), with the description of a new species. J. Zool. Syst. Evol. Res. 56, 364–381. doi: 10.1111/jzs.12205
Liu, Q., Zhu, F., Zhong, G., Wang, Y., Fang, M., Xiao, R., et al. (2015). COI–based barcoding of Chinese vipers (Reptilia: Squamata: Viperidae). Amphibia-Reptilia 36, 361–372. doi: 10.1163/15685381-00003012
Luo, A., Ling, C., Ho, S. Y. W., and Zhu, C.-D. (2018). Comparison of methods for molecular species delimitation across a range of speciation scenarios. Syst. Biol. 67, 830–846. doi: 10.1093/sysbio/syy011
Mendoza, A. M., Torres, M. F., Paz, A., Trujillo-Arias, N., López-Alvarez, D., Sierra, S., et al. (2016). Cryptic diversity revealed by DNA barcoding in Colombian illegally traded bird species. Mol. Ecol. Resour. 16, 862–873. doi: 10.1111/1755-0998.12515
Miralles, A., and Vences, M. (2013). New metrics for comparison of taxonomies reveal striking discrepancies among species delimitation methods in Madascincus lizards. PLoS One 8:e68242. doi: 10.1371/journal.pone.0068242
Moura, C. J., Cunha, M. R., Porteiro, F. M., and Rogers, A. D. (2011). The use of the DNA barcode gene 16S mRNA for the clarification of taxonomic problems within the family Sertulariidae (Cnidaria, Hydrozoa). Zool. Scr. 40, 520–537. doi: 10.1111/j.1463-6409.2011.00489.x
Murphy, R. W., Crawford, A. J., Bauer, A. M., Che, J., Donnellan, S. C., Fritz, U., et al. (2013). Cold code: the global initiative to DNA barcode amphibians and noavian reptiles. Mol. Ecol. Resour. 13, 161–167. doi: 10.1111/1755-0998.12050
Nagy, Z. T., Sonet, G., Glaw, F., and Vences, M. (2012). First large-scale DNA barcoding assessment of reptiles in the biodiversity hotspot of Madagascar, based on newly designed COI primers. PLoS One 7:e34506. doi: 10.1371/journal.pone.0034506
Nicholson, K. E. (2002). Phylogenetic analysis and a test of the current infrageneric classification of Norops (beta Anolis). Herpetol. Monogr. 16, 93–120. doi: 10.1655/0733-1347(2002)016[0093:PAAATO]2.0.CO;2
Nicholson, K. E., Crother, B. I., Guyer, C., and Savage, J. M. (2018). Translating a clade based classification into one that is valid under the international code of zoological nomenclature: the case of the lizards of the family Dactyloidae (Order Squamata). Zootaxa 4461, 573–586. doi: 10.11646/zootaxa.4461.4.7
Padial, J. M., Castroviejo-Fisher, S., Köhler, J., Vilá, C., Chaparro, J. C., and De la Riva, I. (2009). Deciphering the products of evolution at the species level: the need for an integrative taxonomy. Zool. Scr. 38, 431–447. doi: 10.1111/j.1463-6409.2008.00381.x
Poe, S., Nieto-Montes de Oca, A., Torres-Carvajal, O., de Queiroz, K., Velasco, J. A., Truett, B., et al. (2017). A phylogenetic, biogeographic, and taxonomic study of all extant species of Anolis (Squamata, Iguanidae). Syst. Biol. 66, 663–697. doi: 10.1093/sysbio/syx029
Pons, J., Barraclough, T. G., Gomez-Zurita, J., Cardoso, A., Duran, D. P., Hazell, S., et al. (2006). Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Sys. Biol. 55, 595–609. doi: 10.1080/10635150600852011
Puillandre, N., Lambert, A., Brouillet, S., and Achaz, G. (2012). ABGD, automatic barcode gap discovery for primary species delimitation. Mol. Ecol. 21, 864–1877. doi: 10.1111/j.1365-294X.2011.05239.x
R Core Team (2015). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available at: https://www.R-project.org/
Rambaut, A., Suchard, M. A., Xie, D., and Drummond, A. J. (2014). Tracer v.1.6. Available at: http://tree.bio.ed.ac.uk/software/tracer/
Ramirez, J. L., Birindelli, J. L., Carvalho, D. C., Affonso, P. R. A. M., Venere, P. C., Ortega, H., et al. (2017). Revealing hidden diversity of the underestimated neotropical ichthyofauna: DNA barcoding in the recently described genus Megalepoinus (Characiformes: Anostomidae). Front. Genet. 8:149. doi: 10.3389/fgene.2017.00149
Ratnasingham, S., and Hebert, P. D. N. (2007). BOLD: the barcode of life data system (www.barcodinglife.org). Mol. Ecol. Notes 7, 355–364. doi: 10.1111/j.1471-8286.2007.01678.x
Sonet, G., Jordaens, K., Nagy, Z. T., Breman, F. C., De Meyer, M., Backeljau, T., et al. (2013). Adhoc: an R package to calculate ad hoc distance thresholds for DNA barcoding identification. ZooKeys 365, 329–336. doi: 10.3897/zookeys.365.6034
Sunyer, J., and Köhler, G. (2010). “Conservation status of the herpetofauna of Nicaragua”, in Conservation of Mesoamerican Amphibians and Reptiles, eds L. D. Wilson, J. H. Townsend, and J. D. Johnson (Eagle Mountain, UT: Eagle Mountain Publishing LC), 489–509.
Talavera, G., Dincă, V., and Vila, R. (2013). Factors affecting species delimitations with the GMYC model: insights from a butterfly survey. Methods Ecol. Evol. 4, 1101–1110. doi: 10.1111/2041-210X.12107
Tang, C. Q., Humphreys, A. M., Fontaneto, D., and Barraclough, T. G. (2014). Effects of phylogenetic reconstruction method on the robustness of species delimitation using single-locus data. Methods Ecol. Evol. 5, 1086–1094. doi: 10.1111/2041-210X.12246
Thompson, J. D., Higgins, D. G., and Gibson, T. J. (1994). CLUSTALW: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties, and weight matrix choice. Nucleic Acids Res. 22, 4673–4680. doi: 10.1093/nar/22.22.4673
Townsend, J. H. (2009). Morphological variation in Geophis nephodrymus (Squamata: Colubridae), with comments on the conservation of Geophis in eastern Nuclear Central America. Herpetologica 65, 292–302. doi: 10.1655/07-039R2.1
Townsend, J. H., and Wilson, L. D. (2010). “Conservation of the Honduran herpetofauna: issues and imperatives”, in Conservation of Mesoamerican Amphibians and Reptiles, eds L. D. Wilson, J. H. Townsend, and J. D. Johnson (Eagle Mountain, UT: Eagle Mountain Publishing LC), 461–486.
Vasconcelos, R., Montero-Mendieta, S., Simó-Riudalbas, M., Sindaco, R., Santos, X., Fasola, M., et al. (2016). Unexpectedly high levels of cryptic diversity uncovered by a complete DNA barcoding of reptiles of the Socotra Archipelago. PLoS One 11:e0149985. doi: 10.1371/journal.pone.0149985
Wilson, L. D., and Townsend, J. H. (2007). Biogeography and conservation of the herpetofauna of the upland pine-oak forests of honduras. Biota Neotrop. 7, 131–142. doi: 10.1590/S1676-06032007000100018
Xia, X., and Lemey, P. (2009). “Assessing substitution saturation with DAMBE”, in The Phylogenetic Handbook: A Practical Approach to DNA and Protein Phylogeny, 2nd Edn, eds P. Lemey, M. Salemi, and A. M. Vandamme (Cambridge: Cambridge University Press), 615–630. doi: 10.1017/CBO9780511819049.022
Zhang, J., Kapli, P., Pavlidis, P., and Stamatakis, A. (2013). A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29, 2869–2876. doi: 10.1093/bioinformatics/btt499
Keywords: ABGD, Anolis, cytochrome c oxidase subunit I (COI), DNA barcoding, GMYC, Norops, PTP, RESL
Citation: Hofmann EP, Nicholson KE, Luque-Montes IR, Köhler G, Cerrato-Mendoza CA, Medina-Flores M, Wilson LD and Townsend JH (2019) Cryptic Diversity, but to What Extent? Discordance Between Single-Locus Species Delimitation Methods Within Mainland Anoles (Squamata: Dactyloidae) of Northern Central America. Front. Genet. 10:11. doi: 10.3389/fgene.2019.00011
Received: 11 June 2018; Accepted: 11 January 2019;
Published: 11 February 2019.
Edited by:Edward Hollox, University of Leicester, United Kingdom
Reviewed by:Waldir M. Berbel-Filho, Swansea University, United Kingdom
Celia A. May, University of Leicester, United Kingdom
Copyright © 2019 Hofmann, Nicholson, Luque-Montes, Köhler, Cerrato-Mendoza, Medina-Flores, Wilson and Townsend. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Josiah H. Townsend, email@example.com