On the edge of the Shivaliks: An insight into the origin and taxonomic position of Pakistani toads from the Duttaphrynus melanostictus complex (Amphibia, Bufonidae)

Daniel Jablonski; Rafaqat Masroor; Sylvia Hofmann

doi:10.3897/zse.98.79213

Research Article

On the edge of the Shivaliks: An insight into the origin and taxonomic position of Pakistani toads from the Duttaphrynus melanostictus complex (Amphibia, Bufonidae)

Daniel Jablonski^‡, Rafaqat Masroor^§, Sylvia Hofmann^|

‡ Comenius University in Bratislava, Bratislava, Slovakia

§ Zoological Sciences Division, Pakistan Museum of Natural History, Islamabad, Pakistan

| LIB - Leibniz Institute for the Analysis of Biodiversity Change, Museum Koenig, Bonn, Germany

Corresponding author: Daniel Jablonski ( daniel.jablonski@uniba.sk )

Academic editor: Umilaela Arifin

This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Citation: Jablonski D, Masroor R, Hofmann S (2022) On the edge of the Shivaliks: An insight into the origin and taxonomic position of Pakistani toads from the Duttaphrynus melanostictus complex (Amphibia, Bufonidae). Zoosystematics and Evolution 98(2): 275-284. https://doi.org/10.3897/zse.98.79213

ZooBank: urn:lsid:zoobank.org:pub:6623BB8C-693D-4836-91BD-E24FE3AABE2C

Abstract

The common Asian toad Duttaphrynus melanostictus (Schneider, 1799) complex has a wide distribution ranging from western foothills of the Himalaya to the easternmost range of the Wallacea, with the evidence of human-mediated introductions to some other areas. In the entire distribution range, the complex is formed by several evolutionary clades, distributed mostly in South-East Asia with unresolved taxonomy. In the northwestern edge of its distribution (Pakistan), the name D. melanostictus hazarensis (Khan, 2001) has been assigned to local populations but its biological basis remained, so far, understudied and unvalidated. Therefore, we re-evaluated the available genetic data (mitochondrial and nuclear) to show the relationships between Pakistani populations (including the type locality of D. m. hazarensis) and others from across the range. Our results showed that Pakistani populations are associated with one, deeply diverged, well-supported and widely distributed clade (so-called Duttaphrynus sp. 1 according to 16S, or clade B based on tRNAGly-ND3), that has already been detected in previous studies. This clade is further distributed in India, Nepal, Bangladesh, Malaysia, Singapore, and Indonesia and is characterized by a low level of genetic variability. This further suggests that both natural, as well as potential human-mediated dispersal, might have played an important role in setting up the current phylogeographic and distribution pattern of this clade. The clade is deeply divergent from other clades of the complex and represents a taxonomically unresolved entity. We here argue that the clade Duttaphrynus sp. 1/B represents a distinct species for which the name Duttaphrynus bengalensis (Daudin, 1802) comb. nov. is applicable, while the description of D. m. hazarensis does not satisfy the rules of the International Code of Zoological Nomenclature.

Key Words

Bufo, Bufonidae, Indian subcontinent, invasions, phylogeography, taxonomy

Introduction

The distribution of the Duttaphrynus melanostictus (Schneider, 1799) complex (Dubois and Ohler 1999) more or less corresponds to the area of the Oriental zoogeographic realm. It is a widespread group of toads with a geographic range throughout South and South-East Asia, southern China, including islands of Taiwan, Hainan, Sumatra, Java, Borneo, and Sumba (Frost 2021). The group has a high potential to spread into new localities and was artificially introduced to many regions of Asia as well as out of Asia (see Reilly et al. 2017; Vences et al. 2017; Othman et al. 2020; Soorae et al. 2020). These introductions are well-monitored because they can pose a serious threat to the local native fauna (Pettit et al. 2021).

The Indian subcontinent and northern Pakistan form the western border of the D. melanostictus complex natural range where its distribution is limited to the edge of Shivalik Hills, i.e. to several districts of Khyber Pakhtunkhwa, Punjab, and Azad Jammu and Kashmir Provinces (Mertens 1969; Khan 1972, 2001; Masroor 2012). Also, the presence of the complex in Sindh Province cannot be excluded (Masroor 2012) as this toad is known from several localities in Gujarat (India) close to the Pakistani border (Sarkar 1984). Such a distribution pattern signifies that D. melanostictus has been recorded in both the Palearctic and Oriental region, which are the main zoogeographic realms of Pakistan, (see Masroor 2012: fig. 20). Surprisingly, populations in the northwestern part of the range have been described as Duttaphrynus melanostictus hazarensis (Khan, 2001) based on morphological data. According to Khan (2001), specimens of these populations differ from other Asian populations of D. melanostictus by bearing bony ridge, kidney-shaped parotids, subarticular tubercles that are doubled under the penultimate phalanx of all fingers, the absence of rostral ridge, and a light brown dorsal coloration. However, the taxon description and comparison provided by Khan (2001) are very brief. Further, the original description neither provides a voucher nor a catalogue number for the holotype, nor information of the type collection or museum designation. Therefore, it does not fulfil the International Code of Zoological Nomenclature (Art. 16.4.2, ICZN), rendering the subspecific name D. m. hazarensis (Khan, 2001) invalid. In order to simplify the terminology in our text, we will, however, refer to the Pakistani population and related lineage as “hazarensis”.

