Species boundaries, biogeography and evolutionarily significant units in dwarf toads: Duttaphrynus scaber and D. atukoralei (Bufonidae: Adenominae)

Species boundaries and patterns of gene flow in Dwarf toads, Duttaphrynus scaber and D. atukoralei, were assessed using mitochondrial DNA markers. Samples from four populations in Sri Lanka (Mihintale, Ampara, Yala, Galle) were analyzed for three mitochondrial gene fragments (16S rRNA, COI and Cyt b) along with four Genbank sequences of 16S rRNA from Indian samples (Thiruvananthapuram, Maharashtra, Mudigere). Phylogenetic trees and haplotype networks were generated, and morphology was assessed. Analyses suggest a single species (Duttaphrynus scaber) with three major clades: a widespread clade shared between India and Northern Sri Lanka, an Eastern and Southeastern Sri Lankan clade (previously referred to as D. atukoralei, the validity of which, however, our analysis disputes), and a distinct Southern wet-zone clade from Galle (referred previously to as D. atukoralei). Duttaphrynus atukoralei (topotypes from Yala, Sri Lanka) is genetically too close to D. scaber (Indian and northern Sri Lankan clade) to be distinguished as a species; these two clades have a genetic distance of 0.95 – 1.55% for the 16S rRNA fragment. The haplotype networks for the 16S rRNA gene suggest incomplete lineage sorting between the Ampara and Yala populations; COI and Cyt b show complete sorting for all populations analyzed, suggesting strong population structure. All analyses suggest substantially restricted gene flow to the southern wet-zone population (Galle). This population also assumes a basal phylogenetic position, suggesting that D. scaber first evolved in southern Sri Lanka’s wet zone and dispersed across the lowland areas of the island and to India. Here, we provisionally recognize this population (Galle) as an evolutionarily significant unit of D. scaber; future analyses using multiple criteria may indicate this to be a new Dwarf toad species. External morphology is largely uninformative as the Yala, Ampara and Galle populations cannot be distinguished from each other; the morphological distinction between Yala, Ampara, Galle versus Mihintale is restricted to only the shape of the parotid glands – slightly oval versus rounded – a minor difference. Both genetic and morphological evidence so far suggest that there is only a single Dwarf toad species in Sri Lanka, which is also shared with India, namely Duttaphrynus scaber; however, with strong population structure, including an evolutionarily significant unit (Southern wet-zone population).


