Cynoglossus semilaevis used in this study were obtained from a commercial fish farm in Daishan (Zhoushan, China). Five fish (at 2.5 months post-fertilization (mpf), length 5–6 cm) were randomly chosen for this study. All fish sampled were anesthetized with tricaine methanesulfonate (MS-222, 150 mg/L) for subsequent sampling.

The muscle tissue pieces (length 1.5 cm) were dissected out of fish samples at the level of the anal opening and then fixed in 4% paraformaldehyde (PFA) for morphological observation.

Meanwhile, the eyed-side muscle (ESM) and blind-side muscle (BSM) from the same individual were collected for RNA-seq. All samples were frozen in liquid nitrogen, and stored at –80 °C for analysis.

The fixed muscle tissues were dehydrated through a series of graded ethanol and embedded in paraffin wax. Transverse serial wax sections, 4 μm thick, were cut transversely to the long body axis and stained with hematoxylin-eosin (H&E). Sections were mounted on a glass slide with a cover slip and neutral resin, and photographed under light microscopy with a 40× objective lens (Nikon Ni-U, Tokyo, Japan).

Three pictures were taken of each sample, and the muscle thickness in the picture was measured. A vertical dotted line was drawn from the spine toward the epidermis as shown in Fig. 1A. The distance was specified as the subscapularis muscle thickness value. The muscle thickness from two sides was analyzed in all samples. Morphological data of muscle thickness was shown by the mean of three pictures.

Total RNA was extracted using the TRIzol Reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocols (Rio et al. 2010). RNA quality and integrity were evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). Then, total RNA was subjected to next-generation sequencing using the DNBSEQ-T7 platform at the Novogene Bioinformatics Technology Co., Ltd. (Tianjin, China). The cDNA libraries of three ESM and three BSM groups were sequenced using the DNBSEQ-T7 platform with 150 bp paired-end reads.

Raw sequencing reads of six transcriptome datasets were cleaned by removing the adaptors and low-quality reads using Trimmomatic (v. 0.36) (Bolger et al. 2014). The obtained clean data were then quality-controlled by calculating the values of GC-content, Q20, and Q30 values. After quality control of each sample, the clean data were aligned to the reference gene of C. semilaevis (https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000523025.1/) using HISAT2 (v. 2.0.5) (Mortazavi et al. 2008)

The gene expression levels were calculated by fragments per kilobase of transcript per million reads (FPKMs) mapped (Mortazavi et al. 2008). Differentially expressed genes (DEGs) were identified utilizing the DESeq2 R package, with a significance P- value ≤ 0.05 and |log2 (FoldChange)| ≥ 1 (Love et al. 2014). Furthermore, the R package was utilized to perform the analysis and visualization of volcanic diagrams and heat maps of DEGs. Next, these identified DEGs were subjected to functional enrichment analysis, including gene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses via the R package clusterProfiler (v. 3.8.1) (Yu et al. 2012). GO terms and KEGG pathways with a P ≤ 0.05 were considered significantly enriched.

Alternative Splicing (AS) was conducted using rMATs (v. 4.1.0) software, including SE, MXE, A5SS, A3SS, and RI for each sample (Shen et al. 2014). For each AS event, rMATS calculated the percentage of exon inclusion (IncLevel) for each sample across the biological triplicates and detected differential IncLevel (IncLevel Difference) between the two groups. The differential alternative splicing (DAS) events were screened and categorized using summary.py in rMATS based on an IncLevelDifference absolute value greater than 0.1 and false discovery rate (FDR)<0.05.

To evaluate the reliability of the RNA-seq data, 9 DEGs were randomly chosen for RT-qPCR. 18S rRNA was chosen as the internal control. Total RNA was extracted as described above. The cDNA was synthesized from 1 μg of total RNA using an M-MLV reverse transcriptase (RNase H⁻) (TaKaRa Bio Inc., Japan). Primers were designed using Primer 5.0 software (Suppl. material 1). RT-qPCR experiments in a total volume of 12.5 μ contained 0.6 μL of cDNA, 0.25 μL of each primer, 6.25 μL of SYBR Green mix, and 5.15 μL of ddH₂O. The amplification program was 95 °C for 2 min, 40 cycles at 95 °C for 5 s, 60 °C for 30 s, 72 °C for 3s. The expression level of genes was calculated using the 2^{− ΔΔCt} method, and each sample was performed in triplicate (Livak and Schmittgen 2001).

