We aimed to extract mtDNA from type specimens for which molecular data was unavailable. The types of Ptychadena cooperi and P. harenna were not sequenced, as these two species are easily distinguishable morphologically and there is no ambiguity regarding their taxonomic status (Goutte et al. 2021). We obtained authorization from the Museum für Naturkunde Berlin (ZMB) and the Museum d’Histoire Naturelle, Genève (MHNG) to sample a small amount of muscle or liver tissue from the type specimens of Ptychadena neumanni (lectotype, ZMB 26879; paralectotype, ZMB 57183), P. erlangeri (holotype, ZMB 26887), P. largeni (paratypes, MHNG 2513-15 & 2513-56) and P. nana (holotype, ZMB 26878). Tissue sampling did not result in any major visible damage to the vouchers, as most tissue was taken from pre-existing incisions. Dissecting tools were cleaned with a 10% bleach solution before and after tissue extraction.

Although we could not find information about the mode of preservation used at the time of collection, it is likely that the type specimens had been fixed in formalin or some form of spirit, which renders the extraction of DNA challenging and requires a different protocol than when using fresh tissue. We followed the protocol described by Shedlock et al. (1997). In short, approximately ~5 mg of tissue (muscle or liver) was placed in a 2 mL tube and washed in 1X GTE buffer for four days, changing the buffer daily. The tissue was then incubated for four additional days at 65 °C in a 2 mL tube with 500 μL cell lysis buffer, 100 μL proteinase K, and 20 μL 1 mM DTT (dithiothreitol). Proteinase K was added daily until all the tissue was completely digested. A standard potassium acetate (KOA) DNA precipitation protocol was then followed. A detailed protocol for the DNA extraction is provided in the Suppl. material 1. To reduce the possibility of contamination, DNA extraction was carried out using new reagents in a marine biology lab that does not work with amphibian samples, and multiple negative controls (using deionized water) were used along the process.

DNA concentration was measured using a high sensitivity kit in a Qubit fluorometer (Life Technologies) and DNA fragment size distribution and concentration was estimated using a Bioanalyzer 7500 high sensitivity DNA chip (Agilent, Santa Clara, CA, USA). A NEBNext FFPE DNA Repair Mix (New England Biolabs) was used to repair damaged bases prior to library preparation, which was carried out with a NEB library preparation kit. The shredding step was skipped because of the fragmented nature of historical DNA. All libraries were pooled and sequenced on an Illumina NextSeq 550 (75 bp paired-end) at the Genome Core Facility of New York University Abu Dhabi, UAE. The FASTx Toolkit (Gordon and Hannon 2010) was used to remove Illumina adaptors and low-quality reads (mean Phred score < 20). The final average read length was 63 bp after trimming (Suppl. material 2). The program FastQC (Andrews 2010) was then used to observe if base composition was biased towards the end of the raw reads, which is a common phenomenon resulting from deamination (Dabney et al. 2013) when sequencing older DNA, however this was not observed. Summary statistics describing the sequencing data are available in Suppl. material 2: Table S1. All sequences are deposited in GenBank (MW375737–MW375766; Suppl. material 2: Table S2).

Whole mitochondrial genomes of the type specimens were assembled from the Illumina reads using MITObim (Hahn et al. 2013). MITObim uses an iterative baiting method to generate mitochondrial contigs from short reads. First, a published sequence of the mitochondrial genome of Ptychadena mascareniensis (Duméril & Bibron, 1841) (GenBank reference number JX564890) was used as the reference mitogenome, with the default program settings, except for a k-mer length of 21. The analysis was then run again using the resulting contigs from the first MITObim run. An additional 21 mitochondrial genomes from fresh samples collected by us of members of the P. neumanni species complex were also assembled following the same protocol (Suppl. material 2: Table S2). These samples were sequenced for another project on the genomics of the Ptychadena neumanni species group (Hariyani et al. in prep.) following the sampling, tissue handling and DNA extraction protocols of Reyes-Velasco et al. (2018).