Due to the topographic heterogeneity and habitat diversity, it might be plausible that Pakistani amphibian populations of different genera with the Oriental affiliation may hold a different evolutionary history compared to their populations from other parts of South-East Asia. The significance of the far-western part of the Himalaya and the Indo-Gangetic plain for understanding the biogeography of taxa distributed across these geomorphological domains has been recently shown in several species groups (see Allopaa, Chrysopaa, Microhyla, Nanorana, Scutiger, Sphaerotheca; Hofmann et al. 2017, 2019, 2021; Jablonski et al. 2020, 2021). However, except for Khan´s (2001) original work, populations of D. melanostictus from Pakistan, including toads from the type locality where the “hazarensis” population has been described, have not been studied so far in a broader biogeographical context (cf. Hussain et al. 2020; Akram et al. 2021). According to Khan (2001), the type locality has a wider range and is defined as “Ooghi, Mansehra and Datta, Hazara Division, eastern NWFP, Pakistan”, i.e., covering several villages and districts in northern Pakistan where these toads occur. Recent molecular studies focusing on the genus Duttaphrynus were based on dense sampling, covering a wide range of South-East Asia (mostly), and different DNA markers for D. melanostictus (see Hasan et al. 2014; Wogan et al. 2016; Reilly et al. 2017; Vences et al. 2017; Othman et al. 2020; Soorae et al. 2020; Bisht et al. 2021). Because the resulting DNA sequence variation showed deep genetic structuring, we reviewed and reanalysed these data together with new sequence data of D. melanostictus from different localities of Pakistan, including the locality of “hazarensis” (Datta, Khyber Pakhtunkhwa; 34.2962°N, 73.2649°E). Our main goals of the study are to (i) find out the genetic affiliation of the Pakistani population of D. melanostictus in the context of known molecular phylogeny of the complex, (ii) show phylogeographic relationships of the western populations to others inside the detected clade, and (iii) comment on the validity and taxonomy of the Pakistani taxon in light of molecular phylogeny and nomenclature of the D. melanostictus complex.

Materials and methods

Sampling and DNA extraction

Tissue samples of six individuals of the Duttaphrynus melanostictus complex from northern Pakistan were collected from five different localities of the Khyber Pakhtunkhwa Province during the field studies between 2018 and 2019. Individuals were documented photographically by DJ (Nikon D810, and 105 mm macro lens) and sampled for toe or blood before release or being deposited in scientific collections (Suppl. material 1: Table S1). Tissue samples were transferred into 96% ethanol, kept at ambient temperature during the time of the fieldwork, and later stored at -20 °C. DNA was isolated using the DNeasy Tissue Kit (QIAGEN), following the manufacturer’s protocols, and stored at -20 °C.

Amplification and sequencing

