Evaluation of partial 12S rRNA, 16S rRNA, COI and Cytb gene sequence datasets for potential single DNA barcode for hylids (Anura: Hylidae)

Abstract

We evaluated the extent of intraspecific and interspecific genetic distances and the effectiveness of predefined threshold values using the main genes for estimates of biodiversity and specimens’ identification in anurans. Partial sequences of the mitochondrial genes for small (12S) and large (16S) ribosomal subunits, cytochrome c oxidase subunit I (COI) and Cytochrome b (Cytb) of the family Hylidae were downloaded from GenBank and curated for length, coverage, and potential contaminations. We performed analyses for all sequences of each gene and the same species present in these datasets by distance and tree (monophyly)-based evaluations. We also evaluated the ability to identify specimens using these datasets applying “nearest neighbor” (NN), “best close match” (BCM) and “BOLD ID” tests. Genetic distance thresholds were generated by the function ‘threshVal’ and “localMinima” from SPIDER package and traditional threshold values (1%, 3%, 6% and 10%) were also evaluated. Coding genes, especially COI, had a better identification capacity than non-coding genes on barcoding gap and monophyly analysis and NN, BCM, BOLD ID tests. Considering the multiple factors involved in global DNA barcoding evaluations, we present a critical assessment of the use of these genes for biodiversity estimation and specimens’ identification in anurans (e.g. hylids).

Key words
DNA barcoding; identification capacity; mitochondrial genes; molecular identification; SPIDER package

INTRODUCTION

DNA barcoding aims to use large-scale tracing of a reference gene to assign unknown individuals to species and increase the discovery of new species (Hebert et al. 2003). For animals, the use of the cytochrome c oxidase subunit I mitochondrial gene (COI) has been effective for several taxa, such as ground beetles of Germany (Raupach et al. 2018), flies and dragonflies of Brazil (Koroiva et al. 2017, 2018) and anurans in central South America (Koroiva et al. 2020). However, its effectiveness seems to vary depending on the group analyzed (Waugh 2007). In this sense, the use of other genes for identification has been suggested, such as the mitochondrial gene 16S ribosomal subunit (16S) (Vences et al. 2005b). One of the groups in which the best gene to use has remained ambiguous is amphibians.

Within Amphibia, anurans are the richest order and are considered the most threatened vertebrate group in the world; about 32% of their species are at risk of extinction (Stuart et al. 2004). As stated by Lyra et al. (2017), Anura has the potential to greatly benefit from the use of DNA barcoding once this approach provides information on populations that may be cryptic species and assists in taxonomic identification. These procecsses can be quite difficult—even for specialists—because both parallel and convergent evolution have led to an extremely conserved morphology.

Some genetic regions, especially within mitochondrial DNA, have been suggested for the molecular identification of anurans. Hebert et al. (2003), Che et al. (2012), and Murphy et al. (2013) defended the use of COI given its capability of identifying many species. Vences et al. (2005b) discusses that the priming sites of COI are highly variable in anurans and suggested that 16S may represent a better alternative because its universality and robustness of primers. Mitochondrial 12S rRNA (12S) and Cytochrome b (Cytb) genes are standard markers for phylogeny reconstruction in amphibians (Vences et al. 2005a) and have also been used to assist molecular identification of species. Nevertheless, COI and 16S are considered the main genes for the molecular identification of anuran species (Murphy et al. 2013).

Currently, many researchers use 16S as a DNA barcode marker for biodiversity estimates in amphibians. However, the use of COI has been increasing in recent years (e.g. Lehr et al. 2017, Zhao et al. 2017, Gao et al. 2019, Koroiva et al. 2020). As evidenced by Lyra et al. (2017), amphibian taxonomists have used different thresholds for COI and 16S to flag candidate species in the absence of a more resolved taxonomy. Vences et al. (2005a) provided an initial conservative threshold for the identification of potential candidate species from 5% divergence in 16s and 10% divergence in COI. Fouquet et al. (2007) suggested a threshold of 3% for the 16S marker for Neotropical frogs, which has become the standard value applied for the preliminary hypotheses of many candidate species in the amphibian integrative taxonomy (Vieites et al. 2009). Crawford et al. (2010) working in a more restricted geographical area, used an 8.6% limit for COI and more than 2% for 16S. Lyra et al. (2017) suggested 3% for 16S and 6% for COI using their own database of anurans. For 12S, the threshold values at species-level are commomly the same as those used in 16S (eg. Howlader et al. 2016), while in Cytb, the threshold values are those evaluated to large vertebrate datasets such as 5.52 ± 1.34 for sibling species and 10.69 ± 1.34 for species within a genus (Kartavtsev & Lee 2006). Despite the lack of specific studies, they are widely used in phylogeny studies and therefore have a large number of deposited sequences that make it a potential gene on delimiting species based on genetic distances.

Few studies have directly compared the relative efficacy of these molecular markers within a single group of anurans. In general, these studies compare only regional or local databases (e.g. Vences et al. 2005b, Lyra et al. 2017). The presence of a gap between the maximum intraspecific genetic distance value and lowest genetic distance to the nearest neighbor species within a study group—called the barcoding gap—is considered a key factor in the selection of a gene for the rapid identification of unknown samples (Meyer & Paulay, 2005). At regional levels, there are favorable analyses for both COI and 16S (e.g. Vences et al. 2005b, Lyra et al. 2017). Despite this, an overall association between the markers and even between the same species is still non-existent.

Therefore, in order to clarify the effectiveness of these mitochondrial genes for species identification in anurans, the present study aimed to comprehensively sample 12S, 16S COI and Cytb from online repository (GenBank). Here, we evaluate the presence or absence of a global barcoding gap among taxa and the ability to identify them as unique species using thresholds. Considering the pronounced genetic variation within the Anura order, we decided to analyze the Hylidae family. This family is one of the most speciose among the anurans (about 10% of all anurans described; 730 species) and is found worldwide (see Frost 2020). This kind of study is important for assessing the ability to identify anura species with molecular tools as well as their usefulness for biodiversity estimation and species differentiation. The other purpose of the study was to indicate the most favorable choices of genes and primers for taxon identification and high-throughput applications.

MATERIALS AND METHODS