INTRODUCTION
Sri Lanka is a global amphibian hotspot with many unique and threatened taxa, especially of the genus Pseudophilautus (shrub frogs), which have undergone a large endemic radiation resulting in some 100 species (Meegaskumbura et al., 2002). Bossuyt et al. (2004), in an analysis of several vertebrate (shrub frogs, caecilians, earth snakes and freshwater fish) and invertebrate (freshwater crabs and shrimps) groups representing a diversity of lifehistory strategies showed that the wet-adapted taxa are characterized by clade level endemism. For wet-adapted frogs, this notion was reinforced with the recognition of a Sri Lankan endemic genus, Taruga, (Meegaskumbura et al., 2010), delineation of reciprocal endemicity of Hylarana species of India and Sri Lanka (Biju et al., 2014), endemicity of Nannophrys (which exhibits several unique features , and the endemicity of the genus Adenomus Van Bocxlaer et al., 2009). However, it is not known if this pattern is representative of the other anuran taxa, especially for some of the dry-adapted species, which are putatively shared with India.
Members of the family Bufonidae, commonly known as true toads, with over 500 species, show an extensive global distribution (Pramuk et al., 2008). Using molecular data, Van Bocxlaer et al. (2009) allocated some of the South and Southeast Asian bufonids to the subfamily Adenominae Cope, 1861. Although Indian bufonids comprise six genera, Sri Lanka has only two: Adenomus Cope, 1861, endemic to the island, and Duttaphrynus Frost, Grant, Faivovich, et al., 2006, of which some species are shared with India Van Bocxlaer et al., 2009). The Dwarf Toads, a non-taxonomic term used to recognize the small toads inhabiting India and Sri Lanka include Duttaphrynus scaber (Schneider, 1799), a species shared with India, and D. atukoralei, a species thought to be endemic to Sri Lanka. The Dwarf Toads, according to published work (Van Bocxlaer et al., 2009;Meegaskumbura et al., 2015) form a well-supported clade.
In the original description of Duttaphrynus atukoralei, Bogert & Senanayake (1966) considered D. scaber (referred to as Bufo fergusonii Boulenger, 1892, until Dubois andOhler (1999) showed it to be a junior subjective synonym of D. scaber) to be its sister species. The type locality of Duttaphrynus scaber is Thiruvananthapuram, Kerala, South India, while that of D. atukoralei is Yala, southeastern Sri Lanka (Bogert and Senanayake, 1966). Within Sri Lanka, however, the two species are thought to overlap in distribution across the lowland dry zone of the eastern region: D. scaber shows a northwestern, northern and eastern distribution while D. atukoralei appears to be restricted to the eastern and southeastern segments of the island despite the habitat being uniform, without apparent barriers to dispersal (Bogert and Senanayake, 1966).
The two species are morphologically similar, distinguished largely by the shape of the parotid glands: lobulated in D. atukoralei, rounded in D. scaber (Bogert and Senanayake, 1966). Though Bogert and Senanayake (1966) noted that the call of D. atukoralei was different from that of D. scaber (= B. fergusonii), they only sampled the call of D. scaber from Sri Lanka.
The partial sympatry of Duttaphrynus scaber and D. atukoralei, the shared distribution of D. scaber between the Indian mainland and Sri Lanka, and the occurrence of D. atukoralei in the wet zone (Galle), in the context of the clade-level endemicity of wet-adapted species observed by Bossuyt et al. (2004), gives rise to several questions regarding the taxonomy and biogeography of these taxa, such as: is there congruity between molecular and morphological data in the context of species boundaries? Is there clade-level endemism shown by dry zone elements? What are the patterns of haplotype structures for these species?
Here we examine the species boundaries and patterns of gene flow of D. scaber and D. atukoralei in Sri Lanka using three mitochondrial genes (16S ribosomal RNA -16S rRNA, Cytochrome b -Cyt b and Cytochrome c oxidase subunit 1 -COI).

MATERIALS AND METHODS
This study has been cleared by the Ethical Clearance Committee, Postgraduate Institute of Science, University of Peradeniya at its 18th meeting held on 19th May 2015.

Sample localities and tissue collection
Specimens of Duttaphrynus atukoralei and D. scaber were collected from four populations (Mihintale, Galle, Ampara and Yala) between January and December 2014 (Table 1). Three samples each from Mihintale, Galle and Yala, and two from Ampara, were collected during the course of the survey. The specimens from Galle, Yala and Ampara were identified as D. atukoralei on the basis of the shape of their parotid glands (Bogert and Senanayake, 1966) while those from Mihintale were assigned to D. scaber. Tissue samples taken from either thigh muscles or toe tips and were preserved in 100% ethanol. All specimens and tissue samples were deposited in the natural history collection of the Department of Zoology (DZ), University of Peradeniya.

DNA extraction, PCR amplification, DNA sequencing and sequence alignment
DNA was extracted from ethanol-preserved tissues using Qiagen DNeasy tissue extraction kits following manufacturer's protocol. Mitochondrial 16S rRNA, Cyt b and COI gene fragments were amplified using standard PCR conditions and primers; 16Sar and 16Sbr (Palumbi, 1996), which amplified a ~550 bp region of 16S rRNA gene; CB-J-10933 and BSF4 (Bossuyt and Milinkovitch, 2000), which amplified a ~ 590 bp of Cyt b gene; chmf4 and chmr4 (Che et al., 2012), which amplified a ~ 640 bp of COI gene. PCR conditions for 16S rRNA and Cyt b were as follows: denaturation at 95° C for 30 s, annealing at 45° C for 40 s and extension at 72° C for 40 s, 35 cycles, with a final extension of 72° C for 7 min and PCR conditions for COI as; denaturation at 95° C for 1 min, annealing at 48° C for 1 min and extension at 72° C for 1 min, 35 cycles, with a final extension of 72° C for 10 min. PCR Products were purified using QIAquick (Qiagen) PCR purification kit and sequenced on an ABI 3500 automated sequencer following manufacturer's protocols.

