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

The common Asian toad Duttaphrynus melanostictus (Schneider, 1799) complex has a wide distribution ranging from western foot-hills 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.


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 1972Khan , 2001Masroor 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. 2017Hofmann et al. , 2019Hofmann et al. , 2021Jablonski et al. 2020Jablonski et al. , 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 (1999).
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).
For 16S 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).
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 differenti-ated 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.
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%.

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 (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 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 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.  (Jablonski et al. 2020. 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 , 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 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] 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 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.