2016 Volume 91 Issue 4 Pages 217-227
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.
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.
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.
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.
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).
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 (AE), the observed (HO) and expected heterozygosity (HE), 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 FST and Nei’s DA distance (Nei et al., 1983) among the populations. We then constructed multidimensional scaling based on the pairwise FST 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).
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 FST [FST /(1 – FST )] (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 FST 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 r2 values to identify the true outlier.
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.
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.
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 AE, HE 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 FST value was 0.030 and the pairwise FST and Nei’s DA 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).
Abbreviations: N, number of individuals; A, number of alleles; AE, number of effective alleles; HO, observed heterozygosity; HE, expected heterozygosity; F, fixation index.
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.
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.
A moderate IBD pattern was observed within the H. tigerinus population because a plot of the relationship between pairwise linearized FST and geographic distance showed a gradual positive correlation (Fig. 5). Mantel tests detected significant correlations between the geographic distance (r2 = 0.147, P = 0.0017) and the putative barrier of the two major rivers (r2 = 0.142, P = 0.001). Partial Mantel tests also showed significant correlations between geographic distance when controlling the putative barrier (r2 = 0.147, P = 0.005) and the putative barrier when controlling geographic distance (r2 = 0.082, P = 0.0099).
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, r2 = 0.147, thick line) and excluding five putative outlier populations (filled circles only; y = 0.00006x + 0.0571, r2 = 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 r2 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).
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.
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.
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.
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 (FST = 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 (FST = 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.
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.
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.