All data are shown as the mean ± standard error (SE). Statistical analysis was performed with SPSS 26.0 (SPSS, USA). One-way analysis of variance (ANOVA) was used to compare if there were differences between samples. Significant difference was obtained at p < 0.05, which was indicated by an asterisk (*).

Transverse sections of C. semilaevis are shown in Fig. 1A. The muscle cells were arranged neatly, showed an obvious mosaic appearance, and the cell lines were visible. New muscle cell neogenesis (colored dark in Fig. 1A) is produced in the middle parts of the myotomes, neighboring to the horizontal septum, and extends to the end ventrally, which results in a mosaic appearance. In the same visual field, the number of new muscle cells in the ESM was significantly more than that in the BSM (Fig. 1A). Furthermore, the muscle thickness in the ESM (15.30 mm) was significantly greater than that in the BSM (13.33 mm) (Fig. 1B, P < 0.05). These results suggest that the mosaic hyperplastic growth of muscle cells may result in differences in muscle development and growth.

Figure 1. 

A histological transverse section of muscle tissue stained with haematoxylin and eosin in C. semilaevis. The arrow demonstrates the position of small new fibers. SK: skin; SC: Spinal cord; N: notochord; P: permysium; scale bar ​= ​50 ​μm B the muscle thickness of ESM and BSM. Data was shown as mean ± standard error (SE) (N = 5), “*” indicates significant difference (p < 0.05).

In this study, we constructed a total of six cDNA libraries from two groups, each replicated three times, with an average of 51.89 and 48.06 Mb raw reads for ESM and BSM, respectively. After quality filtering and assessment, an average of 50.52 million and 46.38 million clean reads were yielded, respectively. The GC percentage ranged from 51.77% to 52.62%, with Q20 > 98.38%, Q30 > 94.53%, and no more than 0.02% error rate (Table 1). In addition, over 94% of the clean reads (94.12% - 94.69%) were mapped to the reference genome of C. semilaevis. These results indicated that the sequencing data were reliable and could be used for subsequent analysis (Table 1). All raw data obtained in this study were submitted to the NCBI Sequence Read Archive database under the project (PRJNA1196394).

In total, 1177 DEGs were identified by comparing the ESM vs BSM group (Suppl. material 1), of which 291 were up-regulated and 886 were downregulated (Fig. 2A). We constructed a clustered heatmap to depict the expression profile of DEGs between ESM vs BSM group. The clustered heatmap showed that the expression patterns within the two groups were mainly divided into main clades, indicating that some DEGs may perform the same or a similar function in a pathway (Fig. 2B).

Figure 2. 

A volcano plot of DEGs between ESM and BSM groups (red dots represent up-regulated genes, blue dots represent down-regulated genes, and gray dots represent indistinguishable genes) B heatmap plot of DEGs in ESM and BSM groups of C. semilaevis. X-axis represents sample name. Y-axis represents the relative expression of DEGs. The different colors indicated changes in the relative expression of DEGs.

To explore the function of C. semilaevis muscle development-related genes and perform the potential functional annotation of DEGs, GO annotation and KEGG enrichment analyses were performed. The GO enrichment analysis showed that 1177 DEGs were enriched in 631 GO terms, and 43 terms were significantly enriched (Suppl. material 2). Fig. 3 shows the top 30 GO terms, of which 21 are associated with biological process (BP), two are related to cellular component (CC), and seven are related to molecular function (MF). These terms mainly affected peptidase inhibitor activity, peptidase regulator activity, carbohydrate metabolic process, endopeptidase inhibitor activity, etc.

Figure 3. 

The top most significantly enriched GO terms (top 30) of DEGs between ESM and BSM groups in C. semilaevis.

KEGG enrichment analysis indicated that DEGs were annotated into 121 KEGG pathways and 19 pathways were significantly enriched (Suppl. material: Suppl. material 3). As shown in Fig. 4, the main significantly enriched pathways were “cell adhesion molecules”, “Tight junction”, “Amino sugar and nucleotide sugar metabolism”, “ECM-receptor interaction” pathways, etc. Among all the metabolic pathways, glycolysis/gluconeogenesis may be a key regulation pathway in skeletal muscle energy metabolism, including 12 DEGs.

Figure 4. 

