Parapatric speciation found in Limnonectes (Dicroglossidae) species in Thailand

Abstract

Biodiversity stems from speciation, and species identification allows scientists to measure species biodiversity. However, some closely related species, particularly amphibians that choose to inhabit the same environment, are difficult to identify. Juvenile specimens of Limnonectes gyldenstolpei (Lg) and L. taylori (Lt) were used as models to study the speciation process in populations with a parapatric distribution. Samples of the two species were collected from six provinces in the north and northeast of Thailand. Inter-simple sequence repeat (ISSR) data of five loci showed 96 polymorphic bands and species-specific bands of the primers UBC811 (Lg) and UBC824 (Lt). There was limited genetic variation within the same population, 8% in Lg and 2% in Lt, whereas a high genetic variation was observed between populations from different provinces 90% and 97% in Lg and Lt, respectively. An unweighted pair-group method using arithmetic averages (UPGMA) tree was constructed for each population using Jaccard’s similarity efficient. The tree demonstrated genetic differentiation between the two species. Principal coordinates analysis also supported the UPGMA diagram by showing no clear clusters with close genetic distances. The results of the molecular variance analysis within the population were high, indicating that they were genetically similar, and that microhabitat is an important factor influencing genetic distribution. Overall, the present study proves the hypothesis that the distribution of the two species is in a Parapatric form and that the ISSR technique is inappropriate tool to study the differentiation of a species complex.

Introduction

Speciation is an important cause of diversity (Schluter and Pennell, 2017). Speciationmay occur either (1) among species that coexist in the same microhabitat, which includes sympatric and parapatric speciation, or (2) among species with a mutually exclusiv geographic range or allopatric speciation (Cain, 1953; Mayr, 1947; Mayr et al., 1963; Rivas, 1964). Morphologically similar species are difficult to identify, particularly in their larval and juvenile stages. In the present study, we found two similar-looking juveniles of the genus Limnonectes inhabiting the same microhabitats in northern and north-eastern Thailand, which then develop into adults that are easily distinguishable. This makes it interesting to study the genetic variability between L. gyldenstolpei (Lg) and L. taylori (Lt) found in the same areas. It is possible that the genes of these two species are similar.

Dicroglossine frogs of the genus Limnonectes are found across most Asian and Southeast Asian countries. Their dominant characteristic is the odontoid processes called fangs on their lower jaw (Emerson et al., 2000; Inger and Stuart, 2010). L. kuhlii Tschudi (1838) is a complex of cryptic species (McLeod, 2010a). In 2008, a big-headed species of Limnonectes, L. megastomias, phenotypically similar to L. gyldenstolpei in the L. kuhlii complex, was collected from eastern Thailand by David S. McLeod. Two species of the L. kuhlii complex from Thailand were described by Matsui et al. (2010) based on molecular and morphological evidence: L. taylori, a northern Thai lineage, identified as Lineage 12 by McLeod (2010b), and L. jarujini, a southern Thai lineage discussed by Taylor (1962). Amphibians are a critical component of certain ecosystems. In addition, they play a role in the proper functioning of communities that live in their vicinity and are dependent on other animals, for instance, for food. The loss of amphibians from an ecosystem can change primary production, affecting food webs from aquatic insects up to riparian predators (Hocking and Babbitt, 2014; Meredith et al., 2016; Whiles et al., 2006). The permeable skin % of frogs is sensitive to environmental changes, and frogs can act as bio-indicators (Campos et al., 2014; Hopkins, 2007; Lebboroni et al., 2006). Amphibian populations are declining, and some are on the edge of extinction. Small changes in the environment can cause speciation that can be traced back by means of genetic studies. To study genetic variability, many molecular tools are available, including a Technique based on inter-simple sequence repeats (ISSRs). ISSRs are regions in the genome located between microsatellite sequences (Zietkiewicz et al., 1994). To use ISSR markers, prior knowledge of the ISSR target sequences is not required; moreover, the markers are highly reproducible due to their primer length and the high stringency achieved by the annealing temperature, and they have been found to provide highly polymorphic fingerprints (Bornet and Branchard, 2001; Kojima et al., 1998; Zietkiewicz et al., 1994). Here, we examined parapatric speciation in two frog species of the L. kuhlii complex. We randomly sampled the areas where records of the species were found. To study their genetic variability, we then subjected juvenile specimens of L. gyldenstolpei and L. taylori to the ISSR technique.

Material and Methods

Sampling, DNA extraction, and polymerase chain reaction (PCR) amplification using an ISSR marker in total, 64 L. taylori and 59 L. gyldenstolpei samples were collected from six provinces in the north and northeast of Thailand (Table 1 and Figure 1).

Table 1: Population locations and sizes of L. taylori and L. gyldenstolpei