We amplified a fragment of the nuclear (nDNA) rhodopsin (Rho) gene (~316 bp), of the 16S rRNA (~560 bp), and approximately 469 bp of mitochondrial DNA (mtDNA) that encodes part of the cytochrome oxidase c subunit III (COIII), tRNA glycine (tRNAGly), and part of the NADH dehydrogenase subunit 3 (ND3) using the primers and PCR conditions following previous studies (Vences et al. 2017; Hofmann et al. 2019; Bisht et al. 2021). PCR products were Sanger-sequenced in both directions by Macrogen Europe (Amsterdam, The Netherlands; http://www.macrogen-europe.com) and sequence data edited where necessary in Sequencher 5.4 by eye.

Dataset and phylogenetic analyses

We complemented our new mitochondrial and nuclear DNA sequences dataset with published and unpublished sequence data from GenBank (for details see Suppl. material 1: Table S1). We mostly followed datasets provided by Wogan et al. (2016), Vences et al. (2017), and Bisht et al. (2021) as well as others with available sequences for D. melanostictus. To map the phylogeographic pattern of detected clades, we used published georeferenced localities of the sequence origin. If coordinates were not available for sequences, we georeferenced them according to the provided name of the locality or specimen voucher numbers connected to a particular museum collection. In some cases, sequences could not be georeferenced (e.g., “India”; especially for 16S dataset) or with uncertainty only (e.g., “Uttarakhand”, “Kalimantan”; details in Suppl. material 1: Table S1). For an overall view on the distribution of the complex we downloaded data for D. melanostictus (~20,000 records) available in Global Biodiversity Information Facility (GBIF, accessible via www.gbif.org).

New sequences were aligned in MEGA X v.11.0 (Kumar et al. 2018). We followed Vences et al. (2017) and trimmed the mitochondrial sequences, retaining 411 bp corresponding to 71 bp of tRNA-Gly and 340 bp of ND3. We excluded the following three 16S sequences of the data set of Bisht et al. (2021) due to low data quality and potential sequencing errors: EU071741, U52798, EU071765. Also, one sequence (AY458592) was excluded from the data set of Vences et al. (2017) due to uncertain geographic origin. The resulting data sets consisted of 150 sequences for 16S (560 bp), 496 sequences for tRNAGly-ND3 (411 bp), and ten sequences for Rho (300 bp) (Suppl. material 1: Table S1). The latter sequence alignment was used only for the haplotype network construction. All new sequences were deposited in GenBank (see Suppl. material 1: Table S1 for accession numbers).

Bayesian phylogenetic analyses were carried out with MrBayes v.3.2.6 (Ronquist et al. 2012) using the best-fit model of sequence evolution as determined for 16S (GTR+I+G) by jModeltest v.2.1.10 (Darriba et al. 2012) and tRNAGly-ND3 (four partitions: K80, K80+G, HKY+G, GTR+I+G) by PartitionFinder v.1.1.1 (Lanfear et al. 2012). We run MrBayes with a random starting tree for five million generations using four parallel Markov chain Monte Carlo (MCMC) simulations with four chains, sampling trees every 500 generations, until reaching an average standard deviation of split frequencies of <0.01. Stationarity and convergence of the runs were confirmed using the software Tracer v.1.7.2 (Rambaut et al. 2018). The first 25% of the samples of each run was discarded as burn-in.

For reasons of comparison, we also inferred the BI tree for 16S based on the secondary structure of the ribosomal RNA using RNAsalsa 0.8.1 (Stocsits et al. 2009). As the initial input, we used the constraint file based on the secondary structures of Bos taurus 16S rRNA which is provided with the RNAsalsa package. The sequence was then aligned with our 16S data using the MUSCLE algorithm in MEGA X. Starting with this initial alignment and the constraint file, RNAsalsa predicts RNA secondary structure and enhances structural alignment that results in a final multiple sequence alignment together with a consensus structure. Unambiguous stem pairs were then derived based on the consensus structure from RNAsalsa, specified in the MrBayes input file, and analysed under the doublet model (16×16) proposed by Schoninger and von Haeseler (1994).

To assess the evolutionary distance between our samples and the Duttaphrynus sister clades recovered from the phylogenetic analysis, we calculated the uncorrected p-distances for the 16S rRNA and tRNAGly-ND3 data set separately using Mega-X based on 1,000 bootstrap replications, with the pairwise deletion option, and by considering both transitions and transversions. For clades that correspond with the Pakistani population of D. melanostictus, and for the Rho sequence data we created a median-joining haplotype network using PopART (Leigh and Bryant 2015).

Results

The newly generated sequences of 16S, tRNAGly-ND3, and Rho from Pakistan belong to one, well-supported, monophyletic group in all investigated markers (Figs 1–3). This clade is known from previous studies, particularly Wogan et al. (2016), Vences et al. (2017) and Bisht et al. (2021).

Figure 1.

The molecular phylogenetic reconstruction of the Duttaphrynus melanostictus complex on 16S rRNA (A) sensu terminology of Bisht et al. (2021) and tRNAGly–ND3 (C) sensu Vences et al. (2017) together with median-joining networks (B, D) calculated from clades that include sequences of the “hazarensis” toads from Pakistan, now D. bengalensis comb. nov. (Daudin, 1802). The numbers on branches represent Bayesian Posterior Probabilities/bootstraps showing branch support. Sequence details are presented in Suppl. material 1: Table S1. The same colors of clades used for particular mtDNA trees reflect clade affiliation. The pictured individual represents a population from Shahdara, Margalla Hills, Islamabad, Pakistan.

For 16S, populations from Pakistan correspond with the clade called Duttaphrynus sp. 1 (sensu Bisht et al. [2021]) that includes populations from Bangladesh, India, Nepal, and surprisingly, Malaysia (continental as well as Borneo part; Figs 1A, 2A). This clade has a sister position to clade Duttaphrynus sp. 2 (sensu Bisht et al. 2021) including populations from China, Laos, Malaysia, Myanmar, and Vietnam. Overall, seven clades of the D. melanostictus complex were detected in the 16S phylogeny. Uncorrected p-distances between these clades range from 2.3 to 4.9% (Table 1). The genetic distance between Duttaphrynus sp. 2 and Duttaphrynus sp. 1 clades is 2.7%. The clade Duttaphrynus sp. 1 is formed by 12 haplotypes in 24 sequences with low internal sequence variation (π=1.0%). Sequences from Pakistan have a closer relationship (one mutation step) to sequences from peninsular India (KU589218, EU071749, EU071752, EU071754, EU071757, KY000453–KY000456) and Nepal (MF319517; Fig. 1B). Other sequences (haplotypes) differ by two to six mutation steps. Surprisingly, two sequences from Malaysia (including Borneo; AB530648, AB331714) are separated by only two mutation steps from Pakistani sequences, or are even identical to some sequences from India or Nepal (Figs 1B, 2A). Significantly, some sequences from India could not be georeferenced and, thus, not displayed on the map in Fig. 2A (see Suppl. material 1: Table S1). However, the overall phylogeographic pattern suggests a distribution of the clade Duttaphrynus sp. 1 in the northern part of the Indian subcontinent with a potential range limit in Myanmar, where it is replaced by the sister clade Dutaphrynus sp. 2. Populations of this ‘sp. 1’ clade in Western Malaysia and Borneo (Sarawak) might have been artificially introduced (see Discussion below and Fig. 2A).

Table 1.

Download as

CSV

XLSX

Uncorrected genetic distances (%) between investigated 16S clades of the Duttaphrynus melanostictus complex as shown in the Fig. 1A. Given are p-distance values (lower-left) and standard error estimates (upper-right).

Clade p-distances (%)	cf. microtympanum	brevirostris	parientalis	melanostictus ?	melanostictus	Duttaphrynus sp. 2	D. bengalensis comb. nov. (Duttaphrynus sp. 1)
cf. microtympanum	–	0.7	0.7	0.8	0.6	0.9	0.9
brevirostris	2.3	–	0.8	0.9	0.7	1.0	0.9
parientalis	2.7	3.0	–	0.6	0.7	0.9	0.9
melanostictus ?	2.9	3.4	2.3	–	0.6	0.8	0.8
melanostictus	2.8	3.0	2.9	2.5	–	1.0	0.8
Duttaphrynus sp. 2	4.1	4.3	4.7	4.3	4.9	–	0.6
D. bengalensis comb. nov. (Duttaphrynus sp. 1)	4.2	3.7	4.2	3.2	3.9	2.7	–

Figure 2.

Phylogeographic reconstruction and distribution of major clades of the Duttaphrynus melanostictus complex across the species range based on 16S rRNA (A) and tRNAGly–ND3 (B) sequences available in GenBank and from own data (see Suppl. material 1: Table S1). The clade colors follow those of Fig. 1, clade terminology is derived from Bisht et al. (2021; 16S rRNA), and Vences et al. (2017; tRNAGly–ND3). The pictured individual represents the population of Duttaphrynus bengalensis comb. nov. (Daudin, 1802) from Datta, Mansehra, Pakistan. Coloration reflect the main mtDNA clades according to Fig. 1. Distribution records (white circles) have been downloaded from GBIF (gbif.org). The red line (panel A) depicts the area considered as “Bengal”, the type locality of Bufo bengalensis.

For tRNAGly-ND3, Pakistani sequences correspond to the clade B (sensu Vences et al. [2017]) that includes populations from Indonesia (Kalimantan), West Malaysia, and Singapore (Figs 1C, 2B). This clade has a sister relationship with the clade C distributed in Myanmar. Overall, we confirmed seven clades of the D. melanostictus complex in the phylogeny of tRNAGly-ND3. Uncorrected genetic distances between clades range from 2.7 to 9.3% (Table 2). The genetic distance between closely related clade B and C is 5.6%. The clade B includes three haplotypes. Two of them are found in Pakistan and differentiated by one mutation step from each other. The third haplotype differs by four mutation steps from the Pakistani sequences and can be found in Indonesia (KU183404), Malaysia (KU183466, KU183468–KU183470), and Singapore (KU183475, KU183476, KU183484). Internal sequence variation in clade B is very low (π = 0.3%). Insufficient sampling across the Indian subcontinent did not allow for assessing range limits of the clade B in relation to other clades in the resulting phylogeny. However, detection of the clade B in Singapore, West Malaysia, and Indonesia, (Borneo; Kalimantan, Fig. 2B) suggests an artificial, probably human-mediated introduction.

Table 2.

Download as

CSV

XLSX

Uncorrected genetic distances (%) between investigated tRNAGly-ND3 clades of the Duttaphrynus melanostictus complex as shown in the Fig. 1C. Given are p-distance values (lower-left) and standard error estimates (upper-right).

Clade p-distances (%)	A	B	C	D	E	F	G
A	–	3.4	3.4	2.8	2.7	3.2	2.9
B	9.3	–	2.4	2.9	2.1	3.0	2.7
C	9.3	5.6	–	2.6	2.3	2.6	2.4
D	7.9	7.2	6.8	–	1.8	1.9	1.9
E	7.4	4.9	5.4	4.2	–	1.7	1.6
F	8.4	7.2	6.1	4.3	3.6	–	1.2
G	8.0	6.7	5.9	4.7	3.6	2.7	–

The results on rhodopsin (Rho) confirm mitochondrial data. It shows that D. melanostictus from Pakistan represents an independent group (Fig. 3). This group is formed by sequences from Pakistan and differs by three mutation steps from populations from Vietnam (DQ283967; based on the geographic position it probably represents the mtDNA clade Duttaphrynus sp. 2 in 16S, or clade G in tRNAGly-ND3, respectively) and by two to six mutation steps from D. melanostictus from India (AF249097; without precise locality). The Indian sequence cannot be related to any known mitochondrial phylogenetic clade due to insufficient sampling and unknown geographic origin. The uncorrected genetic distances between Rho sequences of the D. melanostictus complex range between 0.7 and 1.3%.

Figure 3.

Median-joining alleles network of sequences from the Duttaphrynus melanostictus complex and closely related members of western Palearctic Bufonidae based on rhodopsin (Rho) gene available in GenBank and from own data. The colors correspond to those used in Fig. 1. Note, the origin of sequence AF249097 from India cannot be assigned due to unknown geographic details. The pictured individual represents a population of Duttaphrynus bengalensis comb. nov. (Daudin, 1802) from Datta, Mansehra, Pakistan.

Discussion

The most prominent results of our study are the following: (i) the “hazarensis” population from Pakistan belongs to a widely distributed and taxonomically so far unresolved clade (see below), (ii) Pakistani mitochondrial sequences represent unique haplotypes with a close relationship to other populations from northern parts of the Indian subcontinent, (iii) comparing the overall phylogeographic pattern, the clade Duttaphrynus sp. 1/clade B has been most likely artificially translocated to South-East Asia, probably from the Indian subcontinent, (iv) the name for the clade Duttaphrynus sp. 1/clade B is available and, thus, should be applied (see Taxonomy). These results are thus important for both, evolution and taxonomy of the D. melanostictus complex. However, they still underline previous conclusions that the complex needs comprehensive biogeographical research with taxonomic revision based on dense DNA sampling, especially from the Indian subcontinent (see Wogan et al. 2016; Vences et al. 2017; Bisht et al. 2021).

Genetic affiliation and origin of Pakistani populations

By re-analysing 16S sequences of the D. melanostictus complex we showed a pattern of three main, geographically widely distributed clades (called as Duttaphrynus melanostictus, Duttaphrynus sp. 1, and sp. 2 sensu Bisht et al. 2021) with potential distribution boundary between them in western Myanmar and Malay Peninsula/Sumatra (Fig. 2). Populations from Pakistan belong to the Duttaphrynus sp. 1 clade that is distributed in the northern part of the Indian subcontinent (Fig. 1; Bisht et al. 2021). Unfortunately, limited sampling in central India does not allow delimitating the range of this clade in the southern parts of India. The Indian subcontinent is, however, absolutely crucial to understand the range of this clade as southern India is inhabited by different evolutionary clades of the complex, i.e. D. brevirostris, D. microtympanum, and D. melanostictus (sensu Bisht et al. 2021; Figs 1, 2). Re-evaluation of tRNAGly-ND3 sequence data of the complex revealed similar results as for 16S, showing two widely distributed clades (A and G) and four geographically disjunct clades assigned to the northern and western parts of Indochina (C, E, D, F; Wogan et al. 2016; Vences et al. 2017; Figs 1, 2). Clade B (sensu Vences et al. 2017 as a part of the “mainland clade” sensu Wogan et al. 2016), where Pakistani toads belong to, was previously reported from Indonesia (Kalimantan, Borneo), Malaysia, and Singapore (Wogan et al. 2016; Vences et al. 2017). However, since these former studies did not focus on this specific clade, no conclusions were drawn about its distribution area. Despite limited data from the Indian subcontinent, the genetic clustering of Pakistani samples within clade B, and the similar phylogenetic placement observed in the 16S tree suggest a wide distribution range of this clade. The low level of intra-clade sequence differentiation and star-like pattern of the haplotype networks implies a very close relationship between geographically distant populations and potentially rapid expansions to wider areas, supplemented by recent human-mediated introductions (see below).

For the clade Duttaphrynus sp. 1/clade B, we revealed a phylogeographic structure resembling patterns detected in other widely distributed amphibians of the Indo-Gangetic plain. These species colonized Pakistan from the Oriental region, namely Microhyla nilphamariensis Howlader, Nair, Gopalan, Merilä, 2015 and Sphaerotheca maskeyi (Schleich & Anders, 1998) (Jablonski et al. 2020, 2021). Both species are wide-ranging, occupying central and northern parts of the Indian subcontinent, from eastern India, through Bangladesh, Nepal, central, northern, and western India to Pakistan. Their distribution more or less follows the Indo-Gangetic plain including the Brahmaputra, Yamuna, Ganga, and Indus rivers. Thus, we assume that the clade Duttaphrynus sp. 1/clade B colonized new areas along these riverine corridors (probably from the northern India and Nepal area; see the structure of the haplotype network Fig. 1B) and that this clade can also be found in central and western India, eventually also southern Pakistan (see Sarkar 1984). However, as shown by the model of the potential distribution of the complex in Asia (see supplementary fig. S2 in Vences et al. 2017), the distribution of these toads might be scattered in central India and denser to the region of the Shivalik Hills (see also the georeferenced records from the GBIF database Fig. 2). The Shivaliks (also known as the Churia Hills) may provide optimal environmental conditions for these toads (or this specific evolutionary clade). It is a mountain range of the outer Himalayas that stretches from Pakistan about 2,400 km eastwards (10–50 km wide) close to the Brahmaputra River, spanning across the northern parts of the Indian subcontinent. As the western edge of these hills is currently inhabited by the clade Duttaphrynus sp. 1/clade B we can assume that such an environmental corridor (the Siwalik Corridor sensu Khan 2006) also affects the distribution of this species complex, enabling migration from the east to the west.

Recently, two basic biogeographic scenarios have been hypothesized for the origin of different amphibian taxa distributed in northern Pakistan: i) dispersal of ancestral lineages during the Miocene from East or Southeast China into the far northwestern part of the Himalaya-Tibet orogen along a warm temperate corridor that existed in significant parts of Paleo-Tibet’s interior (Hofmann et al. 2021), this corridor played probably an important role in the initial diversification process of Himalayan faunal organisms (Schmidt et al. 2012), (ii) dispersion on widely distributed species characterized by the low level of genetic divergence that followed Shivalik Hills and the Indo-Gangetic plain and reached the territory of Pakistan (Jablonski et al. 2020, 2021) as a result of the Pleistocene or Holocene colonisations (see Othman et al. 2020). Although these are presumable biogeographic models, the second scenario could be now applied for D. melanostictus from Pakistan.