The top 20 enriched KEGG pathways of DEGs between ESM and BSM groups in C. semilaevis.

A total of 16728 AS genes were identified from ESM and BSM groups. The SE, A5SS, A3SS, MXE, and RI splicing types accounted for 71.71%, 6.79%, 12.27%, 5.39%, and 3.84% of all splicing events, respectively (Fig. 5A). After filtering (|IncLevelDifference| > 0.1 and FDR < 0.05), 167 differential AS (DAS) were identified between ESM and BSM. The SE, A5SS, A3SS, MXE, and RI splicing types accounted for 56.29%, 10.18%, 15.57%, 8.98%, and 8.98% of all splicing events, respectively (Fig. 5B).

Figure 5. 

Type distribution of AS events (A) and DAS (B) in ESM and BSM of C. semilaevis.

Nine DEGs were randomly selected for qRT-PCR to verify the reliability of the results from the transcriptome sequencing. As shown in Fig. 6, the relative expression patterns of these genes revealed by qPCR were similar to the RNA-seq results, indicating that the RNA-seq results were reliable.

Figure 6. 

Validation of expression patterns of 9 genes from RNA-seq data by RT-qPCR with 18S gene as an internal control. The values represent the gene expression level of ESM compared with the BSM group.

Muscle tissue is the most abundant tissue in the body mass and plays an important role in whole-body metabolism (Zhu et al. 2015). Teleost muscle usually grows in two ways: hypertrophy (increase in fiber size) and hyperplasia (increase in the number of muscle fibers) (Rowlerson and Veggetti 2001). Hyperplasia continues until the fish body grows to a certain stage, after which hypertrophy mainly dominates muscle growth (Johnston et al. 2003). The relative timing of the main hyperplastic growth process concerning the life cycle varies between species (Rowlerson and Veggetti 2001). A previous study found that in the common sole (Solea solea), the first hyperplastic phase starts 12 days post-hatching, and the second hyperplastic phase was already underway at 2.5 months, and by 1 year, the hyperplastic growth of muscle had ceased (Veggetti et al. 1999). In this study, small new fibers were found between the already existing fibers located in the most dorsal and ventral areas of lateral muscle (Fig. 1), which indicated that muscle growth mainly relied on the recruitment of new muscle fibers (hyperplasia) during the juvenile stage of C. semilaevis. The same situation was also observed in other species including, European plaice (Pleuronectes platessa) (Brooks and Johnston 1993), turbot (Scophthalmus maximus) (Gibson and Johnston 1995), blackspot seabream (Pagellus bogaraveo) (Silva et al. 2009). Additionally, the muscle thickness in the ESM was significantly greater than that of the BSM, which may be caused by the hyperplastic growth of muscle fibers on the eyed side.

This is the first time that muscle transcriptome has been compared between the eyed side and blind side of C. semilaevis using transcriptome sequencing to investigate the potential mechanism of muscle development. In this study, a total of 26,550 unique transcripts were generated and 1177 DEGs were identified (Fig. 2). Among these DEGs, some genes are related to muscle development, for example, myosin family genes and ubiquitin family genes. Myosin family protein is the most abundant and important protein in muscle tissue, which participates in a variety of cellular processes, including cytokinesis, cell polarization, signal transduction, intracellular transport, etc. (Furst et al. 1989; Whiting et al. 1989; Hofmann et al. 2009; Hartman et al. 2011). Studies in salmonids, grass fish, and hybrid grouper indicate that myofibrillar component genes concern growth (Devlin et al. 2009; Sun et al. 2016; Lu et al. 2020). In the present study, several genes related to skeletal muscle myofibrillar components, such as actin, myosin-binding protein C (mybp-C), titin (ttn), troponin (tn), tnnil, myoglobin (mb), and astrotactin-2 (astn-2) were highly expressed in the ESM group. Elevated expression levels of these genes were also found in larger rainbow trout (Salem et al. 2012; Kocmarek et al. 2014). Tn, which is related to the calcium signaling pathway in muscle cells, was up-regulated in the BSM, suggesting that this gene controlling calcium homeostasis also regulates muscle growth in fish (Baar and Esser 1999). It’s identical to the results in grass fish, which promote an increase in muscle tissue (Lu et al. 2020). In addition, mb was also found to be significantly up-regulated in ESM, which indicated that myoglobin directly facilitates aerobic exercise along with the growth of muscle through transporting the oxygen consumed during the process of respiration in muscle cell mitochondria (Wittenberg et al. 1975). However, we also observed up-regulation for different copies of myofibrillar components genes such as myosin light chain 4 (mlc4), myosin regulatory light chain 2 (myl2), troponin C (cTnC), alpha- actin (acta), troponin T-2C (tnt-2C) in ESM group, which were found in other fish (Kocmarek et al. 2014).

