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

Nasrin Sultana‡, Takeshi Igawa*‡, Mohammed Mafizul Islam, Mahmudul Hasan, Mohammad Shafiqul Alam, Shohei Komaki, Kensuke Kawamura, Md. Mukhlesur Rahman Khan† and Masayuki Sumida Institute for amphibian Biology, Graduate School of Science, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan Graduate School for International Development and Cooperation, Hiroshima University, 1-5-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8529, Japan Department of Genetics and Fish Breeding, Faculty of Fisheries, Bangabandhu Sheikh Mujibur Rahman Agricultural University, Gazipur 1706, Bangladesh Department of Fisheries Biology and Genetics, Faculty of Fisheries, Bangladesh Agricultural University, Mymensingh 2222, Bangladesh


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.  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 °C freezer until needed.

Sample collection
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. (2010Alam et al. ( , 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 (nonpartitioned) 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-  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.Genetic divergence of Asian tiger frogs 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 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 (codonpartitioned).Bootstrap values below 50% are indicated as '-'.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 2 values to identify the true outlier.

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 codonpartitioned models.The topologies of the resultant NJ and ML trees (-InL = 5,612.10and 5,231.80 in nonpartitioned 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.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.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.
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).

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 2 = 0.147, P = 0.0017) and the putative barrier of the two major rivers (r 2 = 0.142, P = 0.001).Partial Mantel tests also showed significant correlations between geographic distance when controlling the putative barrier (r 2 = 0.147, P = 0.005) and the putative barrier when controlling geographic distance (r 2 = 0.082, P = 0.0099).
Five populations, Vo, BAU, Rang, SK and Fulb, were detected as putative outlier populations by systematic 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 2 = 0.147, thick line) and excluding five putative outlier populations (filled circles only; y = 0.00006x + 0.0571, r 2 = 0.177, thin line).DPR analysis (Fig. 6), and Vo was identified as a true outlier because the highest r 2 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).

DISCUSSION I nter-and intraspecific divers ity 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(Hasan et al., , 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 (codonpartitioned), 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

Fig. 1 .
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 TableS1.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.

Fig. 3 .
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 twodimensional 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. 6 .
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.

Fig. 7 .
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.

Table 2 .
Genotypic data of 21 microsatellite loci for H. tigerinus 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.