It is well known that the D. melanostictus complex is considered an invasive organism, with several evolutionary clades that were introduced out of its native range (Reilly et al. 2017; Vences et al. 2017; Othman et al. 2020; Soorae et al. 2020; Pettit et al. 2021). Evidence for potentially human-mediated introductions is visible from our analysed data, showing large geographic distances between the source and apparently introduced populations (Fig. 2). Consistent with previous studies (Reilly et al. 2017; Vences et al. 2017; Othman et al. 2020; Soorae et al. 2020), probably clades A and G were sources of introduced populations (sensu Vences et al. 2017, Fig. 1). These clades are widely distributed and covered by dense sampling. Moreover, these clades are distributed also in and around populated cities and ports of South-East Asia from where the artificial, human-meditated introduction to other areas is highly evident (even to Madagascar or the United Arab Emirates; Vences et al. 2017; Soorae et al. 2020). On the other hand, limited sampling of the clade B (sensu Vences et al. 2017 which is equal to clade Duttaphrynus sp. 1 sensu Bisht et al. [2021] in the 16S dataset) did not allow to evaluate the origin of this clade in the Malay Peninsula or Borneo. Based on the combination of 16S and tRNAGly-ND3 data and the distribution pattern of both detected clades (Duttaphrynus sp. 1/clade B) we can hypothesize that Malaysian, Singapore, and Indonesian populations are results of the expansion of the range during ice-ages followed by recent introductions due to human activities (Othman et al. 2020) from a hitherto unknown region of the Indian subcontinent (Fig. 2).

