Inter- and intra-specific genetic divergence of Asian tiger frogs (genus Hoplobatrachus), with special reference to the population structure of H. tigerinus in Bangladesh

ABSTRACT

The five frog species of the genus Hoplobatrachus are widely distributed in Asia and Africa, with Asia being considered the genus’ origin. However, the evolutionary relationships of Asian Hoplobatrachus species remain ambiguous. Additionally, genetic diversity and fundamental differentiation processes within species have not been studied. We conducted molecular phylogenetic analysis on Asian Hoplobatrachus frogs and population genetic analysis on H. tigerinus in Bangladesh using the mitochondrial CYTB gene and 21 microsatellite markers. The resultant phylogenetic tree revealed monophyly in each species, notwithstanding the involvement of cryptic species in H. chinensis and H. tigerinus, which are evident from the higher genetic divergence between populations. Bayesian inference of population structure revealed genetic divergence between western and eastern H. tigerinus populations in Bangladesh, suggesting restricted gene flow caused by barriers posed by major rivers. However, genetic distances among populations were generally low. A discrete population is located in the low riverine delta region, which likely reflects long-distance dispersal. These results strongly suggest that the environment specific to this river system has maintained the population structure of H. tigerinus in this region.

INTRODUCTION

The tiger frogs, genus Hoplobatrachus (family Dicroglossidae), comprise large robust frogs with either numerous ridges or warts on the back and extensive webbing between the toes. Individuals are semi-aquatic and live mostly near the edge of ponds, marshes, rivers and flooded rice fields. According to Frost (2014), five species are currently recognized: H. crassus, in south to east India, Sri Lanka, Nepal and Bangladesh; H. occipitalis in western and central Africa; H. rugulosus [= H. chinensis, used by some authors (e.g., Kosuch et al., 2001) as its nomenclature status is unclarified] in Myanmar, southern China, Taiwan, Thailand and peninsular Malaysia; H. tigerinus in east Afghanistan, north Pakistan, India, Sri Lanka, Nepal, Bangladesh and Myanmar; and the recently discovered H. litoralis in the southeastern coastal belt of Bangladesh (Hasan et al., 2012a).

For the evolutionary history of Hoplobatrachus species, Kosuch et al. (2001) suggested a common Hoplobatrachus ancestor of southern Asian origin and dispersal of the H. occipitalis lineage into Africa. However, the phylogenetic relationship among Asian species remains ambiguous. Alam et al. (2008) reported the phylogenetic relationships of three genera distributed in South Asia, including Hoplobatrachus, using partial mitochondrial genes. More recently, Alam et al. (2012) clarified the degree of postmating isolation in six species. Although they outlined genetic divergence among the species, their taxon sampling was not comprehensive for Asian Hoplobatrachus and the evolutionary relationships remain uncertain.

Additionally, within-species genetic diversity has not been fully studied in Hoplobatrachus. As mentioned above, this genus is thought to have originated in South and Southeast Asia, which are recognized as biodiversity hot spots for many organisms including amphibians (Myers et al., 2000). These areas are generally located in developing countries and, thus, potential frog species habitats are decreasing as a result of rapid development before they can be extensively investigated. Three Hoplobatrachus species are distributed in Bangladesh (H. tigerinus, H. litoralis and H. crassus), and their habitats are decreasing because of anthropological activity and rapid environmental changes, such as development of infrastructure, and water pollution by pesticides and agrochemicals (Islam et al., 2012). In terms of species conservation, the existence of intraspecific genetic diversity and hierarchical population structure is crucial for population persistence and evolutionary potential. Moreover, finding environmental factors that affect population structure is also important to understand the forces that drive speciation in biodiversity hot spots. Therefore, population genetic studies are needed before Hoplobatrachus populations decline too severely in Bangladesh. Of the three species found in Bangladesh, the distribution of H. litoralis and H. crassus is restricted to the southeastern corner and southern coastal belt, respectively; on the other hand, H. tigerinus is found all over Bangladesh (Alam et al., 2008; Hasan et al., 2012b), making it an ideal model to investigate population structure and the diversification process.

In the present study, we aimed to elucidate the evolutionary relationships of Asian Hoplobatrachus and the fine-scale population structure and genetic diversity of H. tigerinus in Bangladesh. To do this, we conducted molecular phylogenetic and population genetic analysis using complete nucleotide sequences of the mitochondrial CYTB gene and genotypic data from 21 microsatellite markers. We discussed the inter- and intraspecies genetic diversity of Asian Hoplobatrachus species and the population structure of H. tigerinus related to causal environmental factors.

MATERIALS AND METHODS

Sample collection

