Whole-Transcriptome RNA Sequencing Uncovers the Global Expression Changes and RNA Regulatory Networks in Duck Embryonic Myogenesis

Abstract

Duck meat is pivotal in providing high-quality protein for human nutrition, underscoring the importance of studying duck myogenesis. The regulatory mechanisms governing duck myogenesis involve both coding and non-coding RNAs, yet their specific expression patterns and molecular mechanisms remain elusive. To address this knowledge gap, we performed expression profiling analyses of mRNAs, lncRNAs, circRNAs, and miRNAs involved in duck myogenesis using whole-transcriptome RNA-seq. Our analysis identified 1733 differentially expressed (DE)-mRNAs, 1116 DE-lncRNAs, 54 DE-circRNAs, and 174 DE-miRNAs when comparing myoblasts and myotubes. A GO analysis highlighted the enrichment of DE molecules in the extracellular region, protein binding, and exocyst. A KEGG analysis pinpointed pathways related to ferroptosis, PPAR signaling, nitrogen metabolism, cell cycle, cardiac muscle contraction, glycerolipid metabolism, and actin cytoskeleton. A total of 51 trans-acting lncRNAs, including ENSAPLT00020002101 and ENSAPLT00020012069, were predicted to participate in regulating myoblast proliferation and differentiation. Based on the ceRNAs, we constructed lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA ceRNA networks involving five miRNAs (miR-129-5p, miR-133a-5p, miR-22-3p, miR-27b-3p, and let-7b-5p) that are relevant to myogenesis. Furthermore, the GO and KEGG analyses of the DE-mRNAs within the ceRNA network underscored the significant enrichment of the glycerolipid metabolism pathway. We identified five different DE-mRNAs, specifically ENSAPLG00020001677, ENSAPLG00020002183, ENSAPLG00020005019, ENSAPLG00020010497, and ENSAPLG00020017682, as potential target genes that are crucial for myogenesis in the context of glycerolipid metabolism. These five mRNAs are integral to ceRNA networks, with miR-107_R-2 and miR-1260 emerging as key regulators. In summary, this study provides a valuable resource elucidating the intricate interplay of mRNA-lncRNA-circRNA-miRNA in duck myogenesis, shedding light on the molecular mechanisms that govern this critical biological process.

Keywords: ceRNA regulatory network; duck myogenesis; non-coding RNAs; whole-transcriptome sequencing.

Figure 1

Duck primary embryonic myogenesis. (A) Flow chart of the experimental procedure. (B) Proliferating myoblasts. (C) Differentiating myoblasts on day 4 (myotubes). (D) Immunofluorescence analysis of undifferentiated duck myoblasts was maintained in GM and stained with anti-desmin. (E) Immunofluorescence analysis of differentiating myoblasts on day 4 stained with anti-MyHC. Relative expression of muscle markers, including MYOG (F), MYOD (G), and MYF5 (H), during in vitro differentiation of duck primary myoblasts. RNA was isolated on days 0 (GM), 1 (DM1), 2 (DM2), 3 (DM3), and 4 (DM4) of myoblast differentiation culture. Data (4 biological replicates) are presented as means ± s.e.m. Different letters between two groups represent significant differences (p < 0.05).

Figure 1

Duck primary embryonic myogenesis. (A) Flow chart of the experimental procedure. (B) Proliferating myoblasts. (C) Differentiating myoblasts on day 4 (myotubes). (D) Immunofluorescence analysis of undifferentiated duck myoblasts was maintained in GM and stained with anti-desmin. (E) Immunofluorescence analysis of differentiating myoblasts on day 4 stained with anti-MyHC. Relative expression of muscle markers, including MYOG (F), MYOD (G), and MYF5 (H), during in vitro differentiation of duck primary myoblasts. RNA was isolated on days 0 (GM), 1 (DM1), 2 (DM2), 3 (DM3), and 4 (DM4) of myoblast differentiation culture. Data (4 biological replicates) are presented as means ± s.e.m. Different letters between two groups represent significant differences (p < 0.05).

Figure 2

