Comparative Mitogenome Analyses Uncover Mitogenome Features and Phylogenetic Implications of the Parrotfishes (Perciformes: Scaridae)

Simple Summary Parrotfishes are among the most colorful and diverse inhabitants of the coral reefs and sea grass beds and are ecologically important in these habitats. Here, we presented the complete mitogenome sequences from twelve parrotfish species and conducted comparative analysis of mitogenome features among the seven published species for the first time. The comparative analysis revealed both the conserved and unique characteristics of parrotfish mitogenomes. The mitogenome structure, organization, gene overlaps, putative secondary structures of transfer RNAs, and codon usage were relatively conserved among all the analyzed species. However, the base composition and the intergenic spacers varied largely among species. All of the protein-coding genes were under purifying selection. Phylogenetic analysis revealed that the parrotfishes could be divided into two clades with distinct ecological adaptations. Early divergence of these two clades was probably related to the expansion of sea grass habitat, and later diversifications were likely associated with the geomorphology alternation since the closing of the Tethys Ocean. This work offered fundamental materials for further studies on the evolution and conservation of parrotfishes. Abstract In order to investigate the molecular evolution of mitogenomes among the family Scaridae, the complete mitogenome sequences of twelve parrotfish species were determined and compared with those of seven other parrotfish species. The comparative analysis revealed that the general features and organization of the mitogenome were similar among the 19 parrotfish species. The base composition was similar among the parrotfishes, with the exception of the genus Calotomus, which exhibited an unusual negative AT skew in the whole mitogenome. The PCGs showed similar codon usage, and all of them underwent a strong purifying selection. The gene rearrangement typical of the parrotfishes was detected, with the tRNAMet inserted between the tRNAIle and tRNAGln, and the tRNAGln was followed by a putative tRNAMet pseudogene. The parrotfish mitogenomes displayed conserved gene overlaps and secondary structure in most tRNA genes, while the non-coding intergenic spacers varied among species. Phylogenetic analysis based on the thirteen PCGs and two rRNAs strongly supported the hypothesis that the parrotfishes could be subdivided into two clades with distinct ecological adaptations. The early divergence of the sea grass and coral reef clades occurred in the late Oligocene, probably related to the expansion of sea grass habitat. Later diversification within the coral reef clade could be dated back to the Miocene, likely associated with the geomorphology alternation since the closing of the Tethys Ocean. This work provided fundamental molecular data that will be useful for species identification, conservation, and further studies on the evolution of parrotfishes.


Introduction
The mitochondrial genome (mitogenome) of a vertebrate is a small (16-17 kb), compact, and circular double-stranded molecule, typically encoding 13 protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), two ribosomal RNA genes (rRNAs), and two non-coding regions (the origin of L-strand replication, O L , and control region, CR) [1]. The mitochondrial DNA sequences have been extensively employed in a variety of study areas, from phylogeography, which elucidates the spatial arrangement of genetic variation among populations or closely related species [2][3][4], to phylogenetic studies, which decipher the evolutionary relationships across a wide range of taxa at higher taxonomic levels [5][6][7]. Compared to single or a few mitochondrial gene-based markers, the complete mitogenome sequences generally provide much finer phylogenetic resolution [8]. Moreover, genomelevel characteristics, including nucleotide composition, genome structural arrangement, overlap, and non-coding intergenic spacers between genes, vary largely among different species and might possess evolutionary significance [9][10][11].
Parrotfishes (Scaridae) are among the most colorful and diverse inhabitants of coral reefs and sea grass beds [12]. Currently, a total of 100 species belonging to 10 genera are recognized, with Scarus being the most specious genus (52 species) [13]. These fish are mainly herbivorous, foraging mostly by excavating or scraping surfaces of rocks and carbonate substrate that are encrusted with algae, bacterial mats, and detritus [14]. As such, it is widely recognized that parrotfishes play an important role in marine bioerosion [15,16] and serve as determinants of benthic community structure [17]. For example, parrotfish can exert a top-down control on algal communities to provide more space and resources for coals and promote the attachment and recruitment of coral larvae [18][19][20]. Therefore, it can help to mitigate the competition between coral reefs and macroalgae and increase the resilience of coral reef ecosystems subjected to anthropogenic or natural disturbances [21,22]. In addition, the excavating and scraping species can break the reef framework into sandsized sediments and facilitate the cycling of calcium carbonate on reefs, which are also dispensable agents in reef erosion and sediment production and transport [23,24].
Deciphering mitogenome structures and sequences can provide insights into evolutionary processes and contribute to species delimitation and conservation efforts [25,26]. Despite the fact that parrotfishes play an irreplaceable role in coral reef and sea grass bed habitats due to their unique behavioral and ecological characteristics, only a few studies have addressed their mitogenome characteristics [27][28][29], and comparative analysis is scarce. Although a handful of works have tried to elaborate on the phylogenetic relationships among the parrotfishes [30][31][32], none of them have addressed this question from a mitogenomic perspective. The deficiency of mitogenome data and comparative works hindered us from understanding the evolution of the parrotfish and establishing proper management and conservation decisions. In the present study, we reported twelve parrotfish mitogenomes for the first time and conducted comparative analysis with the published sequences from other seven species to elaborate the detailed features of the parrotfish mitogenomes. Additionally, we also investigated the phylogenetic relationships among these parrotfishes and estimated divergence times using mitogenome data. We hope that our newly generated data and results will provide some insights into the evolution of the parrotfishes as well as contributions towards the identification and conservation of these fishes.

