Comparative mitochondrial genome brings insights to slight variation in gene proportion and large intergenic spacer and phylogenetic relationship of mudskipper species

Fish mitochondrial genome have been largely studied worldwide for evolutionary and other genetic purposes and the structure and gene organization are commonly conservative. However, several studies have demonstrated that this scenario may present variations in some taxa, showing differentiation on the gene rearrangement. In this study, the complete mitogenome of terrestrial fish Boleophthalmus dussumieri was generated and compared with other species of the Exudercidae fishes. The newly complete mitogenome generated is circular and 16,685 bp of length, and it contained 13 protein-coding genes (PCGs), two ribosomal RNA (rRNAs), 22 transfer RNA genes (tRNAs), and one control region (CR), with high conservative structure, like other Mudskippers. Most of the PCG showed similar codon usage bias. The gene length was found to be different specially for the CR, 12S rRNA gene and ND5 gene in some taxon. All the Boleophthalmus species showed a gene duplication in the CR, except for B. dussumieri, and they presented a long intergenic spacer specially on the tRNA-Pro/ OH Tandem duplication/random loss (TDRL) and dimer-mitogenome and nonrandom loss (DMNL) are suitable to explain the mitogenome rearrangement observed in this study. The phylogenetic analysis well supported the monophyly of all mudskipper species and the analysis positioned the Periophthalmus clade as the most basal of the terrestrial fishes. This finding provides basis and brings insights for gene variation, gene rearrangements and replications showing evidence for variety of mitochondrial structure diversity within mudskippers.


Introduction
Mudskippers are amphibious shes with peculiar characteristics, living in muddy areas of mangroves (Gordon et al., 1969;Randall et al., 2004).The genus Boleophthalmus is known to be one of the most terrestrial among mudskippers, exhibiting locomotory, respiratory, vertebral adaptation and specializations that enable overland excursions lasting up to 14 hours (Steppan et al., 2022;Polgar et al., 2017).Despite being a crucial group of sh for evolutionary and ecological studies, only seventeen species are available in literature (Table 2), three of which belong to the genus Boleophthalmus (Zhang et al., 2016;Pan et al., 2021).
Mitogenome studies enable a wide range of investigations, including tandem repeats, phylogenetic analyses, gene rearrangements, gene overlap, analysis short or long intergenic regions (Zhang et  and organization of the mitochondrial genome in many bony shes exhibit a similar structural conformity with 37 genes, including 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs) and 2 ribosomal RNA (rRNAs) genes, and one control region (Li et al., 2018;Cao et al., 2020).However, some organisms present variation in the length of genes (resulting in long gene sequences) and sometimes repetitions may be found in portions of the genes, including in the protein-coding genes (Mar-Silva et al., 2022), or in the control region (Ding et al., 2022).
Several studies have used tools to conduct comparative analysis, aiming to understand the correlations in terms of GC contents percentages, presence of complete or short gene repetition, protein-coding and non-coding regions, as well as the consistency between different regions (Ding et al., 2022, Kabiraj et al., 2022), which provides relevant information on proportions of changes in gene rearrangements within each class or family.These studies have provided valuable insights into the proportions of gene rearrangements within each class or family (Montaña-Lozano et al., 2022;Li et al., 2022).
Variations in mitogenome length and genetic organization are often associated with various evolutionary events, such as nucleotide insertions and deletions in the control region (Lee et al., 1995;Sbisà et al., 1997), as well as duplications followed by deletion in certain regions, particularly on tRNAs (Formenti et al., 2020).This features directly impact functional differences, considering that the distinct gene structure in uences gene function (Gissi et al. 2008).
The analysis of complete mitogenomes offers signi cant advantages in studying vertebrates, and researchers can utilize the available data in public databases, such as the NCBI, to address phylogenetic questions (Montaña-Lozano et al., 2022; Li et al., 2022).
Mudskippers exhibit variations in individual characteristics, encompassing both genetic variations and environmental adaptations (Sakamoto et al., 2002;Tan et al., 2020;Theeranukul et al., 2021).Therefore, through research utilizing these tools, various mechanisms and genetic information can still be uncovered.Studies already carried out with mudskipper mitogenomes generally describe their complete mitogenomes, except for some studies that, in addition to studying phylogenetic relationships, they bring insights on the genetic functioning of the group (Liu et al., 2012;Zhang et al., 2016 ;Qiu et al., 2017;Tan et al., 2020).
In this study, we present for the rst time, the complete mitochondrial genome of the B. dussumieri (Valenciennes 1837), a mudskipper species.The primary objective of the study was rst to describe the complete mitochondrial genome and present the comparative analysis with available mudskipper mitogenomes.Additionally, we aimed to reconstruct the phylogenetic relationships among the seventeen mitogenomes of individuals within the family, based on the thirteen protein-coding genes.Furthermore, we implemented a repeat sequence analysis in the control region to assess the potential presence of gene repetitions as well as duplication in CR and certain gene/regions.