Phylogenetic analyses
Our study included 38 species, representing Adenominae. The data set contained 16 samples of Duttaphrynus atukoralei and D. scaber sampled from the four locations (Table 1). Additionally, 22 sequences of 16S rRNA were obtained from the GenBank (Table 1). Four species of the genus Ansonia (Ansonia spinulifer, A. malayana, A. hanitschi, A. fuiginea and A. leptopus) served as the outgroup; this is the closest sister group of the family Adenominae (Van Bocxlaer et al., 2009). Sequences were aligned using ClustalW as implemented in MEGA v. 5.0 (Tamura et al., 2011), and ambiguous regions were identified. A DNA fragment of 16S rRNA ca. 467 bp was used in the subsequent phylogenetic analyses. Bayesian (BA) and maximum likelihood (ML) analyses were performed to infer relationships among clades and assess node support. The best-fit model for the dataset was chosen using jModelTest v. 2.1.4 (Darriba et al., 2012;Guindon and Gascuel, 2003). Maximum likelihood analysis was performed using the software GARLI (Zwickl, 2006) on the Cipres Science Gateway (Miller et al., 2010), using the model selected by jModeltest (GTR+I+G) with all parameters estimated. The run was repeated twice to ascertain the tree topology. The tree was visualized using FigTree v. 1.3.1 (Rambaut and Drummond, 2010). Bayesian inference as implemented in MrBayes v. 3.1.2 (Huelsenbeck & Ronquist, 2001) was used to assess posterior probability values for each node with the parameters of a GTR model of sequence evolution with gamma-distributed rate variation among sites and a proportion of invariant sites (GTR+I+G) estimated as obtained from the jModelTest. Four Metropolis-Coupled Markov Chain Monte Carlo (MCMCMC) chains were run for ten million generations. Burn-in of 5 million generations was estimated using Tracer v. 1.6 (Rambaut et al., 2014). Bootstrapping was done in a ML framework using RAxML 8 (Stamatakis, 2014). GTRGAMMA model was used with one thousand RAxML searches and 1000 iterations. Clade support was assessed using posterior probability (PP) and ML bootstrapping values.

Haplotype network reconstruction
The TCS 1.21 (Clement et al., 2000) in PopArt (Leigh & Bryant, 2014) was used to reconstruct the haplotype networks for 16S rRNA, Cyt b, and COI gene sequences. TCS is a computer program that implements the estimation of gene genealogies from DNA sequences as described by Templeton et al. (1992).

Molecular diversity indices
Nucleotide diversity (Nd) of Nei (1987), the average number of nucleotide differences per site between sequences, and haplotype diversity (Hd) were calculated for the all the three genes and separately for each population using the program DnaSP, version 5.10.01 (Librado & Rozas, 2009).

Phylogenetic analysis
Sequence data for 467 bp of 16S rRNA were obtained for a total of 38 species of Adenominae, including D. scaber and D. atukoralei. jModelTest results shows that GTR+I+G is the best fit model for the data set. Maximum Likelihood inference was performed using the GTR+I+G model with the following parameters applied: The uncorrected pairwise genetic distance between D. scaber (Indian and Mihintale population) and D. atukoralei (Yala population; topotypes) were between 0.97 -1.55% for the 16S rRNA fragment (Fig. 1C). The highest average pairwise genetic distances (16S rRNA; 2.90% COI; 5.37% Cyt b; 7.68%) were between the population from Galle and D. scaber from Mihintale. The average pairwise genetic distance between Galle and Yala populations was 2.32% for 16S rRNA, 4.72% for COI, 6.64% for Cyt b (Figures 1C  and 2). The maximum genetic distance between samples from Yala (type locality of D. atukoralei) and other D. scaber samples was 1.55% for 16S rRNA, 3.96% for COI and 4.36% for Cyt b. The average pairwise genetic distance between the Indian D. scaber and the Sri Lankan forms are 1.51% (n=48) for 16S rRNA.