Sampling, DNA Extraction, PCR Amplification, and Sequencing
In the present study, we de novo sequenced twelve parrotfish species (with one specimen each): Calotomus carolinus, Cetoscarus bicolor, Hipposcarus longiceps, Scarus Genomic DNA was extracted from either a small piece of flesh or a pelvic fin clip taken from the right side of the specimen using the E.Z.N.A. ® Tissue DNA Kit (OMEGA, Beijing, China) and following the manufacturer's instructions. High-quality DNA samples were randomly broken into fragments with a length of 300~500 bp. Then complete genomic libraries were established using the Illumina TruSeq TM Nano DNA Sample Prep Kit (Illumina, San Diego, CA, USA) following the manufacturer's recommendation. The 150-bp paired-end sequencing was performed on the Illumina HiSeq2500 platform. Library construction and sequencing were performed by the Biozeron Corporation (Shanghai, China).

Sequence Assembly, Annotation, and Analyses
Prior to assembly, raw reads were filtered by Trimmomatic v0.39 [33] in order to remove the reads with adaptors, the reads showing a quality score below 20 (Q < 20), the reads containing a percentage of uncalled bases ("N" characters) equal to or greater than 10%, and the duplicated sequences. GetOrganelle 1.7.5 was used to assemble the mitogenomes [34]. The newly generated mitogenome sequences were deposited in Genbank under the accession numbers OQ349180-OQ349191. Annotation of the mitogenomes (PCGs, tRNAs, rRNAs, and CR) was performed using MITOS [35] and Mitoannotator v3.83 [36]. Transfer RNA (tRNA) genes and their secondary structures were determined by the MITOS webserver [35]. The base composition and codon distributions were analyzed in MEGA 7.0 [37], and the nucleotide composition skewness was calculated using the formulas (A − T)/(A + T) for AT skew and (G − C)/(G + C) for GC skew. Relative synonymous codon usage (RSCU) was calculated using DAMBE 7 [38]. The conserved sequence block domains (CSBs) were determined by comparing them with those of other species [1].