A total of 343 individuals comprising four species from the genus Hoplobatrachus (H. tigerinus, H. litoralis, H. chinensis and H. crassus) were used in the present study (Table 1 and Fig. 1A). Of these, 329 H. tigerinus individuals were collected from 30 localities in Bangladesh and four in India; four H. litoralis individuals were collected from three localities in Cox’s Bazar, Bangladesh; seven H. chinensis individuals were collected from three localities in Laos and Huu Lien in Vietnam; and three H. crassus individuals were collected from two localities in Bangladesh and Assam, India. Sampling was conducted from October 2001 to November 2012. The frogs were mainly captured from paddy fields near water bodies at each site. Their toes were cut with scissors and stored in 100% ethanol. These samples were kept in a –30 ℃ freezer until needed.

Table 1. Specimens used in this study

Species name	Country	Locality	Abbreviation	Number of individuals used in this study	Number of individuals used for mt analyses	Number of individuals used for SSR analyses
H. tigerinus	Bangladesh	Lalmonirhat	Lal	9	2	9
		Dinajpur	Din	20	2	20
		Saidpur, Nilphamari	Syd	21	2	21
		Rangpur	Rang	16	2	16
		Sonamasjid, Chapai Nawabganj	CNSM	2	2	2
		Chapai Nawabganj	CN	7	2	7
		Pabna	Pabna	6	2	6
		Iswardi, Pabna	IS	6	2	6
		Jagannathganj, Jamalpur	Jagan	14	2	14
		Jamalpur	Jamal	16	2	16
		Kishorganj	Kish	1	1	1
		Netrokona	Net	12	2	12
		Sylhet	Syl	25	2	25
		Dirai, Sunamganj	C	20	2	20
		BAU Campus, Mymensingh	BAU	30	3	30
		Fulbaria, Mymemsingh	Fulb	11	0	11
		Phulpur, Mymensingh	Fulp	12	2	12
		Narsingdi	NG	7	2	7
		Khulna	Khul	15	2	15
		Satkhira	SK	11	2	11
		Barguna	Br	15	1	15
		Noakhali	No	10	2	10
		PSTU Campus, Patuakhali	Pkun	5	2	5
		Bhola	Vo	21	2	21
		Sandwip, Chittagong	Snd	1	1	1
		Raojan, Chittagong	Raoj	1	1	1
		Bandarban	Bn	1	1	1
		Ranirhat, Chittagong	Chr	2	2	–
		Churkhai, Mymensingh	Chur	1	1	–
		Shamvuganj, Mymemsingh	Shamb	2	2	–
	India	Bajpe	IBM	3	3	–
		Padil	IPM	2	2	–
		Karnoor	Ikrn	2	2	–
		Shirva	Ishr	2	2	–
H. litoralis	Bangladesh	Ukhia, Cox’s Bazar	lit_Uk	1	1	–
		Teknaf, Cox’s Bazar	lit_Tk	1	1	–
		Cox’s Bazar	lit_Cx	2	2	–
H. chinensis	Laos	Luang Prabang	chi_LaLP	1	1	–
		Vang Vieng	chi_LaVV	2	2	–
		Long Nai Phongsaly	chi_LLNPP	2	2	–
	Vietnam	Huu Lien	chi_VtHL	2	2	–
	Thailand	Chachoengsao*	chi_Chach	3	3	–
		Nong Khai*	chi_Nongk	1	1	–
H. crassus	Bangladesh	Dacope, Khulna	cra_Dac	2	2	–
	India	Assam	cra_IA	1	1	–
Euphlyctis hexadactylus	Bangladesh	Sundarban**		1	1	–
Total				348	81	315

*  data from Alam et al. (2012).

**  data from Alam et al. (2010).

Fig. 1.

Sampling locations, mitochondrial CYTB haplotype network and distribution of H. tigerinus in Bangladesh. A) Map of sampling locations and haplotype distribution of H. tigerinus CYTB in Bangladesh. Size of pie charts represents sample size of each population. Colors of pie chart elements correspond to those displayed in the haplotype network. For detailed distribution of haplotypes, see Supplementary Table S1. B) Haplotype network of mitochondrial CYTB based on the parsimony method. Circle sizes represent haplotype frequencies. Numbers in circles indicate haplotype numbers. Each dot on the lines connecting the haplotypes represents a substitution between the haplotypes.

DNA extraction and CYTB and microsatellite loci genotyping

Total genomic DNA was extracted from the ethanol-preserved clipped toes using a DNeasy Blood and Tissue Kit (Qiagen, CA, USA) according to the manufacturer’s instructions. Extracted DNA was used to genotype the mitochondrial gene and microsatellite loci.

To confirm evolutionary relationships among species and the genetic diversity of H. tigerinus populations in Bangladesh, we defined mitochondrial CYTB gene haplotypes and the genotypes of 21 microsatellite loci. For the CYTB gene, we defined the haplotypes of 76 samples representing every population of the five Hoplobatrachus species. PCR amplification was performed using an EmeraldAmp PCR Master Mix (Takara Bio, Otsu, Japan) using one pair of primers, HoploGluF-18 (forward primer 5′-TAACTCGGACCTGTAGTCTGA-3′) and HoploDloopR-46 (reverse primer 5′-GAGTGTACTTGAAGATATGCTTG-3′), for the complete CYTB sequence. PCR was performed in 10-μl reactions following the manufacturer’s protocol. The amplified fragments were purified by polyethylene glycol precipitation and subjected to DNA sequencing. All fragments were sequenced in both directions using the same primers as for PCR amplification. DNA sequencing was performed using BigDye Terminator ver. 3.1 (Thermo Fisher Scientific, DE, USA) with an ABI 3130xl Genetic Analyzer (Thermo Fisher Scientific) according to the manufacturer’s instructions.