We reconstructed phylogenetic relationships within the Ptychadena neumanni species complex using three different datasets. First, all 13 mitochondrial protein-coding genes and the rRNA 12S and 16S genes were used. No stop codon was found in the protein-coding sequences. We excluded tRNAs because in some cases these were not complete or were difficult to align. Because we did not have the complete mitochondrial genome for all species, we used a subset of genes which was obtained for all species as a second dataset. This dataset included the 12S and 16S rRNA, as well as the protein-coding gene cytochrome c oxidase I (cox1). In the last dataset, we included only the rRNA 16S, in order to be able to include as many individuals as possible. Alignments at each gene were performed in MAFFT v. 7 (Katoh and Standley 2013), and other ambiguously aligned regions were removed using the online server G-Blocks (Castresana 2000) with the least stringent options selected. Geneious v. 9.1.6 (Biomatters Ltd., Auckland, NZ) was used to manually trim any remaining poorly aligned regions and to ensure that protein-coding genes were in the correct reading frame. Our final concatenated datasets consisted of 15,708 bp for the dataset that included all genes and 2522 bp for the reduced dataset (12S, 16S and cox1).

The best-fit model of nucleotide evolution for each gene was selected using the Bayesian Information Criterion (BIC) in PartitionFinder v. 1.1.1 (Lanfear et al. 2012; Suppl. material 2: Table S3). Data was partitioned by gene and by codon position in protein-coding genes. All genes were concatenated using the program Sequence Matrix (Vaidya et al. 2011) and performed Bayesian phylogenetic inference (BI) in MrBayes v. 3.2.2 (Ronquist et al. 2012) on the CIPRES Science Gateway server (Miller et al. 2010).

For both datasets, the BI analyses consisted of four runs of 10⁷ generations, sampling every 1000^th generation, with four chains (three heated and one cold). Convergence between the runs was assessed by visually inspecting overlap in likelihood and parameter estimates between the runs, as well as effective sample sizes and potential scale reduction factor (PSRF) value estimates for each run using Tracer v. 1.6 (Rambaut et al. 2014). Individual runs converged by 10⁵ generations, based on the PSRF, so the first 25% of each run were discarded as burn-in. The runs were combined, and the resulting tree was visualized in FigTree v. 1.4.2 (Rambaut 2014).

In a previous study (Reyes-Velasco et al. 2018), we sequenced thousands of loci from across the genome of all known lineages of Ptychadena inhabiting the Ethiopian highlands (12 putative species, 105 individuals) using double digest restriction-site associated DNA sequencing (ddRAD-seq). Here we briefly describe the methods used as more details were provided in our original article.

Individuals of the Ptychadena neumanni species complex were collected in the highlands of Ethiopia between 2011 and 2018. Tissue samples were extracted and stored in RNAlater or 95% ethanol. Genomic DNA was extracted with one of the following methods: with the use of a DNeasy blood and tissue kit (Qiagen, Valencia, CA), with the use of Serapure beads (Rohland and Reich 2012), or by standard potassium acetate extraction. DNA concentration was measured with a Qubit fluorometer (Life Technologies) so that DNA sample concentrations could be standardized. Genomic DNA was then digested with the enzymes SbfI and MspI (Peterson et al. 2012). Barcoded samples were size-selected between 400 and 550 base pairs using a Pippin Prep (Sage Science, Beverly, MA, USA), and attached to unique Illumina indices (Peterson et al. 2012). Libraries were sequenced on an Illumina HiSeq2500 (100 bp paired-end reads) at the Genome Core Facility of New York University Abu Dhabi, United Arab Emirates.

Ipyrad 0.6.17 (Eaton and Overcast 2020) was used to assemble loci de novo and create SNP datasets. After quality filtering, a total of ~158 million sequencing reads were retained, with a mean coverage of about 1.60 million reads, and a mean of ~11,400 RAD-tags per individual. We obtained between 800 and 2918 polymorphic loci and between 28,000–36,000 SNPs (Reyes-Velasco et al. 2018). The best model of evolution for our concatenated ddRAD-seq dataset was estimated with using BIC in PAUP v. 4.0.a151 (Swofford 1993), which showed the GTR + I + G model as the most supported. Maximum likelihood (ML) was implemented in RAxML v. 8 (Stamatakis 2014) to infer evolutionary relationships between populations and species in this group. RAxML was performed with rapid bootstrapping in the CIPRES portal (Miller et al. 2010).