Phylogenetic Analyses
Prior to the phylogenetic analysis, the method of Xia et al. [39] was used to access substitution saturation of the sequences by comparing the information entropy-based index (I SS ) with critical values (I SS.c ) in DAMBE 7 [38]. If I SS is significantly lower than I SS.c , then sequences have not experienced substitution saturation. The sequence of the control region showed significant substitution saturation (I SS = 1.1897 > I SS.c = 0.7851, p < 0.001) and was thus excluded from further analysis. The phylogenetic relationships were reconstructed using the 13 PCGs and 2 rRNAs of the 19 parrotfish mitogenomes. Three Cheilinus species (C. fasciatus, NC037707; C. oxycephalus, NC061045; C. undulatus, NC013842) were used as outgroup taxa. Multiple sequence alignment was performed using MAFFT [40] implemented in PhyloSuite [41] under default parameters and subsequently checked by eye in SeaView [42]. Our dataset was partitioned by gene and codon position, and then the best-fit nucleotide substitution model for each partition was determined using Modelfinder [43]. Phylogenetic relationships were reconstructed using Bayesian inference (BI) and maximum likelihood (ML) approaches. BI was carried out in Mr. Bayes 3.2.7 [44]. Two independent Markov chains were run with 1 × 10 6 iterations, and 10,000 trees were retained, with the first 25% of the samples discarded as burn-in. ML analysis was conducted in IQTREE [45] under 10,000 ultrafast bootstrap replicates. DNAsp 6 [46] was used to calculate non-synonymous substitution rates (dN), synonymous substitution rates (dS), and the ratio of dN/dS (ω). MCMCTree, implemented in the PAML4.9i software package [47], was used to estimate the divergence time among the parrotfishes. The tree topology generated from BI was calibrated with fossil dates. The information of branch lengths, gradients, and hessian were first estimated with a maximum likelihood method in BsaeML of the PAML package. Then the MCMC approximation was performed with a burn-in period of 50,000 cycles, and a total of 10,000 samples were generated every 50 iterations. Two independent runs were performed. Tracer 1.7 [48] was used to check for effective sample sizes (ESS) of parameters. The ESS larger than 200 were considered to reach convergence.
Two fossil calibration points were used in the divergence time estimation. Calotomus preisli was known from the middle Miocene (~14 Ma) in Austria [49]. We calibrated the minimum age of the split between the sea grass clade and the coral reef clade using this fossil. The fossil elements belonging to the genus Bolbometopon were known from the late Miocene (~5.3 Ma) [49]. These fossils were used to set the minimum age of the separation between Bolbometopon and Cetoscarus. The root age of our phylogeny was set to be lower than 50 Ma, for the oldest known labrid fossil was dated back to 50 Ma from the Monte Bolca in Italy [50].

General Features of Mitochondrial Genomes
The total length of the 12 newly sequenced complete mitogenomes ranged from 16,657 bp in Scarus niger to 16,816 bp in Scarus globiceps. The typical set of 37 genes, including 13 PCGs, two rRNAs, and 22 tRNAs, and a control region, were detected in all the mitogenomes (Table 1, Figure 1, Supplementary File S1: Table S1). All PCGs were encoded on the Heavy (H) strand except for NADH dehydrogenase subunit 6 (ND6), which was located on the Light (L) strand. Eight tRNAs (tRNA Gln , tRNA Ala , tRNA Asn , tRNA Cys , tRNA Tyr , tRNA Ser (UGA) , tRNA Glu and tRNA Pro ) were located on the L-strand, and the remaining 14 tRNAs were on the H-strand ( Figure 1, Table 1). This coding pattern on the H and L-strand was identical among the 19 parrotfish species (Additional File 1: Table S1) and was consistent with most vertebrates [51].

Nucleotide Composition of the Parrotfish Mitogenomes and Unusual AT Skew of Calotomus Species
The nucleotide composition was similar among all of the parrotfish species, with the overall A + T content ranging from 53.0% in Scarus globiceps to 56.4% in Bolbometopon muricatum, and the A + T content was the lowest in ND4L (51.7 ± 3.6%) and the highest in the control region (62.6 ± 4.1%) (Figure 2a, Additional File 1: Table S2). All of the parrotfish mitogenomes exhibited AT bias, with the largest and most positive value observed in the rRNAs and the smallest and most negative value found in ND6 (Figure 2b, Additional File 1: Table S3). Compared with other parrotfishes, species of the genus Calotomus exhibited an unusual AT skew for the whole mitogenomes with a slightly negative value (−0.04 to −0.02), while other species all displayed a positive value (Figure 2b, Additional File 1: Table S3). These results indicated that species of the genus Calotomus displayed an excess of T over A in the whole mitogenome.
Biology 2023, 12, x FOR PEER REVIEW 6 of 16 Figure 1. Organization of the parrotfish mitogenome. Calotomus carolinus was taken as an example. The inner ring indicated GC content.

