Assessing the effectiveness of molecular data for species identification and diversity studies of Nigerian herpetofauna

Abstract

We evaluated the effectiveness of molecular data (mitochondrial cytochrome c oxidase subunit I – COI and 16S ribosomal RNA – 16S) for species identification and diversity studies of Nigerian herpetofauna. We collected 189 specimens consisting of 91 amphibians and 98 reptiles from the largest and under-studied national park in Nigeria – Gashaka Gumti National Park (GGNP), eastern Nigeria. The resulting COI and 16S data allowed us to compare our sequences with others in global databases for species identification. During query searches, we observed poor representation of COI sequences for African herpetofauna in the GenBank and BOLD databases. Neighbour-joining (NJ) trees constructed based on COI and 16S datasets for all specimens sequenced for this study clustered most individuals of the same species. On the other hand, morphologically unidentified amphibian (Leptopelis sp and Ptychadena sp) and reptile (Hemidactylus sp and Lycophidion sp) species formed strongly supported lineages in the NJ trees, and were clearly separated from their sister species. Further, the NJ tree-based analyses recovered at least two well supported divergent lineages within Agama agama and Trachylepis maculilabris. Our study confirms the efficiency of integrating morphological and DNA barcoding approaches in identification and diversity studies of Nigerian herpetofauna. Finally, we recommend development of comprehensive and reliable DNA barcode reference libraries for global amphibians and reptiles.

Keywords:

• Integrative taxonomy
• DNA barcoding
• biodiversity
• amphibians
• reptiles
• Nigeria

1. Introduction

Biodiversity inventory and documentation are one of the leading concepts in conservation biology. However, this is hinged closely to accurate species identification. In last two decades, molecular data (mitochondrial DNA) has gained popularity as a complementary tool for species identification. This involves the use of mitochondrial cytochrome c oxidase subunit I (COI) gene fragment and in some cases, 16S ribosomal RNA (16S) as an alternative gene marker in systematics (Hebert et al. 2003; Vences et al. 2005, 2012, Rockney et al. 2015; Nazarov et al. 2012). This method has proven fast and reliable in identification of species across life stages and unveiling cryptic diversities (Smith et al. 2008; Crawford et al. 2011; Nagy et al. 2012).

Despite the accrued global importance of mitochondrial DNA (COI and 16S) in taxonomy, fewer to no studies have extended this approach in identification of herpetofauna from Nigeria. Nigeria is a country in West Africa with high diversity of amphibians and reptiles (collectively referred to as herpetofauna). Herpetofauna play vital roles in ecosystems as predators, food for certain animal species, herbivores and commensal taxa (Wyman 1998; Hopkins 2007; Bohm et al. 2013). Anthropogenic activities, such as deforestation, bush burning etc., pose threat to its herpetofauna (Nneji et al. 2018). This calls for concerted approaches to document herpetofaunal diversity for initiation of conservation plans. Our study, therefore, aims to evaluate the effectiveness of mitochondrial DNA sequences (COI and 16S genes) in identification of Nigerian herpetofauna and detection of possible hidden species diversity.

2. Materials and methods

2.1. Permits, sample collection and morphological identification

Sample collection followed ethical approval from National Park Service, Nigeria. Herpetofauna were collected from Gashaka Gumti National Park (GGNP) in eastern Nigeria (Supplementary Appendix 1) on two occasions: 21 August–25 September 2017 and 12–16 February 2018. GGNP is located between latitude 06°58^′–08°05^′N and longitude 11°10^′N–12°13^′E (Dunn 1999). We collected 189 specimens including 91 amphibian and 98 reptile specimens (Supplementary Appendix 2) with the aid of visual encounter surveys and opportunistic observations. Tissue samples (toe clips) were collected and preserved in 95% ethanol (Sigma-Aldreich, Germany). Subsequently, some animals were released at the site of collection. About 115 voucher specimens were euthanized humanely, tissues (liver/muscle) sampled and subsequently stored in 95% ethanol. Thereafter, specimens were fixed with 4% formalin (Sigma-Aldreich, Germany) and transferred to 70% ethanol. Voucher specimens were deposited in the museums of Department of Zoology, University of Ibadan and GGNP, Nigeria. Ethanol preserved-tissues were kept frozen for further molecular studies.