To genotype the microsatellite loci, we first genotyped 27 loci previously reported for 315 H. tigerinus individuals, following a published protocol (Sultana et al., 2014). Genotypic data were first checked for the plausible occurrence of null alleles, large allele dropout, and errors in scoring caused by stutter bands or misidentified alleles using Microchecker v.2.2.3 (Van Oosterhout et al., 2004). Tests for Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) were performed for each of the residual 21 microsatellite loci using GENEPOP version 4.2.1 (Rousset, 2008).

Inference of phylogenetic relationships and genetic diversity

The resultant CYTB sequences were aligned using ClustalW (Thompson et al., 1994). We prepared two datasets: (1) a dataset containing all haplotypes from five Hoplobatrachus species in South and Southeast Asia and (2) a dataset containing H. tigerinus haplotypes in Bangladesh. In dataset (1), we also included H. chinensis and Euphlyctis hexadactylus CYTB haplotypes previously reported by Alam et al. (2010, 2012) (accession numbers: AB636597–AB636600 and AP011544). Based on this dataset, we constructed a Hoplobatrachus species phylogenetic tree using the neighbor-joining (NJ) (Saitou and Nei, 1987) and maximum likelihood (ML) methods. An NJ tree was constructed based on the p-distance (Nei and Kumar, 2000) using MEGA 6 (Tamura et al., 2013). We constructed two ML trees by applying different substitution models: either a single model for all codon sites (non-partitioned) or specific models for each codon position (codon-partitioned). The best-fit models were selected by KAKUSAN4 (Tanabe, 2007) and the ML trees were constructed using RAxML ver. 8.1.24 (Stamatakis, 2014). In every analysis, bootstrap values were calculated with 1000 replicates and the resultant trees were rooted with E. hexadactylus as an outgroup. Based on dataset (2), we constructed a minimum-spanning network using TCS 1.21 (Clement et al., 2000), which estimates the parsimony according to Templeton et al. (1992). The maximum number of steps to connect haplotypes parsimoniously was calculated with a 99% limit.

For the H. tigerinus microsatellite genotype data, we examined genetic diversity by calculating the number of alleles (A), the number of effective alleles (A_E), the observed (H_O) and expected heterozygosity (H_E), and the fixation index (F) for each population using GenAlEx 6.5 (Peakall and Smouse, 2012).

For the subsequent population-based analysis, using 20–30 individuals per population is recommended for population genetic studies on wild organisms based on microsatellite data (Pruett and Winker, 2008) because small sample sizes can skew results. However, we could not obtain such numbers in all of the populations because of sampling difficulty. Therefore, we excluded populations with less than 10 individuals (leaving 16 populations) from the analysis. We also essentially confirmed the same results based on a dataset that excluded populations with less than 20 individuals.

To infer the fine-scale population structure of H. tigerinus, we first calculated pairwise F_ST and Nei’s D_A distance (Nei et al., 1983) among the populations. We then constructed multidimensional scaling based on the pairwise F_ST values to visualize genetic differences using the ‘cmdscale’ function in R (R Core Team, 2015). The Bayesian clustering approach was performed in STRUCTURE v.2.3.3 (Pritchard et al., 2000) to assess the level of population genetic structure without allocating individuals to populations and to assign ancestral elements to each individual. STRUCTURE estimates ancestries of individuals under a given number of clusters (K) and identifies migrants and admixed individuals based solely on genetic data. We inferred the most reasonable K (number of clusters or populations) using all of the microsatellite loci to reveal the population genetic structure. We ran 1,000,000 MCMC repetitions after discarding the first 100,000 iterations as burn-in data. Ten simulations were completed for each estimated K (1–8). Analysis was performed using the admixture model with the LOCPRIOR option (Hubisz et al., 2009) for MCMC inference. To estimate a realistic K value, we analyzed our resultant data as described by Evanno et al. (2005) in Structure Harvester v.0.6.93 (Earl and vonHoldt, 2012), in which log-likelihood values and their variances were used to calculate ΔK, an objective alternative to simply choosing a value with the highest log probability (Evanno et al., 2005). CLUMPP v.1.1.2 (Jakobsson and Rosenberg, 2007) was then used to permute clusters across the 10 runs of the selected K. The results of these analyses were used to summarize runs and visualize results using DISTRUCT v.1.1 (Rosenberg, 2004).

Detection of IBD and decomposed pairwise regression analysis