Nucleotide Composition of the Parrotfish Mitogenomes and Unusual AT Skew of Calotomus Species
The nucleotide composition was similar among all of the parrotfish species, with the overall A + T content ranging from 53.0% in Scarus globiceps to 56.4% in Bolbometopon muricatum, and the A + T content was the lowest in ND4L (51.7 ± 3.6%) and the highest in the control region (62.6 ± 4.1%) (Figure 2a, Additional File 1: Table S2). All of the parrotfish mitogenomes exhibited AT bias, with the largest and most positive value observed in the rRNAs and the smallest and most negative value found in ND6 (Figure 2b, Additional File 1: Table S3). Compared with other parrotfishes, species of the genus Calotomus exhibited an unusual AT skew for the whole mitogenomes with a slightly negative value (−0.04 to −0.02), while other species all displayed a positive value (Figure 2b, Additional File 1: Table S3). These results indicated that species of the genus Calotomus displayed an excess of T over A in the whole mitogenome.

Protein-Coding Genes
The total length of PCGs ranged from 11,391 bp to 11,415 bp, with ATP8 being the shortest (168 bp) and ND5 being the longest (1839 bp to 1848 bp). Most genes exhibited the typical start codon ATN. However, COI initiated with GTG in all species, and ATP6 started with GTG in Calotomus japonicus and Scarus oviceps or CTG in Calotomus carolinus (Additional File 1: Table S1). Four types of stop codons were detected, including two canonical (TAA and TAG) and two truncated codons (T-and TA-) (Additional File 1: Table S1). The incomplete stop codons were commonly observed in fish mitogenomes [1] and might be completed by post-transcriptional polyadenylation [52].
For all the parrotfish mitogenomes, Leu (CUN) , Ala, and Thr were the three most frequently translated amnio acids, while Cys was the least used amnio acid (Figure 3a). Moreover, the most frequently used codon was CGA for arginine in all the parrotfish mitogenomes ( Figure 3b). The RSCU revealed that degenerate codons were biased to use more A and T than G and C in the third codon position, which resulted in higher A + T content than G + C content in the third codon position of parrotfish mitogenomes (Figures 2a and 3b, Supplementary File S2: Figure S1).

Protein-Coding Genes
The total length of PCGs ranged from 11,391 bp to 11,415 bp, with ATP8 being the shortest (168 bp) and ND5 being the longest (1839 bp to 1848 bp). Most genes exhibited the typical start codon ATN. However, COI initiated with GTG in all species, and ATP6 started with GTG in Calotomus japonicus and Scarus oviceps or CTG in Calotomus carolinus (Additional File 1: Table S1). Four types of stop codons were detected, including two canonical (TAA and TAG) and two truncated codons (T--and TA-) (Additional File 1: Table  S1). The incomplete stop codons were commonly observed in fish mitogenomes [1] and might be completed by post-transcriptional polyadenylation [52].
For all the parrotfish mitogenomes, Leu (CUN) , Ala, and Thr were the three most frequently translated amnio acids, while Cys was the least used amnio acid (Figure 3a). Moreover, the most frequently used codon was CGA for arginine in all the parrotfish mitogenomes ( Figure 3b). The RSCU revealed that degenerate codons were biased to use more A and T than G and C in the third codon position, which resulted in higher A + T content than G + C content in the third codon position of parrotfish mitogenomes (Figures 3b and 2a, Supplementary File S2: Figure S1).    (Figure 4b), thus not forming the typical clover-leaf structure. A gene rearrangement of the tRNA gene cluster between ND1 and ND2 was detected, with the tRNA Met inserted between the tRNA Ile and tRNA Gln , and the tRNA Gln was followed by a putative tRNA Met pseudogene.

Gene Rearrangement and Secondary Structure of tRNAs
All 22 tRNAs typical of the mitogenomes of vertebrates were found in the parrotfish mitogenomes (Figure 4a). Most tRNAs could be folded into the canonical clover-leaf secondary structure. The secondary structure of tRNAs generally consisted of four domains and a short variable loop: the amino acid acceptor (AA) stem, the dihydrouridine (D) arm (D stem and loop), the anticodon (AC) arm (AC stem and loop), the thymidine (T) arm (T stem and loop), and the variable (V) loop (Figure 4a). However, tRNA Ser(AGN) in Bolbometopon muricatum and Calotomus japonicus possessed only small loop(s) in their D arms (Figure 4b), thus not forming the typical clover-leaf structure. A gene rearrangement of the tRNA gene cluster between ND1 and ND2 was detected, with the tRNA Met inserted between the tRNA Ile and tRNA Gln , and the tRNA Gln was followed by a putative tRNA Met pseudogene.

