DNA Barcode of Barred Mudskipper (Periophthalmus argentilineatus Valenciennes, 1837) from Tekolok Estuary (West Nusa Tenggara, Indonesia) and Their Phylogenetic Relationship with Other Indonesian Barred Mudskippers

Barred mudskipper (Periophthalmus argentilineatus) has a potency to be developed as protein for human consumption and ornamental fish. The fish also has an important role in mangrove ecosystems. Nevertheless, many barred mudskippers have been considered a cryptic species. Therefore, accurate identification is needed to clarify species identification of the barred mudskipper using DNA barcoding.  This research aimed to identify barred mudskippers from Tekolok Estuary (East Lombok, West Nusa Tenggara, Indonesia) using COI mitochondrial gene as a DNA barcode and analyze genetic relationship with other barred mudskippers from several regions of Indonesia recorded in GenBank. This study used a PCR method with universal primers FishF2 and FishR2.  The data was then analysed using DNASTAR, BLAST, Mesquite, MEGA, DnaSP, BEAST, GenAlEx, and NETWORK. The results revealed that barred mudskipper from Tekolok Estuary has been verified as Periophthalmus argentilineatus. The results also exhibited that P. argentilineatus from Tekolok Estuary has a close genetic relationship to P.argentilineatus from Tukad Bilukpoh (Jembrana, Bali).  In addition, phylogenetic analysis showed that P.argentilineatus from Indonesia consisted of two clades with a genetic distance of approximately 6.64%. This analysis revealed evidence of the cryptic diversity of P.argentilineatus from Indonesia. Further detailed studies are needed to clarify whether Indonesian P.argentilineatus should be categorized into more than one species or single species with several subspecies.


INTRODUCTION
Barred mudskipper (Periophthalmus argentilineatus) belongs to suborder Gobioidei, family Gobiidae and subfamily Oxudercinae. Family Gobiidae consists of 130 genera with 1,866 species (Nelson 2016;Eschmeyer 2019). One of the most widely distributed of this family in the Indo-Pacific region is the genus Periophthalmus and according to Froese and Pauly (2004), the genus can be bent to function like arms that can be used to crawl or jump on mud and perch on mangrove roots (Elviana & Sunarni 2018). Barred mudskipper spends mostly outside of the water and prefers muddy habitat especially in mangrove areas (Vanhove et al. 2012). The mangrove habitat is a suitable habitat for barred mudskipper survival life due to lots of detritus which is used for the fish food. In addition, barred mudskipper is also a carnivorous opportunist feeder eating small crabs, small fish, shrimps, and other small arthropods that are commonly found in mangrove habitat. Barred mudskipper can also inhabit extreme conditions (Lee et al. 2005).
Previous studies have reported a common phenomenon of cryptic species in the members of the suborder Gobioidei including barred mudskippers which have a high similar morphology but genetically different. As a consequence of this phenomenon is taxonomic confusion. This is due to inaccurate identification of a species indicated by a cryptic species may lead to erroneous conclusions in naming species (Leys et al. 2016). This phenomenon can cause synonymous problems where there is a double name in the same species or cause an increase in cryptic diversity defined as two or more different species and classified as a single species due to morphological similarities (Tronteli & Fisher 2009). Double species name due to misidentification will have an impact on conservation efforts and sustainable use.
The discovery of cryptic species was reported from members of family Rhyacichthyidae, Odontobutidae, Xenisthmidae, Eleotridae Gobiidae, Kraemeriidae, Schindleriidae, Microdesmidae, and Ptereleotridae using mitochondrial genes: ND1, ND2 and COI (Thacker 2003;Thacker & Hardman 2005), genus Schindleria (family Schindleriidae) using 16S mitochondrial gene (Kon et al. 2007), genus Economidichthys, genus Knipowitschia, and genus Pomatoschistus using mitochondrial genes 12S and 16S (Vanhove et al. 2012). Comprehensive research was conducted by Agorreta et al. (2013) in members of the suborder Gobioidae using four mitochondrial genes (12S, 16S, COI, Cyt b) and two nuclear genes (Zicf and RAG-1), but this study did not include members of the genus Periophthalmus. The cryptic diversity of P. argentilineatus and P. kalolo was also investigated by Polgar et al. (2014) using two mitochondrial genes namely D-Loop and 16S as well as one nuclear gene namely RAG-1. Research on the molecular identification of mudskippers belonging to the genus Periophthalmus using the COI mitochondrial gene as a DNA barcoding marker in Indonesia is still very limited. Dahruddin et al. (2017) identified three species of members of the genus Periophthalmus namely P. novemradiatus, P.kalolo, and P. argentilineatus in Java (Pandeglang, Cilacap, Tulung Agung) and Bali (Jembrana). Furthermore, Arisuryanti et al. (2018) successfully identified two species namely P. argentilineatus and P. kalolo in Bogowonto Lagoon (Yogyakarta). This study also revealed that P. argentilineatus from Bogowonto Lagoon has a close genetic relationship with P. argentilineatus from Pandeglang. However, no genetic information of barred mudskipper has been conducted from Tekolok Estuary (West Nusa Tenggara). Therefore, this study used COI mitochondrial gene as a DNA barcoding marker to identify barred mudskippers from Tekolok Estuary and analyze their genetic relationship with other Indonesian barred mudskippers recorded in GenBank. The use of COI mtDNA as a DNA barcoding marker for a molecular-based approach to the taxonomic investigation such as identification of closely related species, revealing cryptic speciation and phylogeographic structures within a species has been conducted for over recent years (Healey et al. 2018;Shelley et al. 2018;Rathipriya et al. 2019). This marker is a molecular standard approach developed by Hebert et al. (2003) which is applied to identify species and support species delimitations.

