Population variations in the Fungoid Frog Hylarana malabarica ( Anura : Ranidae ) from northern Western Ghats of India

Widely distributed species often show interpopulation variation. Studying such variations can be helpful in understanding contributing factors and distinguishing widespread species and species complexes. We studied six populations of Hylarana malabarica distributed along the northern Western Ghats of India using morphometric and genetic analysis. Of 24 size-adjusted morphometric characters, 14 were significantly different among populations. Hierarchical clustering and discriminant analysis of morphometric characters suggested that the six populations form at least four distinct clusters. Analysis of morphometric data was supported by genetic polymorphism data obtained by the Randomly Amplified Polymorphic DNA (RAPD) method. Since the similarity and variation observed among populations was independent of their spatial distribution, it is possible that this widely-distributed species may be a species complex.


INTRODUCTION
The high level of endemism among vertebrate and plant species has led the Western Ghats of India and Sri Lanka to be considered a hotspot of global biological diversity (Myers et al. 2000).The Western Ghats are rich in amphibian fauna, and while the first species was discovered in the early 1800s the discovery trend for Western Ghats amphibians has yet to reach a plateau (Aravind et al. 2007).A recent record of a new frog family from the Western Ghats (Biju & Bossuyt 2003) reflects the limitation of our knowledge of the amphibian diversity of this important biogeographic region (Hedges 2003), as do recent descriptions of new amphibian species (Gururaja et al. 2007;Kuramoto et al. 2007;Biju & Bossuyt 2009;Zachariah et al. 2011).
Recent investigations of Western Ghats amphibians have shown that several populations contain cryptic species revealed by in-depth study (Kuramoto et al. 2007).The broad application of molecular techniques to phylogenetic reconstruction has been used effectively to unveil haplotypes or morphologically cryptic species (Kuramoto et al. 2007;Biju et al. 2009).Discovering these cryptic species is important for our understanding of species richness and essential for the design and implementation of conservation action plans (Bickford et al. 2007;McLeod 2010).This is especially true for species which show wide distribution, which may turn out to be species complexes (Inger et al. 2009).
The Fungoid Frog Hylarana malabarica (Tschudi, 1838) is widely distributed in peninsular India (Daniel 2000;Padhye & Ghate 2002), Assam and Meghalaya (Dutta 1997).Because of its wide distribution the species is categorized under Least Concern in the IUCN Red List

Sampling
Six study sites were identified (Table 1, Image 1).Field visits were conducted during the breeding season (July to September) in 2009 and 2010.While a number of individuals were collected for morphometry, only two to three individuals were brought back to the laboratory; remaining individuals were released back in the same habitat.Individuals brought back to the laboratory were etherized and fixed in 100% ethanol.The tissues (liver and muscles) from these specimens were used for DNA extraction.Amboli and Dhamapur specimens were collected from the same congregation, while Velneshwar, Tamhini, Kolvan and Ghatghar specimens were collected from different places in the same area.Fifteen specimens collected during the of Threatened Species (Biju et al. 2004).Despite its widespread distribution, H. malabarica shows patchy distribution in the northern Western Ghats.In our initial studies we observed variation among the individuals from different populations of H. malabarica from the northern Western Ghats, and in the current work we have studied morphological and genetic variation among six populations of H. malabarica collected from six isolated locations.Our analysis, based on morphological and genetic studies, indicates that the six populations of H. malabarica form at least four separate clusters, raising the possibility that H. malabarica could be a species complex.19, 20, 21, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 35 and 36.