When populations are maintained under gene drift-flow equilibrium, population structure is expected to correlate with geographic distance, which is known as isolation by distance (IBD) (Wright, 1943). The most common pattern shows greater genetic divergence among geographically distant populations as a result of higher effect of genetic drift, and greater genetic affinity among geographically close populations as a result of more gene flow. However, environmental factors may also affect migration between populations; these are typically geographic barriers, which prevent migration between populations.

To confirm IBD patterns and the existence of geographic barriers, we conducted Mantel and partial Mantel tests (Mantel, 1967) using the ‘mantel’ and ‘mantel.partial’ functions in the ‘vegan’ package (Oksanen, 2011) in R with 9,999 randomizations. We used a linearized F_ST [F_ST /(1 – F_ST )] (Rousset, 1997) matrix as a dependent variable and Euclidian distance and putative barrier matrices as predictor variables. Taking the STRUCTURE results, the putative barrier matrix was constructed based on the two major rivers that flow in a north-south direction in this region (the Jamuna and Meghna Rivers) (1 = crossing the rivers, 0 = not).

Moreover, specific information on individual populations was obscured on all of the pairwise plots of genetic distance against geographic distance. We therefore applied decomposed pairwise regression (DPR) analysis (Koizumi et al., 2006) to determine the characteristics of individual H. tigerinus populations in Bangladesh. In brief, this analysis first detected a putative outlier population, which was not in accordance with the regression lines from the pairwise plot of linearized F_ST and geographic distances (for details, see Koizumi et al., 2006). In the original method, true outlier populations are then identified by choosing the best model based on the corrected Akaike information criterion (AIC). However, using AIC in DPR analysis is technically invalid (Reynolds and Fitzpatrick, 2013; I. Koizumi, personal communication), so we used r² values to identify the true outlier.

RESULTS

Mitochondrial haplotypes and phylogenetic relationships

We sequenced 1,146 bp of the mtDNA CYTB region, and 41, 4, 7 and 2 haplotypes were found in H. tigerinus, H. litoralis, H. chinensis and H. crassus (54 haplotypes in total), respectivley. The nucleotide sequences of these haplotypes were deposited in DDBJ (accession numbers LC120953–LC121023). For the ML based on dataset (1), GTR + G was selected as the best fit model for both non-partitioned and every codon site of codon-partitioned models. The topologies of the resultant NJ and ML trees (–InL = 5,612.10 and 5,231.80 in non-partitioned and codon-partitioned models, respectively) were almost identical in terms of the relationships between major clades, and thus the NJ tree is shown in Fig. 2 with bootstrap values. The resultant tree based on dataset (1) exhibited clear monophyletic clades for each species with high bootstrap values (> 95%), except for H. chinensis (Fig. 2). Within the ingroup (E. hexadactylus was an outgroup), although the bootstrap support values were low, H. crassus was the first Asian Hoplobatrachus species to diverge. Hoplobatrachus chinensis then diverged from a common ancestor of H. tigerinus and H. litoralis, which is in agreement with a previous study (Alam et al., 2008). Within the H. tigerinus clade, the Bangladeshi and Indian populations diverged with a mean p-distance of 0.093.

Fig. 2.

Neighbor-joining tree of Asian Hoplobatrachus frog mitochondrial CYTB haplotypes based on p-distance. The numbers at each node are bootstrap values (≥ 50%) in order for NJ, ML (non-partitioned) and ML (codon-partitioned). Bootstrap values below 50% are indicated as ‘–’.

We found 33 haplotypes in 76 H. tigerinus individuals (dataset (2), see Supplementary Table S1 for details). Genetic divergence within Bangladesh was very small, showing a mean p-distance of 0.002 and a total nucleotide diversity of 0.00456. A H. tigerinus haplotype network tree formed two major groups (haplogroups) having their centers at H6 and H20 (Fig. 1B). H6 was the most frequent haplotype (9 of 76) and was common to eight of the populations.

Microsatellite genotypes and H. tigerinus population structure in Bangladesh

Of the 27 loci examined in this study, we could not repeat the scoring of six loci (Htigr2250, Htigr1262, Htigr1163, Htigr156, Htigr76 and Htigr1284) in more than one population. These loci were excluded from subsequent analysis. For the remaining 21 loci, the A ranged from 4 to 39 with a mean of 17.1 per locus (Supplementary Tables S2 and S3). After calculating the Bonferroni correction, significant LD was detected in Jamal (Htigr961/Htigr823) and SK (Htigr243/Htigr735) (adjusted P < 0.00001035), and significant deviation from HWE was detected in Lal (Htigr409), Jagan (Htigr828), Jamal (Htigr1984), BAU (Htigr409), Fulp (Htigr828) and Vo (Htigr409) (adjusted P < 0.00010352). Because none of these loci exhibited significant LD in the same combinations in other populations, we did not exclude them from the analysis. The A_E, H_E and F of all loci for each population varied from 1.571 to 5.353 (mean 3.842), from 0.286 to 0.718 (mean 0.615) and from −0.175 to 0.151 (excluding populations consisting of a single individual) (mean 0.028) in each population, respectively (Table 2). The overall H. tigerinus F_ST value was 0.030 and the pairwise F_ST and Nei’s D_A among H. tigerinus populations (excluding populations with a single individual) are summarized in Supplementary Table S4. Multidimensional scaling shows the presence of two regional clusters (western and eastern) and two outlier populations (Vo and No) (Fig. 3). For the STRUCTURE analysis, the mean log-likelihood value was highest in K = 3 (Supplementary Fig. S1A), while ΔK was highest in K = 2 (Supplementary Fig. S1B). At K = 2, two ancestral elements were split between the east and west in the northern region (Fig. 4).