Overlaps and Non-Coding Intergenic Spacers
A total of four gene overlaps were detected in the mitogenome of Calotomus carolinus and five were observed in the mitogenomes of other parrotfishes (Table 1, Additional File 1: Table S1). The longest overlap was found between ATP8 and ATP6, with highly conserved 10-bp motifs of "ATGGCACTAA" or "ATGACACTAA" detected in most parrotfish mitogenomes except for that of the genus Calotomus. The latter genus showed 16-bp overlaps of "CTGACCTTGGCACTAG" or "GTGGCCCTGGCACTAG". Apart from that,

Overlaps and Non-Coding Intergenic Spacers
A total of four gene overlaps were detected in the mitogenome of Calotomus carolinus and five were observed in the mitogenomes of other parrotfishes (Table 1, Additional File 1: Table S1). The longest overlap was found between ATP8 and ATP6, with highly conserved 10-bp motifs of "ATGGCACTAA" or "ATGACACTAA" detected in most parrotfish mitogenomes except for that of the genus Calotomus. The latter genus showed 16-bp overlaps of "CTGACCTTGGCACTAG" or "GTGGCCCTGGCACTAG". Apart from that, a 7-bp overlap was observed between ND4L and ND4 in all parrotfish mitogenomes with highly conserved sequences of "ATGCTAA" or "ATGTTAA".
Two long intergenic spacers (IGS; tRNA Gln -ND2 and O L ) were found in all the parrotfish mitogenomes. Moreover, another long IGS between tRNA Glu and Cyt b was also found in the mitogenomes of the genus Calotomus. As mentioned above, the IGS between tRNA Gln and ND2 was assumed to be a pseudogene of tRNA Met . O L is located within the five tRNA gene cluster (WANCY), and its secondary structure showed a stable stem-loop hairpin, which is strengthened by 9 to 10 G-C base pairs (Figure 5). The G-C base pairs on the stem were highly conserved, while the loop varied in its base composition, with T being scarce.
Biology 2023, 12, x FOR PEER REVIEW 9 of 16 a 7-bp overlap was observed between ND4L and ND4 in all parrotfish mitogenomes with highly conserved sequences of "ATGCTAA" or "ATGTTAA". Two long intergenic spacers (IGS; tRNA Gln -ND2 and OL) were found in all the parrotfish mitogenomes. Moreover, another long IGS between tRNA Glu and Cyt b was also found in the mitogenomes of the genus Calotomus. As mentioned above, the IGS between tRNA Gln and ND2 was assumed to be a pseudogene of tRNA Met . OL is located within the five tRNA gene cluster (WANCY), and its secondary structure showed a stable stem-loop hairpin, which is strengthened by 9 to 10 G-C base pairs (Figure 5). The G-C base pairs on the stem were highly conserved, while the loop varied in its base composition, with T being scarce. The control region, located between tRNA Pro and tRNA Phe , was the most variable region and constituted the majority of the length variation of the parrotfish mitogenomes (Additional File 1: Table S1). Only three conserved sequence blocks (CSB-D, CSB-I, and CSB-II) were detected ( Figure 6), with CSB-III completely missing in all the parrotfish mitogenomes. The base composition was extremely unique to each CSB, with CSB-D being T rich, CSB-I being AT rich, and CSB-II being C rich ( Table 2).