Haplotype Networks
The four populations from Sri Lanka, for the most part, did not share haplotypes. Haplotype sharing only occur for the 16S rRNA gene fragment between the Ampara and Yala populations.
( Figures 1B and 3A -C). Haplotype networks for the COI and Cyt b genes showed an absence of haplotype sharing (sequence data for D. scaber from India, however, were not available for these two genes: Figure 3B,C). All genetic markers also indicated high intraspecific genetic diversity between different localities (Figures 1B and 3). Duttaphrynus scaber had 7, 2 and 3 unique haplotypes, while D. atukoralei had 4, 5 and 6 unique haplotypes for the 16S rRNA, COI and Cyt b mitochondrial markers, respectively ( Figures 1B and 3B,C). In all three markers, D. 'atukoralei' samples from Galle were deeply divergent from their conspecifics from other locations (Figures 1B and 3B,C). The molecular differences between D. scaber and D. 'atukoralei' in 16S rRNA were low (maximum number of changes = 11) compared to intraspecific differences among D. 'atukoralei' from different localities (maximum number of changes = 11) ( Figure 1B).

DISCUSSION
Our molecular phylogeny suggests a single species of Dwarf toad (D. scaber) with three closely related clades: a widespread clade shared between India and Northern dry zone of Sri Lanka, a clade restricted to the Eastern and Southeastern dry zone of Sri Lankan (previously referred to as D. atukoralei; referred to as Ampara-Yala population), and a clade known from the southern wet zone of Sri Lanka (Galle population), which is the most distinct form from all currently recognized taxa. The genertic distance between Duttaphrynus scaber and 'D. atukoralei' is sufficiently low (Figs. 1C, 2) for the latter to be considered a synonym of the former. For the examined populations, the phylogeny and haplotype networks for the 16S rRNA gene suggests incomplete lineage sorting between two populations (Ampara and Yala); COI and Cyt b show complete sorting. The incomplete lineage sorting of 16S rRNA and the close genetic distance suggests that D. atukoralei, for which the type specimen is from Yala (Southeastern dry zone of Sri Lanka) to be synonymous with D. scaber; complete lineage sorting of COI and Cyt b suggests that this species is a distinct population of D. scaber. However, all analyses suggest restricted gene flow to the Southern wet-zone population ( Figure 1A).
The clades from Ampara-Yala are resolved as the sister to the Indian and Northern Sri Lankan D. scaber clade rather than the clade from Galle. Although the locality data for the D. 'atukoralei' from Sri Lanka for Van Bocxlaer et al.'s (2009) analysis was not provided, it is 100% similar to the Galle population.
In our analysis, the uncorrected pairwise genetic distances calculated for Mihintale and Galle were 2.51-2.90% for 16S rRNA, 4.75-5.37% for COI and 6.91-7.68% for Cyt b. These are the highest pairwise genetic distances observed among all populations studied by us. Earlier studies have proposed 3% and 4-8% as threshold values for species discrimination using 16S rRNA and COI, respectively (Fouquet et al., 2007;Smith et al., 2008). Indeed, Fouquet et al. (2007) show that the levels of divergence between lineages, populations and even most sister species in temperate areas for the 16S rRNA gene are well below the 3% threshold limit. For Cyt b, species-level benchmarks have not been suggested so far. Therefore, these benchmarks only come close to the distances between Mihintale and Galle populations. Given that the Galle population is morphologically indistinguishable from the Mihintale and Ampara-Yala populations, at present we recognize the Galle population also as D. scaber.
The pairwise genetic distance calculated for the two clades (Ampara + Yala) and (Mihintale + India) were 0.19-1.55% for 16S rRNA, 0.32-4.08% for COI and 0.21-4.36% for Cyt b. The highest genetic distance within the Mihintale + India D. scaber cluster is 1.16% for 16S rRNA (COI and Cyt b gene sequence data for Indian D. scaber were unavailable in the GenBank; we therefore used only 16S rRNA data); the Mihintale samples are morphologically and genetically consistent with D. scaber.
A rough molecular divergence estimate based on the largest genetic distances of the Cyt b gene between suggests that northern and southern populations of D. scaber diverged about 4.7-7.6 million years ago in the Pliocene. When two rates of change of 0.69% (Martin and Palumbi, 1993) and 1.0% (Tan and Wake, 1995) are applied to the Cyt b molecular distances for the Galle and Mihinthale populations, a range of 4.7-7.6 million years since divergence is given. Though we do not have Cyt b material from an Indian population to compare with the Mihintale population, the patterns that we see from 16S rRNA suggests that their separation is more recent. This suggest recent gene exchange between India and Sri Lanka's northern population. The more basal position of the Galle population in the context of the phylogeny, suggests that Dwarf toads evolved in Sri Lanka and moved across the Sri Lanka's lowland areas to India, however with later migrations (genetic interchanges). For a better understanding of the biogeography and phylogeography of these toads, dense population sampling and gene sampling is needed, especially from India.
The morphology of the Dwarf toads portrays a more complicated situation. The samples we collected from the type locality of D. atukoralei (Yala, southeastern Sri Lanka) appear morphologically similar to D. atukoralei (by the presence of a lobulated parotid gland), though genetically they form a clade with D. scaber, whose type locality is Thiruvananthapuram, Kerala State, India. The population from Galle also show a lobulated parotid gland, though it is genetically distinct from the Yala population.
Despite the small sample sizes, high mitochondrial haplotype diversity was present among all three markers. Interestingly, none of the 16S rRNA haplotypes of D. scaber was shared between the Indian and Sri Lankan populations. All three mitochondrial markers showed similar patterns of genetic structure though COI and Cyt b were similar to each other. Haplotype networks of 16S rRNA, COI and Cyt b corroborate the patterns of the phylogenetic analysis. According to our Cyt b haplotype network, no haplotypesharing was observed between the four populations (Mihintale, Ampara, Yala and Galle). In contrast, haplotype sharing was observed only between two samples of the Ampara population and two of the Yala population. The COI TCS haplotype network is comparable with the Cyt b TCS network. When considering both the COI and the Cyt b gene sequences, the Cyt b sequences seems to be more suitable for constructing haplotype networks because of the greater number of variations.
16S rRNA haplotype networks indicated that there is no haplotype sharing between Mihintale (yellow), Thiruvananthapuram (pink), Maharashtra (brown), Galle (red) and Mudigere (India); (grey) samples. Ampara (blue) and Yala (green) samples share the same haplotype and all the Galle samples within the population share the same haplotype ( Figure 1B). Although Indian and Sri Lankan populations do not share any haplotypes, there are only a few mutational steps among them. Bogert & Senanayake (1966), considering the external morphology, also noted that these toads in peninsular India are larger but identical with those in northern Sri Lanka despite the seawater barrier that prevents gene exchange between toads in the two areas. In both the Cyt b and 16S rRNA haplotype networks, haplotypes from Yala show only a few mutational steps with D. scaber haplotypes, though the Galle haplotypes show significant differences from other D. scaber haplotypes.
Within Sri Lanka, there are no apparent extrinsic barriers to gene flow and D. scaber populations, except for the Galle population which is isolated by a climatic barrier. According to our results, the intra-population genetic diversity (16S rRNA and COI) of D. scaber ('atukoralei') is low for the Galle population, suggesting a small population persisting in relative isolation. Hence our results suggest that the Galle population is an evolutionarily significant unit (ESU) worthy of special conservation attention.

CONCLUSIONS
Our molecular analysis suggests that there is only one Dwarf toad species in Sri Lanka, which is shared with India -Duttaphrynus scaber. However, there is strong population structure within the populations analyzed. The analysis of the 16S rRNA gene fragment shows incomplete lineage sorting, suggestive of a recent divergence event; the other two gene fragments analyzed, COI and Cyt b, show complete lineage sorting. The wet-adapted population from Galle is here considered to be D. scaber, which stands out as a distinct clade that we recognize as an evolutionarily significant unit that does not share any haplotypes with any other population. Deeper analyses using multiple criteria are necessary to ascertain if the Galle population of D. scaber can be considered a distinct species.