Table 2. Genotypic data of 21 microsatellite loci for H. tigerinus populations in Bangladesh

Population	N	A	AE	HO	HE	F
Lal	9	5.8	4.113	0.636	0.645	–0.002
Din	20	8.0	4.705	0.605	0.657	0.075
Syd	21	7.7	4.425	0.594	0.644	0.093
Rang	16	7.9	4.760	0.661	0.689	0.051
CNSM	2	2.4	2.289	0.524	0.440	–0.175
CN	7	5.2	3.763	0.612	0.640	0.044
Pabna	6	5.0	3.724	0.698	0.648	–0.087
IS	6	4.9	3.683	0.654	0.638	–0.035
Jagan	14	6.9	4.318	0.585	0.689	0.151
Jamal	16	8.4	5.353	0.653	0.718	0.091
Kish	1	1.7	1.667	0.667	0.333	–
Net	12	6.6	4.225	0.674	0.670	0.030
Syl	25	7.2	4.630	0.696	0.711	0.019
C	20	7.8	4.629	0.668	0.695	0.041
BAU	30	8.8	4.562	0.650	0.709	0.076
Fulb	11	6.0	4.198	0.672	0.671	0.015
Fulp	12	6.8	4.092	0.615	0.667	0.094
NG	7	5.2	4.081	0.652	0.662	0.040
Khul	15	7.7	4.838	0.669	0.717	0.097
SK	11	6.3	4.106	0.630	0.665	0.058
Br	15	6.7	4.289	0.673	0.674	–0.005
No	10	5.8	4.227	0.671	0.665	–0.010
Pkun	5	4.5	3.466	0.669	0.611	–0.074
Vo	21	7.6	4.541	0.662	0.713	0.064
Snd	1	1.6	1.571	0.571	0.286	–
Raoj	1	1.8	1.810	0.810	0.405	–
Bn	1	1.7	1.667	0.667	0.333	–
Mean	11.667	5.8	3.842	0.650	0.615	0.028

Abbreviations: N, number of individuals; A, number of alleles; A_E, number of effective alleles; H_O, observed heterozygosity; H_E, expected heterozygosity; F, fixation index.

Fig. 3.

Multidimensional scaling plots of the populations based on linearized genetic distances among H. tigerinus populations in Bangladesh. The populations were plotted using a two-dimensional scale, and, except for the No and Vo populations, those located in western (Din, Syd, Rang, Khul and SK) and eastern (Jagan, Jamal, Net, Syl, C, BAU, Fulb, Fulp) areas were clustered in the lower right and upper right, respectively.

Fig. 4.

Assigned ancestral elements from STRUCTURE for H. tigerinus populations in Bangladesh. Individuals were assigned two ancestral elements using STRUCTURE 2.3.3 (Hubisz et al., 2009) at K = 2 under the admixture model with the LOCPRIOR option. Pie charts show the total ratios of assigned ancestral elements of individuals in each population. Sizes of pie charts represent sample sizes in each population.

Detection of IBD, barrier effects and the outlier population

A moderate IBD pattern was observed within the H. tigerinus population because a plot of the relationship between pairwise linearized F_ST and geographic distance showed a gradual positive correlation (Fig. 5). Mantel tests detected significant correlations between the geographic distance (r² = 0.147, P = 0.0017) and the putative barrier of the two major rivers (r² = 0.142, P = 0.001). Partial Mantel tests also showed significant correlations between geographic distance when controlling the putative barrier (r² = 0.147, P = 0.005) and the putative barrier when controlling geographic distance (r² = 0.082, P = 0.0099).

Fig. 5.

Relationship between genetic and geographic distances for all 16 populations that consisted of ten or more individuals (open and filled circles; y = 0.00006x + 0.0545, r² = 0.147, thick line) and excluding five putative outlier populations (filled circles only; y = 0.00006x + 0.0571, r² = 0.177, thin line).

Five populations, Vo, BAU, Rang, SK and Fulb, were detected as putative outlier populations by systematic DPR analysis (Fig. 6), and Vo was identified as a true outlier because the highest r² value (0.185) was observed when this population was excluded (Table 3). The DPR of Vo gave lower slope and intersect values, which means lower genetic differentiation even among distant populations (Fig. 7), and is thus recognized as Pattern 4, i.e., gene flow > genetic drift (Koizumi et al., 2006).

Fig. 6.