Morphometric analysis
Morphometry was done in the field with the help of digital Vernier calipers (Areospace®, least count 0.01mm).Morphometric measurements were taken for snout-vent length and 24 different characters (Table 2).Since females were poorly represented in our samples we considered only males for morphometric analysis to avoid any effect of sexual dimorphism in the morphological analysis.
Since different individuals differed in their snout to vent length and all other morphometric measurements were strongly correlated with the snout-vent length, we first adjusted the morphological measurements of all the individuals for mean snout to vent length so as to remove the size and shape effect.We used allometric adjustment (Lleonart et al. 2000) given by the formula, M adj = M obs (SL 0 /SL) b , where, M adj is the size adjusted length of a morphological character, M obs is the observed length of the character, SL 0 is the mean snout to vent length of all the individuals, SL is the snout to vent length of the given individual, and, b is the allometric exponent of power function relation between the character and the standard length of all individuals (in other words it is the slope of the line between logM obs versus logSL).Efficiency of size adjustment was assessed by testing the significance of correlation between transformed variables and standard length, which was not significantly different from zero.This corrected morphometric data of 24 characters was used for further statistical analysis.Since allometric adjustment nullifies both the size and shape bias the resultant data was independent of the ontogenic variations among the individuals.
We performed ANOVA (Analysis of Variance) to understand which standardized morphometric characters differed among the populations.Since multiple tests were performed on the same data we applied sequential Bonferroni correction to the a values wherever applicable.Variables which where significant after sequential Bonferroni correction were further analyzed using unpaired t test to find out differences between populations.We performed hierarchical cluster analysis based on mean values of significantly different standardized characters for a given population to understand the general pattern in similarity among the populations.Euclidian distances between populations were calculated and Ward's method was used for clustering.
Discriminant analysis (DA) was performed on the

DNA Extraction and Randomly Amplified Polymorphic DNA (RAPD) analysis
The tissue was digested at 50°C for two hours using the extraction buffer (0.1M NaCl, 0.05 M Tris-HCl, 0.01M EDTA, 1%SDS) with 15µl Proteinase K (20mg/ ml).DNA was then extracted using the conventional phenol-chloroform method (Sambrook et al. 1989).Polymerase chain reaction was performed to randomly amplify the polymorphic DNA.Primers used for the study were based on Wei et al. (2001).Primers that gave consistent results under repeated experiments are given in the Table 3.The PCR amplifications were conducted in 20µl reaction volume containing 2µl of template DNA, 2µl of 10X reaction buffer (100 mM Tris pH9.0, 500 mM KCl, 15 mM MgCl2, 0.1% Gelatin), 1µl of 25mM MgCl 2, 1µl of 10mM dNTPs, 1µl of primer, 0.8µl Taq polymerase and 11.2µl sterile distilled water.The cycling profile used was 5 min at 95 0 C, and 35 cycles of 1 min at 95 0 C, 1 min at 30 0 C and 2 min at 70 0 C, followed by 10 min at 72 0 C. Amplified DNA fragments were checked using 1% agarose gel electrophoresis and further analyzed.
Presence and absence of a given fragment amplified in RAPD was represented by '1' and '0' characters respectively.Only clear and reproducible bands were recorded as '1'.No or non-reproducible bands were recorded as '0'.We used Nei and Li (NL) coefficient for comparison between the RAPD patterns between diffident individuals (Lamboy 1994).The NL coefficient, which denotes a value of the similarity between two samples (Nei & Li 1979), is given by the formula, NL = 2a/(b+c), where a is the number of similar bands from two samples, and b and c are the total numbers of bands from each sample.Based on the NL similarity coefficient, which ranged between 0 and 1, we performed cluster analysis using five methods, viz., single linkage, complete linkage, flexible linkage, unweight pair-group average and weight pair-group average (Sneath & Sokal 1973) and the most consistent tree topology was chosen for plotting.