Representatives of the 12S, 16S, COI and Cytb sequences were analyzed separately in order to evaluate each gene. All sequences were downloaded on April 09, 2020 from GenBank by sending the results (in FASTA format) of the following search queries to separate files: “txid8418 [Organism:exp] 12S [product]” and “txid8418 [Organism:exp] 12S ribosomal RNA [product]” for 12S rRNA gene; “txid8418 [Organism:exp] 16S [product]” and “txid8418 [Organism:exp] 16S ribosomal RNA [product]” for 16S rRNA gene; “txid8418 [Organism:exp] cytochrome oxidase subunit I [product]”, “txid8418 [Organism:exp] co1 [product]”, “txid8418 [Organism:exp] coi [product]”, for COI gene; “txid8418 [Organism:exp] cytb [product]” and “txid8418 [Organism:exp] cytochrome b [product]” for Cytochrome B gene. The duplicate records were removed. This resulted in 2,766, 5,534, 3,667 and 3,758 individual sequences for 12S, 16S, COI and Cytb, respectively (Table I).

Table I
Curation and alignment size of genetic markers in each dataset.

Gene	Download sequences	Number of removed sequences					Final dataset	Aligment size (bp)
		Duplicates	No Hylidae species (e.g Phyllomedusidae) and other genetic regions	Imprecise taxonomic labels	Less than 400 bp	Singletons
12s	5484	2718	323	194	453	106	1690	969
16s	10986	5452	1651	818	240	113	2712	658
COI	5558	1891	1284	501	98	47	1737	658
Cytb	6484	2726	312	28	981	12	2425	651

The datasets were further manipulated using Geneious v.9.0.5 (Kearse et al. 2012) in combination with Microsoft Excel™ (Microsoft Corporation, Redmond, WA). To enable comparisons of verified data and to ensure intraspecific and interspecific comparisons, all undeterminated sequences were removed; i.e., sequences with imprecise taxonomic labels (e.g. “sp.”, “cf”, “gr”, “aff”) were excluded from the alignements to be studied. Note that these instances only composed a small fraction of the full dataset (Table I). Valid names and taxonomic synonyms were checked against a standardized list of Frost (2020) on April 13, 2020.

Potential contaminants were then controlled by a BLASTn search of the genes’ sequences against the non-redundant sequence database on GenBank (Koroiva & Kvist 2018). In addition to the evaluation of all sequences (“All sequences” in the Results section), we also analyzed only the species that had sequences for 12S, 16S, COI and Cytb (“Species in common” in the Results section). Considering the limited number of species present in the Cytb dataset, we also performed comparing the database of the other three genes (“Species in common without Cytb” in the Results section), which have more three times the number of species in common. When comparing only datasets of common species between 12S and 16S, the results are close to the analysis of “All sequences” datasets (data not shown).

Target region and directionality

First, the 5’–3’ direction of the sequences was confirmed by comparison to sequences that were previously determined to be in the correct direction (KY385762 Dryophytes japonicus for COI, FJ784413 Smilisca phaeota for 16S, MG282247 Dryophytes immaculatus for 12S and AF205095 Dryophytes japonicus for Cytb). This was controlled through direction plots using MAFFT ver. 7.309 (Katoh & Standley 2013) implemented in Geneious v.9.0.5. Also, in order to increase robustness in the homology statement and elevate matrix occupancy, long sequences were truncated to cover only commonly used regions for DNA barcoding and phylogenetic studies. For hylids (Vences et al. 2012), the COI region can be amplified using the “Folmer” primer pair (HCO2198 and LCO1490) and the 16S rRNA region using the “Palumbi” primer pair (16SA-L and 16SB-H). The 12S and Cytb region were truncated to cover the region from the primer pairs t-Phe-frog/MVZ59 and 12Sb-H (Goebel et al. 1999) and L14850 and H15502 (Tanaka-Ueno et al. 1998), respectively. Several articles used these primers to amplify Hylidae species (e.g. Vences et al. 2012, Jeong et al. 2013, Orrico et al. 2017). The truncation was carried out using Geneious v.9.0.5, following the positioning of the aforementioned primers. Sequences of the four genes with a length below 400 basepairs (bp) were removed.

Genetic analysis

MAFFT ver. 7.309 was also used to align and realign sequences within each of the datasets separately by employing the FFT-NS-1 strategy with a gap opening penalty of 3.0, the 200PAM/K = 2 scoring matrix, and an offset value of 0.0. These alignments were inspected and manually refined. We use genetic distance to assess the ability to correctly identify and the effectiveness of using the proposed thresholds for the identification of specimens. First, the identification capacity followed the design of Yodphaka et al. (2018), where it was assessed whether there was a sufficient gap between intraspecific and interspecific distances, called “barcoding gap” and whether the closest neighbor was co-specific (the test of the closest neighbor; Meier et al. 2006). The genetic distance between individual sequences was calculated using the p-distance model in the “dist.dna” function in the “SPIDER” package (Brown et al. 2012) of the R software (R Development Core team 2020). Traditionally, many studies use the K2P distance model in barcode analyzes, however, this has been challenged and the p-distance has been proposed to be a better model (Collins et al. 2012, Srivathsan & Meier 2012).