Taxonomy

All the hereby presented markers and overall phylogenetic pattern of the clade to which sequences of D. melanostictus from Pakistan belong to (see Figs 1–3) support previous studies discussing the unsettled taxonomy of this toad complex in its entire distribution (Wogan et al. 2016; Bisht et al. 2021). The taxonomy of this group is difficult to resolve due to missing data from different regions, especially from the western part of the range, and similar morphology. In our study, we provide, for the first time, the evidence that the “hazarensis” population from Pakistan belongs to the widely distributed clade specified as Duttaphrynus sp. 1 (16S, sensu Bisht et al. 2021), or clade B (tRNAGly-ND3; sensu Vences et al. 2017), respectively. It implies that Khan (2001) correctly observed morphological differences on these toads (see Introduction) but the lack of broader comparative material probably did not allow to compare specimens from Pakistan with other specimens, especially from the Indian subcontinent. The description, moreover, is formally invalid, since the holotype lacks information about the name and location of the repository where it has been deposited (Art. 16.4, ICZN 1999). Khan (2001) compared the Pakistani population with island populations from Taiwan (which, genetically, probably belong to Duttaphrynus sp. 2 clade in the 16S tree, and clade G in the tRNAGly-ND3 tree, respectively) and Sri Lanka from where genetic data are lacking.