Means and 95% confidence intervals of residuals from the least-squares regression (calculated from Fig. 5, all populations included). The populations in bold type (Vo, BAU, Rang, SK and Fulb) were detected as putative outliers.

Table 3. Fit of alternative model with and without putative outlier populations

Population excluded	Intercept	Slope*102	r2	P value (Mantel test)
Vo	0.055	0.007	0.185	0.001
Vo, BAU, Rang, SK, Fulb	0.057	0.006	0.177	0.006
Vo, BAU, Rang	0.061	0.006	0.155	0.002
none	0.055	0.006	0.147	0.005
Vo, BAU	0.060	0.005	0.123	0.007
Vo, BAU, Rang, SK	0.061	0.005	0.132	0.013

Fig. 7.

Decomposed pairwise regression of genetic and geographic distances for each of the 16 populations. Each of the five putative outlier populations (labeled; colored lines) was regressed with the 11 non-outlier populations (black lines), while each of the 11 non-outliers was regressed with the other ten populations. Solid and dashed lines represent statistically significant and nonsignificant regressions, respectively.

DISCUSSION

Inter- and intraspecific diversity in Asian Hoplobatrachus species

Our phylogenetic analysis based on the CYTB gene indicated a sister relationship between H. tigerinus and H. litoralis, in agreement with previous studies (Hasan et al., 2012a, 2012b). However, the relationship at deeper nodes was unclear. Alam et al. (2008) also reported similar results based on shorter partial fragments of three genes (CYTB and 12S and 16S rRNAs). In this study, we used a complete CYTB gene sequence and more samples, but our results eventually failed to uncover the evolutionary relationships. Specifically, H. crassus was located at the most basal position of Asian Hoplobatrachus and this was supported by a moderate bootstrap value in the ML tree (codon-partitioned), but the generally lower bootstrap values imply alternative plausible relationships. The uncertainty of the relationships may be because of a lack of nucleotide substitution accumulation, which represents phylogenetic information. The stem lengths of common ancestral Asian Hoplobatrachus lineages were comparatively short, which suggests rapid speciation in these lineages. To clarify the true relationships, further studies using more molecular markers for multiple nuclear genes are needed.

Our phylogenetic tree also clearly indicates the presence of a cryptic species (putative new taxon) in the H. chinensis specimens examined. Additionally, H. tigerinus in India and Bangladesh may be different species because the haplotypes in these populations accumulated a genetic difference between them comparable to that of H. litoralis. Further phylogeographic studies using multiple nuclear markers are needed to elucidate the complete speciation process and species status of Asian Hoplobatrachus.

Population structure of H. tigerinus in Bangladesh and habitat effects

Although genetic differences among H. tigerinus populations in Bangladesh were generally low based on the mitochondrial haplotype data, the distribution of haplogroups exhibited a weak division between southeastern and northwestern regions (Fig. 1), which suggests the existence of isolated populations in the past that resulted in these two major haplogroups. Although the genetic distance was still not high (F_ST = 0.030) in the microsatellite data, the population structure became apparent in our analysis. As shown in Fig. 5, although the slope of the regression line was shallow, genetic distances were correlated with geographic distances and thus the genetic diversity of H. tigerinus in Bangladesh represents a weak IBD pattern. The topography of Bangladesh, i.e., low flatlands in a region with enormous wetland areas (Bhuiyan and Baky, 2014), and the semi-aquatic habitat of H. tigerinus, including marshland and agricultural paddy, should maintain constant gene flow among populations and explain the weak IBD pattern. While the multidimensional scaling and STRUCTURE results revealed longitudinal divergence between eastern and western regions delineated by the Jamuna and Meghna Rivers (Fig. 4), the Mantel and partial Mantel tests supported a negative effect of these rivers on transversal gene flow across them. Vast quantities of water flow into the Jamuna and Meghna Rivers, which are 8–10 km wide (Bhuiyan and Baky, 2014). Even though H. tigerinus individuals can live in the water body, they cannot traverse these rivers easily by swimming. Although further landscape genetic analyses using other environmental factors are needed before firm conclusions can be drawn, the rivers are therefore liklely to restrict migration and to have a negative effect on gene flow as exemplified in other frogs (Angelone et al., 2011; Fouquet et al., 2012). However, because the absolute value of genetic distances was not very high (F_ST = 0.03), the magnitude of this barrier effect is probably either low or weakened by occasional migration and/or artificial transfer. Alternatively, H. tigerinus in Bangladesh may have experienced a recent rapid population expansion. In the southern region, the two ancestral elements inferred from STRUCTURE were intermixed (Fig. 4). In Bangladesh, the boundaries of the rivers in the southern delta region change seasonally and floods occur (Mirza et al., 2003). The bifurcated population structure observed in the northern region may also occur in the southern area in association with disturbance of the landscape from seasonal changes in the river system. However, ecological and ethological studies, including direct observation, are needed to verify this.