The quantification of the barcoding gap was calculated from the difference between the shortest interspecific distance and the longest intraspecific distance. These distances were determined with the “nonConDist” and “maxInDist” functions in the “SPIDER” package (Brown et al. 2012). The graphs are plotted as violin plots using PlotsOfData (https://huygens.science.uva.nl/PlotsOfData/). The Kruskal-Wallis and Mann-Whitney tests were used to determine differences in the mean number of genetic distances by gene (P<0.01) using PAST v.3.18 software (Hammer et al. 2001). The percentage of correct identification was calculated from the number of sequences with a specific closest neighbor, divided by the total number of sequences. The test was performed with the “nearNeighbor” function on the “SPIDER” package (Brown et al. 2012). We also evaluated which species had barcoding gap. A marker with high discriminatory power must have a high percentage of correct identifications from the nearest neighbor test and a positive value for the barcode gap.

For threshold analysis, we evaluated the quality of the data set by simulating a specimen identification scenario based on sequence-sequence using R software with the APE and SPIDER packages (Brown et al. 2012, Popescu et al. 2012). The identification was provided following two different criteria: Best closed correspondence (BCM) and barcode of all species (BOLD ID). The BCM criterion assigns identifications to the closest match with a distance below a defined limit. The BOLD ID criterion through the “threshID” function simulates the BOLD ID mechanism (http://www.boldsystems.org/index.php/IDS_OpenIdEngine) applying the predetermined limit and consulting all the sequences of a single species below it. The results are reported as correct when they correspond to previous morphological identifications, otherwise, a result is considered incorrect. A query can provide ambiguous results if the sequences of different species are those that correspond below the limit (BCM) and if the divergences of sequences of different species are below the limit (BOLD ID). A query that results in “no ID” does not match below the defined threshold.

We used six different threshold values: the value obtained by using the “ThreshVal” function in the “SPIDER” package (from now on “ThreshVal”), which minimizes cumulative identification errors, that is, the sum of the false positive (without specific matches within the query limit) and the negative-negative (sequences of several species within the limit). The value obtained by the ‘localMinima’ function (SPIDER package), which indicates the minimum in density, corresponds to the transition between intra and interspecific distances. The value of 1%, which is the default used by the BOLD ID mechanism. The value of 3%, which is used in studies of species delimitation with the 12s and 16s genes in anurans (e.g. Howlader et al. 2016), and the value of 6%, which is the standard suggested by Lyra et al. (2017) for the delimitation of species for COI and close to the value suggested by Kartavtsev & Lee (2006) for Cytb. Finally, we use the value of 10%, as suggested by Vences et al. (2005a) and Kartavtsev & Lee (2016) for COI and to congeneric species for Cytb, respectively.

We perform the counting of monophyletic groups in each database using the “Monophyly” function in the “SPIDER” package. This metric was used only and simply to determine the ability of a marker to recover monophyly among sequences of the same species. As presented by Knowles & Carstens (2007), the use of a single gene does not fully qualify for species identification. In this sense, the high number of species with reciprocal monophyly in a gene may be indicative that the time needed to coalesce has elapsed, these being preferably chosen in multiloci studies.

RESULTS

Datasets

The “All sequences” datasets for 12S, 16S, COI and Cytb consisted of a total of 1,690, 2,712, 1,737 and 2,425 individual sequences from 259, 276, 156 and 53 species, respectively (Fig. 1a; nucleotide alignments provide in Supplementary Appendix S1-S4). The “Species in common” datasets for 12S, 16S, COI and Cytb included 444, 499, 504 and 1,449 individual sequences from 31 species, respectively. Finally, the “species in common without Cytb gene” datasets for 12S, 16S and COI included 1,036, 1,341 and 1,113 individual sequences from 100 species, respectively. Detailed information about the sequence curation and alignment size are presented in Table I.

Figure 1
Comparison among the studied barcoding markers for Hylidae from “All sequences”, “Species in common” and “Species in common without Cytb” datasets. (a) Number of sequences, species and genera per marker. (b) Distribution of barcoding gap, as defined by the difference between the minimum non-conspecific distance and the maximum conspecific distance. Different letters represents a statistically significant difference in means based on the Mann-Whitney test (P<0.01). (c) The percentage of correct identifications from the nearest neighbor test.

Distance-based evaluations

In “All sequences” datasets, on average (mean ± SD), the distribution of the DNA barcoding gap showed that the interspecific distances were greater than the intraspecific distances in 12S (0.00±0.03), COI (0.03±0.05) and Cytb (0.04±0.06), while the intraspecific distances were greater than the interspecific distances in 16S (negative barcoding gap; -0.01±0.06, Fig. 1b). Considering each gene dataset, the percentage of species with barcoding gap was 72.43% for COI, 68.72% for 12S, 63.40% for 16S and 62.26% for Cytb (see full list of species that presented barcoding gap in Supplementary Material - Table SI).

When comparing the same species (“Species in common” datasets), only COI (0.01±0.04) and Cytb (0.05±0.05) had a positive gap in DNA barcoding (Fig. 1b). When comparing only 12S, 16S and COI (“Species in common without Cytb” datasets), all genes had a positive gap with a higher average for the COI gene (0.04±0.06) (Fig. 1b). The barcoding gaps varied significantly between the studied markers (Kruskal-Wallis test, P<0.01). In “All sequences” and “Species in common” datasets, the nearest neighbor test showed the highest percentage of more specific neighbors in Cytb (98.56% and 98.89%) and the lowest percentage in 12S (94.32%) and 16S (95.10%), respectively (Fig. 1c). In “Species in common without CytB” datasets, the highest percentages were observed in COI gene, 98.74%.

Threshold values

In “All sequences” datasets (Table II), for the BCM approach, the correct identifications ranged from 73.49% in 12S (“ThreshVal” value of 0.3%) to 97.47% in COI (threshold value of 10%). Sequences “no ID” ranged from 0.00% in Cytb (threshold value of 10%) to 20.36% in 12S (“ThreshVal” value of 0.3%). The incorrect and ambiguous identifications ranged from 0.11% to 4.87%. For the BOLD ID approach, the correct identifications ranged from 4.73% in 12S (threshold value of 10%) to 85.26% in COI (“ThreshVal” value of 1.15%), incorrect and ambiguous identifications ranged from 0.12% to 94.91%, and “no ID” identification had values from 0.00% in Cytb (threshold value of 10%) to 20.36% in 12S (“ThreshVal” value of 0.3%).

Table II
Results of the identification simulations from “All sequences” datasets using Best Close Match (BCM) and BOLD ID criteria based on SPIDER and tree-based comparison of efficiency among the studied barcoding markers using the percentage of monophyletic groups recovered from the neighbor-joining phylogenetic reconstructions.

			Distance-based evaluation																								Tree-based evaluation
			threshVal				LocalMinima				1.00%				3.00%				6.00%					10.00			Monophyly
All sequences
	All sequences			Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)
		12s (259 sp; 1690 spm)	Threshold value	0.30%				0.66%				1.00%				3.00%				6.00%				10.00%				71.43	28.57
			BCM	73.49	2.13	4.02	20.36	81.24	2.78	4.32	11.66	84.50	3.02	4.38	8.11	89.41	3.85	4.44	2.31	90.59	4.32	4.44	0.65	90.77	4.73	4.44	0.06
			BOLD ID	65.56	1.54	12.54	20.36	64.32	1.24	22.78	11.66	62.31	1.18	28.40	8.11	38.17	0.95	58.58	2.31	15.44	0.65	83.25	0.65	4.73	0.30	94.91	0.06
		16s (276 sp; 2712 spm)	ThreshID value	0.50%				0.37%				1.00%				3.00%				6.00%				10.00%				68.12	31.88
			BCM	82.45	1.92	4.54	11.10	80.64	1.70	4.42	13.24	87.72	2.36	4.61	5.31	90.63	2.69	4.65	2.03	91.78	2.91	4.65	0.66	92.00	3.20	4.65	0.15
			BOLD ID	76.51	1.11	11.28	11.10	76.00	0.92	9.85	13.24	68.73	1.18	24.78	5.31	49.37	0.44	48.16	2.03	24.19	0.37	74.78	0.66	5.53	0.29	94.03	0.15
		COI (156 sp; 1737 spm)	ThreshID value	1.15%				0.91%				1.00%				3.00%				6.00%				10.00%				75.00	25.00
			BCM	90.39	0.86	0.12	8.64	88.72	0.86	0.12	10.31	89.23	0.86	0.12	9.79	94.13	1.09	0.12	4.66	96.03	1.21	0.12	2.65	97.47	1.21	0.11	1.21
			BOLD ID	85.26	0.46	5.64	8.64	84.80	0.46	4.43	10.31	85.20	0.46	4.55	9.79	81.81	0.46	13.07	4.66	77.14	0.46	19.75	2.65	69.77	0.29	28.73	1.21
		CytB (53 sp; 2425 spm)	ThreshID value	0.70%				1.25%				1.00%				3.00%				6.00%				10.00%				67.92	32.08
			BCM	96.08	0.83	3.01	0.08	94.72	0.66	3.01	1.61	94.56	0.58	3.01	1.85	95.71	0.78	3.01	0.50	96.04	0.78	3.01	0.17	96.12	4.87	3.01	0.00
BOLD ID			62.43	0.17	37.32	0.08	82.84	0.29	15.26	1.61	83.05	0.25	14.85	1.85	81.77	0.25	17.48	0.50	68.48	0.16	31.19	0.17	54.64	0.12	45.24	0.00

In “Species in common” datasets (Table III), for the BCM approach, the correct identifications ranged from 74.77% in 12S (“ThreshVal” value of 0.15%) to 96.62% in Cytb (threshold value of 10%). Sequences “no ID” ranged from 0.00% (threshold values ​​of 3% and 6%) to 18.47% in 12S (“ThreshVal” value of 0.15%). The incorrect and ambiguous identifications ranged from 0.00% to 10.47%. For BOLD ID approach, the correct identifications ranged from 2.70% in 12S (threshold value of 10%) to 95.38% in Cytb (“ThreshVal” value of 0.7%), incorrect and ambiguous identifications ranged from 0.00% to 97.39% and “no ID” identification had values ​​of 0.00% (threshold value of 3%, 6% and 10%) to 18.47% in 12S (“ThreshVal” value of 0.15%).

Table III
Results of the identification simulations from “Species in common” datasets using Best Close Match (BCM) and BOLD ID criteria based on SPIDER and tree-based comparison of efficiency among the studied barcoding markers using the percentage of monophyletic groups recovered from the neighbor-joining phylogenetic reconstructions.

			Distance-based evaluation																								Tree-based evaluation
			threshVal				LocalMinima				1.00%				3.00%				6.00%					10.00			Monophyly
Species in common (31 sp)			Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)
	12s (444 spm)	Threshold value	0.15%				0.86%				1.00%				3.00%				6.00%				10.00%				58.06	41.94
		BCM	74.77	1.58	5.18	18.47	88.06	2.48	6.08	3.38	88.74	2.48	6.08	2.70	90.77	2.93	6.31	0.00	90.77	2.93	6.31	0.00	90.77	2.93	6.30	0.00
		BOLD ID	67.79	12.39	1.35	18.47	39.64	0.68	56.31	3.38	38.06	0.68	58.56	2.70	20.95	0.23	78.83	0.00	5.86	0.00	94.14	0.00	2.70	0.00	97.30	0.00
	16s (499 spm)	ThreshID value	0.45%				1.02%				1.00%				3.00%				6.00%				10.00%				51.61	48.39
		BCM	77.51	2.90	10.02	9.58	82.18	3.34	10.24	4.23	82.18	3.34	10.24	4.23	85.52	3.79	10.47	0.22	85.52	3.79	10.47	0.22	85.52	3.79	10.46	0.23
		BOLD ID	63.92	1.34	25.17	9.58	58.80	1.34	35.63	4.23	58.80	1.34	35.63	4.23	32.29	0.22	67.26	0.22	15.59	0.22	83.96	0.22	7.57	0.00	92.20	0.23
	COI (504 spm)	ThreshID value	0.65%				7.72%				1.00%				3.00%				6.00%				10.00%				70.97	29.03
		BCM	92.86	0.20	0.00	6.94	98.41	0.99	0.00	0.60	93.45	0.40	0.00	6.15	96.23	0.99	0.00	2.78	97.42	0.99	0.00	1.59	98.61	0.99	0.00	0.40
		BOLD ID	89.88	0.20	2.98	6.94	77.38	0.60	21.43	0.60	86.31	0.40	7.14	6.15	77.38	0.60	19.25	2.78	77.38	0.60	20.44	1.59	75.99	0.40	23.21	0.40
	Cytb (1449 spm)	ThreshID value	0.70%				1.14%				1.00%				3.00%				6.00%				10.00%				77.42	22.58
		BCM	96.34	0.35	0.69	2.62	96.89	0.35	0.69	2.07	96.75	0.35	0.69	2.21	98.14	0.48	0.69	0.69	98.55	0.55	0.69	0.21	98.62	0.69	0.69	0.00
		BOLD ID	95.38	0.07	1.93	2.62	93.79	0.00	4.14	2.07	94.00	0.00	3.79	2.21	93.30	0.14	5.87	0.69	89.85	0.14	9.80	0.21	60.25	0.07	39.68	0.00