Preliminary identification of specimens was based on morphology. Further, we referred to taxonomic papers and books on African amphibians and reptiles (Romer 1953; Schi⊘tz 1963; Walker 1967; Dunger 1967a, 1967b; 1968; Walker 1968; Perret and Amiet 1971; Hoogmoed 1974; Perret 1977; Hughes 1983; Butler and Reid 1990; Reid et al. 1990; Schi⊘tz 1999; Rödel 2000; Mediannikov et al. 2012) for clarification of species-level assignments. Specimens with uncertain species-level identification were referred to its genus. Species nomenclatures were in accordance with the Amphibian Web (2018) and The Reptile Database (Uetz et al. 2018).

2.2. Genomic DNA extraction, PCR amplification and sequencing

Total genomic DNA was isolated using standard phenol-chloroform method (Sambrook and Russell 2001). Primer sets used for Polymerase Chain Reaction (PCR) amplification and Sanger sequencing (Bossuyt and Milinkovitch, 2000; Che et al. 2012) are shown in Supplementary Appendix 3. The amplification was performed in 25 µl volume reaction that contained 1.5 µl of genomic working DNA, 18.5 µl of PCR water, 2.5 µl of PCR buffer, 2 µl of dNTP, 1 µl of each of the forward and reverse primers (10 pm/µl) and 0.30 µl of rTaq polymerase. The PCR cycle profile consisted of 5 min initial denaturation at 94 °C, followed by 35 cycles of 1 min at 94 °C, annealing for 45 s at 46 °C (COI) or 55 °C (16S), extension for 1 min at 72 °C; final extension for 10 m at 72 °C. The PCR products were checked on 1.2% agarose gel. DNA sequencing of the purified PCR products was performed in both directions using the primers described previously with an automated DNA sequencer (ABI 3730, DNA Analyzer; Applied Biosystems, Foster City, CA).

2.3. Genetic analyses

The nucleotide sequences were viewed and confirmed by eye using SeqManTMII (DNASTAR Lasergene 7) and aligned in MEGA 6.0 using ClustalW (Tamura et al. 2013) with default parameters. The aligned COI sequences were translated to amino acids to check for premature stop codons and to confirm that the open reading frame was maintained in the protein-coding locus.

We attempted to use variety of methods for analyzing our molecular data, including BLAST searches for sequence similarity, the Barcode Index Number (BIN) system, neighbour-joining (NJ) trees and percentage sequence divergence. To confirm the identity of amplified sequences, we conducted BLAST searches on GenBank. However, the unavailability of COI sequences in the database limited its use for species identification. For instance, BLAST searches in GenBank (November 2018) for ‘anurans, Africa, COI’ and ‘reptiles, Africa, COI’ revealed 263 (128 of which were from Xenopus) and 351 (182 of which were from Acontias meleagris) COI records for amphibians and reptiles, respectively. Thus, we failed to use BIN system for our COI data for the species identification. Given the large availability of 16S sequences for African herpetofauna (2905 for anuran and 3022 reptile sequences), we used the data for molecular identification of specimens. After the BLAST search, we downloaded all the related 16S sequences in the Genbank and conducted NJ tree searches for each of the morphologically identified species as implemented in MEGA 6.0 based on the Kimura 2 parameter distance (K2P) (Kimura, 1980). Our identification using 16S NJ trees were based on criteria set by Deichmann et al. (2017) as follows: (1) specimens were identified as the same species if our sequences clustered with same individuals from GenBank; (2) specimens were treated as geographic variation of the same species if specimens were morphologically identified as the same but in different clades (population divergence); (3) specimens were treated as different species if they do not fit morphological description of that species and also showed substantial genetic variation; and (4) specimens were treated as different species even when identified morphologically as the same species but in different clades, sister to other species. Following our specimen identification, we used MEGA v. 6.0 to create NJ trees based on the K2P distance (K2P) and estimated inter- and intraspecific sequence divergences for 16S and COI DNA barcode data for all specimens sequenced for this study.

3. Results and discussion