Mitogenome organization, composition, and skewness
The complete mitochondrial genome of B. dussumieri is (GenBank accession XX) 16,685pb of length, which is similar to other mudskippers mitogenomes (16,470 − 17,243pb) (Table 1 and Table 2).The mitochondrial genome and the structure were also typical of the mudskipper species and with highly conservative sites, comprising 37 mitochondrial genes (13 PCGs, 22 tRNAs, and 2 rRNAs) and one control region (CR) (Figs. 1 and 2).The mitogenome composition and features are presented in Table 1 and Table 2, highlighting the characteristics of the mitogenome, including gene content and speci c regions.The base composition of the complete mitochondrial genome of B. dussumieri was determined to be A = 29,68%; C = 27,90%; T = 26,87% and G = 15,53%, respectively.In the mitogenome composition, we observed the GC and AT content, as well as the GC skew + and GC skew region.The overall GC and A + T content for the entire genome were calculated to be 43.43% and 56.54%, respectively, which were similar to the concatenated composition of the PCGs, tRNA and rRNA.Furthermore, the nucleotide composition of all mitochondrial genes exhibited a bias, with A + T content notably higher (56.54%) than the GC content (43.43%).This pattern of nucleotide composition is consistent with other mudskippers mitogenomes, and no signi cant differences were observed among them (Table 2).Furthermore, to access weather there were biased or not in the nucleotide composition, the GC and AT skew were measured for the whole mitogenome, as well as for the PCGs, ribosomal genes, and the control region.
The GC and AT skew trends of all analyzed genes and control region are presented in (Table 2).Generally, most of the sh mitogenome present a clear bias in the nucleotide composition, which is the case of this study, where the GC skew was notably negative as most of the teleost sh (Sharma et  The GC and AT skew followed the same conventional preference of most mitogenomes, with AT skew positive, except for PCGs, and GC skew negative, except for tRNAs and CR, that varies from negative to positive in some mitogenomes.The AT skews of whole mitogenome ranged from − 0.0 (Periophthalmus novaeguineaensis) to 0.065 (Boleophthalmus.sp.JZ-2015), rRNAs 0.22 (P.minutus) to 0.254 (Boleophthalmus pectinirostris) whereas AT skew of the tRNAs range from 0.009 (Periophthalmus novemradiatus) to 0.039 (Scartelaos gigas), respectively (Fig. 3).This preferences in the sequenced mitogenomes are like other teleost mitogenomes (Zhang et al.,2016).

PCGs, Amino acid composition and codon usage
The PCG size of each mudskipper mitogenome ranged from 11,399 pb (B.dussumieri) to 11.713pb (P.novaeguineaensis) (Supplementary Table S1).We compared seventeen mudskippers mitogenomes, and among these almost all PCG encoded in the H-strand, except for ND6 that encoded in the L-strand, as typical of most vertebrate mitogenomes (Zhang et al., 2021).Overall, all PCGs showed the same con guration with ATG as start codon with exception of the COI gene that started with GTG codon.The preferences for codon termination of ve PCGs (NADH-1, NADH-2, COI, NADH-4L, and NADH-5) were TAA and three (ATP8, NADH-3 and NADH-6) terminated with TAG (Table 2).Some protein coding genes had an incomplete start and stop codon as in most sh mitogenomes (Zhao et al., 2015;Sun et al., 2021).For example, NADH-4, COII, CYT B gene, had an incomplete nucleotide T-whereas ATP6 and COIII had incomplete TA-codon as stop codon (Table 1).These features are in accordance with most mudskipper shes (Jin eat al., 2012).However, we found one unusual start and stop codon for most shes (CCT/AAC) in ATP8 gene of Periophthalmus novemradiatus species.The CCT codon, codes the Proaline amino acid whereas the AAC codes Asparagine amino acid.
Curiously, P. novaeguineaensis presented very peculiar codon preferences, totally different from others with only ten (ATG) as start codon, common GTG for COI, and three non-common codons (ATA for ATP6 a, TTA for ATP6 b, TTC and ACT for NADH5-0 and NADH5-1, respectively).In contrast, the codon termination was also different from most mudskippers, presenting only four genes with TAA, as stop codon in the (NADH-2, COI, NADH-4L and NADH-6) genes and six different stop codons, AGG for (NADH-1), a common TAG stop codon for (ATP8 and NADH-3), ACC for ATP 6 a, GTT for ATP 6 b, CTT and CTC for NADH5-0 and NADH5-1, respectively.
The relative Synonyms codon usage (RSCU) is shown in (Fig. 4).The use of RSCU was biased for both two and six-fold degenerate codons.In overall, 5,556 codons were analysed with exclusion of stop codons for B. dussumieri mitochondrial genome.

Ribosomal and Transfer RNAs genes
Gene size showed similar values among all mudskippers.The 12S rRNAs gene presented a size of 948 pb in a group that presented values between 947 and 956 pb, while for the 16S rRNAs gene, the size of the genes varied from 1,683pb belonging to (Parapocryptes serperaster) to 1,699 pb (B.dussumieri) (Supplementary Table S1).In general, the tRNAs genes were very conserved compared to some sh species (Zhao et al., 2015;Ruan et al., 2020).We found 14 tRNAs out of 22 encoding the H-strand and eight (trnA, trnC, trnE, trnN, trnP, trnQ, trnS2, and trnY) in L-strand.
The structure of almost all tRNAs were typical clover structure.Although most of them presented a secondary structure similar to other species with four arms and a central loop.Some of those tRNAs (trnM, trnS1, trnE, trnH, trnW, trnT, trnN, trnF, trnV, trnL) presented a variation in certain arms with presence of small loops (Fig. 6).The entire length of tRNA of the B. dussumieri was 1,556 pb and the length of 22 individual tRNAs gene ranged from 65 to 76 bp.This value was between the length range of the other mudskippers analysed.(Supplementary Table S1).

Gene Rearrangement, Repeated Region, and Concerted Evolution
The organization of almost all vertebrate animals generally have their mitogenome structurally conserved with one set gene copy, no intron or long intergenic spacer and single CR (Kurabayashi et al., 2013;Wolstenholme, 1992).Nevertheless, gene rearrangement in mitogenome of vertebrates are deeply explained using different models including shu ing, translocation, and inversion (Gong et al., 2013), but most common is duplication and deletion model approaches (Kong et al., 2009;Gong et al., 2020).The CR of this study followed typical rearrangement located between the tRNA-Pro and tRNA-Phe gene, as most vertebrates mitogenomes (Chandhini et al., 2021;Zhang et al., 2021).
We analysed seventeen mudskippers' species and nine of them presented duplication in their mitogenome specially in the control region.The duplication was found in three different regions, most of them, in a non-coding gene, the control region (Boleophtalmus boddarti, B. pectinirostris, Boleophthalmus sp JZ-2015, Oxuderces dentatus, S. gigas); three in the PCGs, ATP6, NADH-5 gene (P.novaeguineaensis), ATP8 (P.novemradiatus and P. minutus); and one ribosomal gene, the 12S rRNA gene duplication in the Periophthalmodon schlosseri.Our study shows two species with the OH region duplicated, P. serperaster with four OH region, with the sequence length ranging from 78 bp to 320 bp, and S. gigas with one duplication involving 31 bp and 316 bp, respectively (Table 4 and Fig. 3).Furthermore, our results also presented one species (P.serperaster) with both CR duplication and one coding-region duplication, the NADH-1 gene.

Gene duplication in the Boleophthalmus genus
The genus Boleopthalmus was represented by four species, B. dussumieri (this study), B. boddarti, Boleopthalmus pectinirostris and Boleophthalmus.sp.The rst was considered the basal species and thus being the rst to be analysed and further identi ed the motif fragment in their mitogenome, however no duplication was found.On the other hand, the B. boddarti presented 131 pb of the motif sequence, which was duplicated 2.4 times, whereas the other two species (B.pectinirostris and Boleophthalmus.sp) had this region repeated 5 times, meaning that there was a systematic increasing of the region (Fig. 3).This phenomenon is common in most sh mitogenomes (Sun et al.,2021), however there are other different gene organization the sh mitochondrial genome (Zhang et al., 2021).All the repeats were found in the 5´ to 3´ directions.
The alignment of the motif region of all Boleophthalmus was clearly similar showing that this duplication evolved in concert (Melters et al., 2013).In addition, studies show that, when the motif region is considered homologues, the phenomena is therefore considered concerted evolution ( In the newly sequenced mitogenome of B. dussumieri there were 17 intergenic spacer regions with a total sequence length of 899 bp (Table 3).Four of these regions were considered long IGS with the sequence length ranging from 14 bp, located between COX III and tRNA-Gly, 51 bp between tRNA-Val and r16S; 117 bp between ATP6/ COX III on the H-strand and one longest IGS of 495 bp located between tRNA-Pro/ OH on the L-strand.Large intergenic spacer has been reported in many vertebrate mitogenomes (Yoon et al., 2011;Li et al., 2021).Besides those, all of IGS had a length ranging from (1bp to 8bp).To explain this phenomenon, two main evolutionary mechanisms of IGS origin in mitogenomes is discussed, the tandem duplication/random loss (TDRL) model and slipped-strand mispairing (Shi et al.,2013;Li et al., 2015;Du et al.,2017).
The comparative analysis of the IGS and overlapping regions presented several variations in numbers, length, and locations within mudskippers.The number of IGS regions ranged from 13 to 17 with the total length varying from 418 (P.schloseri) to 1,207 (Boleophthalmus sp.) which have the largest IGS (706 bp) located in the OH region ( Table 3).Most of large IGS were located between tRNA-Pro and OH region, but only one (385 bp), belonging to P. novaeguineaensis, was located between tRNA-Leu1 and NADH5_0.
In overall, the overlapping length varied from 1b to 7bp in most mitogenomes.However, two species of the genus Periophthalmus presented unusual overlapping, involving 20 bp (P.minutus and P. novaeguineaensiss).This may be due to the split/duplication of the ATP8 gene in these species.The study also presents some peculiarities involving the IGS of the OH regions.We found large IGS of the OH region ranging from 61 bp (Periophthalmus modestus), to 185 bp (B.dussumieri), and very large varying from 314 bp (B.boddarti) to 701 bp (B.pectinirostris) (Table 3).

Phylogenetic analysis
The phylogenetic relationship within mudskippers were determined using a combination of 13 proteincoding genes (PCGs) and 2 ribosomal genes from a total of 17 mudskipper species.Additionally, three species of gobies (Tridentiger bifasciatus, Tridentiger barbatus and Stiphodon alcedo) were included as outgroups.The evolutionary tree was constructed using the ML method.The resulting phylogenetic tree for the mudskippers positioned them as monophyletic group with a high support value (ML bootstrap = 100) (Fig. 7-A).Several other studies have provided evidence for the monophyly of mudskippers using complete mitogenome analysis and analysis of the COI (Chen et al., 2014).These studies have identi ed two major groups within the mudskippers.The rst consists of all Periophthalmus species, which form a sister group of Periotphthalmodon, with high support values (> 90).The second group is comprised of two subclades.One subclade includes Boleophthalmus, which is sister group of Scatelaos, with a support value of > = 90.The other subclade consists of Apocryptodon, Paraprocryptes, and Oxuderces, with bootstrap value > 50.Scartelaos is grouped together with the Periophthalmus clade.The Periophthalmus clade is considered monophyletic and there is high similarity (> 99%) of Scartelaos histophorus with Periophthalmus magnuspinnatus (GenBank accession numbers KT284931.1,KT357639.1),so the presence of S. histophorus (accession number JQ654459.1) in the Periophthalmus clade should be considered misidenti cation (Sokefun et al., 2022).Therefore, the specie revalidation in the Genbank should be considered to avoid further taxonomic problems.
The topology position of P. barbarus presented a different clustering in both trees.In the ML tree the taxon clustered together with P. minutus, P. modestus, P. novaeguineaensis and P. argentilineatus as a sister group, whereas in the IB the species was positioned out of the sister genus.The results presented in this study brings a clari cation in relation to the separations of the groups within mudskippers.

Estimation of Divergence Times
The target species of the study (B.dussumieri) was co-related with other species of its genus and the genus Periophthalmus.The diversi cation of mudskippers from the family Gobiidae family was estimated to have occurred approximately 28.72 million years (Ma) (Fig. 7-B).Within the mudskippers, the Periophthalmus and Periophthalmodon genus were the rst to diverge and diversify, dating back to around 19.52 Ma in the Early Miocene.(Fig. 7-B).This suggests that these genera represent the basal lineages within terrestrial gobies.The divergence between the clades of Boleophthalmus and Periophthalmus occurred around 24.42 Ma in the early Miocene (Fig. 7-B).Boleophthalmus dussumieri and its sister species within the genus are estimated to have diverged around 8.

Conclusion
In this study we present the complete mitochondrial genome of the mudskipper species B. dusumieri with 16,685 bp of length and the rst comparative mitogenome gathering all available mudskipper mitogenomes.Generally, the organization followed almost the same con guration as most sh mitogenomes.However, we found different peculiarities regarding gene rearrangement, long intergenic spacer, and duplication of gene fragments mostly on the CR but also in a coding gene and one ribosomal gene in the analysed mitogenomes.The phylogenetic reconstruction using 13 concatenated PCG positioned all mudskipper species as monophyletic group within the family Oxudercidae.The phylogenetic clustering recovered in this study are in accordance with other previous study.The time scale estimation demonstrated that within mudskippers, the Periophthalmus genus was the rst to diverge in early Miocene and thus considered basal group.The complete mitogenome generated in this study is a source of genetic information for further molecular studies including conservation strategies.
After sequencing the fragments, the Bowtie tool, implemented in Geneious 9, was used to map the generated sequences against the reference and generate the consensus sequence.For mapping the consensus sequence, the species B. boddarti was used as a reference.The complete mitogenome was annotated using the Mitos 2 (http://mitos.bioinf.uni-leipzig.de)(Bernt et al., 2013) and later con rmed in MitosZ 3.6 (https://github.com/linzhi2013/MitoZ)with reference to Metazoa 63.The analysis of the repetition identi cation regions was performed using the nder tinder repeat tool (https://academic.oup.com/nar/article/27/2/573/1061099?login=true).The base composition, amino acid calculation as well as relative synonymous codon usage (RSCU) of PCGs were estimated using MEGA 11 (Tamura et al., 2021).

Phylogenetic analysis
To investigate the phylogenetic relationship among all mudskippers, a phylogenetic tree was constructed using a nal dataset of 20 goby species.This dataset was based on the combinations of 13 PCG and two RNAs genes.
To accomplish that, we performed Maximum Likelihood (ML) analysis in the IQTree v 2.2.0 program.The best partitioning scheme for the database were adopted considering the position of the codons of each gene, as well as the evolutionary models for the respective partitioning scheme.Subsequently the ML analysis was run with 1000 bootstrap pseudo replicates.

Divergence time analysis
Divergence time analyses were conducted using BEAST v.
Gene order of seventeen mudskippers species and three outgroups (Gobiidae).

Figures Figure 1
Figures

Figure 4 Relative
Figure 4
Intergenic spacers (IGS) are non-coding regions that are found between genes and are typically found in Metazoa and most vertebrate animals and in most cases serve as a transcription promotor mitogenomes(Saton etal., 2016; Yang et al., 2018) and, they have in uence on the growth rate of some invertebrate animals (Gorokhova et al., 2002).They are important signal for evolutionary studies and even species delimitation (Xi et al., 2015; Xu et al., 2022; Xiao et al., 2022).Besides IGS event, there is also segments of gene overlapping in mitochondrial genome (Chandhini et al., 2021).

Table 3
Intergenic spacer and overlap regions.

Table 4
22 Ma during the late Miocene (Fig. 7-B).Comparing the divergence time of B. dussumieri with another sister genus, they were roughly 14.88 Ma for Exuderces, Apocriptodon and Paraprocryptes and 12.82 Ma for Scartelaos.These values are within the values recovered from other studies with closer families which are sister family of Mudskippers (Oxudercidae), the sub family Amblyopinae Günther, 1861 Gobionellinae Bleeker, 1874 and Sicydiinae T.N.Gill, 1860 (Herler et al., 2009; Lu et al., 2023).The Minimal and Maximum, with 95% Highest Posterior Density (HPD) of individual and group node ages are shown in Table 4. Con dence intervals of the diverge time scales.
(Drummond et al., 2012)2012).The uncorrelated relaxed clock(Drummond et al., 2006)was employed, and the Yule process served as the prior model for the tree(Gernhard, 2008).Calibration points were incorporated to estimate divergence time, including the separation between B. boddarti and B. pectinirostris being calibrated to 3.8 Ma, based on Chen et al., (2014), while the separation between P. schlosseri and Periophthalmus genus was calibrated to 20.09 Ma, according to Zhanget al., (2017).The analysis was run for 100 million generations, with log parameters recorded every 100 registered at each 5000 generations.Trees were summarized in TreeAnnotator v.1.8.4(Drummond et al., 2012).A burnin of 20% was used and the run was considered satisfactory when all ESS values checked in Tracer v.1.6(Rambautet al., 2014)were equal to or larger than 200.