Finally, in “Species in common without Cytb” datasets (Table IV), for the BCM approach, the correct identifications ranged from 71.81% in 12s (“ThreshVal” value of 0.15%) to 96.86% in COI at the threshold value of 6%. Sequences “no ID” ranged from 0.30% (threshold value of 6%) to 10.95% in 16S (“ThreshVal” value of 0.50%). The incorrect and ambiguous identifications ranged from 0.09% to 4.84%. For the BOLD ID approach, the correct identifications ranged from 21.43% in 12S (threshold value of 6%) to 85.62% in COI (“ThreshVal” value of 1.00%), incorrect and ambiguous identifications ranged from 0.27% to 77.51% and “no ID” identification had values ​​of 0.30% in 16S (threshold value of 6%) to 24.13% in 12S (“ThreshVal” value of 0.15%).

Table IV
Results of the identification simulations from “Species in common without Cytb” datasets using Best Close Match (BCM) and BOLD ID criteria based on SPIDER and tree-based comparison of efficiency among the studied barcoding markers using the percentage of monophyletic groups recovered from the neighbor-joining phylogenetic reconstructions.

			Distance-based evaluation																				Tree-based evaluation
			threshVal				LocalMinima				1.00%				3.00%				6.00%				Monophyly
Species in common without Cytb (100 sp)			Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)	Ambiguous (%)	No ID (%)	Correct (%)	Incorrect (%)
	12s (1036 spm)	Threshold value	0.15%				0.78%				1.00%				3.00%				6.00%				65.00	35.00
		BCM	71.81	1.06	2.99	24.13	86.20	1.83	3.57	8.40	88.32	2.03	3.57	6.08	92.37	2.51	3.67	1.45	92.76	2.80	3.67	0.77
		BOLD	68.82	0.87	6.18	24.13	61.97	0.77	28.86	8.40	60.33	0.97	32.63	6.08	46.81	0.58	51.16	1.45	21.43	0.29	77.51	0.77
	16s (1342 spm)	ThreshID value	0.50%				2.60%				1.00%				3.00%				6.00%				67.00	33.00
		BCM	83.01	1.42	4.62	10.95	91.43	2.24	4.84	1.49	88.30	2.01	4.77	4.92	91.73	2.31	4.84	1.12	92.47	2.38	4.84	0.30
		BOLD	75.04	0.67	13.34	10.95	61.77	0.37	36.36	1.49	71.09	0.82	23.17	4.92	58.57	0.22	40.09	1.12	34.43	64.98	0.30	0.30
	COI (1113 spm)	ThreshID value	1.00%				0.89%				1.00%				3.00%				6.00%				80.00	20.00
		BCM	89.04	0.36	0.09	10.51	87.96	0.36	0.09	11.59	89.04	0.36	0.09	10.51	94.61	0.63	0.09	4.67	96.86	0.63	0.09	2.43
		BOLD	85.62	0.27	3.59	10.51	84.64	0.27	3.50	11.59	85.62	0.27	3.59	10.51	84.91	0.36	10.06	4.67	81.58	0.36	15.63	2.43