Morphometric analysis
Out of total 24 size-adjusted morphological characters, 14 were significantly different for six populations, while another three were different but could not qualify as significant after Bonferroni correction (Table 2).Differences in these 14 characters between six populations are given in Table 4. Tamhini Population differed from all the other populations in three characters namely, toe 5 length, hind limb length and inter narial distance.Even though, there was no particular character in which the Dhamapur population differed from all others, this population could be distinguished by a combination of characters: inter-narial distance, width of upper eyelid, tympanum diameter vertical, tympanum diameter horizontal, tympanum to eye distance, forelimb length, hindlimb length, tibia length, toe 4 length and toe 5 length.A maximum number of characters were different between the Tamhini and Dhamapur populations, followed by Tamhini and Velneshwar.Least differences were seen between Kolvan and Ghatghar and Velneshwar and Amboli.Image 2 shows individuls from different populations used in the analysis.
A dendrogram based on the mean values of the standardized characters for a population is shown in Fig. 1.Tamhini and Dhamapur populations were distinctly different from other populations, while Kolvan and Ghatghar, and Amboli and Velneshwar shared more similarity with each other.
Discriminant Analysis extracted five factors, out of which first three factors explained around 86.19% of the total variation in the data.The means vectors of the six populations were significantly different (Pillai's trace = 3.439, F 70,220 = 6.923,P < 0.0001), indicating that the six populations formed six significantly different clusters (Fig. 2).Higher values of variables such as inter-narial distance followed by toe 3 length and eye diameter separated Tamhini from other populations, while higher values of variables like width of upper eyelid, tympanum diameter horizontal and forelimb length separated Dhamapur population from all the rest of the populations on the first two axes (Table 5).Among the remaining four populations, Kolvan and Ghatghar had negative factor loading on the third canonical axis, while Amboli and Velneshwar had positive factor loading (Fig. 2a).Fisher's distances between centroids of all six populations were significant indicating that these six populations formed different clusters (Fig. 2b).However, Fisher's distance between  centroids of Amboli and Velneswar was the lowest (3.205) while Fisher's distances between Tamhini and all other populations were very high (Fig. 2b).
Standardized factor coefficient (Table 5) suggests that on the third axis head width, tympanum diameter vertical, tympanum diameter horizontal had high positive factor loading, while hindlimb length, forelimb length, inter-orbital distance and nostril to snout distance had negative factor loading with high magnitude.Thus these characters separate Amboli and Velneshwar from, Kolvan and Ghatghar populations.

RAPD analysis
Out of the five methods of cladistic analysis used: single linkage, complete linkage, flexible linkage, unweight pair-group average and weight pair-group average, the last three methods gave consistent tree topology while the first two gave two different tree topologies.A consensus tree based on NL coefficient and flexible linkage, unweight pair-group average and weight pair-group average methods is shown in Image 3.
RAPD analysis revealed four different clusters.Tamhini population had the least genetic similarity with all other populations, while Dhamapur formed a separate cluster.Among the remaining four populations, Kolvan and Ghatghar formed one cluster while Amboli and Velneshwar formed another cluster (Image 3).This pattern is similar to the pattern depicted in the morphometric data (Fig. 1) and hence supports the results of morphometric analyses.

DISCUSSION
Both morphological and genetic analysis revealed that the six populations in the current study lie in at least four different clusters: (1) Tamhini, (2) Dhamapur, (3) Kolvan and Ghatghar, and (4) Amboli and Velneshwar.The morphological differences between the populations are likely to be independent of sexual dimorphism, as we have considered only males in the population.Further, we nullified the effect of ontogenetic allometry as we applied allometric adjustments to the data to nullify any size and shape bias as suggested by Lleonart et al. (2000).Further, morphological as well as genetic similarity and differences among the six isolated populations were not dependent on their geographical distances (Table 1).Kolvan and Ghatghar populations shared more similarity (Fig. 1 any similarity (Fig. 1, Image 3) yet are 20km apart.Similarly, Amboli population shared more similarity with Velneshwar population than with Dhamapur population, even though Amboli and Dhamapur are just 44km apart while Amboli and Velneshwar are 182km apart.Further, Amboli and Velneshwar have a large difference in altitude (Table 1), as Velneshwar is on the coastline while Amboli is on the crest line of the Western Ghats (which form a geographical barrier between these two populations).There is also a difference of 2 0 in the latitudinal distribution of these two populations (Table1).Such kind of pattern suggests possibility of more than one species that are together considered as Hylarana malabarica.
Recent trends in the discoveries of new species of anurans (Gururaja et al. 2007;Kuramoto et al. 2007;Biju & Bossuyt 2009;Biju et al. 2009;Zachariah et al. 2011) suggests that there are likely to be many more species of anurans still waiting to be described from the Western Ghats of India.The continual increase in the discovery trend of Western Ghats' amphibians (Aravind et al. 2007), further bolsters this fact.It is possible that several of these species could be cryptic with no apparent easily distinguishable morphological differences (Biju et al. 2009).Such species will require detailed morphological and molecular phylogenetic studies for establishing their taxonomic status.For example, Kuromoto et al. ( 2007) described four cryptic species of anuran genus Fejervarya from central Western Ghats.These four species of Fejervarya are not easily distinguishable morphologically but show their distinctness in both detailed morphometric analysis and molecular analysis based on DNA sequencing (Kuromoto et al. 2007;Meenakshi et al. 2010).Our finding of morphological and genetic  variations in different populations of Hylarana malabarica suggests the possibility of recent events in speciation and presence of cryptic species.However, since genetic studies using RAPD suffers from low reproducibility of results, to bolster our arguments it is essential to study the molecular markers in different populations of H. malabarica.Molecular phylogeny of these different populations will be able to help us in separating the population level variations from the species level variations and help us in studying the monophyly of H. malabarica.Méndez et al. (2004) also suggested similar strategy to reveal the presence of recent speciation in Bufo spinulosus, who conducted similar study on this species.Additionally, ecological studies on these populations using niche modeling method (Raxworthy et al. 2007) will also be interesting and may probably help in resolving the phylogeny of H. malabarica.Rissler & Apodaca (2007) have shown the application of ecological niche modeling in defining cryptic species in Black Salamander Aneides flavipunctatus.
Recent trends in the amphibian taxonomy have revealed several lineages of cryptic species (Stuart et al. 2006;Elmer et al. 2007) especially in the wide spread species (Inger et al. 2009;McLead 2010).Understanding this cryptic diversity is essential for species management and conservation (Bickford et al. 2007;McLead 2010).Hylarana malabarica is assessed as Least Concern (Biju et al. 2004) owing to its widespread distribution in India with no major widespread threats.However, if our assertion that the Hylarana malabarica is a species complex with several cryptic species is true, then it is possible that some of the cryptic species might have more restricted distribution and may require immediate conservation attention.
in the collection of Department of Zoology, Abasaheb Garware College, Pune, under the accession numbers AGCZRL Amphibia

Figure 1 .
Figure 1.Dendrogram based on mean values of 24 standardized morphological characters for six populations.

Figure 2 .
Figure 2. Discriminant Analysis of size adjusted morphological data for six populations.(a) Factor loading of individuals on the first three discriminant axes and (b) Fisher's distances between centroids of the six populations (blue coloured cells) and P values for Fisher's distances (red coloured cells).
analysis.(a) Consensus tree based on NL coefficient and flexible linkage, unweight pair-group average and weight pair-group average methods (b) RAPD pattern for RM15 and (c) RAPD pattern for RM17.RAPD patterns for populations not considered in the present study are deleted from the right hand side of the gel pictures.Code in the parenthesis is specimen number.M is marker DNA.

Table 1 . Sampling localities. Table 2. ANOVA depicting the size adjusted characters that are significantly different in six populations. Significant P values are depicted in bold.
(Legendre & Legendre 1998)ations form distinct clusters as well as to identify the discriminating characters(Legendre & Legendre 1998).DA supposes priory groups, which in the current case populations situated in different geographical locations.DA explicitly attempts to model the difference between the classes of data by extracting factors that maximize inter class variation and minimize intra class variations.DA is therefore, more appropriate choice than principle component analysis (PCA), which gives equal weight to all the available variables because of which it cannot reveal the differences among closely related clusters in less number of dimensions.However, since DA considers prior groups, to test whether our analysis is biased by this grouping we tested for intra-group homogeneity by two methods.First, the null hypothesis, which states that the mean vectors of the six populations are equal, (Harris 2001)cant after sequential Bonferroni correction significantly different characters was tested using Pillai's trace(Harris 2001).Second, we calculated Mahalanobis distances(Harris 2001)among the individuals and computed Fisher's distances between six populations as the distance between the centroids of the clusters, divided by the sum of their standard deviations to check if the clusters formed by six populations are significantly different.Statistical analysis was performed in Microsoft EXCEL ® and Systat 12 ® .

Table 4 . Population wise differences in the size adjusted morphometric characters which were significant in ANOVA after Bonferroni correction. P values of unpaitred t test are provided. Significant P values after Bonferroni correction are depicted in bold.
* Locality abbreviation as perTable 1; ** Character abbreviation as per Table 2.