MATERIALS AND METHODS
Sample collection, DNA extraction, amplification, and sequencing of COI mitochondrial gene Ten barred mudskippers (code MSL) were collected from Tekolok Estuary, East Lombok, West Nusa Tenggara, Indonesia (8 o 20'30.0"S 116 o 42'31.0"E) ( Figure 1). Before dissecting the tissue sample, the barred mudskipper was documented ( Figure 2). The tissue sample contained 50-100 mg of muscle tissue. Each tissue sample was dissected with a sterilized surgical scissor, placed into a 1.5 ml tube, and preserved in 99% ethanol in the field, and stored at -20°C in the laboratory for further analysis.  Total Genomic DNA was extracted from muscle tissue of each specimen using DNeasy blood and tissue kit (QIAGEN, Valencia, CA, USA) according to the manufacture's protocols. The partial COI mitochondrial gene was amplified using primers FishF2 (5'-TCGACTAATCATAAAGATATCGGCAC-3') and FishR2 (5'-ACTTCAGGGTGACC GAAGAATCAGAA-3') (Ward et al. 2005). The MyTaq HS Red Mix PCR kit (Bioline) was used for the polymerase chain reaction (PCR) and PCR amplifications were conducted in 50µL reaction volumes containing 10-100 ng of genomic DNA, 25µL MyTaq HS Red Mix PCR, 1mM MgCl 2 , 0.6µM of each primer and 11µL double distilled water (ddH2O). Negative control was set up by omitting template DNA from the reaction mixture to evaluate the reliability of the DNA amplification. The amplification process was conducted using the following cycling conditions as initially conducted by Arisuryanti et al. (2019). The PCR products were then run on 1% agarose gel buffered with Tris-acetate EDTA (TAE), stained with Gel Red (Bioline), and visualized under UV light. All amplification samples were then sent to First Base Sdn Bhd (Malaysia) through P.T. Genetika Science (Jakarta) for purification and sequencing in both forward and reverse directions using the Big Dye Terminator (Applied Biosystems) and the ABI 3730xl Genetic Analyzer (Applied Biosystems).

Genetic Analysis
Data obtained from DNA sequencing results were edited with SeqMan and EditSeq on the DNASTAR program (DNASTAR Inc. Madison, USA). For each individual, sequencing reactions were performed using both forward and reverse primers. Chromatograms were inspected manually to check ambiguous bases and stop codons. Next, the COI sequence data is converted to FASTA format through the Mesquite v.3.5 program (Maddison & Maddison 2018). COI sequences of barred mudskipper (P. argentilineatus) were then aligned using Opal in Mesquite v.3.5 (Maddison & Maddison 2018) and ClustalW on the MEGAX program (Kumar et al. 2018). In this study, 17 COI mitochondrial sequence data of P. argentilineatus (10 samples from Tekolok Estuary and 7 samples from several regions in Indonesia recorded in GenBank with accession number KU692743, KU692744, KU692745, KU692746, MT439598, MT439599, MT439600) were used to analyze the intraspecific genetic relationship and their phylogenetic relationship. The composition of the COI nucleotides was calculated through the EditSeq program in the DNA Statistics menu using DNASTAR software. The composition of C+G is validated using the DnaSP program v. 6 (Rozas et al. 2017). Intraspecific diversity was estimated as the number of haplotypes, the number of polymorphic sites, haplotype diversity, and nucleotide diversity using the software DnaSP v.6 (Rozas et al. 2017). Next, intraspecific genetic distance was analysed using the Kimura-2 Parameter (K2P) model and was summarised in a Neighbour-Joining (NJ) tree as this is the standard methodology used in barcoding studies with bootstrap 1,000 replicates (Hebert et al. 2003). Phylogenetic relationships were estimated using a Bayesian approach. The Akaike Information Criterion (AIC) implemented in jModelTest ver 2.1.10 (Darriba et al. 2012) was used to determine the best fit evolutionary model. Bayesian inference (BI) was performed with BEAST ver 1.10 (Suchard et al. 2018) under the best-fit model. Two simultaneous Markov chain analyses (MCMC) were run for 10 6 generations to estimate the posterior probabilities distribution with a sampling frequency set to every 1,000 generations. The analysis was done until the standard deviation of split frequencies was below 0.01. The analysis used a relative burn-in of 25% for diagnostics. Consensus trees were visualised in FigTree 1.4.4 (Rambaut 2019). COI sequences of eight other species included in the genus Periophthalmus (accession number KU692751, KU692755, KU692756, KU692759, KP966189, KP966192, KT368097, and KT368106) and Pseudapocryptes lanceolatus (JX260938) were included in this analysis as an outgroup and to construct phylogeny tree. Patterns of divergence between haplotypes were analyzed using Principal Coordinate Analysis (PCoA) based on the genetic distance of COI mitochondrial gene conducted in GenAlEx version 6.51 (Peakall & Smouse 2012) and NETWORK version 10.1 (https://www.fluxus-engineering.com).

PCR amplification result and DNA barcode of barred mudskipper from Tekolok Estuary
The results revealed that the amplification of COI mitochondrial gene from all 10 barred mudskipper fish from Tekolok Estuary produced fragment lengths between 600 and 700 bp (Figure 3). After editing on chromatograms, the resulting consensus sequence was 678 bp which can be translated into 226 amino acids. Furthermore, the results of the analysis using Nucleotide BLAST, 10 barred mudskipper samples from Tekolok Estuary investigated in this study had a similarity of 99.23-100% with Periophthalmus argentilineatus with accession number KU692743. Therefore, all of the ten barred mudskippers from Tekolok Estuary were verified as Periophthalmus argentilineatus. All of the COI mitochondrial genes of P. argentilineatus from Tekolok Estuary have been submitted to GenBank with accession number MW514015-MW514024.

Composition of mtDNA COI nucleotide
In this study, 17 samples of P. argentilineatus consisted of 10 samples of P. argentilineatus from Tekolok Estuary and 7 samples of P. argentilineatus from several regions in Indonesia recorded in GenBank were analyzed. After alignment on 17 samples of P. argentilineatus and obtained fragment length 519 bp. We can arrange the composition of mitochondrial COI nucleotide (Table 1). No similar value can be seen from all of the samples. From Table  1, the differentiation of the composition of nucleotide T, C, A, and G were 0.0-0.20%, 0-0.46%, 0-0.20%, and 0-39% respectively. In addition, P. argentilineatus from Tekolok Estuary had the highest frequency of nucleotide C (25.70%) while P. argentilineatus from Pandeglang (KU692745 and KU692746) had similar nucleotide T, A, C, G with P. argentilineatus from Bogowonto Lagoon (MT439599 and MT439600). The result revealed the divergence of T, C, A, and G nucleotides in COI mitochondrial genes among P. argentilineatus samples analyzed in this study. This data revealed polymorphism in the COI mitochondrial gene of P. argentilineatus in Indonesia which indicated the intraspecific genetic variation of P. argentilineatus.

COI mitochondrial sequence variation and phylogenetic relationship
Seventeen COI sequences that were obtained from P. argentilineatus individuals from Tekolok Estuary and other Indonesia regions recorded at the GenBank database (Table 2) represented 10 distinct haplotypes. The 519bp alignment contained 42 polymorphic sites (8.09%) with 34 being parsimony-informative ( Figure 4) and no gaps and nonsense codon was found among the 17 sequences after the COI sequences were translated using the vertebrae mitochondrial code. From 42 polymorphic sites, 35 sites were detected transitions whereas 7 sites were transversions (Figure 4).
The optimal model of nucleotide sequence substitution for the COI gene including the outgroup samples was the HKY model with the gammadistributed rate as inferred by jModelTest2 under the Akaike information criterion (AIC). The Bayesian analysis of COI mitochondrial data, together with additional sequences of eight species that belong to the genus Periophthalmus and Pseudapocryptes lanceolatus taken from GenBank, revealed that the barred mudskippers from Indonesia fall into two distinct clades (Clade A and B), with the associated nodes supported by posterior probability 1.00 ( Figure 5). Similarly, the neighbour-joining (NJ) tree exhibited two distinct clades (Clade A and B) with a bootstrap value of 99% ( Figure 6). The two clades or OTUs generally showed strong geographic patterns of differentiation. The clade A consisted P. argentilineatus from eastern part of Indonesian region (Tekolok Estuary, West Nusa Tenggara and Tukad Bilukpoh, Jembrana, Bali) while the clade B comprised P. argentilineatus from western (Pandeglang) and central part of Indonesian region (Bogowonto Lagoon, Yogyakarta). The geographic distribution within clades showed considerable overlap between haplotypes such as haplotype A 1 and B 2 (Figure 7).
The genetic distance based on Kimura-2 Parameter (K2P) model among all COI haplotypes within P. argentilineatus was quite variable ranging from 0.0 to 8.24% (mean= 3.43%) ( Table 3). The average nucleotide sequence divergence between haplotypes within clade A and clade B is        0.402% (range 0-0-0.97%) and 0.452% (0-1.17%) respectively. The genetic divergence between clade A and B is quite high ranging from 6.09 to 8.24% (mean=6.64%). The two distinct clades of P. argentilineatus from Indonesia supported by genetic distance (6.64%) show that P. argentilineatus is a cryptic species. According to Zemlak et al. (2009), if the genetic distance between the two populations or clade exceeds 3.5% then both are considered different species. The phenomenon of cryptic species is commonly found in many freshwater invertebrates such as Penaeus monodon (Yudhistira & Arisuryanti 2019) and Etheria elliptica (Elderkin et al. 2018) as well as in freshwater fish such as Monopterus albus (Arisuryanti et al. 2016) and Cirripectes alboapicalis (Delrieu-Trottin et al. 2018). The discovery of cryptic species of P. argentilineatus in Indonesia with the use of DNA barcoding is a novel finding. Therefore, a more detailed taxonomic revision of P. argentilineatus in Indonesia needs to be justified whether it consists of different species or remains one species with two subspecies. The correct and clear taxonomic status of P. argentilineatus can be applied for the conservation of this fish species in its habitat.
Levels of differentiation within and between P. argentilineatus clades are summarised by Principal Coordinate Analysis (PCoA) in Figure 8. This analysis indicates that the degree of differentiation between clades is much greater than within clades. From Figure 8, the clades remain distinct over significant geographic space, which includes a zone of geographic overlap (Figure 7), provide support for speciation of P. argentilineatus in Indonesia.

Variation within clades
Six COI haplotypes were found in clade A from 11 individuals and four haplotypes were detected in clade B from six individuals (Table 2). Within clade A the level divergence among haplotypes was between 1-4 bp. From 519 bp, seven polymorphic sites were detected with three parsimony informative sites. The seven polymorphic sites contained six transitions and one transversion. In addition, haplotype diversity and nucleotide diversity were 0.727 and 0.0040 respectively. Next, within clade B the level divergence among haplotypes was slightly higher than clade A which was between 1-6 bp consisted of seven variables sites without parsimony informative sites. The seven polymorphic sites comprised four transitions and three transversions. The haplotype diversity and nucleotide diversity were 0.800 and 0.0045 respectively. The finding of variation within clades exhibited polymorphism of the COI mitochondrial gene. This gene polymorphism indicates the intrapopulation genetic variation of P. argentilineatus within clades.

CONCLUSION
Barred mudskippers from Tekolok Estuary are successfully identified as a Periophthalmus argentilineatus using COI mtDNA as a DNA barcoding marker. In addition, P. argentilineatus from Tekolok Estuary (West Nusa Tenggara, Indonesia) has a close genetic relationship with P. argentilineatus from Tukad Bilukpoh (Jembrana, Bali). Furthermore, DNA barcoding using COI mitochondrial gene has provided strong evidence that P. argentilineatus from Indonesia contains cryptic species. In this analysis, the Indonesian P. argentilineatus fell into two distinct groups, which are genetically distinct from each other. Further investigations are needed to clarify the two genetically distinct groups of Indonesian P. argentilineatus.

AUTHORS CONTRIBUTION
T.A. designed the research and supervised all the processes, F.A.R collected and analyzed the data, D.J.A and L.H provided and documented the samples, and helped F.A.R. to analyze the data. All authors wrote and revised the manuscript. All authors read and approved the final paper.