Following the recommendation of Dubois and Ohler (1999), i.e. “If further studies demonstrate that a distinct species, close to Bufo melanostictus, occurs in Bengal, the name Bufo bengalensis Daudin, 1802 is available for it”, and based on the 16S phylogeny, the geographic distribution of Duttaphrynus sp. 1 clade (sensu Bisht et al. 2021), and the Principle of priority rule (ICZN 1999), we suggest that the name Duttaphrynus bengalensis (Daudin, 1802) comb. nov. (type locality “Bengale”, neotype MNHN 4967; for description see Dubois and Ohler 1999) is applicable for this clade. “Bengale” is equal to the area that is divided between Bangladesh and several states in neighbouring India (states West Bengal, Jharkhand, Bihar and the Cachar region of Assam matching the distribution of the Duttaphrynus sp. 1 clade (Fig. 2A). According to our data, the clade is native to Bangladesh, India, Nepal and Pakistan, and apparently introduced to Indonesia, Malaysia, and Singapore. (see above and Fig. 2). Morphologically, the descriptions of D. m. hazarensis (Khan, 2001) and D. bengalensis (Daudin, 1802) comb. nov. are almost identical in several aspects: the interorbital space is concave in both “hazarensis” and D. bengalensis comb. nov., the space is larger than the width of the upper eyelid and internarial distance. Tympanic oblique, its diameter in both these species is about two thirds of the eye diameter; the tympanic-eye distance is about less than one half of its distance. Moreover, the larger size of “hazarensis” (max. SVL of female up to 111.0 mm) and D. bengalensis comb. nov. (max. SVL of female neotype up to 133.2 mm) differentiates this species from D. melanosctictus (Schneider, 1799) specimens from India whose maximum SVL is 71.6 mm (lectotype ZMB 3462; Dubois and Ohler 1999). The snout length in D. bengalensis comb. nov. and “hazarensis” is longer than the horizontal diameter of eye (vs. shorter in D. melanostictus). The canthus rostralis in D. bengalensis comb. nov. and “hazarensis” is sharp, with a ridge in comparison to D. melanostictus whose canthus rostralis is angular. The tympanum diameter in D. bengalensis comb. nov. and “hazarensis” is oval, oblique, its diameter about two-thirds of the eye diameter, whereas it is vertical and its diameter is more than half of the eye diameter in D. melanostictus. The tympanum-eye distance in D. bengalensis comb. nov. and “hazarensis” is about less than one half of its diameter, whereas it is less than one sixth of its diameter in D. melanostictus. Thus we recommend to assign the name D. bengalensis (Daudin, 1802) for the clade Duttaphrynus sp. 1/clade B for future studies.

Acknowledgements

We would like to thank our colleagues, friends and local people for their help during the field work in Pakistan. For technical and laboratory support we thank Martina Lawson and Jana Poláková. For recommendations regarding the taxonomy we thank Frank Tillack. We are also grateful to S. R. Chandramouli and one anonymous reviewer for their useful comments on the manuscript. This work was supported to DJ by the Slovak Research and Development Agency under the contract APVV-19-0076, and by the German Research Foundation (DFG, grant no. HO 3792/8-1) to SH.

References

Akram A, Rais M, Lopez-Hervas K, Tarvin RD, Saeed M, Bolnick DI, Cannatella DC (2021) An insight into molecular taxonomy of bufonids, microhylids, and dicroglossid frogs: First genetic records from Pakistan. Ecology and Evolution 11(20): 14175–14216. https://doi.org/10.1002/ece3.8134

Bisht K, Garg S, Akalabya Sarmah AND, Sengupta S, Biju SD (2021) Lost, forgotten, and overlooked: Systematic reassessment of two lesser-known toad species (Anura, Bufonidae) from Peninsular India and another wide-ranging northern species. Zoosystematics and Evolution 97(2): 451–470. https://doi.org/10.3897/zse.97.61770

Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: More models, new heuristics and parallel computing. Nature Methods 9(8): 772–772. https://doi.org/10.1038/nmeth.2109

Dubois A, Ohler A (1999) Asian and Oriental toads of the Bufo melanostictus, Bufo scaber and Bufo stejnegeri groups (Amphibia, Anura): A list of available and valid names and redescription of some name-bearing types. Journal of South Asian Natural History, Colombo 4: 133–180.

Frost DR (2021) Amphibian Species of the World: an Online Reference. Version 6.1 (23 November 2021). Electronic Database accessible at https://amphibiansoftheworld.amnh.org/index.php. American Museum of Natural History, New York, USA.

GBIF.org (4 December 2021) GBIF Occurrence Download: https://doi.org/10.15468/dl.tdbg7m

Hasan M, Islam MM, Khan MMR, Igawa T, Alam MS, Djong HT, Kurabayashi A, Joshy H, Sen YH, Belabut DM, Kurabayashi A, Kuramoto M, Sumida M (2014) Genetic divergences of South and Southeast Asian frogs: a case study of several taxa based on 16S ribosomal RNA gene data with notes on the generic name Fejervarya. Turkish Journal of Zoology 38(4): 389–411. https://doi.org/10.3906/zoo-1308-36

Hofmann S, Stöck M, Zheng Y, Ficetola FG, Li JT, Scheidt U, Schmidt J (2017) Molecular phylogenies indicate a paleo-Tibetan origin of Himalayan lazy toads (Scutiger). Scientific Reports 7: 3308. https://doi.org/10.1038/s41598-017-03395-4

Hofmann S, Baniya CB, Litvinchuk SN, Miehe G, Li JT, Schmidt J (2019) Phylogeny of spiny frogs Nanorana (Anura: Dicroglossidae) supports a Tibetan origin of a Himalayan species group. Ecology and Evolution 9(24): 14498–14511. https://doi.org/10.1002/ece3.5909

Hofmann S, Jablonski D, Litvinchuk SN, Masroor R, Schmidt J (2021) Relict groups of spiny frogs indicate Late Paleogene-Early Neogene trans-Tibet dispersal of thermophile faunal elements. PeerJ 9: e11793. https://doi.org/10.7717/peerj.11793

Hussain S, Bukhari SM, Javid A, Hussain A, Rashid M, Ali W (2020) Molecular identification of genus Duttaphrynus from Punjab, Pakistan. Mitochondrial DNA. Part B, Resources 5(3): 3218–3220. https://doi.org/10.1080/23802359.2020.1810143

ICZN [International Commission on Zoological Nomenclature] (1999) International Code of Zoological Nomenclature, 4^th Edn. London, International Trust for Zoological Nomenclature, [xxix +] 306 pp.

Jablonski D, Khan MA, Masroor R (2020) The genus Microhyla (Anura: Microhylidae) in Pakistan: species status and origins. Zootaxa 4845(2): 293–296. https://doi.org/10.11646/zootaxa.4845.2.11

Jablonski D, Masroor R, Hofmann S (2021) Revisited molecular phylogeny of the genus Sphaerotheca (Anura: Dicroglossidae): the biogeographic status of northernmost populations and further taxonomic changes. Diversity (Basel) 13(5): e216. https://doi.org/10.3390/d13050216

Khan MS (1972) The “commonest toad” of West Pakistan and a note on Bufo melanostictus Schneider. Biologia (Lahore) 18: 131–133.

Khan MS (2001) Notes on cranial-ridged toads of Pakistan and description of a new subspecies (Amphibia: Bufonidae). Pakistan Journal of Zoology 33: 293–298.

Khan MS (2006) Amphibians and reptiles of Pakistan. Krieger Publishing Company, Malabar, Florida, 311 pp.

Kumar S, Stecher G, Li M, Knyaz C, Tamura K (2018) MEGA X: Molecular evolutionary genetics analysis across computing platforms. Molecular Biology and Evolution 35(6): 1547–1549. https://doi.org/10.1093/molbev/msy096

Lanfear R, Calcott B, Ho SY, Guindon S (2012) Partitionfinder: Combined selection of partitioning schemes and substitution models for phylogenetic analyses. Molecular Biology and Evolution 29(6): 1695–1701. https://doi.org/10.1093/molbev/mss020

Leigh JW, Bryant D (2015) PopART: Full-feature software for haplotype network construction. Methods in Ecology and Evolution 6(9): 1110–1116. https://doi.org/10.1111/2041-210X.12410

Masroor R (2012) A contribution to the herpetology of northern Pakistan: The amphibians and reptiles of Margalla Hills National Park and surrounding regions. Society for the Study of Amphibians and Reptiles (SSAR), Ithaca, New York and Chimaira Buchhandelsgesellschaft mbH, Edition, Frankfurt, 217 pp.

Mertens R (1969) Die Amphibien und Reptilien West-Pakistans. Stuttgarter Beiträge zur Naturkunde 197: 1–96.

Othman SN, Chen Y-H, Chuang M-F, Andersen D, Jang Y, Borzée A (2020) Impact of the Mid-Pleistocene revolution and anthropogenic factors on the dispersion of Asian Black-spined toads (Duttaphrynus melanostictus). Animals (Basel) 10(7): e1157. https://doi.org/10.3390/ani10071157

Pettit L, Somaweera R, Kaiser S, Ward-Fear G, Shine R (2021) The impact of invasive toads (Bufonidae) on monitor lizards (Varanidae): An overview and prospectus. The Quarterly Review of Biology 96(2): 106–125. https://doi.org/10.1086/714483

Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA (2018) Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Systematic Biology 67(5): 901–904. https://doi.org/10.1093/sysbio/syy032

Reilly SB, Wogan GOU, Stubbs AL, Arida E, Iskandar DT, McGuire JA (2017) Toxic toad invasion of Wallacea: A biodiversity hotspot characterized by extraordinary endemism. Global Change Biology 23(12): 5029–5031. https://doi.org/10.1111/gcb.13877

Ronquist F, Teslenko M, Van der Mark P, Ayres DL, Darling A, Hohna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61(3): 539–542. https://doi.org/10.1093/sysbio/sys029

Sarkar AK (1984) Ecological studies on the amphibians of Gujarat. Bulletin Zoological Society India 6: 87–93.

Schmidt J, Opgenoorth L, Holl S, Bastrop R (2012) Into the Himalayan exile: The phylogeography of the ground beetle Ethira clade supports the Tibetan origin of forest-dwelling Himalayan species groups. PLoS ONE 7(9): e45482. https://doi.org/10.1371/journal.pone.0045482

Schoninger M, von Haeseler A (1994) A stochastic model and the evolution of autocorrelated DNA sequences. Molecular Phylogenetics and Evolution 3: 2–240. https://doi.org/10.1006/mpev.1994.1026

Soorae PS, Frankham GJ, Mohamed AA (2020) The first record of the Asian common toad Duttaphrynus melanostictus Schneider, 1799 in Abu Dhabi, United Arab Emirates. BioInvasions Records 9(2): 434–443. https://doi.org/10.3391/bir.2020.9.2.29

Stocsits RR, Letsch H, Hertel J, Misof B, Stadler PF (2009) Accurate and efficient reconstruction of deep phylogenies from structured RNAs. Nucleic Acids Research 37(18): 6184–6193. https://doi.org/10.1093/nar/gkp600

Vences M, Brown JL, Lathrop A, Rosa GM, Cameron A, Crottini A, Dolch R, Edmonds D, Freeman KLM, Glaw F, Grismer LL, Litvinchuk S, Milne MG, Moore M, Solofo JF, Noël J, Nguyen TQ, Ohler A, Randrianantoandro C, Raselimanana AP, van Leeuwen P, Wogan GOU, Ziegler T, Andreone F, Murphy RW (2017) Tracing a toad invasion: Lack of mitochondrial DNA variation, haplotype origins, and potential distribution of introduced Duttaphrynus melanostictus in Madagascar. Amphibia-Reptilia 38(2): 197–207. https://doi.org/10.1163/15685381-00003104

Wogan GO, Stuart BL, Iskandar DT, McGuire JA (2016) Deep genetic structure and ecological divergence in a widespread human commensal toad. Biology Letters 12(1): e20150807. https://doi.org/10.1098/rsbl.2015.0807

Supplementary material

Supplementary material 1

Table S1

Daniel Jablonski, Rafaqat Masroor, Sylvia Hofmann

Data type: excel file

Explanation note: Dataset of sequences used for the study.

This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.

Download file (62.25 kb)