We provide in Supplementary Appendix 4, the sequence information of 16S and COI gene fragments. Two species of reptiles (Trachylepis affinis and Boaedon lineatus) failed to amplify for COI fragment even after repeated attempts. The amplified COI sequences were without gaps or stop codons. We deposited nucleotide sequences in GenBank with accession number MH708899-MH709081 (16S sequences) and MH700787-MH700946 (COI sequences; Supplementary Appendix 2). The associated 16S NJ trees of our specimens alongside others in GenBank can be found in Supplementary Appendix 5–6.

Using GenBank BLAST tool to confirm identity of the amplified sequences, 81 (83.50%) 16S sequences showed sequence similarity of ≥97% corresponding to 10 morphologically identified amphibian species (66.70%). For reptiles, 74 (80.40%) 16S sequences belonging to 11 species (78.60%) showed sequence similarity of ≥97%. COI sequences yielded low sequence query successes due to poor representation of COI sequences of African herpetofauna in global databases. Given to this, we failed to use BIN system in BOLD database for the species identification. Our study, therefore, shows that GenBank’s BLAST tool for comparing DNA sequence similarity supports the use of 16S gene marker for confirming identities of amplified sequences for Nigerian herpetofauna due to its wide availability in GenBank. It is obvious that the use of molecular data in species identification would not be effective if most African herpetofauna are not well documented in sequence libraries. Thus, efforts in building comprehensive DNA barcode libraries for Nigerian as well as African herpetofauna are recommended.

Anurans morphologically identified as Hoplobatrachus occipitalis, Amnirana galamensis, Sclerophrys maculata, Sclerophrys regularis, Kassina senegalensis, Ptychadena bibroni, and Xenopus fischbergi formed 16S clade with others of same species in the GenBank (Supplementary Appendix 5). Following this, we maintain the species names as valid for the specimens. Our NJ tree analysis showed that morphologically unidentified Sclerophrys sp clustered with S. maculata from GenBank including our newly collected S. maculata from GGNP. Therefore, we refer our specimen (Sclerophrys sp) as S. maculata. Individuals morphologically identified as Kassina cassinoides formed a distinct clade and we observed no records of K. cassinoides in the Genbank, thus, refer our specimen as K. cassinoides pending further morphological comparisons.

The phenetic analyses (Supplementary Appendix 5) showed that individual morphologically identified as Hyperolius concolor formed a sister clade to Hyperolius chlorosteus (GU443986 and FJ594076) and Hyperolius laurenti (GU443897). Given this, we refer our specimen as H. cf. concolor pending further taxonomical identification. Further, morphologically unidentified Hyperolius sp from GGNP clustered with Hyperolius nitidulus and Hyperolius viridiflavus nitidulus but are likely part of Hyperolius viridiflavus complex. Pending further revision, we refer our specimen as H. cf. viridiflavus.

Additionally, individuals morphologically identified as Leptopelis aubryi and Ptychadena sp formed sister clades to Leptopelis viridis (KT976110-KT976112) and Ptychadena sp. Mikumi (DQ525944) respectively (Supplementary Appendix 5). The exact species-level identification of our specimen (L. aubryi and Ptychadena sp) cannot be ascertained and thus, we refer it to the genus level as Leptopelis sp and Ptychadena sp pending further morphological and genetic identification.

Morphologically identified Ptychadena longirostris formed a sub-clade and sister clade within P. longirostris group that comprised other individuals from Cote D'Ivoire (Supplementary Appendix 5). For this reason, we refer our specimen as P. longirostris and highlight possibility of population divergence within this species group. Individuals identified as Ptychadena pumilio formed a clade within P. pumilio complex comprising Ptychadena taenioscelis and P. pumilio. Our data, consistent with Deichmann et al. (2017), suggest that P. pumilio is a species complex in need for further revision, and pending this, we identify our species as P. pumilio.

Our NJ tree analyses showed that reptile species morphologically identified as Agama agama, Agama parafricana, Psammophis lineatus and Varanus niloticus formed 16S clade with others of same species in the GenBank, and thus, we maintain species names as valid (Supplementary Appendix 6). On the other hand, individuals morphologically identified as Boaedon lineatus clustered with Boaedon paralineatus in the GenBank, but are likely part of Boaedon lineatus complex. However, we refer our specimen as ‘B. paralineatus’ pending further work is done on this group. Individuals identified as Boaedon capensis clustered with Boaedon perisilvestris from Congo and Chad. The type locality of B. perisilvestris is Congo, thus, we refer our individuals as B. perisilvestris.