Moreover, in contrast to the negative effect of the rivers on migration and gene flow across them, the water current itself could enhance directional migration from upper to lower regions (north to south) by transporting individuals. The Vo population, which we identified as an outlier population in the DPR analysis, may reflect such migration. As mentioned in the original DPR analysis (Koizumi et al., 2006), a population showing a stronger gene flow effect may be the result of long-distance dispersal by water currents from upper regions. This is also applicable to the Vo population because it is located in the middle of the delta area neighboring the Meghna River.

As for intrapopulation genetic diversity regarding H. tigerinus conservation, every population retained the same level of diversity and no population bottlenecks were observed. This is also explained by low genetic difference and higher gene flow among populations. Nonetheless, Bangladesh is a developing country and some populations close to urban areas may be suffering from anthropological effects and rapid environmental changes. Therefore, continuous monitoring is needed and our present data will serve as a foundation for the sustainable conservation of H. tigerinus in Bangladesh.

CONCLUSION

We have elucidated the evolutionary relationships of Asian Hoplobatrachus species and the fine-scale population structure of H. tigerinus in Bangladesh. Our phylogenetic tree based on the mitochondrial CYTB gene suggested the presence of a cryptic species within H. chinensis. Within H. tigerinus, our microsatellite data revealed genetic divergence between western and eastern areas separated by two major rivers, which suggests that the rivers can restrict gene flow. The changing river system and annual flooding should disturb the population structure and long-distance dispersal of populations in the delta region. We believe that our results will also be valuable for the sustainable conservation of H. tigerinus in Bangladesh.

ACKNOWLEDGMENTS

We thank Dr. A. Dubois, Muséum National d’Histoire Naturelle, Paris, France, for supplying the samples from Laos, Vietnam and Thailand, and Dr. S. H. Joshy, St. Aloysius College, Mangalore, India and Dr. M. Kuramoto, Professor Emeritus, Fukuoka University of Teacher Education, Fukuoka, Japan, for supplying the samples from India. We thank Dr. Itsuro Koizumi, Hokkaido University, Sapporo, Japan, for providing us with the DPR analysis program and valuable comments. We also thank two anonymous reviewers for their critical review of early versions of the manuscript and the Wildlife Management and Nature Conservation Division, Bangladesh for permission to collect and transport frogs from Bangladesh to Japan. This research was supported by a Grant-in-Aid for Scientific Research (B) (No. 24310173) to M. S. and a Grant-in-Aid for Young Scientists (B) (No. 23710282) to T. I. from the Japan Society for the Promotion of Science.

REFERENCES