A GPS location for each specimen was recorded. Fragments of the liver, muscle, or toe were preserved in 95% ethanol for DNA extraction, which was performed with a Genomic DNA Extraction Kit (RBC, Taiwan) following the manufacturer’s instructions.

Spectrophotometry (BioDrop DUO, UK) and was used to determine the concentration of the DNA obtained, and its quality was assessed using electrophoresis in agarose/TAE (1.5%) gel. A total of 20 primers of the UBC primer set (University of British Columbia) were tested for the ISSR analysis (Table 2). PCR amplifications were performed in a 15 μl reaction volume containing 10 ng of template DNA, 0.8 μM primer, 0.2 mM dNTP, 1x buffer containing 2 μM Mg²⁺, and 1U Taq polymerase. The cycling condition was as follows: denaturation step at 94°C for 5 min, 35 cycles of denaturation at 94°C for 1 min, annealing temperature (Table 2) for 1 min, and an extension at 72°C for 1.3 min, followed by a final extension at 72°C for 5 min (Machkour-M'Rabet et al., 2009).

Table 2: Primers and their annealing temperature, polymorphic band and size range

After amplification, samples comprising 7 μL of the PCR reaction mixture were subjected to electrophoresis on 1.5% agarose gels and 3× GelRed (Biotium, USA). Electrophoresis was performed in 1× TAE buffer at 130 V for 34 min. The bands were visualized under ultraviolet light, and the gel was photographed for analysis.

Data analysis

The fragments of ISSR were scored as binary characters: present (1) or absent (0) for each individual and each primer. To describe the genetic structure and variability among and between populations, non-parametric analysis of molecular variance (AMOVA) (Excoffier et al., 1992) was performed and a scatter plot of principal coordinates was constructed using the program principal coordinates analysis (PCoA) and GENALEX version 6.502. Nei’s genetic distances were used to perform an unweighted pair group method using arithmetic averages (UPGMA) cluster analysis, and phylogenetic trees were constructed using the POPTREE2 program (Takezaki et al., 2009). The relationship between genetic isolation and geographic1distance was analyzed using Isolation by Distance (IBD) Web Service version 3.23.

Results

ISSR PCR is a technique that uses microsatellite sequences and randomly amplifies them in a similar manner as in RAPD. Thus, ISSR PCR is sensitive to reaction parameters. A total of 20 ISSR primers were screened with the 123 samples as mentioned above. After comparing the effects of the template DNA concentration and temperature during the annealing stage of amplification, only five primers were able to generate clear and reproducible DNA fragments (Table 3). These were then selected for further analysis. A total of 96 alleles were revealed by five ISSR loci across 123 individuals with a size ranging between 100–2,000 base pairs (Figure 2). The primers UBC845 and UBC811 showed the highest and lowest number of ISSR bands, respectively. We found UBC811 to be species-specific with the size of the amplified DNA fragment being 300 base pairs in Lg and 1500 base pairs in Lt.

Table 3: AMOVA results of L. gyldenstopei and L. taylori based on ISSR

The dendrogram obtained from the UPGMA cluster analysis of the ISSR marker is shown in Figure 3. The UPGMA dendrogram revealed a close association between L. gyldenstolpei and L. taylori. According to the tree, L. gyldenstolpei divided from L. taylori at IV. The similarity index was 0.385–1.260 (Table 4) and 0.000–2.585 (Table 5) for L. gyldenstolpei and L. taylori, respectively. The sub-clusters correlated with the height of the sampling locations.

Table 4: Genetic distance of L. gyldenstopei accession using ISSR markers

Population structure

AMOVA was used to partition the genetic variation using hierarchical analysis from the distance matrix. AMOVA showed a high percentage of variation within populations; 90% and 97% for Lg and Lt, respectively. There was little genetic difference between the populations collected from the six different regions (provinces); Table 5. Genetic distance of L. taylori accession using ISSR markerss 8% and 2% in Lg and Lt, respectively. Likewise, the genetic difference among populations within the same region was low; only 2% and 1% in Lg and Lt, respectively (Table 3). The PCoA analysis did not clearly separate the two-species based on the first two principal component scores, which accounted for 4.62% and 3.84% of the total variations (Figure 4). To study the relationships within the populations of L. gyldenstolpei and L. taylori, PCoA was performed. According to the result, the first group (red), which accounted for 4.62%, did not clearly separate from the second group (gray), which accounted for 3.84% (Figure 4a).

Table 5: Genetic distance of L. taylori accession using ISSR markers

Isolation by distance

The box plots that resulted from height geography daHta analysis provided a clear graphical representation of the distribution of L. gyldenstolpei and L. taylori (Figure 4b). It was observed that the L. taylori population has a wider distribution than that of L. gyldenstolpei. The mean elevation was found to be 709.88 m and 891.45 m for L. gyldenstolpei and L. taylori, respectively, and a t-test showed no significant difference (t=1.025, p-value=0.3194) in the elevation at which the two species were present. The IBD test showed a relationship between the genetic distance (Fst) and geographic distance. For L. gyldenstolpei, there was a non-significant correlation (r=−0.0963 and p-value=0.6330) between the genetic and geographical distances (Figure 4c). However, for L. taylori, the correlation was significant (r=0.4257, p-value=0.0020) (p<0.05 is the level of significance) (Figure 4d).

Discussion

In the present study, we investigated the genetic variability of L. gyldenstolpei and L. Taylori in the north and northeast of Thailand using ISSR markers. We found a clear genetic difference between the two species from the UPGMA analysis (data from the individual is grouped into the population), but the genetic distance between the two-frog species was very low, showing that they were closely related. The UPGMA analysis presented here provides an interesting result, particularly when considering the population of L. taylori in Ban Pa Sak Ngam, Chaing Mai (latitude=18°59¢16.0¢¢, longitude=99°06¢52.4¢¢, elevation=402 km).

The reason for this population being distinct from the others could be due to the presence of the Mae Kuang Udom Thara dam, which serves as a geographical barrier for frogs. We propose that the population of L. taylori in Ban Pa Sak Ngam has an allopatric distribution, similar to that of the Chinese wood frog (Rana chensinensis) whose genetic variation is caused by a mountain ridge barrier (Zhan et al. (2009). For both frog species, our data provides a high genetic variance within the same populations based on the AMOVA results. The PCoA plot showed an overlapped distribution of the two populations suggesting that there might share the genetic information from deep coalescence as the two populations are treated as the member of the L. khulii complex. In the present study, we found no significant correlation between elevation and genetic distance. The height at which the frogs inhabit does not significantly affect their genetic distance.

Among the populations of L. taylori, where inbreeding and a high degree of population subdivision are observed, the IBD test indicated a significant effect of geographical distance on the animals’ genetic distance. This effect, however, was not significant among the specimens of L. gyldenstolpei. The AMOVA and PCoA analyses support the phylogenetics of the two species as proposed by Funk et al. (2008). The genetic distance between the two species is very low, indicating their close genetic relationship. A very low genetic variation within the genus Limnonectes suggests that there is widespread inbreeding within this genus, which increases its risk of extinction. Inbreeding can be explained by the fact that amphibians are conservative when choosing their spawning ponds and have a low migration activity (Chernyshov and Truveller, 2006; Zhan et al., 2009). The inbreeding caused by geographical barrier was also reported in previous researches showing that mountain ridges act as barriers to gene flow in Columbia spotted frogs (Funk et al., 2005) and the isolation-by-distance factor applies to the plains populations (Palo et al., 2003). Our results could be interpreted as being consistent with the speciation theory that depends on distribution ranges and geographic barriers, which could be classified as allopatric, sympatric, or parapatric speciation (Coyne and Orr, 2004; Mayr, 1942; Schluter, 2001). The study of species co-occurrence and parapatric distributions can provide an evolutionary history of amphibian species, for example, in tree frogs (Hyla molleri and H. meridionalis) and midwife toads (Alytes obstericans and A. cisternasii). Their combined fuzzy intersection and niche comparison metrics showed the effect of the environment or niche on the survival until segregation (Reino et al., 2017). Similarly, Posso-Terranova and Andrés (2016) showed the relative contribution of geographical and environmental factors to the diversification of poison frogs (Oophaga) using gene genealogies and the phylogenetic method. The study found a parapatric distribution of sister taxa and no significant effect of elevation on the genetic distance. High levels of genetic variability and differentiation are typical for amphibians in general and exceed that for other groups of vertebrates (Mezhzherin, 1992).

Conclusion

In conclusion, juvenile L. gyldenstolpei and L. taylori inhabit geographically different habitats, but they have a similar morphology. Our study shows that both internal and external dynamics affect spatial genetic subdivision, which could potentially lead to Parapatric and allopatric speciations. The present study demonstrates that the ISSR technique is an appropriate tool to study the differentiation of a species complex. L. gyldenstolpei and L. taylori are genetically distinct. The genetic variation and distance point to the recent separation of the two species of Limnonectes. Further studies in other organisms using techniques mentioned in this paper should be conducted in the future.

Acknowledgments

We would like to express our sincere gratitude to Mr. Somchai Nimuan and Ms. Patcharawadee Thongkumkoon for their precious advice. This project was funded by the Graduate School, Chiang Mai University, Thailand. The Institute of Animals for Scientific Purpose Development (IAD) provided permission for the use of animals in this study (No. U1-01205-2558).