Genetic Identification of Two Mudskipper Species (Oxudercidae: Periophthalmus) from Kulon Progo, Special Region of Yogyakarta, Indonesia

Mudskippers are commonly cryptic species, making identification based solely on morphological characteristics challenging. This study used the DNA bar-coding method to identify mudskipper species based on the COI mitochondrial gene. The analysis revealed two distinct species, P. kalolo (20 samples) and P. argentilineatus (3 samples) with high GC contents ranging from 42.94-45.2%. The genetic divergence analysis of P. kalolo showed that they divided into two clades, while P. argentilineatus is divided into three clades with two of the clades (C and D) still conspecific groups, and those two clades with clade E exhibit a genetic distance greater than 3.5%, suggesting the presence of cryptic species. These findings provide valuable insights into the intraspecies genetic diversity of mudskippers in Indonesia, which could have essential implications for conservation efforts and highlight the potential of DNA barcoding as a powerful tool for the identification of cryptic species. Further research combining molecular and morphological identification could lead to a more comprehensive understanding of species identification and help address the challenges posed by cryptic species.

terrestrial locomotion via pectoral fins, and greater immunological protection against pathogens (You et al. 2018). The versatility of the fish is further reflected in the variety of pelvic fin morphologies, which suggest their adaptability to both terrestrial and aquatic environments (Hidayat et al. 2022). In addition to their unique adaptations, mudskipper fish have the potential as a source of food, as they are rich in carbohydrates, lipids, proteins, and well-balanced amino acids (Duggal et al. 2020;Mahadevan et al. 2021). Moreover, they play an important role as filter feeders and potential indicators of environmental contamination (Sangur et al. 2021).
The taxonomic classification of members within the family Oxudercidae (previously known as Gobiidae) is often hindered by their morphological similarities, presenting significant difficulties in species identification (Thacker 2003). This phenomenon has been observed in mudskippers of the genus Boleophthalmus and Periophthalmus, with numerous reports of cryptic species (Chen et al. 2014;Polgar et al. 2014;Arisuryanti et al. 2018;Aji & Arisuryanti 2021;Rha'ifa et al. 2021). Using traditional morphological methods alone has proven to be insufficient in accurately identifying Periophthalmus species, as evidenced by several misidentifications such as Periophthalmus walailakae with Periophthalmodon schlosseri (1), Periophthalmus takita, Jaafar & Larson, 2008with Periophthalmus novaeguineaensis Eggert, 1935, and Periophthalmus variabilis Eggert, 1935 with Periophthalmus novemradiatus Hamilton, 1822 (3) (Murdy 1989;Jaafar et al. 2006;Jafaar & Larson 2008;Jaafar et al. 2009). To overcome these limitations, the incorporation of molecular approaches, such as DNA barcoding, into the taxonomic classification process is crucial. DNA barcoding is a fast and accurate technique that utilizes short DNA sequences to identify species. The cytochrome oxidase subunit I (COI) mitochondrial gene is a commonly used marker in animal DNA barcoding due to its rapid evolutionary rate, lack of introns, and high copy number, making it an effective tool for identifying numerous animal taxa and cryptic species (Hebert et al. 2003b;Imtiaz et al. 2017;Roesma et al. 2020).
Pasir Mendit Beach in Kulon Progo, Special Region of Yogyakarta, Indonesia is a unique and valuable ecosystem, encompassing both a tourism and conservation area with dense mangroves and muddy terrain. It is located near the coastline of Yogyakarta, close to Muara Bogowonto and Baros Beach, the habitats for several species of mudskipper (Arisuryanti et al. 2018;Aji & Arisuryanti 2021). Despite its significance, the genetic information of the mudskipper populations in Pasir Mendit Beach is yet to be explored. Therefore, this study aimed to identify mudskippers from Pasir Mendit Beach in Kulon Progo, Special Region of Yogyakarta, Indonesia using the mitochondrial gene of COI as a DNA barcoding marker.

MATERIALS AND METHODS Sample Collection
Twenty-three samples of mudskippers (code MSP) were collected from Pasir Mendit Beach, Special Region of Yogyakarta, using a hand net in November 2019 (Figure 1). The collected specimens were thoroughly cleaned and documented (as shown in Figure 2). The samples were placed in a ziplock bag and stored in an icebox, and transported to the Laboratory of Genetics and Breeding at the Faculty of Biology, Universitas Gadjah Mada. The samples were preserved in 99% ethanol and stored at -20°C until utilized for DNA extraction. Given the difficulties in accurately identifying mudskipper species based on morphological characteristics, this study exclusively relied on molecular identification methods. Specifically, Nucleotide BLAST analysis (https://blast.ncbi.nlm.nih.gov) and BOLD (https://www.boldsystems.org/index.php/ IDS_OpenIdEngine) of the COI gene sequences for each sample code to verify the species of mudskipper.

DNA Extraction
The complete DNA genomic samples were isolated from muscle tissue on the pectoral fins above the back area of the gills using the kit (DNeasy Blood and Tissue kit from QIAGEN, Valencia, CA, USA) following the manufacturer's protocols.

DNA Amplification, Electrophoresis and Sequencing
A partial COI mitochondrial gene was amplified using a pair of primers, FishF2 (5'-TCGACTAATCATAAAGATATCGGCAC-3') and FishR2 (5'-ACTTCAGGGTGACC GAAGAATCAGAA-3') (Ward et al. 2005), in a Biorad T100 Thermal-Cycler. The amplification reaction was carried out in a 50 µL volume, containing 10-100 ng of genomic DNA, 25 µL of MyTaq HS Red Mix kit (Bioline), 1 mM of MgCl 2 , 0.6 µM of each primer (FishR2 and FishR2), and 11 µL of ddH 2 O. To test the effectiveness of DNA amplification, a negative control was added by removing the DNA template from the reaction mixture. The PCR was set up in 35 cycles of reaction for 2 minutes of pre-denaturation at 95°C, 15 seconds of denaturation at 95°C, 30 seconds of annealing at 50°C, 30 seconds of extension at 72°C and 5 minutes of final extension at 72°C (Arisuryanti et al. 2020).
The PCR products were run on a 1% agarose gel electrophoresis at 100 volts for 25 minutes, with Florosafe (Bioline) for staining and Trisacetate EDTA (TAE) 1X as buffer. The visualization was carried out using UV light. Each amplification sample was delivered to First Base (Malaysia) via PT Genetika Science in Jakarta for purification and sequencing in forward and reverse directions utilizing the Big Dye Terminator (Applied Biosystems) and the Genetic Analyzer ABI 3730xl.

Sequence Editing
The consensus sequence was constructed and edited using GeneStudio software and validated using the SeqMan and EditSeq menus in the DNASTAR software (DNASTAR Inc., Madison, USA). During this process, the chromatograms were carefully inspected for any ambiguous bases or stop codons.

Sequence Identification
The sequence of mudskipper species was identified using the Nucleotide BLAST analysis available on the NCBI website (https:// blast.ncbi.nlm.nih.gov/Blast.cgi) and the BOLD Identification Engine menu on the BOLD website (https://v3.boldsystems.org/index.php/ IDS_OpenIdEngine). The similarity of the sample and the query cover were used to estimate the species of the mudskipper fish compared to the data available in GenBank and BOLD.

COI Sequence Alignment
The mudskipper COI sequences were aligned using the Mesquite v.3.70 software on the Opal menu (Maddison & Maddison 2021) and the MEGAX software on the ClustalW menu (Kumar et al. 2018).

Intrapopulation and Intraspecies Analysis
The focus of the intrapopulation analysis was to determine the nucleotide composition of each species of mudskipper within the population. In addition, the intraspecies analysis expands the investigation by comparing the nucleotide composition, genetic distance, and phylogenetic relationship across different populations.

Nucleotide Composition & Genetic Distance
The MEGAX program was used to compute the COI nucleotide composition. The genetic distance was examined using the Kimura-2 parameter model and summarized in a Neighbor-Joining (NJ) tree, which is a commonly used methodology in DNA barcoding research (Hebert et al. 2003a)

Phylogenetic Relationship
The MEGAX software (Kumar et al. 2018) was used to assess the phylogenetic tree reconstruction using the Neighbor-Joining and Maximum Likelihood methods with 1,000 bootstraps, as well as the BEAST program for Bayesian Inference method (Suchard et al. 2018). The Bayesian Information Criterion (BIC) developed in jModelTest 2.1.10 was used to find the optimal evolutionary model (Darriba et al. 2012). The best sequence substitution model on the Bayesian Inference Criterion is HKY with gamma (HKY + G). The posterior probabilities distribution was estimated using the method of Markov chain Monte Carlo (MCMC) over 10 7 generations using a 1,000-generation sampling frequency.

RESULTS AND DISCUSSION Result Species Verification
We successfully amplified the COI mitochondrial gene from twenty-three mudskipper fish samples collected from Pasir Mendit Beach, yielding a fragment length of 700 bp (Figure 3). Sequence editing using GeneStudio software resulted in consensus sequences ranging from 675 to 702 bp with translated amino acids numbering from 225 to 234. Our analysis of the genetic similarity of the samples to reference species in GenBank and BOLD showed that 20 out of 23 samples had a similarity range of 99.27 to 100% to P. kalolo, while the remaining three samples had a similarity range of 99.69 to 99.85% to P. argentilineatus. The genetic similarity of the samples to their conspecific references is presented in Table 1. The COI sequences of P. argentilineatus has been deposited in GenBank under accession number OQ804359-OQ804361 whereas the COI sequences of P. kalolo has been deposited in GenBank under accession number OQ804373-OQ804392.
We performed alignment of the COI mitochondrial gene of P. kalolo and P. argentilineatus collected from Pasir Mendit Beach which produced consensus sequences of 666 bp for each species. These consensus se-quences were then used to perform intrapopulation analysis for each species. In addition, we also conducted intraspecies analysis, using COI sequences of mudskippers from Pasir Mendit Beach and other Indonesian regions recorded in GenBank and BOLD, yielding fragment lengths of 624 bp for P. kalolo and 579 bp for P. argentilineatus. To understand the evolutionary relationships of the two species, we aligned the sequences from Pasir Mendit Beach with those from other locations in Indonesia documented in GenBank, as well as with the COI sequence of Boleophthalmus boddarti (accession number KU692377) as an outgroup. The resulting aligned sequences (579 bp) were then used to construct a phylogenetic tree.
Our intraspecies analysis revealed variations in the composition of nucleotides T, C, A, and G within the P. argentilineatus. These variations ranged from 0-1.21% for T, 0.17-1.21% for C, 0-0.52% for A, and 0-0.52% for G. Interestingly, P. argentilineatus displayed a higher A+T composition than G+C, with A+T ranging from 56.99 to 57.86% and G+C ranging from 42.14 to 43.01%. The highest A+T composition was observed in samples collected from Pasir Mendit Beach (MSP-03), Baros Beach (MZ606685), and Bogowonto (MT439598), while the highest G+C composition was recorded in samples from Baros Beach (MZ606680, MZ606682). In addition, we also identified the highest number of similarities in nucleotide composition among specimens sourced from various locations, such as Pandeglang (KU692745, KU692746), Bogowonto (MT439599), Baros Beach (MZ606684, MZ606686).
The resulting phylogenetic tree revealed that P. kalolo species was divided into two separate clades, supported by high bootstrap values of 100% (NJ), 99% (ML), and a posterior probability of 1 (BI). Meanwhile, P. argentilineatus was separated into three clades with values of bootstrap 99% (NJ), 95% (ML), and 1 in the posterior probability (BI) value for the separation of clade C and clade D. The bootstrap values between clade D and clade E were 100% (NJ), 99% (ML), and 1 (BI) for posterior probability. The phylogenetic tree reconstruction is presented in Figure 4.
Our genetic distance analysis of P. kalolo samples from Pasir Mendit Beach compared to other P. kalolo documented in BOLD and GenBank from other Indonesian regions yielded an average genetic distance of 1.90% (ranging from 1.30 to 2.63%) between Clade A and Clade B ( 2.84%), while the highest was between clade E and clade C at 6.01% (5.62 -6.38%) ( Table 3). The average genetic distance between clade E and clade D was 5.53% (5.24-5.81%).

Discussion
DNA barcoding has emerged as a robust tool in the field of taxonomy, enabling the rapid and accurate identification of species. By combining genetic distance and phylogenetic tree approaches, this method has been extensively validated across a wide range of taxa (Ude et al. 2020;Wang et al. 2021;Sachithanandam et al. 2022). The present study applied DNA barcoding to identify of mudskipper species, utilizing the COI mtDNA gene. We compared our sequences with those from NCBI GenBank and BOLD database and discovered that the 23 mudskipper samples in this study belonged to two distinct species, P. kalolo (20 samples) and P. argentilineatus (3 samples). According to Triandiza and Maduppa (2018), the most similar GenBank sequence is characterized by the same maximum and total score, query coverage close to 100%, E value close to 0, and percentage of identity or similarity close to 100%. Based on that criterion, all mudskipper samples in this study showed greater similarity to the sequence of both species in GenBank.
We observed that the nucleotide composition analysis of the COI gene in both species of mudskippers revealed GC contents ranging from 42.94-45.2%. This result is consistent with the study by Ward et al. (2005) which discovered a total GC content COI ranging from 42.2-47.1%. Additionally, we found a T>C>A>G nucleotide composition pattern of the COI gene that was comparable to studies in other fish families (Bingpeng et al. 2018;Kombong & Arisuryanti 2018;Wu et al. 2018;Linh et al. 2018;Aji & Arisuryanti 2021).
Our phylogenetic tree reconstruction, based on mitochondrial COI sequences from various locations, including Pasir Mendit Beach, revealed the existence of two distinct clades in P. kalolo and three distinct clades within P. argentilineatus. We found that P. kalolo samples from Bogowonto were in the same clade as samples from Baros Beach  Arisuryanti et al. (2018) and Aji and Arisuryanti (2021). Additionally, we also observed that P. argentilineatus from Baros Beach (MZ606684-86) were in the same clade as samples from Pandeglang (KU692744-46) and samples from Bogowonto (MT4359598-600), separated from the clade with P. argentilineatus from Baros Beach (MZ606687) and a sample from Tukad Bilukpoh, Jembrana Bali (KU692743) (Dahruddin et al. 2017;Aji & Arisuryanti 2021). We evaluated the robustness of our phylogenetic tree by assessing bootstrap and posterior probability values. We found that the formation of the clades in both P. kalolo and P. argentilineatus had high values of bootstrap and posterior probability, demonstrating the robustness of these clades. According to Apriliyanto and Sembiring (2016), if the bootstrap value in methods of Neighbor-Joining and Maximum Likelihood is closer to 100% and the posterior probability value in the Bayesian Inference method is close to 1, then the formation of the clades is considered robust. Our results meet this criterion, further supporting the reliability and accuracy of our phylogenetic tree reconstruction.
Our analysis revealed a genetic divergence of 1.90% (ranging from 1.30 to 2.63%) in two distinct clades of P. kalolo. This value is fall within the intraspecies genetic distance threshold of 3.5% proposed by Zemlak et al. (2009), indicating that P. kalolo samples from different Indonesian locations belong to the same species (conspecific). This finding aligns with previous studies by Arisuryanti et al. (2018) and Aji and Arisuryanti (2021), which also found no significant genetic differentiation among P. kalolo populations in Indonesia. In contrast, the genetic distance analysis of P. argentilineatus revealed a more complex picture. Our finding indicates that P. argentilineatus is divided into two conspecific groups, clades C and D, while one group (clades E) exhibited a genetic distance greater than 3.5%. This result is consistent with previous studies (Arisuryanti et al. 2018;Rha'ifa et al. 2021;Aji & Arisuryanti 2021) that have suggested the presence of cryptic species in P. argentilineatus in Indonesia. These results support the clade formation observed in the phylogenetic tree.
Interestingly, we did not observe a clear clades separation by geographical region in the phylogenetic tree, suggesting a lack of geographical structuring in the populations of P. kalolo and P. argentilineatus. This lack of geographical structure was particularly evident in clade E, where the sample MSP-03 from Pasir Mendit Beach was grouped with a sample from Baros Beach (MZ606687) and a sample from Tukad Bilukpoh, Jembrana, Bali (KU692743). The distance between Pasir Mendit Beach to Tukad Bilukpoh is more than 400 km, and we assume that there should be little or no gene flow between them. We hypothesize that this result could be attributed to geological events that allowed for the dispersion and gene flow between the remote populations, such as variations in sea level or oceanic currents (Miller et al. 2005;García-De León et al. 2018). Additionally, it is important to note that the lack of clear separation by geographical region in the phylogenetic tree may be influenced by the small sample size of some populations. Future studies with larger sample sizes may be needed to better understand the patterns of genetic diversity and geographic structuring in these two species.

CONCLUSIONS
Our study demonstrates the effectiveness of DNA barcoding in identifying and differentiating mudskipper species. Through analysis of the COI mtDNA gene, we were able to identify two distinct species of mudskippers, P. kalolo and P. argentilineatus, and assess their genetic diversity and relationships. The high GC content and nucleotide composition pattern observed in the COI gene were consistent with other fish families. Our phylogenetic tree analysis revealed the existence of two distinct clades in P. kalolo and three in P. argentilineatus, with the robust formation of these clades. Genetic divergence analysis indicated that P. kalolo samples from different Indonesian locations are conspecific, while P. argentilineatus is divided into two conspecific groups (clades C and D) and one group with greater genetic distance (clade E), suggesting the presence of cryptic species. These findings provide valuable insights into the evolutionary history of mudskippers and highlight the potential of DNA barcoding as a powerful tool for species identification and conservation.

AUTHOR CONTRIBUTION
T.A. was in charge of planning, designing, supervising the entire research process and reviewed the initial manuscript and revised the final manuscript. D.F. worked in a laboratory (samples collections, DNA extraction, amplification using a PCR technique, electrophoresis using agarose gel, analysis of data, and writing of the manuscript). KWA and DSP guide the laboratory work, data analysis and writing the manuscript.

ACKNOWLEDGMENTS
We would like to thank the Head of the Laboratory of Genetics and