Monophyly analysis

COI recovered the highest percentage of monophyletic groups in the reconstructions of both “All sequences” (75.00%) and “Species in common without Cytb” (80.00%) datasets. In “Species in common” dataset, Cytb produced the largest number of monophyletic groups (77.42%). All markers recovered a percentage of monophyletic groups above 50% (Tables II-IV).

DISCUSSION

In this work, the efficiencies of four barcode markers species identification within hylids were evaluated. The coding genes, COI and Cytb, showed greater discriminatory power than the ribosomal genes. However, the advantages and disadvantages of each marker must be considered.

Despite the curation of the sequences, an intrinsic limitation of this study is the use of the public database. Genbank is the leading public sequence database with more than 216 million sequences deposited (April 2020). It uses basic checks such as vector contamination analysis, adequate translation of the contamination regions, and genetic annotation. However, unlike the main DNA barcoding database, the BOLD system, GenBank does not store chromatograms or collection metadata, which increases the risk of errors regarding each sequence. Note that GenBank and BOLD system exchange sequences.

Regardless sequences deposition, the generation and submission of incorrect sequences must occur mainly due to incorrect identification of the original material. Murphy et al. (2013) indicated that about 5.00% of the specimens deposited in reference collections of anurans and reptiles have some identification error. Other factors should also be considered for analyzes with pre-defined thresholds such as variability or heterogeneity in the rates of evolution and the presence of cryptic species. (e.g. Hoskin et al. 2005, Murphy et al. 2013, Caminer & Ron 2014). In addition, we cannot rule out the potential for heterotachy, heteroplasmy, pseudogenes, and other changes in evolution rates. We emphasize that, unlike other groups (see Summerbell et al. 2007), there are few sequences (less than 0.10%, April 2020) deposited with an earlier date of 2000 for Hylidae and Anura, which rules out the hypothesis of errors caused by artifacts of sequencing technologies of the 1980–1990s.

In the proposal for standardized barcode regions for animals, COI was indicated as the main marker (Hebert et al. 2003). Among the main attributes for the use of the COI gene were universality and high rate of substitution in the third codon position (Frati et al. 2000). The great development of public genetic databases for this gene since the publication of Hebert et al. (2003) has provided an important framework for its use in specimen identification. Recent evaluations have advocated its use compared to other genes (e.g. Lyra et al. 2017). The development of new primers and traceable databases, such as the BOLD platform, has been supporting its effectiveness and choice. However, the COI gene is still highly questioned in its application capacity, especially in high-throughput applications, which has been defended in recent years due to the choice of new primers and methodologies for this type of sequencing (see Andújar et al. 2018, Pierce 2019).

