Histomorphological features of the hindlimb thigh muscles of O. tormota at different developmental stages
From G26 to G36, the hindlimb thigh muscles consisted of myogenous cells (Fig. 1a). At G36, 11 groups of pre-myogenic masses of cells and a small number of long cylindrical multinuclear myotubes had formed (Fig. 1b). From G37 to G41, the hindlimb thigh muscles consisted mainly of myotubes (Fig. 1c, d) and the nuclei of multinucleated myotubes moved from the center to the edge until G41 (Fig. 1d). At G41, 11 muscle groups had formed (Fig. 1e). From G42 to 10 months old, the hindlimb thigh muscles consisted mainly of myotubes and myofibers. At G42, six muscle groups in the females (vastus medialis, gracilis minor, gracilis major, adductor magnus, musculus semimembranosus, and sartorius) and four muscle groups in the males (vastus medialis, gracilis minor, gracilis major, and adductor magnus) had differentiated into primary myofibers (Fig. 1f). At G45, a few secondary myofibers had formed around the primary myofibers (Fig. 1g). At G46, the male and female froglets both had a new muscle group that had differentiated into myofibers (biceps femoris for males; vastus lateralis for females) and the fusion of primary and secondary myofibers to form myofibers with two nuclei began to promote myofiber hypertrophy (Fig. 1h). The number of muscle groups with differentiated myofibers in one-month-old froglets after metamorphosis did not change compared with G46. These myofibers had 1 or 2 nuclei (Fig. 2a). Myofibers with three and four nuclei appeared at the 3-month-old stage (Fig. 2b) and these myofibers had increased hypertrophy. Exception for semitendinosus (in the form of tendons) and adductor longus, in the 5-month-old froglets, the other nine muscle groups differentiated into myofibers (Fig. 2c). In the 10-month-old froglets hindlimb thigh muscles, six muscle groups in the females (biceps femoris muscle, gracilis minor, gracilis major, adductor magnus muscle, musculi semimenbranosus and sartorius) and five muscle groups in the males (gracilis minor, gracilis major, adductor magnus muscle, musculi semimenbranosus and sartorius) completely developed into myofibers (Fig. 2d). At 14 months old, all these muscle groups existed as multinucleated myofibers (Fig. 2e–h). From the 14-month to 4-year-old stages, these myofibers underwent further hypertrophy. Comparative analysis of the number of muscle fibers in each muscle group at the 3-month-old (pre-hibernation) and 5-month-old (period of hibernation) stages showed that the number of muscle fibers decreased in three muscle groups in the females (gracilis major, musculus semimembranosus, and gracilis minor) and six muscle groups in the males (vastus lateralis, musculus semimembranosus, gracilis minor, gracilis major, sartorius, and rectus femoris) (Fig. 3). Histomorphological features for other developmental stages were presented in Fig. S1.
Transcriptome sequencing, de novo assembly, and gene annotation
Samples with obvious differences in histomorphological features and unique life history (metamorphosis and hibernation) were selected to screen for genes and pathways related to development and adaptation to the unique life history. The stages for selected samples were as following, G36 (myogenic cell state), G40 (myotube state), G42 (appearing primary myofibers, early metamorphosis), G45 (appearing secondary myofibers, late metamorphosis), 3-month-old (appearing 3 and 4 nuclear myofibers, pre-overwinter), 5-month-old (9 muscle groups differentiating into myofibers, overwinter), 14-month-old (all these muscle groups completely developing into multinucleated myofibers), and 2-year-old (myofibers undergoing further hypertrophy, adult). The correlation analysis of gene expression levels in the 18 hindlimb thigh muscle samples showed that the correlation of biological duplication was very high (R = 0.96–1.00) (Fig. S2). After filtering out the low-quality reads and trimming the adaptor sequences, we obtained 846,522,958 valid reads from the 18 libraries (Table S1). Quality control checks for each library showed that the clean sequences were of high quality, with Q30 values > 92.65% for all libraries. A total of 35,974 genes were identified with lengths of 201–74,670 bp (average length 403 bp) and N50 of 1151 bp. All genes were annotated against the GO, KEGG, Pfam, Swiss-Prot, eggNOG, and NR databases, which returned 9668 (26.87%), 9104 (25.31%), 8884 (24.7%), 9092 (25.27%), 10,582 (29.42%), and 12,388 (34.44%) matches, respectively.
Time series analysis of gene expression
Samples at eight different developmental stages formed seven time series profiles of gene expressions (Fig. 4). The GO functional enrichment analysis revealed significant differences among the different profiles (Table S2). Profile 2 contained many terms related to muscle growth, development, and movement, including structural constituent of muscle, myofibril assembly, skeletal muscle thin filament assembly, muscle fiber development, locomotory behavior and muscle contraction. The different profiles may provide new molecular characteristics related to muscle development and its adaptation to life history.
Analysis of differentially expressed genes (DEGs)
A total of 11,482, 8352, 1254, 1911, 6865, 9561, and 8032 DEGs were detected in G36 vs G40, G40 vs G42, G42 vs G45, G45 vs 3- month -old, 3-month-old vs 5-month-old, 5-month-old vs 14-month-old, and 14-month-old vs 2-year-old adult comparisons, respectively (fold change ≥2; adjusted P value < 0.05) (Table S3).
Among the DEGs in the G36 vs G40 comparison, 280 had significantly enriched GO terms (P < 0.05), including cell proliferation, nervous system development, immune system process, skeletal muscle fiber development, cell division, hematopoietic progenitor cell differentiation, and chondrocyte differentiation (Table S4), among which 18 terms were related directly to muscle growth and development (Table S5). We screened out 141 DEGs that may be related to muscle growth and exercise, including MyoD, MyoG, MYF6, KLHL40, KLHL41, ALDOA, OBSCN, and TPM3 (Table S6). 11,482 DEGs were significantly enriched in 49 KEGG pathways (P < 0.05), including tight junction, focal adhesion, notch signaling pathway, oxytocin signaling pathway, and ECM–receptor interaction. (Fig. 5a).
Among the DEGs in the G40 vs G42 comparison, 181 had significantly enriched GO terms (P < 0.05), including skeletal muscle contraction, myosin filament, motor activity, regulation of myoblast differentiation, glycogen metabolic process, and glycogen biosynthetic process (Table S4), among which 21 terms were related directly to muscle growth, development, and exercise (Table S5). We screened out 149 DEGs that may be related to muscle growth and exercise, including MYF6, MEF2D, MEF2C, MYH8, MYH4, MYH3, MAIA, TNNI1, and MYOZ3 (Table S6), and 36 DEGs that may be related to glycogen metabolism, glycolysis process, and energy metabolism, including PGM5, PYGB, PYGM, PRKAG2, PRKAG3, AMPD1, and OXCT1. Thirty-six of these DEGs were up-regulated in G42 (Table S7). 8352 DEGs were significantly enriched in 46 KEGG pathways (P < 0.05), including tight junction, insulin signaling pathway, MAPK signaling pathway, ECM–receptor interaction, focal adhesion, and hedgehog signaling pathway (Fig. 5b).
Among the DEGs in the G42 vs G45 comparison, 178 had significantly enriched GO terms (P < 0.05), including myosin filament, cardiac muscle hypertrophy, skeletal muscle thin filament assembly, negative regulation of myoblast differentiation, and cardiac myofibril assembly (Table S4), among which 40 terms were related directly to muscle growth and development (Table S5). We screened out 55 DEGs that may be related to muscle growth and exercise, including MYH2, MYH4, LMOD3, MYH7B, MYH1, MYH8, MYL4, and TTN (Table S6).
Among the DEGs in the G45 vs 3-month-old comparison, 226 had significantly enriched GO terms (P < 0.05), including glucocorticoid receptor binding, glycogen metabolic process, skeletal muscle thin filament assembly, skeletal muscle myosin thick filament assembly, and skeletal muscle fiber development (Table S4), among which 38 terms were directly related to muscle growth and development (Table S5). We screened out 49 DEGs that may be related to muscle growth and exercise, including MA1A, TTN, ACTA1, MYH4, MYLK2, PRM, SLC8A3, and XK70A (Table S6), and 10 DEGs related to lipid metabolism, namely CEBPB, NOCT, PNOLA2, PLIN2, GM2A, PDK4, ACACB, MLYCD, PRKAG2, and HPGDS. Except for HPGDS, which was up-regulated at 3-month-old, all the other genes were up-regulated at G45 (Table S8). We also screened out 10 DEGs related to glucose metabolism, namely PRKAG2, PPPLR3C, PPPLR3CB, CAV3, PDK4, SLC8B1, SLC2A4, NCOA5, BAD, and PCK2. Except for NCOA5, which was down-regulated at 3-month-old, all the other genes were up-regulated at G45 (Table S8). 1911 DEGs were significantly enriched in 37 KEGG pathways (P < 0.05), including thyroid hormone signaling pathway, P53 signaling pathway, insulin signaling pathway, focal adhesion, FoxO signaling pathway, PI3K-Akt signaling pathway, MAPK signaling pathway, PPAR signaling pathway, and adipocytokine signaling pathway (Fig. 5c).
Among the DEGs in the 3-month-old vs 5-month-old comparison, 255 had significantly enriched GO terms (P < 0.05), including ATP metabolic process, muscle cell development, lysosome, autophagy, ubiquitin ligase complex, glycogen catabolic process, glucose metabolic process, and extrinsic apoptotic signaling pathway (Table S4), among which 47 terms were related to directly muscle growth and development (Table S5). We screened out 144 DEGs that may be related to muscle growth and exercise, including MYH8, MYH4, MYH6, TWF2-A, SCO1, FLOT2, TMOD4, MSTN, and PDGFRA (Table S6). We detected 21 DEGs related to glycogen synthesis and glucose metabolism that were down-regulated at 5-month-old (Table S9) and 76 DEGs related to mitochondria and energy metabolism, 74 of which were up-regulated at 3-month-old; the exceptions were KDM5C and HSD17B14, which were up-regulated in 5-month-old (Table S9). We also found 90 DEGs related to protein degradation in lysosome and ubiquitination pathways, most of which were up-regulated at 5-month-old; the exceptions were CTSK, CAT-1, RNF7, ASB15, ANAPC16, and UBE2C, which were up-regulated at 3-month-old (Table S9). 6865 DEGs were significantly enriched in 30 KEGG pathways (P < 0.05), including ECM–receptor interaction, focal adhesion, apoptosis, tight junction, and PI3K-Akt signaling pathway (Fig. 5d).
Among the DEGs in the 5-month-old vs 14-month-old comparison, 375 had significantly enriched GO terms (P < 0.05), including skeletal muscle contraction, skeletal muscle cell differentiation, skeletal muscle thin filament assembly, and positive regulation of myoblast fusion (Table S4), among which 44 terms may be related to muscle growth and development (Table S5). We screened out 215 DEGs that may be related to muscle development and exercise, including MSTN, MYBPC2, IGFBP4, IGFBPL1, MYOG, MYF6, MYF5, and PDLIM3 (Table S6). 9561 DEGs were significantly enriched in 34 KEGG pathways (P < 0.05), including ECM–receptor interaction, focal adhesion, tight junction, PI3K-Akt signaling pathway, and p53 signaling pathway (Fig. 5e).
Among the DEGs in the 14-month-old vs 2-year-old adult comparison, 429 had significantly enriched GO terms (P < 0.05), including muscle contraction, skeletal muscle myosin thick filament assembly, myofibril assembly, adult walking behavior, and muscle fiber development (Table S4), among which 55 terms were related directly to muscle growth, development, and exercise (Table S5). We screened out 192 DEGs that may be related to muscle development and exercise, including CAPZA1, ACTN4, SYNE2, SlC4A7, MYL12A, MEF2D, FBXO40, and MSTN (Table S6). 8032 DEGs were significantly enriched in 32 KEGG pathways (P < 0.05), including insulin signaling pathway, glycerophospholipid metabolism, PPAR signaling pathway, FoxO signaling pathway, focal adhesion, and hippo signaling pathway-fly (Fig. 5f).
RT-qPCR verification of randomly selected DEGs
To confirm the RNA-seq results, we randomly selected five DEGs (MYH4, MYH8, MYLPF, MYF6, and MYF5) related to muscle development and growth and one DEG (DNMT1) related to DNA methylation for verification by RT-qPCR. The RT- qPCR and RNA-seq data showed the selected genes had similar expression trends. Relative expression changes of qRT-PCR data were highly (r = 0.87–0.99) correlated with RNA-Seq data (Fig. 6), suggesting the reliability of the RNA-Seq method.