Abstract: Background and Objective: Clarias batrachus is native to Asia, Sumatra was one of the islands that were divided due to the division of Sundaland because of a transgression event. The Bukit Barisan Mountains divide Sumatra into two parts: West and East. Therefore, it is important to research C. batrachus phylogeography. One of the mtDNA genes that can be used for phylogeographic is the CO1 gene. This study examined the association between geographic isolation in Sumatra's Eastern and Western watersheds and the phylogenetic relationships between C. batrachus populations in Sumatra and other Southeast Asia using the CO1 gene. Materials and Methods: Clarias batrachus samples were collected from five populations 7 (Mungo River, Sinamar River, Agam River, Anai River and Lubuak Paraku River) from West Sumatra. The method of this research used the molecular method. The DNA isolation based on the COI gene follows the Invitrogen PureLinkTM Genomic DNA Mini Kit protocol. Haplotype network using software Popart v.1.7 and molecular evolution analysis using software MEGA. Results: Based on a total of 28 sequences (including comparisons and outgroups) with 504 base pairs of the CO1 gene from C. batrachus, three haplogroups and 19 haplotypes were identified using the haplotype analysis. Speciation events between haplogroup 1 estimates are 0.41 Mya haplogroup 2 estimates are 0.50 Mya. Conclusion: The separation of Sundaland and the separation in Sumatra’s Eastern and Western watersheds affects the kinship relationship of these populations with variations in subspecies level.
INTRODUCTION
Siluriformes include 3093 species in 478 genera and 36 families. One of the biggest families in the order Siluriformes, which has 57 species and 16 genera, is the Clariidae1. Characteristics of Clarias have a cylindrical body that resembles an eel, a hard, flat head, a broad mouth with four pairs of barbels, a long dorsal fin and a second organ (the suprabranchial organ) that enables them to survive in waters with little oxygen2,3. Currently, 16 species of Clarias have been found in Indonesia1,4. One of them is Clarias batrachus Linnaeus (1758), which has an elongated body character, no adipose fins, serrated pectoral-fin spines, gray head and back, pale gray belly and white spot on the sides of the body5,6. In addition, there are morphological differences in the shape and width of the supraoccipital process of C. batrachus from populations in Java, Malaysia and Thailand7. According to Teugels et al.8 Clarias can be identified by the shape of the supraoccipital process.
Sumatra is one of the distribution locations for C. batrachus, according to Ferraris and Kottelat et al.1,2. Kottelat6 later proved that C. batrachus originates in Sundaland. During the last glacial maximum, a region known as Sundaland united the Malay Peninsula, Sumatra, Java and Borneo. Sumatra Island is one of the islands created by sea level rise in the Andaman, South China and Java seas9,10. Evolution can occur over a wide range and geographical differences can lead to evolution. Adaptation to the environment causes the accumulated genes to differ from one population to another11. Diversity and specific traits, both morphometrically and genetically, arise as part of the process by which an organism adapts to the environment12. Therefore, it is important to conduct a phylogeographic study of C. batrachus.
Phylogeography combines phylogenetic and biogeographic studies, which aim to determine kinship relationships based on geological history13. Phylogenetics is a method used to model the proximity of a species to other species14. For example, many phylogenetic studies have been carried out, Segherloo et al.15. They reported that the Garra typhlops and Garra lorestanensis exist in sympatry in a single subterranean habitat following a colonization event from the Karun-Dez-Karkheh basins in the Zagros Mountains of Iran. Biogeography is a science that brings together concepts and information from ecology, evolutionary biology, taxonomy, geology, physical geography, paleontology and climatology16. Molecular markers that can be used for phylogeographic studies are genes in mitochondrial DNA, because they are maternal, non-recombinant, have a fast evolutionary rate and have wide intraspecific polymorphism12. One of the mitochondrial DNA genes that can be used as a marker for phylogeographic studies is the CO1 gene.
The CO1 gene is one of the mtDNA genes with slow evolution, with very few deletions and insertions in the sequence and many parts that are conserved17,18. Therefore, the CO1 gene is effective for identifying phylogenetic and evolutionary relationships to the species level, for example, identification in several animal taxa such as insects, birds and fish18-20. Molecular identification studies based on the CO1 gene have been widely used in the phylogenetic and phylogeographic studies in Sumatra. For example, the CO1 gene confirmed Puntius cf. binotatus of the lake of Gunung Tujuh as Barbodes banksi21. Also, by phylogenetic analysis based on the CO1 gene, Rinuak and Badar fish (local name for small fish like Gobies) are the same species22.
Phylogeographic studies with molecular markers have been carried out on several Clarias. The phylogeography of C. gariepinus in the Congo basin, Africa, which is divided into four clades, is supported by data that there is a population divergence of C. gariepinus at 400 kya23. Research using Cyt-b and D-loop genes reported that the phylogeography of C. macrocephalus in Malaysia had not been significant genetic differentiation from the three populations24.
Based on the research results that have been reported, it is known that phylogeographic studies can determine relationships based on geological history. The geographic history of Sundaland allows for the evolution of C. batrachus populations. However, a phylogeographic study of C. batrachus has never been done. Thus, research on phylogeographic studies of C. batrachus on the island of Sumatra and other populations in Southeast Asia based on the CO1 gene needs to be carried out.
MATERIALS AND METHODS
Study area: All the experiments in the studies were performed from November, 2021 to November, 2022 in the Genetic and Biomolecular Laboratory of the Department of Biology, Faculty of Mathematics and Natural Sciences Andalas University, Padang, Indonesia.
Sample collection: The samples used in this study consisted of liver samples from nine individuals C. batrachus. The samples collected were from five populations (Mungo River, Sinamar River, Agam River, Anai River and Lubuak Paraku River) in West Sumatra.
DNA isolation, amplification and sequencing: The DNA isolation steps follow the Invitrogen PureLinkTM Genomic DNA Mini Kit protocol. The results of DNA isolation were then electrophoresed using 1.2% agarose gel. Furthermore, DNA amplification using a PCR reaction mixture consisting of 11 μL Supermix Bioline, 9 μL Nuclease-free water, 1 μL forward primer, 1 μL reverse primer and 3 μL DNA isolation results. Primers for amplification used FISH F1 and FISH R1 referring to Ward et al.20. The PCR was carried out for 35 cycles. The PCR process begins with predenaturation at 95°C for 2 min, followed by denaturation at 94°C for 30 sec, annealing at 53°C for 30 sec and elongation at 72°C for 1 min. The final extension was carried out in the last cycle at 72°C for 10 min. The PCR product was electrophoresed using a 2% agarose gel. The PCR products were sent for sequencing at 1st Base Malaysia.
Haplotype network analysis: The contiguous DNA sequences of C. batrachus were then BLAST and compared with 19 comparison sequences downloaded from the NCBI GenBank. The DNA sequence alignment, haplotype networks using software Popart v.1.725. Before running the application, prepare the alignment of the sequence to be used and the trait that contains the origin of the sequence area. Then the data is imported and select the network menu then select the method of using NP joining, wait until the haplotype network is made. Making haplotype networks using sequence data from comparison species in NCBI.
Molecular clock analysis: Molecular clock using MEGA (Molecular Evolutionary Genetics Analysis) software26. The first step is to enter the sequence and create a phylogenetic tree that has been saved in the nexus format. Then select molecular clock, then input the sequence file and tree file in nexus format. Then, enter the calibration data. For calibration first, enter the year of sample collection and enter the calibration between two different species, this data can be seen in time tree resources. Then enter the analysis settings and then click launch analysis then click ok. After waiting a few moments, the molecular clock display will appear.
RESULTS
Blast analysis: The BLAST results of the CO1 gene sequences of C. batrachus West Sumatra showed a similarity value of 98.5%-99.6% with the sequences found in Genbank. The length of the CO1 gene that was aligned with the comparison sequence was 504 bp. The results of the alignment of C. batrachus are in the 5635-6139 base sequence of the complete mitochondrial genome and the 166-668 base sequence of the complete CO1 gene, consisting of 334 bp (66.2%) conserved sites, 170 bp (33.7%) variable sites, 108 bp (21.4%) parsimony sites and 62 bp (12.3%) singleton sites. The composition of the nucleotide bases in the 504 bp CO1 gene of C. batrachus was 27.5% A (adenine), 15.8% G (guanine), 29.1% T (thymine) and 27.5% C (cytosine). The average percentage of nucleotide bases of adenine+ thymine (A+T) was 56.7%, while guanine+cytosine (G+C) was 43.3%.
The amino acids produced by C. batrachus individuals have two types of mutations that occurred in C. batrachus individuals’ silent mutations and missense mutations. Among the 28 sequences of the COI gene, 22 of them had mutations in their amino acids. One of the silent mutations occurred in the Anai River population at the 17th amino acid that changed the CCC nucleotide bases to CCA nucleotide bases. However, they still produced the amino acid proline (P), the same as other populations. There are seven missense mutations, the 2nd amino acid sequence is the first amino acid formed due to a mutation. Based on the analysis results, almost all C. batrachus in Sumatra population have amino acid mutated, except for the Agam River population. There are four amino acids mutated in the Sumatra populations. One of the mutations occurs in the Mungo River, Sinamar River, and Lubuak Paraku river populations, where the amino acids change from Isoleusine (I) to Valine (V) because the base change from ATT to GTT.
Haplotype network: Haplotype network analysis of C. batrachus CO1 gene sequences with a length of 504 bp formed three haplogroups. In the Fig. 1, the black dot is the median vector that describes unsampled sequences or sequences not used in the study or depicts an extinct ancestor sequence. The hatch sign represents the number of mutations in the nucleotide base. Haplogroup 1 is a haplotype consisting of populations of the Mungo River, Sinamar River, Agam River, Aceh and Malaysia. The C. batrachus, Mungo River, Sinamar River populations belong to haplotype 1(H1). The C. batrachus individuals collected from the Agam River belong to haplotype 2 (H2).
| Fig. 1: | Visualization of haplotype network of C. batrachus based on the CO1 gene using PopArt software |
The C. batrachus individuals collected from the Lubuak Paraku River belong to haplotype 4 (H4). The Clarias batrachus population of Aceh is classified into haplotype 5 (H5) and the population of Malaysia is classified into haplotypes 8 and 9 (H8 and H9). Based on these haplotypes, there is a high variation in haplogroup 1.
Haplogroup 2 is a haplotype consisting of populations of Anai River and Java. The population of Anai River belongs to haplotype 3 (H3) and the population of Java belongs to haplotypes 6 and 7 (H6 and H7). Haplogroup 3 is a haplotype consisting of populations from Thailand and the Philippines. The Thailand population belongs to haplotypes 10 and 11 (H10 and H11) and the Philippines population belongs to haplotypes 12 and 13 (H12 and H13). Haplotype describes a high variation according to the variation that occurs in haplogroup 2.
Molecular clock: Based on molecular clock analysis, it can be seen that the separation of taxa between Clarias batrachus and other Clarias is estimated to have occurred 20 million years ago in the early Miocene. The results in Fig. 2, were in accordance with the haplotype network, the divergence between C. batrachus Philippines-Thailand (haplogroup 3) and C. batrachus Sumatra, Java and Malaysia (haplogroups 1 and 2) happened about 6 Mya in the late Miocene. Then, Clarias batrachus Sumatra, Java and Malaysia experienced divergence of about 5.2 Mya in the early Pliocene, the haplotype network was haplogroup 1 and haplogroup 2. Speciation events between C. batrachus Aceh, Malaysia, Mungo River, Sinamar River and Lubuak Paraku River with Agam River is forecast at 0.41 Mya. While in C. batrachus Java with Anai River, it is forecast at 0.50 Mya.
| Fig. 2: | Visualization of molecular clock of C. batrachus based on the CO1 gene using MEGA 11 software |
DISCUSSION
The conserved site value indicated that the nucleotide bases remained or did not change in all species in the same genus27. Therefore, the value of the site variable shows different nucleotide bases in each species so that these bases can be a differentiator between species in the same genus. In line with the research on nucleotide composition C. batrachus was 27.4% A (adenine), 28.4% T (thymine), 17.9% G (guanine) and 26.2% C (cytosine)28. In vertebrates, the number of bases (A+T) is higher than that of bases (G+C) because the double bonds of the nucleotide bases of adenine and thymine are unstable and more prone to mutations29.
Silent mutations only change the base composition but do not change the amino acids produced, while missense mutations change the base and change amino acids to other amino acids30. For example, C. batrachus Sumatra has a difference of 46 bases from C. batrachus Java, Aceh and Southeast Asia, however, seven amino acids experience missense mutations. Suppose it is connected with the separation of Sundaland around 12,000 years ago. In that case, it can be interpreted that the changing of the amino acid gene of C. batrachus CO1 took place slowly. The interpretation follows the opinion of Lynch and Jarrell31, who stated that in the CO1 gene, amino acid changes occur slowly.
Haplotype 1 (H1) illustrates the occurrence of haplotype sharing in C. batrachus from Mungo River and Sinamar River. The two populations could be in one haplotype because they come from a common ancestor32. In this study, the haplotype diversity (Hd) value of C. batrachus was 0.92992. According to the haplotype diversity of an individual has a range of values from 0 to 130. The haplotype diversity is high if the value is >0.5 and low if the value is 0.5. Based on the range of values, the diversity of haplotypes of C. batrachus in this study was high. The high haplotype diversity in C. batrachus was thought to be because the samples used were from distant populations.
Calibration information is primarily based on the fossil record, but geological events from that period are also available33. The molecular clock calibration used tMRCA used data, C. gabonensis-C. pachynema divergence around 12.2-17.0 Mya and C. fuscus-C. macrocephalus around 13.8-20.2 Mya34,35. The separation of taxa between C. batrachus and other Clarias is estimated to have occurred 20 million years ago in the early Miocene. Clariidae based gene cyt b originated about 50 Mya (million years ago) for Asian Clarias species. However, recent lineages colonized Africa and South-East Asia from a central Asian origin about 15 Mya independently36.
The divergence between C. batrachus Philippines-Thailand and C. batrachus Sumatra, Java and Malaysia happened about 6 Mya in the late Miocene. Then, C. batrachus Sumatra, Java and Malaysia experienced divergence of about 5.2 Mya in the early Pliocene. Hall37 stated that the 5 Mya collision of the Philippines arc and the Eurasian continental margin occurred in Taiwan. Speciation events between C. batrachus Aceh, Malaysia, Mungo River, Sinamar River and Lubuak Paraku River with Agam River is forecast at 0.41 Mya. While in C. batrachus Java with Anai, it is forecast at 0.50 Mya. There was an overlap between Java and Sumatra at 20-10 Mya (pre-middle Miocene) when Java experienced a counter-clock wise 30°, at that time Sumatra and Java were closer37.
During the sea level low-stand, Borneo, Sumatra and Java were connected to the mainland and Rivers got connected through the development of four main palaeodrainages straddling islands38. Sumatra (Anai River) and Java (Banten) are connected to the north Sunda River system. Meanwhile, Sumatra (Mungo River, Sinamar River, Agam River, Lubuak Paraku River and Aceh) with Malaysia are connected to the Malacca straits River system. Other research, found the similarity in the genetic diversity of the populations from western Borneo and Southeastern Sumatra of the catfish Hemibagrus nemurus were a consequence of past drainage connectivity during sea level low-stands39. So, it can be concluded that the individuals found in haplogroups 1 and 2 results from past Rivers connectivity. Further accurate morphological examination of these clades will hopefully enable us to state their new species status by the definition of diagnostic characters.
CONCLUSION
The result showed that based on haplotype and molecular clock analysis in this study. It is known that based on phylogeographic studies on C. batrachus in Sumatra and other populations in Southeast Asia. The separation of Sundaland affects the kinship relationship of these populations with variations in subspecies level. Speciation events between haplogroup one estimates are 0.41 Mya and haplogroup two estimates are 0.50 Mya.
SIGNIFICANCE STATEMENT
This research was conducted to determine the association between geographic isolation and genetic variation, especially the phylogeography of Sumatra's Eastern and Western watersheds and the phylogenetic relationships between C. batrachus populations in Sumatra and other Southeast Asia using the CO1 gene.
ACKNOWLEDGMENTS
The research was supported by Andalas University under the RPT Unand scheme (No.T/47/UN.16.17/PT.01.03/ IS-RPT/2022). The authors would like to express our sincere thanks to the Head of the Biology Department of Andalas University for the Laboratory work permit. We also would like to thank our team of Geneticists who helped us with sample collection and laboratory works in the Genetic and Biomolecular Laboratory, Faculty of Mathematics and Sciences Andalas University, Padang, Indonesia.