DNA was successfully extracted from the type specimens of Ptychadena neumanni, P. erlangeri, P. largeni and P. nana. After quality filtering, a total of ~1.1 billion reads were retained, with highly variable coverage across individuals (49–359 million reads; Suppl. material 2: Table S1). The complete mitochondrial genome was recovered for five out of the six type specimens and a partial sequence was obtained for the holotype of P. nana. The total number of reads from the type specimens that mapped to the reference mitochondrial genome ranged between 802 (P. nana holotype) to > 900,000 (P. erlangeri holotype; Suppl. material 2: Table S1). We found no correlation between the amount of DNA recovered and the final number of reads for each specimen.

As the assignment of species names to genetic lineages was based on mitochondrial sequences, we first verified that the assignment of individuals to genetic lineages using mitochondrial markers was congruent with that obtained using genome-wide loci from ddRAD-seq (Fig. 2). The clustering of 105 individuals based on nuclear SNPs and on mitochondrial sequences was perfectly congruent (Figs 2, 3a, b); thus, we considered that, at least in this particular case, mitochondrial DNA was sufficient to assign samples to the genetic clusters identified by previous authors (Freilich et al. 2014; Smith et al. 2017a; Reyes-Velasco et al. 2018; Goutte et al. 2021).

Figure 2. 

Comparisons of the topologies of the mitochondrial rRNA 16S (left) and ddRAD-seq (right) for members of the Ptychadena neumanni species complex. Type specimens are indicated in red in the 16S phylogeny. Red lines indicate clades that differ in their placement between 16S and ddRAD-seq, however, the assignment of individuals to a particular species is identical between datasets. Numbers at nodes represent posterior support (pp), while black dots represent nodes with posterior support of 1.

Our phylogenetic analysis based on three mitochondrial genes (Fig. 3c) recovered 11 mitochondrial clusters corresponding to 11 of the 12 genetic lineages, with high support. The only exception was P. delphina, for which mitochondrial sequences did not form a clade (Fig. 3c). The lectotype and paralectotype of P. neumanni (ZMB 26879 and ZMB 57183) were nested with strong support (PP = 1) within the clade that had been referred to as P. erlangeri by multiple authors (Freilich et al. 2014; Smith et al. 2017a; Reyes-Velasco et al. 2018). The holotype of P. erlangeri (ZMB26887) and the paratypes of P. largeni (MHNG 2513–15 and 2513–56) all grouped with strong support (PP = 1) with individuals previously assigned to either P. cf. neumanni 2 (Freilich et al. 2014) or P. largeni (Smith et al. 2017a). This result demonstrates that P. largeni represents a junior synonym of P. erlangeri.

Figure 3. 

Phylogenetic inference of members of the Ptychadena neumanni species complex based on mtDNA and ddRAD-seq data A unrooted UPGMA tree of members of the P. neumanni species complex based on 2182 SNPs obtained with ddRAD-sequencing B unrooted Bayesian phylogenetic inference based on the complete mitochondrial genomes of members of the group C bayesian phylogenetic inference based on the concatenated sequences of the 12S and 16S rRNA and cox1. Black circles represent nodes with a posterior support of 1. Names in bold indicate type specimens, while stars indicate historical type specimens sequenced here and are color-coded as in Figure 1. Photographs represent members of the P. neumanni species complex; P. erlangeri (top), P. neumanni (bottom).

The holotype of Ptychadena nana (ZMB26878) did not group with individuals from the Bale Mountains identified as P. nana by previous authors (Freilich et al. 2014; Smith et al. 2017a; Reyes-Velasco et al. 2018; Fig. 3c) but instead grouped with individuals from the Arussi Plateau (= Didda Plateau; Fig. 1), which corresponds to the type locality for the species (Perret 1980). The holotypes for the three species described by Smith et al. (2017a), P. amharensis, P. goweri and P. levenorum, clustered with genetic lineages that did not include sequences derived from historical type specimens (Figs 2, 3) and thus constitute valid species. Finally, four lineages did not include any type specimen of species described by Smith et al. (2017a) or prior and the corresponding species were recently described by Goutte et al. (2021; P. beka, P. delphina, P. doro and P. robeensis). Notably, the lineage corresponding to P. cf. neumanni 1 of Freilich et al. (2014), which was previously suggested to be conspecific with P. neumanni by Smith et al. (2017a), was in fact genetically distinct and corresponded to the recently described P. beka.