Non-Synonymous and Synonymous Substitutions
To better understand the role of selective pressure and the evolutionary patterns of the protein coding genes, the dN/dS value (ω) of each PCG was calculated (Figure 7). All of the PCGs were subject to purifying selection, with a dN/dS value lower than 1 (ω < 1). The control region, located between tRNA Pro and tRNA Phe , was the most variable region and constituted the majority of the length variation of the parrotfish mitogenomes (Additional File 1: Table S1). Only three conserved sequence blocks (CSB-D, CSB-I, and CSB-II) were detected ( Figure 6), with CSB-III completely missing in all the parrotfish mitogenomes. The base composition was extremely unique to each CSB, with CSB-D being T rich, CSB-I being AT rich, and CSB-II being C rich ( Table 2).
Biology 2023, 12, x FOR PEER REVIEW 9 of 16 a 7-bp overlap was observed between ND4L and ND4 in all parrotfish mitogenomes with highly conserved sequences of "ATGCTAA" or "ATGTTAA". Two long intergenic spacers (IGS; tRNA Gln -ND2 and OL) were found in all the parrotfish mitogenomes. Moreover, another long IGS between tRNA Glu and Cyt b was also found in the mitogenomes of the genus Calotomus. As mentioned above, the IGS between tRNA Gln and ND2 was assumed to be a pseudogene of tRNA Met . OL is located within the five tRNA gene cluster (WANCY), and its secondary structure showed a stable stem-loop hairpin, which is strengthened by 9 to 10 G-C base pairs (Figure 5). The G-C base pairs on the stem were highly conserved, while the loop varied in its base composition, with T being scarce. The control region, located between tRNA Pro and tRNA Phe , was the most variable region and constituted the majority of the length variation of the parrotfish mitogenomes (Additional File 1: Table S1). Only three conserved sequence blocks (CSB-D, CSB-I, and CSB-II) were detected ( Figure 6), with CSB-III completely missing in all the parrotfish mitogenomes. The base composition was extremely unique to each CSB, with CSB-D being T rich, CSB-I being AT rich, and CSB-II being C rich ( Table 2).

Non-Synonymous and Synonymous Substitutions
To better understand the role of selective pressure and the evolutionary patterns of the protein coding genes, the dN/dS value (ω) of each PCG was calculated (Figure 7). All of the PCGs were subject to purifying selection, with a dN/dS value lower than 1 (ω < 1).

Non-Synonymous and Synonymous Substitutions
To better understand the role of selective pressure and the evolutionary patterns of the protein coding genes, the dN/dS value (ω) of each PCG was calculated (Figure 7). All of the PCGs were subject to purifying selection, with a dN/dS value lower than 1 (ω < 1). Among which, ATP8 and COI presented the highest and lowest ω values (ω = 0.300 and 0.016), respectively.