Among the most used genes for identification in high-throughput applications is 16s. Vences et al. (2005b) in their seminal article indicated the use of 16S as a true DNA barcode for Amphibia (and even for vertebrates) evidenced its use for pooled samples. At this time until today, as we found out in our database analysis for Hylidae, this gene has the largest number of sequences and species deposited. It is also the main gene used in anuran studies of species description, phylogeny and phylogeography (e.g. Neves et al. 2017, Mângia et al. 2018, 2020, Andrade et al. 2019). Another advantage of this gene is the numerous metabarcoding protocols. Several examples already allow and validate its use in the identification of anurans (e.g. Bálint et al. 2018). However, the use of 16S for species identification has important limitations for interpretation due to its hypervariable domains. The choice of the alignment method should result in different levels of genetic divergence estimates, which can result in a strong effect in analyzes downstream of phylogenetic inference (see Lyra et al. 2017).

The other two genes evaluated in this work have received attention for their use in species identification. 12S has been used in the identification by metabarcoding (e.g. Lopes et al. 2017), besides being historically used in phylogenetic and phylogeographic analyzes (e.g. Garda & Cannatella 2007, Andrade et al. 2019). Similar to 16S, 12S has a problem with hypervariable regions. In the case of Cytb, the use of DNA barcoding for vertebrates has been advocated for nearly two decades (e.g. Bradley & Baker 2001). Cytb is considered more variable than the COI (Vences et al. 2005a) and has been used even in population analyzes (Recuero et al. 2006). However, for other taxonomic groups, the large number of highly conserved gene regions may limit the use of this gene in species identification (e.g. Satoh et al. 2016). Furthermore, for Cytb in hylids, the number of species present in the public databases is still far from those found for the other evaluated genes.

The high efficiency in the identification of COI and Cytb samples was corroborated by extent of the so-called barcoding gap and our simulations using BCM and BOLD ID with different tolerance levels. As evidenced by Chapple & Ritchie (2013), the presence of barcoding gap is essential to species delimitation and is the basis for specimen identification (local gap) and species discovery (global gap). For all datasets, the barcoding gap value in coding genes values were significantly higher than those of non-coding genes. Regarding the threshold values, this work did not aim to define the best threshold value of molecular divergence for hylids. The accuracy of species delimitation using this methodology is influenced by factors such as the quality of the reference database, the geographic extent of the sample, the representativeness of the intraspecific variation, nucleotide replacement rate of the species, among others (see Chapple & Ritchie 2013), which were not evaluated in this work. The threshold values evaluated here are the most commonly used in taxonomic studies and allowed us to simulate their results to sequence-based identification process using our broad genetic databases. The use of these threshold values also allows us to highlight low interspecific distances (“ambiguous” entries in Tables II-III) and to exclude correspondences with high distances (“no ID” entries in Tables II-III). The Bold ID is considered the most rigorous, which does not provide correct identification if all sequences of a given species are not below the proposed limit (see Raupach et al. 2015). In these evaluations, the ability to identify the coding genes was even more prominent compared to the “near neighbor”. Finally, the monophyly analysis showed the COI gene as the best indication of the taxa, with the other genes showing similar capacities.

Therefore, our results indicate that the coding genes, especially COI, have a better identification capacity than the non-coding genes, 16S and 12S, for hylids. It is noted that there is still a limitation in the number of species of hylids with coding gene sequences in relation to non-coding genes. This is a fundamental topic to defend the continued use of non-coding genes for species identification. In addition, as evidenced in DNA barcoding analysis of global databases, the use of local/regional databases should always be prioritized considering their best identification precision (see Bergsten et al. 2012), which can change the efficiency in the use of each gene. Besides, these results provide a guide to the use of a single gene in studies to assess biodiversity. Considering that amphibians are facing a series of severe threats and the limitation of taxonomists and financier resources for conservation, we provide in this study an indicative for a better allocation of resources in genetic studies (e.g. construction of DNA barcoding reference databases, specimen identification work) for this representative group.

ACKNOWLEDGMENTS

RK received a Post-doctoral scholarship (PNPD/CAPES) from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES). DJS received research fellowship (proc. 311492/2017-7) from the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code 001. We also thank the Editor-in-Chief, Dr. Igor Luis Kaefer, and an anonymous reviewer for helpful comments.

SUPPLEMENTARY MATERIAL

https://minio.scielo.br/documentstore/1678-2690/fbmw4RfDv34jvL6WP3SVRNj/fb9c30ccdf3b8c516b00fa23bfe31656e0861203.pdf https://minio.scielo.br/documentstore/1678-2690/fbmw4RfDv34jvL6WP3SVRNj/bfeea1d55861f06c0d75ffd512d790446022bd42.pdf https://minio.scielo.br/documentstore/1678-2690/fbmw4RfDv34jvL6WP3SVRNj/33e99e9fa2e623dcdbe7f65e0941dc3f307eb6b5.pdf https://minio.scielo.br/documentstore/1678-2690/fbmw4RfDv34jvL6WP3SVRNj/c5d241a9b8d1c50f0b7062cb57f1e79cd0b88874.pdf https://minio.scielo.br/documentstore/1678-2690/fbmw4RfDv34jvL6WP3SVRNj/8e6cb5deed098e2abfe06b86c47f49e1b23afbf7.pdf

Table SI, Appendix S1, S2, S3, S4

REFERENCES