Further, individuals morphologically identified as Chamaeleo gracilis formed a distinct sub-clade in a clade comprising C. gracilis, Chamaeleo necasi and Chamaeleo laevigatus. Although we refer our specimen as C. gracilis, we note the need for a detailed comparison of C. gracilis, C. necasi and C. laevigatus including the type material, to determine true identity and species identification of these individuals.

Further, individuals morphologically identified as Trachylepis affinis, Trachylepis perrotetii, Trachylepis maculilabris and Trachylepis quinquetianata showed significant structures (sub-clades) within clades comprising respective same species from GenBank (Supplementary Appendix 6). We, therefore, retain species names assigned to the individuals as valid and highlight the need for further investigation of processes driving the intraspecific population divergence. We also note that individuals morphologically identified as Hemidactylus sp and Lycophidion sp formed distinct 16S clades. Thus, we refer these individuals as Hemidactylus sp and Lycophidion sp, and highlight the need for detailed genetic and molecular comparisons to ascertain if these unknown individuals represent new species or unidentified species.

In sum, our study precludes the use of either morphology or molecular-based approach alone in identification of herpetofauna. Combined use of both approaches would enable accurate species identification.

From the NJ tree analyses (Figures 1–2) for all specimens sequenced for this study, individuals of the same species clustered together. The clustering pattern is in accordance with morphological and genetic identification, enabling the differentiation of individual species. Most of the species clusters in 16S and COI NJ trees were strongly supported (Figures 1–2). Likewise, morphologically unidentified amphibian (Leptopelis sp and Ptychadena sp) and reptile (Hemidactylus sp and Lycophidion sp) species formed strongly supported lineages in the NJ trees, and were clearly separated from their sister species. Further, our NJ tree-based analyses provide evidence of at least two strongly supported divergent lineages within A. agama (Figures 1–2). In contrast to 16S NJ tree, our COI phenetic analysis revealed strongly supported divergent lineages within T. maculilabris (Figure 2). This finding suggests possible hidden diversities within A. agama and T. maculilabris. Using morphology alone, this diversity remained undetected due to morphological similarity within individuals belonging to these groups.

Figure 1. Unrooted neighbour-joining tree of Kimura-2-parameter distances based on (a) 16S (b) COI gene sequences of the amphibian species. Bootstrap values of greater than 50% (1000 pseudo-replicates) are given above the nodes. Nodal support less than 50% were not shown in the tree.

Figure 2. Unrooted neighbour-joining tree of Kimura-2-parameter distances based on (a) 16S (b) COI gene sequences of the reptile species. Bootstrap values of greater than 50% (1000 pseudo-replicates) are given above the nodes. Nodal support less than 50% were not shown in the tree.

The K2P interspecific divergence obtained for amphibians was in a range of 5.70 to 31.70% in 16S gene and 19.60 to 32.80% in COI gene (Supplementary Appendix 7). We observed high intraspecific divergences for COI (0.409%) and 16S (0.4175%) genes in K. cassinoides (Supplementary Appendix 8). In reptiles, we obtained K2P interspecific divergence ranged from 4.50 to 36.80% in 16S and 8.20 to 44.20% in COI respectively (Supplementary Appendix 9). However, the highest intraspecific divergences for COI (1.90%) and 16S (0.47%) genes were observed in A. agama (Supplementary Appendix 10).

Finally, our study illustrates the relevance of molecular data in species identification and unveiling hidden diversity within Nigerian herpetofauna. Our results largely support the integration of morphological and molecular data in taxonomy. In an ongoing effort to build comprehensive regional sequence libraries, this study adds new sequences to the global databases that would aid in identification of herpetofauna. Further, we recommend validation of sequences deposited in global databases with respect to its voucher specimen for accurate species identification.

Acknowledgements

We thank the management of Gashaka Gumti National Park, Taraba, Nigeria for providing technical assistance during our field survey.

Disclosure statement

The authors declare that they have no competing interests.

Additional information

Funding

This study was funded by Rufford Small Grant for Nature Conservation (22507-1) and IDEA WILD Nature Conservation Grant.