Features of duck lncRNA, mRNA, and circRNA. (A) The positional classification and proportion of lncRNAs. A transcript with one of the class codes, ‘i, j, o, u, and x’, was defined as a lncRNA transcript. (i) Class ‘i’ refers to a transcribed fragment, which could be either in the sense or anti-sense orientation, that is entirely contained within a reference intron. (j) Class ‘j’ pertains to transcripts that exhibit at least one splicing junction shared with the reference transcript. (u) Class ‘u’ designates a lncRNA as an intergenic transcript of unknown function. (o) Class ‘o’ is applied when an ordinary exon of a predicted lncRNA partially overlaps with a reference transcript. (B) The type and proportion of circRNAs. (circRNA) exon-derived circular RNA. (ciRNA) intron-derived circular RNA. (intergenic) intergenic-derived circular RNA. (C) Distribution of exon numbers for lncRNA and mRNA. (D) Transcript lengths of lncRNA and mRNA. Lengths of ORFs for lncRNA (E) and mRNA (F). nt: nucleotides; aa: amino acids.

Figure 3

Differential expression analysis of mRNAs, circRNAs, and lncRNAs during duck primary myoblast differentiation. Volcano plots of gene expression levels of all mRNAs (A), lncRNAs (B), and circRNAs (C). The vertical dotted lines indicate |log2FC| = 1, and the horizontal dotted lines indicate p value = 0.05. GO (D) and KEGG (G) analyses of DE-mRNAs. GO (E) and KEGG (H) analyses of DE-lncRNAs. GO (F) and KEGG (I) analyses of DE-circRNAs.

Figure 4

A regulatory network of trans-acting DE-lncRNAs and DE-mRNAs. Teal circles represent lncRNAs, and orange octagons represent mRNAs. The thickness of lines indicates the magnitude of free energy; the smaller lines indicate less free energy, and the thicker lines indicate more free energy. Light yellow lines indicate a positive correlation, and light-teal lines indicate a negative correlation.

Figure 5

Expression profiling and differential expression of miRNAs in duck myogenesis. (A) The type and proportion of detected small RNAs. (B) Length distribution of miRNA reads. (C) Venn diagrams of detected miRNAs between DM and GM. (D) The number of DE-miRNAs in DM and GM. (E) Volcano plot for miRNA expression. The vertical dotted lines indicate |log2FC| = 1, and the horizontal dotted lines indicate p value = 0.05. (F) Top GO terms of DE-miRNAs. (G) Top 10 KEGG pathways of DE-miRNAs. nt: nucleotides.

Figure 6

ceRNA co-regulatory network of DE-mRNAs, DE-lncRNAs, DE-circRNAs, and DE-miRNAs. (A) Co-regulatory network of lncRNA-miRNA-mRNA. (B) Co-regulatory network of circRNA-miRNA-mRNA. The triangle represents miRNAs. Rhombus represents mRNAs. The circle represents lncRNAs or circRNAs. Up-regulation is shown in red and down-regulation is shown in blue.

Figure 7

Functional enrichment analysis of mRNAs involved in the ceRNA network and the candidate ceRNA co-regulation network. (A) Top 10 GO terms of DE-mRNAs involved in the ceRNA network. (B) Top 10 KEGG pathways of DE-mRNAs involved in the ceRNA network. (C) LncRNA/circRNA-miRNA-mRNA pathway regulatory network. The triangle represents miRNAs. Rhombus represents lncRNAs. The square represents mRNAs. The circle represents circRNAs. The hexagon represents a pathway. Red nodes indicate up-regulation, and blue nodes indicate down-regulation.

Figure 8

qRT-PCR validation of DE-mRNAs (A), DE-circRNAs (B) DE-lncRNAs (C), and DE-miRNAs (D). FPKM: fragment per kilobase of transcript per million mapped reads. SRPBM: spliced reads per billion mapping. NORM: global mean normalizers for miRNA quantification. Data (3 biological replicates) are presented as means ± s.e.m. * p < 0.05, ** p < 0.01, and *** p < 0.001.