De novo Transcriptome Analysis and the Phylogenetic Position of Acanthopagrus schlegelii

Yilin Wang1, Jianfeng Guan2, Dongxiao He3 and Xiangbin Meng4*

1International Education College, Nanjing Forestry University, Nanjing 210037, China

2Rizhao Marine Communication Management and Service Center, Rizhao 276800, China

3Rizhao Ocean and Fisheries Research Institute, Rizhao 276800, China

ABSTRACT

The black porgy (Acanthopagrus schlegelii), a protandrous hermaphroditic marine fish, holds significant economic and ecological value. However, its reproductive biology and phylogenetic position remain poorly understood at the molecular level. This study provides essential transcriptomic resources to explore the regulatory mechanisms involved in sex reversal and to clarify the phylogenetic relationships of A. schlegelii. Using Illumina RNA-seq technology, we generated a comprehensive transcriptome dataset from various tissues, yielding 43,848,600 high-quality clean reads. The de novo assembly produced 36,322 unigenes, of which 30,131 (82.59%) were successfully annotated in multiple protein databases. Phylogenetic analysis, based on transcriptome data from nine species within Perciformes, revealed that A. schlegelii is closely related to Epinephelus fuscoguttatus, with a divergence time estimated at 49.02 to 59.51 million years ago. The transcriptomic data and phylogenetic insights provided by this study contribute to a deeper understanding of the evolutionary biology and reproductive ecology of A. schlegelii, offering valuable resources for future investigations into non-model fish species.

Article Information

Received 17 March 2024

Revised 22 October 2024

Accepted 31 October 2024

Available online 23 June 2025

(early access)

Published 12 March 2026

Authors’ Contribution

YW: Writing manuscript, reviewing and editing, resources, formal analysis. JG: Reviewing and editing, writing manuscript. DH: Visualization, writing review and editing. XM: Supervision, writing review and editing.

Key words

Acanthopagrus schlegelii, Transcriptome analysis, RNA-seq, Phylogenetics, Sex reversal, Marine fish biology

DOI: https://dx.doi.org/10.17582/journal.pjz/20240317050911

* Corresponding author: [email protected]

0030-9923/2026/0003-1009 $ 9.00/0

Copyright 2026 by the authors. Licensee Zoological Society of Pakistan.

This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Introduction

The black porgy, Acanthopagrus schlegelii, is belongs to the Sparidae of Perciforms, which is widely distributed along the West Pacific coasts from Japan and Korea to the East China Sea (Wu et al., 2005; Liu et al., 2007). As a warm-water demersal fish, A. schlegelii typically inhabits the depth from 40 to 60 m with muddy fine sand bottoms and usually appears in the shadow of rocks. In addition, A. schlegelii is a eurythermic fish, the survival temperature ranges from 4 to 35 °C, with the suitable survival temperature ranges from 17 to 25℃. A. schlegelii also can tolerate wide range of salinity varying from 4.09 to 35.0%, and the suitable survival salinity ranges from 10.0 to 30.0% (Zhu, 2017). Although this species have stronger environmental tolerance and supports an important commercial fishery for China, growing climate change pressures and human activities having detrimental and irreversible consequences for the wild resources (Chang and Yueh, 1990; Iwatsuki and Carpenter, 2015). In fact, the intensive cultivation of A. schlegelii has become a topical issue in the past years. However, A. schlegelii is a marine protandrous hermaphrodite fish and have sexually reversal peculiarity, which ultimately might have brought about negative effects on aquaculture industry (Chang et al., 1994). This is probably because a higher proportion of milters were existed in mating phase, albeit there are no significant sex-specific differences in the body size. Therefore, investigating the regulatory mechanisms of sexually reversal may be essential to the breeding of A. schlegelii. In addition, it is worth noting that accurate construction of phylogenetic relationships is a fundamental step toward unraveling the evolutionary processes of sexually reversal of A. schlegelii (Yang and Rannala, 2012).