Phylogenetic Relationships of the Parrotfishes
The ML and BI trees based on the thirteen PCGs and two rRNAs yielded identical gene tree topologies (Figure 8

Phylogenetic Relationships of the Parrotfishes
The ML and BI trees based on the thirteen PCGs and two rRNAs yielded identical gene tree topologies (Figure 8), which congruently revealed two main clades. The first clade (clade A), located at the basal part of the tree, includes species of the genus Calotomus.

Phylogenetic Relationships of the Parrotfishes
The ML and BI trees based on the thirteen PCGs and two rRNAs yielded identical gene tree topologies (Figure 8

Divergence Time Estimation
The estimated divergence time and 95% credible intervals (CIs) are shown in Figure 9. The split between clade A and clade B occurred at 26

Divergence Time Estimation
The estimated divergence time and 95% credible intervals (CIs) are shown in Figure  9. The split between clade A and clade B occurred at 26.9 Ma (95% CI 16.0~36.0 Ma) during the late Oligocene. The Bolbometopon-Cetoscarus clade differentiated at 15.9 Ma (95% CI 9.3~21.2 Ma) during the middle Miocene. Hipposcarus diverged at 14.3 Ma (95% CI 8.3~18.8 Ma). The split between Chlorurus and Scarus was dated back to 8.6 Ma (95% CI 5.2~11.6 Ma). The Scarus species diverged relatively recently, ranging from 0.2 Ma to 7.0 Ma. Figure 9. Divergence time estimation of the parrotfishes derived from MCMCTree analysis. Numbers at nodes indicate estimated age. Blue bars represent 95% credible age intervals for each node.

Discussion
The comparative analysis revealed that the mitogenome structure, organization, codon usage, and putative secondary structures of tRNAs were highly similar among all the analyzed parrotfish species. The gene rearrangement of the tRNA gene cluster between ND1 and ND2 was detected, which is typical of parrotfish [27][28][29]. In parrotfish mitogenomes, the tRNA Met was located between the tRNA Ile and tRNA Gln , then a putative tRNA Met pseudogene was located after the tRNA Gln . The gene rearrangements had been proposed to occur with tandem duplication of gene regions as a consequence of slipped-strand mispairing, followed by deletions of redundant genes [53]. The tRNA Met pseudogene was believed to function as punctuation marks for mitochondrial ND2 mRNA processing [27].
Previous studies on insects suggested that the intergenic spacers were important for transcription and might be associated with gene rearrangement [54][55][56]. Our results showed significant variance in IGS among the parrotfish mitogenomes, especially for the longest IGS, the control region. Despite the great length variations found in the CR of the parrotfish mitogenomes, three conserved sequence blocks could still be detected. Compared with most fish species [1], the CSB-III cannot be observed in the CR of the parrotfish mitogenomes. The lack of CSB-III was also reported in other vertebrates [57]. Up until now, the functions of the CSBs were still not clear, however, the common existence of CSB-D and CSB-I in vertebrate mitogenomes suggested that they were vital in the replication and transcription of the genome [1].
CR was commonly used as genetic markers in phylogenetic and population genetic analysis due to its high variability among populations and closely related species [58,59]. However, our analysis suggested that the CR of parrotfishes experienced significant substitution saturation. Substitution saturation reduces the amount of phylogenetic signals to the point that sequence similarities could probably be the consequence of chance alone rather than homology. Therefore, phylogenetic signals are lost, and the sequences are no longer informative about the underlaying evolutionary processes that generate them if substitution saturation is reached [60]. For example, the mitochondrial markers COI and

Discussion
The comparative analysis revealed that the mitogenome structure, organization, codon usage, and putative secondary structures of tRNAs were highly similar among all the analyzed parrotfish species. The gene rearrangement of the tRNA gene cluster between ND1 and ND2 was detected, which is typical of parrotfish [27][28][29]. In parrotfish mitogenomes, the tRNA Met was located between the tRNA Ile and tRNA Gln , then a putative tRNA Met pseudogene was located after the tRNA Gln . The gene rearrangements had been proposed to occur with tandem duplication of gene regions as a consequence of slipped-strand mispairing, followed by deletions of redundant genes [53]. The tRNA Met pseudogene was believed to function as punctuation marks for mitochondrial ND2 mRNA processing [27].
Previous studies on insects suggested that the intergenic spacers were important for transcription and might be associated with gene rearrangement [54][55][56]. Our results showed significant variance in IGS among the parrotfish mitogenomes, especially for the longest IGS, the control region. Despite the great length variations found in the CR of the parrotfish mitogenomes, three conserved sequence blocks could still be detected. Compared with most fish species [1], the CSB-III cannot be observed in the CR of the parrotfish mitogenomes. The lack of CSB-III was also reported in other vertebrates [57]. Up until now, the functions of the CSBs were still not clear, however, the common existence of CSB-D and CSB-I in vertebrate mitogenomes suggested that they were vital in the replication and transcription of the genome [1].
CR was commonly used as genetic markers in phylogenetic and population genetic analysis due to its high variability among populations and closely related species [58,59]. However, our analysis suggested that the CR of parrotfishes experienced significant substitution saturation. Substitution saturation reduces the amount of phylogenetic signals to the point that sequence similarities could probably be the consequence of chance alone rather than homology. Therefore, phylogenetic signals are lost, and the sequences are no longer informative about the underlaying evolutionary processes that generate them if substitution saturation is reached [60]. For example, the mitochondrial markers COI and ND3 that are commonly used in phylogenetic studies and DNA barcoding were proven to be subjected to significant substitution saturation in Caryophyllidean cestodes. Therefore, arbitrary application of these markers to the phylogenetic inference of this group of cestodes would jeopardize the well-supported phylogenetic estimates and evolutionary relationships [61]. In our case, the CR sequences have never been employed to infer the phylogenetic relationships among parrotfishes so far [30][31][32]. Future studies should avoid using the CR sequences when it comes to phylogenetic relationship inference or identification via DNA barcoding of the parrotfishes.
The RSCU revealed that degenerate codons were prone to use more A and T than G and C in the third codon position, therefore higher A + T content than G + C content was observed in the third codon position of parrotfish mitogenomes (Figures 2a and 3b, Additional File 2: Figure S1). This phenomenon was frequently observed in other teleosts [1] and might be related to genome bias, optimal selection of tRNA, or DNA repair efficiency [62].
Compared with other parrotfish species, the Calotomus species displayed an unusual AT skew for the whole mitogenome with a slightly negative skewness, while other species all showed positive values. Nucleotide skewness might be related to the balance between mutational and selective pressures during replication [63][64][65]. Some previous studies had indicated that the preference for certain nucleotides might be associated with selection rather than mutation [66]. For example, Sinorhodeus microlepis, a bitterling species with highly specialized ecological and behavioral preferences [67], also exhibited an unusual AT skew in its mitogenome [68] and this unique AT skewness was believed to be associated with unique selective forces [68]. Compared with other parrotfish species, the Calotomus species possessed some unique ecological aspects, such as the browsing feeding behavior and the lack of breeding territories [30]. It is suspected that distinct selective pressures or processes might lead to the preference of T in their mitogenomes. However, what and how the selective processes account for the unusual AT skew in the Calotomus species needs further investigation.
All of the PCGs were evolved under the purifying selection (ω < 1). The lower ω value on the whole suggested a prevalent signature of strong functional restrictions across the mitogenome, which was largely in agreement with the functional importance of mitochondria as a respiration chain necessary for OXPHOS and electron transport [69]. Furthermore, the lower ω value indicated fewer variations in the amino acids; therefore, COI and Cyt b could serve as potential barcoding markers for the identification of parrotfish.
The phylogenetic relationships among the parrotfish genus based on thirteen PCGs and two rRNAs of the mitogenome indicated two distinct clades (A and B), which were identical with previous studies based on concatenated data of both mitochondrial and nuclear markers [30][31][32]. These two clades recovered by the phylogenetic analysis correspond to two distinct groups with different aspects of ecological adaptation, which had been defined as the sea grass clade and the coral reef clade [30]. The sea grass clade, as represented by Calotomus in this study, exhibited some less modified morphological characteristics (e.g., discrete teeth without cementation) [70] and showed some distinct ecological and behavioral aspects (e.g., browsing, no breeding territories, and no harem) [30]. These features differed greatly from the coral reef clades. In addition to the phylogenetic analysis, our results also revealed some unique features of the mitogenome composition and organization in the Calotomus mitogenomes (e.g., unusual AT skewness in the mitogenome and additional IGSs), indicating the evolutionary distinctiveness of the sea grass clade.
The first split of the parrotfish was estimated to have occurred in the late Oligocene (26.9 Ma, 95% CI 16.0~36.0 Ma), separating the sea grass clade from the coral reef clade. Geological evidence suggests that tectonic movements in the Indo-West Pacific region during the late Oligocene and early Miocene resulted in the formation of vast areas of shallow-water habitat between Australia and Indonesia [71], facilitating the expansion of sea grass habitat [72]. Our divergence time estimation was largely in congruence with the timing of the large-scale development of the sea grass habitat. This result probably indicated that the ecological differences between these two habitats acted as the major driving force in the early diversification of the parrotfishes. The differentiation within the coral reef clade had been initiated since the middle Miocene (about 15.9 Ma), which is well consistent with the closure of the Tethys Ocean [73]. Alterations in geomorphologies such as sea levels, sea surface temperatures, and ocean circulations exerted a great impact on coral reefs [74][75][76], likely functioning as the driving forces behind the rapid radiation of coral reef species. Previous studies indicated that the extensive diversification of coral reef taxa occurred during this period and was likely associated with the geomorphological reconfiguration of the marine realm [77]. In addition, natural and sexual selections might have also contributed to the diversification of parrotfishes. Some studies suggested that the protogynous mating system of parrotfishes might function as a possible driving force of speciation [30]. Though some research has suggested that ecological and selection may operate in tandem in the speciation processes [31], the function mechanisms and their relative roles still require further investigation.

Conclusions
In the present study, comparative analysis revealed both the conserved and unique characteristics of parrotfish mitogenomes. The mitogenome structure, organization, gene overlaps, putative secondary structures of tRNAs, and codon usage were relatively conserved among all the analyzed species. However, the base composition and the intergenic spacers varied largely among species. All of the PCGs were under purifying selection. Phylogenetic analysis revealed that the parrotfishes could be divided into two clades with distinct ecological adaptations. Early divergence of the sea grass and coral reef clades occurred in the late Oligocene, probably related to the expansion of sea grass habitat. Later diversification within the coral reef clade could be dated back to the Miocene, likely associated with the geomorphology alternation since the closing of the Tethys Ocean. This study offered fundamental molecular materials for further studies on the evolution and diversification of the parrotfishes and would contribute to their identification and conservation.