In addition, genes involved in the cell cycle (such as murine double minute 2 and cyclin A2) were also significantly up-regulated in the ESM. The organs’ growth depends on an increase both in biomass of individual cells and cell number (Sablowski and Carnier Dornelas 2014). Usually, the increase in cell number is accomplished by cell division (cell cycle progression), which plays an important role in organism development (Goranov et al. 2009). Cyclins are important regulatory proteins, and which bind to cyclin-dependent kinases (CDKs) to regulate the cell cycle progression (Darzynkiewicz et al. 1996). In different animals, cyclin has been demonstrated to play a role in the cellular context of terminally differentiated muscle (Rao and Kohtz 1995). In this study, cyclin A2 was significantly up-regulated in ESM. Li et al. (2024a) showed that cyclin was involved in Yangtze sturgeon (Acipenser dabryanus) growth difference as a potential key hub gene through WGCNA analysis. Therefore, the results suggested that the upregulation of genes related to cell cycle might be another endogenous factor that is conducive to different muscle thicknesses between the eyed-side and blind-side muscle in C. semilaevis.

Muscle growth reflects the balance between protein synthesis and degradation. The ubiquitin Proteasome Pathway (UPP) is the main system for degrading unnecessary or damaged proteins (Lecker et al. 2006), and the E3-ubiquitin ligases are an important component of this system, which conjugate ubiquitin to the target protein, and then induce protein degradation (David et al. 2011). Furthermore, a large number of ubiquitin system genes is proven to be regulated in atrophying muscles (Lecker et al. 2006). In our study, 22 ubiquitin family genes were identified to be DEGs, 17 genes of which were up-regulated in BSM and 5 were up-regulated in ESM (Suppl. material 1). These results are in accordance with the finding that E3-ubiquitin ligases activate myogenesis in zebrafish. This suggests that ubiquitin system genes could affect muscle development through the ubiquitin-proteasome pathway to promote the degradation of protein in C. semilaevis.

Glycolysis is an important metabolic pathway in which glucose is converted into pyruvate producing ATP and NADH for cellular metabolic activity (Hers 1983; Muirhead and Watson 1992). This pathway produces some important precursor metabolites, and also plays a vital role in muscle tissue (Wang et al. 2021). Several studies have shown that glycolysis promotes muscle growth in zebrafish, pigs, and chickens (Lametsch et al. 2006; Teltathum and Mekchay 2009; Wang et al. 2021). Regarding C. semilaevis, a previous study indicated that glycolysis/gluconeogenesis contributed to flatfish size dimorphism (Wang et al. 2021). In our study, two important genes in glycolysis (acetyl-CoA synthetase and glucose 6-phosphatase) were more highly expressed in eyed-side muscle tissue than in blind-side muscle tissue. These two DEGs, acs and g6pase, are the key enzymes playing important roles in the regulation of glycolytic/gluconeogenic pathways (Werve et al. 2011; James et al. 2016). The high expression levels of acs and g6pase in the ESM may contribute to the greater muscle thickness in the eyed side of C. semilaevis. However, the precise molecular mechanisms of glycolysis/gluconeogenesis and how they participate in muscle development in C. semilaevis require further study by the functional analysis of DEGs in this pathway.

The authors have declared that no competing interests exist.

No ethical statement was reported.

This study was supported by “the Fundamental Research Funds for Zhejiang Provincial Universities and Research Institutes” (2024J005-2), and Natural Science Foundation of China (42171069).

Zhenming Lü and Liqin Liu planned and designed the research, and afforded the funding. Yuzhen Wang, Jing Yu, Yijing Yang, and An Xu collected samples, performed experiments, and analyzed the data. Yuzhen Wang and Liqin Liu wrote the original draft. Li Gong, Jing Liu, and Fenghui Li reviewed and edited the manuscript.

Jing Yu https://orcid.org/0009-0004-4611-6102

All of the data that support the findings of this study are available in the main text or Supplementary Information.