A growing number of studies had characterized the sexually reversal of A. schlegeliis based on histology, endocrine and molecular genetics (Huang et al., 2002; Wu et al., 2010). For example, Lee et al. (2001) considered that female A. schlegelii usually spawn during from January to March and they begin to sexually reversal to males after the third year. However, limited evidence was used to evaluate the regulatory mechanisms associated with sexually reversal of A. schlegelii. Additionally, although Ma et al. (2016) sequenced the complete mitochondrial genome of A. schlegelii and reconstructed the phylogenetic position, Knowledge of the phylogeny of A. schlegelii also is far from complete. Currently established high-throughput sequencing technologies, such as RNA-seq, has been widely used in the functional gene expression profiling and pathways study. This technology allows simultaneous analyses of all of the processes that are regulated at the transcription level and therefore it has been provide the opportunity to investigate the specific biological process (Tirosh et al., 2006). Recently, RNA were applied to identify the regulatory mechanisms associated with sexually reversal of some hermaphrodite fishes, such as Syngnathus scovelli, Acipenser gueldenstaedtii and others (Rose et al., 2015; Hagihara et al., 2014). Furthermore, RNA-seq enable systematists to analyze phylogenetic analysis of fishes based on hundreds to thousands of loci (Hughes et al., 2018). Therefore, the RNA-seq technology has clear advantages over existing approaches, and it can be used to accurate reveal the regulatory mechanisms and evolutionary processes of sexually reversal of A. schlegelii.

Materials and Methods

Samples collection, RNA extraction and illumina sequencing

A. schlegelii is not an endangered or protected species in China or other countries. In addition, frost anesthesia was made to minimize suffering of all A. schlegelii. Five healthy A. schlegeliis were obtained from an aquaculture farm in Zhoushan (China) in 2018. Then, muscle, liver, gill, heart, kidney, swim bladder and sexual gland of each individual were rapidly sampled, snap-frozen in liquid nitrogen and stored at -80℃ prior to the RNA extraction. Total RNA of every tissues of 5 fishes were extracted using a standard Trizol Reagent Kit following the manufacturer’s protocol and then pooled in equal amounts. The detection of degradation degree, purity (OD260/280), concentration and integrality of total RNA based on the agarose gel electrophoresis, Nanodrop, Qubit and Agilent 2100, respectively. RNA purification beads with Oligo (dT) were used to purify mRNA from total RNA. Then, mRNA fragmentation by using fragmentation buffer. First-strand cDNA synthesis was performed by combined fragmented mRNA and random hexamers. For second strand synthesis, we added buffer solution, dNTPs, DNA polymerase I and RNase H into 1st Strand cDNA. Then, purified double-stranded cDNAs by using AMPure XP beads were further to suffered end repair, added A-tailing and adapter ligation. The fragment size of cDNA was selected by using AMPure XP beads. PCR reaction was carry out and PCR products were further purified to obtain the sequencing libraries. The quantitative evaluation of libraries by using Qubit and then were diluted to 1.5 ng/ul. Finally, Agilent 2100 and Q-PCR was applied to detect the insert size and concentration of libraries, respectively. Then the library was sequenced on the Illumina HiSeqTM 2500 platform and 150bp paired-end reads were generated.

Transcriptome de novo assembly and annotation