• ANDRADE FS DE, SILVA LA DA, KOROIVA R, FADEL RM & SANTANA DJ. 2019. A New Species of Pseudopaludicola Miranda-Ribeiro, 1926 (Anura: Leptodactylidae: Leiuperinae) from an Amazonia-Cerrado Transitional Zone, State of Tocantins, Brazil. J Herpetol 53: 68-80.
• ANDÚJAR C, ARRIBAS P, YU DW, VOGLER AP & EMERSON BC. 2018. Why the COI barcode should be the community DNA metabarcode for the metazoa. Mol Ecol 27: 3968-3975.
• BÁLINT M, NOWAK C, MÁRTON O, PAULS SU, WITTWER C, ARAMAYO JL, SCHULZE A, CHAMBERT T, COCCHIARARO B & JANSEN M. 2018. Accuracy, limitations and cost efficiency of eDNA-based community survey in tropical frogs. Mol Ecol Resour 18: 1415-1426.
• BERGSTEN J ET AL. 2012. The Effect of Geographical Scale of Sampling on DNA Barcoding. Syst Biol 61: 851-869.
• BRADLEY RD & BAKER RJ. 2001. A test of the genetic species concept: Cythrochrome-b sequences and Mammals. J Mammal 82: 960-973.
• BROWN SDJ, COLLINS RA, BOYER S, LEFORT M-C, MALUMBRES-OLARTE J, VINK CJ & CRUICKSHANK RH. 2012. Spider: an R package for the analysis of species identity and evolution, with particular reference to DNA barcoding. Mol Ecol Resour 12: 562-565.
• CAMINER M & RON S. 2014. Systematics of treefrogs of the Hypsiboas calcaratus and Hypsiboas fasciatus species complex (Anura, Hylidae) with the description of four new species. Zookeys 370: 1-68.
• CHAPPLE DG & RITCHIE PA. 2013. A Retrospective Approach to Testing the DNA Barcoding Method. PLoS ONE 8: e77882.
• CHE J, CHEN HM, YANG JX, JIN JQ, JIANG K, YUAN ZY, MURPHY RW & ZHANG YP. 2012. Universal COI primers for DNA barcoding amphibians. Mol Ecol Resour 12: 247-258.
• COLLINS RA, BOYKIN LM, CRUICKSHANK RH & ARMSTRONG KF. 2012. Barcoding’s next top model: An evaluation of nucleotide substitution models for specimen identification. Methods Ecol Evol 3: 457-465.
• CRAWFORD AJ, LIPS KR & BERMINGHAM E. 2010. Epidemic disease decimates amphibian abundance, species diversity, and evolutionary history in the highlands of central Panama. Proc Natl Acad Sci 107: 13777-13782.
• FOUQUET A, GILLES A, VENCES M, MARTY C, BLANC M & GEMMELL NJ. 2007. Underestimation of Species Richness in Neotropical Frogs Revealed by mtDNA Analyses. PLoS ONE 2: e1109.
• FRATI F, DELL’AMPIO E, CASASANTA S, CARAPELLI A & PAOLO FANCIULLI P. 2000. Large Amounts of Genetic Divergence among Italian Species of the Genus Orchesella (Insecta, Collembola) and the Relationships of Two New Species. Mol Phylogenet Evol 17: 456-461.
• FROST DR. 2020. Amphibian Species of the World: an online reference. Version 6.1. Available from: https://amphibiansoftheworld.amnh.org/index.php
  » https://amphibiansoftheworld.amnh.org/index.php