•  Alam,  M. S.,  Igawa,  T.,  Khan,  M. M. R.,  Islam,  M. M.,  Kuramoto,  M.,  Matsui,  M.,  Kurabayashi,  A., and  Sumida,  M. (2008) Genetic divergence and evolutionary relationships in six species of genera Hoplobatrachus and Euphlyctis (Amphibia: Anura) from Bangladesh and other Asian countries revealed by mitochondrial gene sequences. Mol. Phylogenet. Evol. 48, 515–527.
•  Alam,  M. S.,  Kurabayashi,  A.,  Hayashi,  Y.,  Sano,  N.,  Khan,  M. M. R.,  Fujii,  T., and  Sumida,  M. (2010) Complete mitochondrial genomes and novel gene rearrangements in two dicroglossid frogs, Hoplobatrachus tigerinus and Euphlyctis hexadactylus, from Bangladesh. Genes Genet. Syst. 85, 219–232.
•  Alam,  M. S.,  Islam,  M. M.,  Khan,  M. R.,  Hasan,  M.,  Wanichanon,  R., and  Sumida,  M. (2012) Postmating isolation in six species of three genera (Hoplobatrachus, Euphlyctis and Fejervarya) from family Dicroglossidae (anura), with special reference to spontaneous production of allotriploids. Zool. Sci. 29, 743–752.
•  Angelone,  S.,  Kienast,  F., and  Holderegger,  R. (2011) Where movement happens: scale-dependent landscape effects on genetic differentiation in the European tree frog. Ecography (Cop.) 34, 714–722.
•  Bhuiyan,  S. R., and  Baky,  A. Al (2014) Digital elevation based flood hazard and vulnerability study at various return periods in Sirajganj Sadar Upazila, Bangladesh. Int. J. Disaster Risk Reduct. 10, 48–58.
•  Clement,  M.,  Posada,  D., and  Crandall,  K. A. (2000) TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9, 1657–1659.
•  Earl,  D. A., and  vonHoldt,  B. M. (2012) STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361.
•  Evanno,  G.,  Regnaut,  S., and  Goudet,  J. (2005) Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620.
•  Fouquet,  A.,  Ledoux,  J. B.,  Dubut,  V.,  Noonan,  B. P., and  Scotti,  I. (2012) The interplay of dispersal limitation, rivers, and historical events shapes the genetic structure of an Amazonian frog. Biol. J. Linn. Soc. 106, 356–373.
•  Frost,  D. R. (2014) Amphibian Species of the World: an online reference. Version 6.
•  Hasan,  M.,  Islam,  M. M.,  Khan,  M. M. R.,  Alam,  M. S.,  Kurabayashi,  A.,  Igawa,  T.,  Kuramoto,  M., and  Sumida,  M. (2012a) Cryptic Anuran Biodiversity in Bangladesh Revealed by Mitochondrial 16S rRNA Gene Sequences. Zool. Sci. 29, 162–172.
•  Hasan,  M.,  Kuramnoto,  M.,  Islam,  M.,  Alam,  H.,  Khan,  M., and  Sumida,  M. (2012b) A new species of genus Hoplobatrachus (Anura, Dicroglossidae) from the coastal belt of Bangladesh. Zootaxa 3321, 45–58.
•  Hubisz,  M. J.,  Falush,  D.,  Stephens,  M., and  Pritchard,  J. K. (2009) Inferring weak population structure with the assistance of sample group information. Mol. Ecol. Resour. 9, 1322–1332.
•  Islam,  M. N.,  Shaikat,  A. H.,  Khan,  S. A.,  Hassan,  M. M., and  Kriger,  K. M. (2012) Concern over destroying frog habitat on the occasion of save the frogs day in Bangladesh. FrogLog 20, 28.
•  Jakobsson,  M., and  Rosenberg,  N. A. (2007) CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics 23, 1801–1806.
•  Koizumi,  I.,  Yamamoto,  S., and  Maekawa,  K. (2006) Decomposed pairwise regression analysis of genetic and geographic distances reveals a metapopulation structure of stream-dwelling Dolly Varden charr. Mol. Ecol. 15, 3175–3189.
•  Kosuch,  J.,  Vences,  M.,  Dubois,  A.,  Ohler,  A., and  Böhme,  W. (2001) Out of Asia: mitochondrial DNA evidence for an Oriental origin of tiger frogs, genus Hoplobatrachus. Mol. Phylogenet. Evol. 21, 398–407.
•  Mantel,  N. (1967) The detection of disease clustering and a generalized regression approach. Cancer Res. 27, 209–220.
•  Mirza,  M. M. Q.,  Warrick,  R. A., and  Ericksen,  N. J. (2003) The implications of climate change on floods of the Ganges, Brahmaputra and Meghna rivers in Bangladesh. Clim. Change 57, 287–318.
•  Myers,  N.,  Mittermeier,  C. G.,  de Fonseca,  G. A. B., and  Kent,  J. (2000) Biodiversity hotspots for conservation priorities. Nature 403, 853–858.
•  Nei,  M., and  Kumar,  S. (2000) Molecular Evolution and Phylogenetics Oxford Univ.
•  Nei,  M.,  Tajima,  F., and  Tateno,  Y. (1983) Accuracy of estimated phylogenetic trees from molecular data. J. Mol. Evol. 19, 153–170.
•  Oksanen,  J. (2011) Multivariate analysis of ecological communities in R: vegan tutorial.
•  Van Oosterhout,  C.,  Hutchinson,  W. F.,  Wills,  D. P. M., and  Shipley,  P. (2004) Micro-Checker: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538.
•  Peakall,  R., and  Smouse,  P. E. (2012) GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics 28, 2537–2539.
•  Pritchard,  J. K.,  Stephens,  M., and  Donnelly,  P. (2000) Inference of population structure using multilocus genotype data. Genetics 155, 945–959.
• R Core Team (2015) R: A language and environment for statistical computing. R Foundation for Statistical Computing.
•  Reynolds,  R. G., and  Fitzpatrick,  B. M. (2013) Tests of two methods for identifying founder effects in metapopulations reveal substantial type II error. Genetica 141, 119–131.
•  Rosenberg,  N. A. (2003) Distruct: a Program for the Graphical Display of Population Structure. Mol. Ecol. Notes 4, 137–138.
•  Rousset,  F. (1997) Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145, 1219.
•  Saitou,  N., and  Nei,  M. (1987) The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4, 406–425.
•  Stamatakis,  A. (2014) RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313.
•  Sultana,  N.,  Igawa,  T.,  Nozawa,  M.,  Islam,  M. M.,  Hasan,  M.,  Alam,  M. S.,  Khan,  M. M. R., and  Sumida,  M. (2014) Development and characterization of 27 new microsatellite markers of Indian Bullfrog, Hoplobatrachus tigerinus and its congeneric species. Genes Genet. Syst. 89, 137–141.
•  Tamura,  K.,  Stecher,  G.,  Peterson,  D.,  Filipski,  A., and  Kumar,  S. (2013) MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol. Biol. Evol. 30, 2725–2729.
•  Tanabe,  A. S. (2007) KAKUSAN: a computer program to automate the selection of a nucleotide substitution model and the configuration of a mixed model on multilocus data. Mol. Ecol. Notes 7, 962–964.
•  Templeton,  A. R.,  Crandall,  K. A., and  Sing,  C. F. (1992) Cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132, 619–633.
•  Wright,  S. (1943) Isolation by distance. Genetics 28, 114–318.