FastQC software was applied to evaluate the quality of all raw reads in FASTQ format. Then, clean reads were obtained by removing reads with sequencing adaptors, unknown nucleotides (N ratio > 10%) and low quality (quality scores <= 20) based on Trimmomatic 0.36 (Bolger et al., 2014). The remaining high-quality clean reads were de novo assembled to obtained the reference transcripts based on Trinity 2.4.0 (Grabherr et al., 2011) with the parameter: --min_kmer_cov 3. Then the redundancy transcripts were removed using Corset 1.05 with default parameters and further spliced into the longest unigenes for further analyses. Homology searches were acquired by comparing all unigenes against the NR, NT, Pfam, KOG and Swiss-prot databases. On the protein annotation information of NR and Pfam database, Blast2GO 2.5 (Götz et al., 2008) was applied to analyze the biological processes, molecular function and cellular environment of gene products, e-value= 1e-6 was used as the filtering thresholds to determine the significance of GO terms. Furthermore, the biochemical metabolic pathways of gene products were predicted based on KEGG pathway annotation using KAAS (KEGG Automatic Annotation Server; http://www.genome.jp/kaas-bin/kaas_main), e-value = 1e-10 was used as the filtering thresholds to determine the significance of KEGG pathways (Moriya et al., 2007).

Predict the gene structure

Coding sequences (CDS) were extracted by comparing all unigenes with the protein databases and sequentially followed a fixed order of Nr and Swissprot. Then, ESTscan 3.0.3 was used to predict the CDS of unigenes that matched failure in the databases. All CDS will be translated into protein sequences according to the standard codes. Simple sequence repeats (SSRs) of all unigenes were analyzed by using MISA 1.0 with default parameters. SSRs were identified according to the minimum repeats of every unit size, the minimum repeats of 1bp, 2bp, 3bp, 4bp, 5bp and 6bp was 10, 6, 5, 5, 5 and 5, respectively. Then, Primer 3 was applied to the design of primers of SSRs.

Phylogenetic analysis for A. schlegelii and other Perciformes species

Except for the transcriptome data of A. schlegelii, we also downloaded the transcriptome data of other 8 Perciformes species (Anoplopoma fimbria, Trematomus bernacchii, Epinephelus fuscoguttatus, Sebastes caurinus, Perca fluviatilis, Chionodraco hamatus, Gymnodraco acuticeps, Lepidonotothen nudifrons) from NCBI. Firstly, CDS were extracted from every transcriptome data using TransDecoder 3.0.1 with the parameter: -m 200 (the length of protein sequences larger than 200bp). Then, all-against-all reciprocal blastp search of all protein sequences was carry out to obtain the single-copy orthologous genes among 9 species using OrthoMCL (Li et al., 2003). Multiple alignments of single-copy orthologous genes were performed using MAFFT software with default settings and then 9 protein sequence super-matrices were generated by concatenating aligned single-copy orthologous genes (Katoh and Standley, 2013). Conserved sequences were extracted from each super-matrice using Gblocks with parameter: -t=p (Castresana, 2000). The optimal amino acid substitution model JTT+G+F+I was calculated using ProTest software. Finally, the maximum-likelihood method implemented in the RAxML package with the amino acid substitution model were used to construct the phylogenetic tree based on the 11 protein sequence super-alignment (Stamatakis, 2014). iTOL software was applied to draw the phylogenetic tree. The divergence time was estimated using the r8s software and a molecular clock data from the divergence time between A. fimbria and T. bernacchii from the TimeTree database (Hedges et al., 2015).

Results and Discussion

Transcriptome de novo assembly

A total of 44,036,042 raw reads in FASTQ format were obtained from A. schlegelii. We further to evaluate the quality of all raw reads, results showed that the base valid ratio, Q20 %, Q30 % and GC % was 99.98%, 97.43%, 92.82% and 49.39%, respectively. Therefore, higher-quality transcriptome data was obtained in this study. After reads filtering, 43,848,600 high-quality clean reads were generated, corresponding to 6.58 G clean bases. All high-quality clean reads were de novo assembled to produce 59,603 transcripts, corresponding to 115,340,632 nucleotides. The min length, max length, mean length, median length, N50, N90 of transcripts was 301bp, 71,080bp, 1,935bp, 1,148bp, 3,411bp and 805bp, respectively. Then the redundancy transcripts were removed and further spliced into 36,322 longest unigenes, corresponding to 59,019,412 nucleotides. The min length, max length, mean length, median length, N50, N90 of unigenes was 301bp, 71,080bp, 1,625bp, 929bp, 2,791bp and 632bp, respectively. The length distribution of transcripts and unigenes is shown in Figure 1.

 

Unigenes annotation information

In order to obtain the gene functional information associated with sexually reversal, all unigenes were applied to compare with the sequences of seven databases. Results showed that homology sequences of 30,131 (82.59%) unigenes were acquired in at least on database and 7,309 (20.12%) unigenes were annotated in seven databases. Of all homology searches, 23,812 (65.55%), 28,339 (78.02%), 15,544 (42.79%), 20,776 (57.19%), 18,920 (52.08%), 18,920 (52.08%), 10,062 (27.70%) unigenes had significant matches with sequences in the NR, NT, KO, Swiss-prot, Pfam, GO and KOG databases, respectively. The number of unigenes successfully matched in the NR, NT, KOG, GO and Pfam databases were shown through a Venn diagram (Fig. 2).

We can evaluate the genetic sequence similarity of A. schlegelii and other species by comparing all unigenes against the NR protein database. Result showed that 23,812 unigenes of A. schlegelii were matched with protein sequences of 169 species (Fig. 3A). It is worth noting that a large proportion of unigenes of A. schlegelii have a strong similarity with the protein sequences from Larimichthys crocea (7,317, 33.14%), Lates calcarifer (3,896, 17.64%), Seriola dumerili (2,119, 9.60%), Stegastes partitus (1,022, 4.62%) and Labrus bergylta (1,016, 4.60%). It is worth noting that all of these species belong to Percoidei. Based on evolutionary relationship, homologous genes of different species are divided into different ortholog unigenes in KOG database. In the

 

present study, 10,062 unigenes were categorized into 26 KOG categories (Fig. 3B). Among these categories, the first three largest groups were T category (2,015 unigenes), R category (1,584 unigenes) and O category (974 unigenes), representing signal transduction mechanisms, general function prediction only and posttranslational modification, protein turnover, chaperones, respectively. GO classification can exactly define gene characteristics and then help us understanding what gene function might associated with sexually reversal. GO classification result showed that the terms of cellular process and binding were dominant in biological process and molecular function, respectively. In addition, cell and cell process were dominant in cellular component (Fig. 3C). KEGG analysis was used to analyze the metabolic pathways and function of gene products in cells. In this study, a total of 15,544 unigenes were assigned to 5 terms according to metabolic pathways: Cellular processes, environmental information processing, genetic information processing, metabolism and organismal systems (Fig. 3D). Additionally, the most significant pathways in 5 terms were cellular commiunity, folding, sorting and degradation, signal transduction, lipid metabolism and endocrine system, respectively.

 

Gene structure

We further analyzed the structure of all unigenes in the present study. Firstly, we extracted 17,298 CDS (Direction of the sequences is 5’→3’) by blast comparing all unigenes with the Nr and Swissprot databases. Then, we predicted the 15,636 CDS for these unigenes that matched failure in the database by ESTScan software. Additionally, 30,051 SSRs were identified and these were containing in 14,365 sequences. Among them, 6,869 sequences containing more than 1 SSR and 4,173 SSRs present in compound formation. It is worth noting that the relative abundances of specific repeat motifs were highly variable among the repeats (Fig. 4). Results showed that mono-nucleotide repeats was the most frequent repeat among all repeat classes, followed by di-nucleotide repeats and tri-nucleotide repeats. Moreover, AC, AGG, AAAC, AAAAC and ACCTCC repeats was the most frequent di-nucleotide, tri-nucleotide, tetra-nucleotide, penta-nucleotide and hexa-nucleotide motifs, respectively.

 

In this study, we identified a set of 303 single-copy, conserved exons (>60% identity among taxa) longer than 200 bp based on OrthoMCL software. The concatenated alignment of 303 sequences produced a data matrix with 81,948 amino acids for 9 species. Then, the data matrix was applied to construct phylogenetic tree. Phylogenetic analysis showed that A. schlegelii was clustered together with E. fuscoguttatus. Additionally, the divergence time of A. schlegelii and E. fuscoguttatus around 49.02-59.51 million years ago.

Conclusions

In this study, we successfully captured a significant portion of the A. schlegelii transcriptome based on RNA-seq. This was the first transcriptome study of A. schlegelii, and the results provided a fundamental understanding for research on the regulatory mechanisms associated with biological characteristics. Additionally, transcriptomic data of 9 Perciformes species were applied to reconstruct the phylogenomic relationships. We found that the divergence time of A. schlegelii and E. fuscoguttatus around 49.02-59.51 million years ago. Our phylogeny analysis of A. schlegelii based on significantly more genetic loci relative to previous studies. Therefore, we believe that the transcriptome data and analysis method in this work provided a valuable resource for research on the phylogenetic and adaption investigation of A. schlegelii and other fishes.

Declarations

Funding

The study received no external funding.

Statement of conflicts of interest

The authors have declared no conflict of interest.

References

Bolger, A.M., Lohse, M. and Usadel, B., 2014. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics, 30: 2114-2120. https://doi.org/10.1093/bioinformatics/btu170

Cao, Y., 2016. Genetic diversity of 3 Sparid species in coastal waters of China based on mitochondrial control region sequences. Dessertation Jinan University, Jinan.

Castresana, J., 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol., 17: 540-552. https://doi.org/10.1093/oxfordjournals.molbev.a026334

Chang, C.F., Lee, M.F. and Chen, G.R., 1994. Estradiol-17β associated with the sex reversal in protandrous black porgy, Acanthopagrus schlegeli. J. exp. Zool., 268: 53-58. https://doi.org/10.1002/jez.1402680107

Chang, C.F. and Yueh, W.S., 1990. Annual cycle of gonadal histology and steroid profiles in the juvenile males and adult females of the protandrous black porgy, Acanthopagrus schlegeli. Aquaculture, 91: 179-196. https://doi.org/10.1016/0044-8486(90)90187-R

Götz, S., García-Gómez, J.M., Terol, J., Williams, T.D., Nagaraj, S.H., Nueda, M.J., Robles, M., Talón, M., Dopazo, J. and Conesa, A., 2008. High-throughput functional annotation and data mining with the Blast2GO suite. Nucl. Acids Res., 36: 3420-3435. https://doi.org/10.1093/nar/gkn176

Grabherr, M.G., Haas, B.J., Yassour, M., Levin, J.Z., Thompson, D.A., Amit, I., Adiconis, X., Fan, L., Raychowdhury, R., Zeng, Q.D., Chen, Z.H., Mauceli, E., Hacohen, N., Gnirke, A., Rhind, N., di Palma, F., Birren, B.W., Nusbaum, C., Lindblad-Toh, K., Friedman, N. and Regev, A., 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol., 29: 644-652. https://doi.org/10.1038/nbt.1883

Hagihara, S., Yamashita, R., Yamamoto, R., Ishihara, S., Abe, M., Ijiri, T. and Adachi, S., 2014. Identification of genes involved in gonadal sex differentiation and the dimorphic expression pattern in undifferentiated gonads of Russian sturgeon Acipenser gueldenstaedtii Brandt and Ratzeburg, 1833. J. appl. Ichthyol., 30: 1557-1564. https://doi.org/10.1111/jai.12588

Hedges, S.B., Marin, J., Suleski, M., Paymer, M. and Kumar, S., 2015. Tree of life reveals clock-like 418 speciation and diversification. Mol. Biol. Evol., 32: 835-845. https://doi.org/10.1093/molbev/msv037

Huang, J.D., Lee, M.F. and Chang, C.F., 2002. The morphology of gonadal tissue and male germ cell in the protandrous black porgy, Acanthopagrus schlegeli. Zool. Stud., 41: 216-227.

Hughes, L.C., Ortí, G., Huang, Y., Sun, Y., Baldwin, C.C., Thompson, A.W., Arcila, D., Betancur-R, R., Li, C., Becker, L., Bellora, N., Zhao, X., Li, X., Wang, M., Fang, C., Xie, B., Zhou, Z., Huang, H., Chen, S., Venkatesh, B. and Shi, Q., 2018. Comprehensive phylogeny of ray-finned fishes (Actinopterygii) based on transcriptomic and genomic data. Proc. natl. Acad. Sci. USA, 115: 6249-6254. https://doi.org/10.1073/pnas.1719358115

Iwatsuki, Y. and Carpenter, K., 2015. Acanthopagrus sivicolus. The IUCN red list of threatened species. Available at http://www.iucnredlist.org/details/170199/0/

Katoh, K. and Standley, D.M., 2013. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol., 30: 772-780. https://doi.org/10.1093/molbev/mst010

Lau, E.L., Lee, M.F. and Chang, C.F., 2013. Conserved sex-specific timing of meiotic initiation during sex differentiation in the protandrous black porgy Acanthopagrus schlegelii. Biol. Reprod., 88: 1-13. https://doi.org/10.1095/biolreprod.112.105536

Lee, Y.H., Du, J.L., Yueh, W.S., Lin, B.Y., Huang, J.D., Lee, C.Y., Lee, M.F., Lau, E.L., Lee, F.Y., Morrey, C., Nagahama, Y. and Chang, C.F., 2001. Sex change in the protandrous black porgy, Acanthopagrus schlegeli: A review in gonadal development, estradiol, estrogen receptor, aromatase activity and gonadotropin. J. exp. Zool., 290: 715-726. https://doi.org/10.1002/jez.1122

Li, L., Stoeckert, C.J. and Roos, D.S., 2003. OrthoMCL: Identification of ortholog groups for eukaryotic genomes. Genome Res., 13: 2178-2189. https://doi.org/10.1101/gr.1224503

Liu, Y.G., Liu, L.X., Wu, Z.X., Lin, H., Li, B.F. and Sun, X.Q., 2007. Isolation and characterization of polymorphic microsatellite loci in black sea bream (Acanthopagrus schlegeli) by cross-species amplification with six species of the Sparidae family. Aquat. Living Resour., 20: 257-262. https://doi.org/10.1051/alr:2007038

Ma, X.L., Xie, Z.Z., Zhou, L.B. and Chen, Y.Z., 2016. The complete mitochondrial genome of Acanthopagrus schlegelii (Perciformes: Sparidae) with phylogenetic consideration. Mitochond. DNA B, 1: 348-349. https://doi.org/10.1080/23802359.2015.1137839

Moriya, Y., Itoh, M., Okuda, S., Yoshizawa, A.C. and Kanehisa, M., 2007. KAAS: An automatic genome annotation and pathway reconstruction server. Nucl. Acids Res., 35: W182-W185. https://doi.org/10.1093/nar/gkm321

Nadia, M.D. and Alicia, O., 2014. Corset: Enabling differential gene expression analysis for de novo assembled transcriptomes. Gene Biol., 15: 410. https://doi.org/10.1186/PREACCEPT-2088857056122054

Rose, E., Flanagan, S.P. and Jones, A.G., 2015. The effects of synthetic estrogen exposure on the sexually dimorphic liver transcriptome of the sexrole-reversed Gulf Pipefish. PLoS One, 10: e0139401. https://doi.org/10.1371/journal.pone.0139401

Stamatakis, A., 2014. RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics, 30: 1312-1313. https://doi.org/10.1093/bioinformatics/btu033

Tirosh, I., Weinberger, A., Carmi, M. and Barkai, N., 2006. A genetic signature of interspecies variations in gene expression. Nat. Genet., 38: 830–834. https://doi.org/10.1038/ng1819

Wu, G.C., Du, J.L., Lee, Y.H., Lee, M.F. and Chang, C.F., 2005. Current status of genetic and endocrine factors in the sex change of protandrous black porgy, Acanthopagrus schlegeli (Teleostean). Annls N.Y. Acad. Sci., 1040: 206-214. https://doi.org/10.1196/annals.1327.026

Wu, G.C., Tomy, S., Lee, M.F., Lee, Y.H., Yueh, W.S., Lin, C.J., Lan, E.L. and Chang, C.F., 2010. Sex differentiation and sex change in the protandrous black porgy, Acanthopagrus schlegeli. Gen. Comp. Endocrinol., 167: 417-421. https://doi.org/10.1016/j.ygcen.2009.11.003

Yang, Z. and Rannala, B., 2012. Molecular phylogenetics: Principles and practice. Nat. Rev. Genet., 13: 303-314. https://doi.org/10.1038/nrg3186

Zhu, B.Q., 2017. Differences of morphological phenotypes, dorsal musculoskeletal textures and related physiological characteristics on juvebile Sparus macrocephalus with different anti-flowing property. Zhejiang Ocean University, (In Chinese).