• GAO XY, DONG BJ, LI JT, WANG G, JIANG JP, YANG BT & WANG B. 2019. Phylogeographic investigation on the spiny frog Quasipaa shini (Amphibia: Anura: Dicroglossidae) using mitochondrial DNA: cryptic species and species complex. Mitochondrial DNA Part B 4(1): 1479-1483.
• GARDA AA & CANNATELLA DC. 2007. Phylogeny and biogeography of paradoxical frogs (Anura, Hylidae, Pseudae) inferred from 12S and 16S mitochondrial DNA. Mol Phylogenet Evol 44: 104-114.
• GOEBEL AM, DONNELLY JM & ATZ ME. 1999. PCR Primers and Amplification Methods for 12S Ribosomal DNA, the Control Region, Cytochrome Oxidase I, and Cytochromebin Bufonids and Other Frogs, and an Overview of PCR Primers which Have Amplified DNA in Amphibians Successfully. Mol Phylogenet Evol 11: 163-199.
• HAMMER Ø, HARPER DAT & RYAN PD. 2001. PAST: Paleontological Statistics Software Package for Education and Data Analysis. Palaeontol Electron 4: 1-9.
• HEBERT PDN, CYWINSKA A, BALL SL & DEWAARD JR. 2003. Biological identifications through DNA barcodes. Proc R Soc B Biol Sci 270: 313-321.
• HOSKIN CJ, HIGGIE M, MCDONALD KR & MORITZ C. 2005. Reinforcement drives rapid allopatric speciation. Nature 437: 1353-1356.
• HOWLADER MSA, NAIR A & MERILÄ J. 2016. A New Species of Frog (Anura: Dicroglossidae) Discovered from the Mega City of Dhaka. PLoS ONE 11: e0149597.
• JEONG TJ, JUN J, HAN S, KIM HT, OH K & KWAK M. 2013. DNA barcode reference data for the Korean herpetofauna and their applications. Mol Ecol Resour 13: 1019-1032.
• KARTAVTSEV YP & LEE J-S. 2006. Analysis of nucleotide diversity at the cytochrome b and cytochrome oxidase 1 genes at the population, species, and genus levels. Russ J Genet 42: 341-362.
• KATOH K & STANDLEY DM. 2013. MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability. Mol Biol Evol 30(4): 772-780.
• KEARSE M ET AL. 2012. Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28: 1647-1649.
• KNOWLES LL & CARSTENS BC. 2007. Delimiting Species without Monophyletic Gene Trees. Syst Biol 56: 887-895.
• KOROIVA R, DE SOUZA MS, ROQUE F DE O & PEPINELLI M. 2018. DNA Barcodes for Forensically Important Fly Species in Brazil. J Med Entomol 55: 1055-1061.
• KOROIVA R & KVIST S. 2018. Estimating the barcoding gap in a global dataset of cox1 sequences for Odonata: close, but no cigar. Mitochondrial DNA Part A DNA Mapping, Seq Anal Jul 28: 765-771.
• KOROIVA R, PEPINELLI M, RODRIGUES ME, ROQUE FDO, LORENZ-LEMKE AP & KVIST S. 2017. DNA barcoding of odonates from the Upper Plata basin: Database creation and genetic diversity estimation. PLoS ONE 12: e0182283.
• KOROIVA R, RODRIGUES LRR & SANTANA DJ. 2020. DNA barcoding for identification of anuran species in the central region of South America. PeerJ 8: e10189.
• LEHR E, MAY R VON, MORAVEC J & CUSI JC. 2017. Three new species of Pristimantis (Amphibia, Anura, Craugastoridae) from upper montane forests and high Andean grasslands of the Pui Pui Protected Forest in central Peru. Zootaxa 4299: 301-336.
• LOPES CM, SASSO T, VALENTINI A, DEJEAN T, MARTINS M, ZAMUDIO KR & HADDAD CFB. 2017. EDNA metabarcoding: a promising method for anuran surveys in highly diverse tropical forests. Mol Ecol Resour 17: 904-914.
• LYRA ML, HADDAD CFB & DE AZEREDO-ESPIN AML. 2017. Meeting the challenge of DNA barcoding Neotropical amphibians: polymerase chain reaction optimization and new COI primers. Mol Ecol Resour 17: 966-980.
• MÂNGIA S, KOROIVA R, NUNES PMS, ROBERTO IJ, ÁVILA RW, SANT’ANNA AC, SANTANA DJ & GARDA AA. 2018. A New Species of Proceratophrys (Amphibia: Anura: Odontophrynidae) from the Araripe Plateau, Ceará State, Northeastern Brazil. Herpetologica 74: 255-268.
• MÂNGIA S, OLIVEIRA EF, SANTANA DJ, KOROIVA R, PAIVA F & GARDA AA. 2020. Revising the taxonomy of Proceratophrys Miranda-Ribeiro, 1920 (Anura: Odontophrynidae) from the Brazilian semiarid Caatinga: Morphology, calls and molecules support a single widespread species. J Zool Syst Evol Res 58(4): 1151-1172.
• MEIER R, SHIYANG K, VAIDYA G & NG PKL 2006. DNA Barcoding and Taxonomy in Diptera: A Tale of High Intraspecific Variability and Low Identification Success. Syst Biol 55: 715-728
• MEYER CP & PAULAY G. 2005. DNA barcoding: Error rates based on comprehensive sampling. PLoS Biol 3: 1-10.
• MURPHY RW ET AL. 2013. Cold Code: The global initiative to DNA barcode amphibians and nonavian reptiles. Mol Ecol Resour 13: 161-167.
• NEVES MO, DA SILVA LA, AKIEDA PS, CABRERA R, KOROIVA R & SANTANA DJ. 2017. A new species of poison frog, genus Ameerega (Anura: Dendrobatidae), from the southern amazonian rain forest. Salamandra 53: 485-493.
• ORRICO VGD ET AL. 2017. Integrative taxonomy supports the existence of two distinct species within Hypsiboas crepitans (Anura: Hylidae). Salamandra 53: 99-113.
• PIERCE MP. 2019. Filling in the Gaps: Adopting Ultraconserved Elements Alongside COI to Strengthen Metabarcoding Studies. Front Ecol Evol 7: 469.
• POPESCU AA, HUBER KT & PARADIS E. 2012. Ape 3.0: New tools for distance-based phylogenetics and evolutionary analysis in R. Bioinformatics 28: 1536-1537.
• R DEVELOPMENT CORE TEAM. 2020. R A Language and Environment for Statistical Computing. Vienna, Austria: 2020.
• RAUPACH MJ ET AL. 2015. The Application of DNA Barcodes for the Identification of Marine Crustaceans from the North Sea and Adjacent Regions. PLoS ONE 10: e0139421.
• RAUPACH MJ, HANNIG K, MORINIÈRE J & HENDRICH L. 2018. A DNA barcode library for ground beetles of Germany: the genus Amara Bonelli, 1810 (Insecta, Coleoptera, Carabidae). Zookeys 759: 57-80.
• RECUERO E, MARTÍNEZ-SOLANO Í, PARRA-OLEA G & GARCÍA-PARÍS M. 2006. Phylogeography of Pseudacris regilla (Anura: Hylidae) in western North America, with a proposal for a new taxonomic rearrangement. Mol Phylogenet Evol 39: 293-304.
• SATOH TP, MIYA M, MABUCHI K & NISHIDA M. 2016. Structure and variation of the mitochondrial genome of fishes. BMC Genomics 17: 719.
• SRIVATHSAN A & MEIER R. 2012. On the inappropriate use of Kimura-2-parameter (K2P) divergences in the DNA-barcoding literature. Cladistics 28: 190-194.
• STUART SN, CHANSON JS, COX NA, YOUNG BE, RODRIGUES ASL, FISCHMAN DL & WALLER RW. 2004. Status and Trends of Amphibian Declines and Extinctions Worldwide. Science 306: 1783-1786.
• SUMMERBELL RC, MOORE MK, STARINK-WILLEMSE M & VAN IPEREN A. 2007. ITS barcodes for Trichophyton tonsurans and T. equinum. Med Mycol 45: 193-200.
• TANAKA-UENO T, MATSUI M, SATO T, TAKENAKA S & TAKENAKA O. 1998. Phylogenetic Relationships of Brown Frogs with 24 Chromosomes from Far East Russia and Hokkaido Assessed by Mitochondrial Cytochrome b Gene Sequences (Rana: Ranidae). Zoolog Sci 15: 289-294.
• VENCES M, NAGY ZT, SONET G & VERHEYEN E. 2012. DNA barcoding amphibians and reptiles. Methods Mol Biol 858: 79-107.
• VENCES M, THOMAS M, BONETT RM & VIEITES DR. 2005a. Deciphering amphibian diversity through DNA barcoding: chances and challenges. Philos Trans R Soc Lond B Biol Sci 360: 1859-1868.
• VENCES M, THOMAS M, VAN DER MEIJDEN A, CHIARI Y & VIEITES DR. 2005b. Comparative performance of the 16S rRNA gene in DNA barcoding of amphibians. Front Zool 2: 5.
• VIEITES DR, WOLLENBERG KC, ANDREONE F, KOHLER J, GLAW F & VENCES M. 2009. Vast underestimation of Madagascar’s biodiversity evidenced by an integrative amphibian inventory. Proc Natl Acad Sci 106: 8267-8272.
• WAUGH J. 2007. DNA barcoding in animal species: progress, potential and pitfalls. BioEssays 29: 188-197.
• YODPHAKA S, BOONPRAGOB K, LUMBSCH HT & KRAICHAK E. 2018. Evaluation of six regions for their potential as DNA barcodes in epiphyllous liverworts from Thailand. Appl Plant Sci 6: e01174.
• ZHAO H, YANG J, WANG C, LI P, MURPHY RW, CHE J & YUAN Z. 2017. A new species of the genus Rana from Henan, central China (Anura, Ranidae). Zookeys 694: 95-108.

Publication Dates

• Publication in this collection
  05 Dec 2022
• Date of issue
  2022

History

• Received
  5 June 2020
• Accepted
  19 May 2021