Abstract

Background

   Circular RNAs (circRNAs) function as essential regulatory elements with
   pivotal roles in various biological processes. However, their
   expression profiles and functional regulation during the
   differentiation of goat myoblasts have not been thoroughly explored.
   This study conducts an analysis of circRNA expression profiles during
   the proliferation phase (cultured in growth medium, GM) and
   differentiation phase (cultured in differentiation medium, DM1/DM5) of
   skeletal muscle satellite cells (MuSCs) in goats.

Results

   A total of 2,094 circRNAs were identified, among which 84 were
   differentially expressed as determined by pairwise comparisons across
   three distinct groups. Validation of the expression levels of six
   randomly selected circRNAs was performed using reverse transcription
   PCR (RT-PCR) and quantitative RT-PCR (qRT-PCR), with confirmation of
   their back-splicing junction sites. Enrichment analysis of the host
   genes associated with differentially expressed circRNAs (DEcircRNAs)
   indicated significant involvement in biological processes such as
   muscle contraction, muscle hypertrophy, and muscle tissue development.
   Additionally, these host genes were implicated in key signaling
   pathways, including Hippo, TGF-beta, and MAPK pathways. Subsequently,
   employing Cytoscape, we developed a circRNA-miRNA interaction network
   to elucidate the complex regulatory mechanisms underlying goat muscle
   development, encompassing 21 circRNAs and 47 miRNAs. Functional assays
   demonstrated that circTGFβ2 enhances myogenic differentiation in goats,
   potentially through a miRNA sponge mechanism.

Conclusion

   In conclusion, we identified the genome-wide expression profiles of
   circRNAs in goat MuSCs during both proliferation and differentiation
   phases, and established that circTGFβ2 plays a role in the regulation
   of myogenesis. This study offers a significant resource for the
   advanced exploration of the biological functions and mechanisms of
   circRNAs in the myogenesis of goats.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12864-024-11008-4.

   Keywords: Goat, circRNA sequencing, Differential expression, circTGFβ2,
   Differentiation

Introduction

   Skeletal muscle constitutes the primary component of food animals,
   directly influencing both the quantity and quality of meat
   [[42]1–[43]3]. Consequently, the investigation of skeletal muscle
   growth and development holds significant importance for breeding
   practices within the field of animal husbandry. Nonetheless, the
   process of skeletal muscle growth and development is intricate,
   encompassing the activation of skeletal muscle satellite cells, the
   proliferation and differentiation of myoblasts, the fusion of myoblasts
   into myofibers, and the expression of a diverse array of genes [[44]4,
   [45]5]. The introduction of high-throughput sequencing technologies has
   facilitated the identification of an expanding array of noncoding RNAs
   involved in the regulation of muscle growth and development, including
   microRNAs (miRNAs), long noncoding RNAs (lncRNAs), and circular RNAs
   (circRNAs) [[46]6–[47]9]. Importantly, these noncoding RNAs have the
   potential to exert both direct and indirect effects on the process of
   myogenesis [[48]10, [49]11], although the specific regulatory
   mechanisms remain largely undefined. As a result, a considerable number
   of noncoding RNAs require further investigation.

   CircRNAs represent a recently identified class of single-stranded
   non-coding RNAs characterized by a covalently closed loop structure
   generated through back-splicing [[50]12, [51]13]. These circRNAs can
   originate from various genomic regions, including exons, introns,
   untranslated regions, non-coding RNA loci, as well as intergenic and
   antisense transcripts [[52]14]. Compared to their linear RNA
   counterparts, circRNAs exhibit enhanced structural stability and
   resistance to ribonuclease R (RNase R) [[53]14]. CircRNAs are
   ubiquitously present in cells and tissues and play a crucial role in
   regulating gene expression through diverse mechanisms. These mechanisms
   include functioning as miRNA sponges, modulating RNA-binding proteins,
   and serving as templates for protein translation [[54]15–[55]17].
   Emerging evidence underscores the significant role of circRNAs in the
   regulation of muscle growth and development, primarily through their
   function as molecular sponges for miRNAs [[56]18, [57]19]. For example,
   circSGCB has been shown to inhibit myoblast proliferation and promote
   differentiation by sequestering bta-miR-27a-3p [[58]20]. Similarly,
   circNDST1 acts as a competitive endogenous RNA, modulating the
   proliferation and differentiation of bovine myoblasts via the
   miR-411a/Smad4 axis [[59]21]. Additionally, circLRRFIP1 serves as a
   sponge for the miR-15 family, facilitating the proliferation and
   differentiation of chicken skeletal muscle satellite cells [[60]22]. In
   addition to the well-documented sponge effect, alternative regulatory
   mechanisms of circRNAs in myogenesis have been proposed. For instance,
   the protein encoded by circKANSL1L has been shown to interact with Akt,
   thereby enhancing FoxO3 transcriptional activity to regulate skeletal
   myogenesis [[61]23]. Furthermore, another study suggests that the
   binding of PUR proteins to circSamd4 facilitates muscle development by
   inhibiting MHC transcription [[62]24]. Despite these advances, the
   dynamic expression of circRNAs throughout the formation of goat
   skeletal muscle are still unclear, and their functions in myogenesis
   deserve further exploration.

   In this study, we performed circRNA sequencing on skeletal muscle
   satellite cells cultured under both proliferative and differentiative
   conditions to explore the roles of circRNAs in the development of goat
   skeletal muscle. We identified differentially expressed circRNAs and
   subsequently constructed a DEcircRNA-miRNA interaction network to
   elucidate the miRNA-sponge functions of circRNAs. Notably, we
   discovered a novel circRNA, circTGFβ2, which exhibited a skeletal
   muscle-specific expression pattern and facilitated the differentiation
   of goat skeletal muscle satellite cells. These findings will facilitate
   the compilation of a comprehensive catalog of circRNAs involved in goat
   muscle development and enhance our understanding of the regulatory
   mechanisms underlying the development of goat skeletal muscle.

Results

Overview of circRNA sequencing data

   To identify potential circRNAs involved in the regulation of skeletal
   myogenesis in goats, we conducted RNA sequencing to profile circRNA
   expression in skeletal muscle satellite cells (MuSCs) across different
   stages (GM, DM1 and DM5). The flowchart for this study is shown in
   Fig. [63]1A. A total of 721,950,292 raw reads were generated from nine
   samples. Post quality filtering, which involved the removal of
   low-quality reads and adaptor sequences, yielded 720,121,930 clean
   reads that were subsequently aligned to the goat genome (Supplementary
   Table [64]1). Quality control metrics indicated that each sample
   achieved a Q20 score of ≥ 97.82% and a Q30 score of ≥ 93.74%, with an
   average GC content of 50.29% (Supplementary Table [65]1). Alignment to
   the goat reference genome ARS1 revealed that over 90% of the reads were
   uniquely mapped. A total of 2,094 circRNAs were identified, including
   4, 9, and 15 circRNAs specifically expressed in GM, DM1, and DM5,
   respectively, along with 1,465 circRNAs common to all three phases
   (Fig. [66]1B) (Supplementary Table [67]2). To ensure the
   reproducibility of the sequencing data, three replicates of each sample
   were collected. Principal Component Analysis (PCA) revealed high sample
   similarity (Fig. [68]1C), indicating the reliability of the sequencing
   data. Subsequent statistical analysis demonstrated that the majority of
   circRNA types were exon circRNAs, with only a minor proportion
   comprising intron, antisense, and intergenic region circRNAs
   (Fig. [69]1D). Taken together, these results confirmed that the sample
   and sequencing data generated in this study were of good quality and
   were reliable for subsequent bioinformatic analysis.

Fig. 1.

   [70]Fig. 1
   [71]Open in a new tab

   Identification and characteristics of circRNAs in goat skeletal muscle
   satellite cells (MuSCs). (A) Schematic illustration of the experimental
   design. (B) Venn analysis of circRNAs detected at each time point. (C)
   Principal component analysis based on the expression of circRNAs. (D)
   Types of all circRNAs

Identification of differentially expressed circRNAs

   To elucidate the differences in circRNA expression patterns during
   myogenic differentiation in goats, the edgeR package was employed to
   identify DEcircRNAs across various sample pairs. A total of 84
   DEcircRNAs were identified across the three stages, with 25, 55, and 22
   DEcircRNAs detected in the GM-vs-DM1, GM-vs-DM5, and DM1-vs-DM5
   comparison groups, respectively (Fig. [72]2A) (Supplementary Table
   [73]3). Further analysis using a Venn diagram revealed 14, 38, and 14
   stage-specific DEcircRNAs in the GM-vs-DM1, GM-vs-DM5, and DM1-vs-DM5
   comparisons, respectively (Fig. [74]2B). Notably, no DEcircRNAs were
   common across all three comparisons. Finally, hierarchical clustering
   analysis was used to assess the expression patterns of DEcircRNAs to
   explore the similarities and to compare the relationships between the
   different libraries, and it was found that samples within the same
   group were well clustered together (Fig. [75]2C).

Fig. 2.

   [76]Fig. 2
   [77]Open in a new tab

   Analysis of DEcircRNAs in goat MuSCs across different stages. (A)
   Numbers of up-regulated and down-regulated circRNAs in goat MuSCs
   across three stages. (B) Venn diagram showing the DEcircRNAs at the
   three comparisons. (C) Hierarchical clustering heat map of all
   DEcircRNAs. Data are expressed as RPM. Red represents relatively high
   expression, blue represents relatively low expression

Functional enrichment analysis of DEcircRNA host genes

   To investigate the potential functions of the identified DEcircRNAs, GO
   and KEGG pathway enrichment analyses were conducted for the host genes
   of DEcircRNAs across the DM1-vs-DM5, GM-vs-DM1, and GM-vs-DM5
   comparison groups. The results indicated significant enrichment of
   DEcircRNAs in 603, 375, and 543 GO terms, respectively, encompassing a
   wide array of biological processes (Supplementary Table [78]4).
   Specifically, in the DM1-vs-DM5 comparison group, the top 10 enriched
   biological process terms, as illustrated in Fig. [79]3A, revealed that
   the host gene functions were predominantly associated with
   muscle-related biological processes, including the regulation of
   striated muscle contraction, regulation of skeletal muscle contraction,
   regulation of cell adhesion, muscle tissue development, and regulation
   of muscle system processes. In the GM-vs-DM1 group, the top 10 terms
   for biological processes included nucleic acid metabolic process,
   developmental process, and striated muscle contraction (Fig. [80]3B).
   Conversely, in the GM-vs-DM5 group, the identified host genes
   predominantly exhibited established roles in muscle biology,
   encompassing skeletal muscle contraction, musculoskeletal movement,
   striated muscle hypertrophy, muscle system processes, and cell
   differentiation (Fig. [81]3C).

Fig. 3.

   [82]Fig. 3
   [83]Open in a new tab

   Functional enrichment analysis for host genes of differentially
   expressed circRNAs. (A) The top 10 significance terms of biological
   process in DM1 vs. DM5. (B) The top 10 significance terms of biological
   process in GM vs. DM1. (C) The top 10 significance terms of biological
   process in GM vs. DM5. (D) The top 10 significance pathways in DM1 vs.
   DM5. (E) The top 10 significance pathways in GM vs. DM1. (F) The top 10
   significance pathways in GM vs. DM5

   To elucidate the functions of significantly enriched parental genes, a
   pathway analysis was performed using the KEGG pathway database
   (Supplementary Table [84]5). In the DM1-vs-DM5 comparison group, a
   substantial proportion of source genes for differentially expressed
   circRNAs (DEcircRNAs) exhibited significant enrichment in the Hippo
   signaling pathway and the MAPK signaling pathway (Fig. [85]3D), both of
   which are intimately linked to muscle development. Furthermore,
   pathways such as the adherens junction, TGF-beta signaling pathway, and
   FoxO signaling pathway, which rank among the top ten identified
   pathways, also play essential roles in muscle development
   (Fig. [86]3D). In a similar manner, within the GM-vs-DM1 group, several
   of the top ten pathways, such as the Hippo signaling pathway, Rap1
   signaling pathway, TGF-beta signaling pathway, adherens junction, and
   tight junction, are associated with muscle development (Fig. [87]3E).
   Additionally, the parental genes of DEcircRNAs in the GM-vs-DM5 group
   were significantly enriched in the thyroid hormone, TGF-beta, FoxO, and
   Hippo signaling pathways (Fig. [88]3F).

Construction of the circRNA-miRNA network

   More and more evidence suggests that circRNAs can act as miRNA sponges
   to competitively bind miRNA and regulate the expression of miRNA target
   genes, thereby participating in related biological processes
   [[89]25–[90]29]. To predict interactions between circRNAs and miRNAs,
   we utilized the miRanda and TargetScan algorithms. By integrating the
   interaction data with the differential expression profiles of circRNAs,
   we constructed interaction network data files, which were subsequently
   imported into Cytoscape software for visualization and analysis. Our
   analysis identified a total of 21 DEcircRNAs and 47 miRNAs, culminating
   in 71 interaction pairs between miRNAs and circRNAs (Fig. [91]4)
   (Supplementary Table [92]6). Notably, novel_circ_003559 was found to
   interact with 18 miRNAs, whereas novel_circ_002145 interacted with 14
   miRNAs. Moreover, specific miRNAs, such as chi-miR-7-5p,
   chi-miR-487a-5p, and chi-miR-130b-5p, were found to interact with
   multiple circRNAs. In addition, our analysis identified several miRNAs
   implicated in skeletal muscle development within this interaction
   network, including chi-miR-206, chi-miR-381, chi-miR-380-5p,
   chi-miR-16b-5p, and chi-miR-15b-5p [[93]30–[94]32]. These observations
   indicate that these circRNAs may exert a regulatory influence on
   skeletal muscle development through their binding to miRNAs.

Fig. 4.

   [95]Fig. 4
   [96]Open in a new tab

   The circRNA-miRNA interaction networks, including 21 DEcircRNAs and 47
   miRNAs. Blue represent DEcircRNAs and yellow represent miRNAs

Validation of DEcircRNAs

   To validate the reliability of RNA-Seq data, six highly expressed
   DEcircRNAs (novel_circ_000360, novel_circ_000805, novel_circ_002075,
   novel_circ_001349, novel_circ_003418, and novel_circ_004171) were
   selected for validation experiments, and divergent primers
   corresponding to circRNAs were designed (Supplementary Table [97]7).
   The PCR product was examined by sanger sequencing to confirm the
   back-splicing junction site (Fig. [98]5A). Subsequently, we performed
   reverse transcription polymerase chain reaction (RT-PCR) to evaluate
   the exonuclease resistance of circRNAs in comparison to linear RNA.
   Post RNase R digestion, circRNA remained detectable via RT-PCR, whereas
   the linear RNA, GAPDH, was rendered undetectable. These findings
   suggest that, relative to the control group, the six circRNAs
   demonstrated resistance to RNase R digestion as illustrated in
   Fig. [99]5B. Furthermore, divergent primers specific to each circRNA
   successfully amplified the anticipated fragments from complementary DNA
   (cDNA), but failed to produce PCR products when genomic DNA (gDNA) was
   used as a template (Fig. [100]5B), thereby confirming the presence of
   back-site junctions. Furthermore, we analyzed the expression patterns
   of these circRNAs across three phases of MuSCs using qRT-PCR. The
   results demonstrated that the expression profiles were consistent with
   the trends observed in the RNA-seq data (Fig. [101]5C). This
   concordance confirms the reliability of the RNA-seq data and the
   expression of the circRNAs identified in this study.

Fig. 5.

   [102]Fig. 5
   [103]Open in a new tab

   Experimental validation of circRNAs. (A) Sanger sequencing shows
   head-to-tail junctions. Blue arrow shows back spliced junctions. (B)
   Agarose gel electrophoresis test for PCR product of six circRNAs using
   different templates. gDNA represent genomic DNA; R + indicates the
   samples were treated with exonuclease. (C) Expression profiles of the
   six randomly selected circRNAs based on RNA-Seq and qRT-PCR

The expression and characteristics of circTGFβ2

   From the aforementioned network, we selected novel_circ_000975
   (designated as circTGFβ2) for further investigation, hypothesizing its
   potential role in regulating myoblast differentiation through the
   predicted circTGFβ2-chi-miR-206/chi-miR-211 pathway. CircTGFβ2 is an
   circRNA derived from exons 4 to 7 of the TGFβ2 gene. The back-splice
   junction sequence of circTGFβ2 was validated via Sanger sequencing
   (Fig. [104]6A). Subsequently, we performed RT-PCR using cDNA
   synthesized with either random primers or Oligo(dT) primers. The
   absence of bands in the Oligo(dT) group indicated that circTGFβ2 lacks
   a 3’ poly-A tail (Fig. [105]6B). Divergent primers successfully
   amplified the predicted fragments utilizing cDNA, but failed to
   generate PCR products when gDNA was employed as a template
   (Fig. [106]6B), thereby corroborating the presence of back-site
   junctions. Additionally, we assessed the exonuclease resistance of
   circTGFβ2 using RT-PCR and qRT-PCR methodologies. Treatment with RNase
   R demonstrated that circTGFβ2 exhibited resistance to RNase R-mediated
   degradation, in contrast to linear TGFβ2 mRNA, which was predominantly
   degraded (Fig. [107]6B and C). The expression of circTGFβ2 was detected
   in various tissues, with the highest expression levels observed in
   skeletal muscle tissues, followed by the heart and liver
   (Fig. [108]6D). We performed a comprehensive analysis of circTGFβ2
   expression during the differentiation of MuSCs, employing growth medium
   for 1 and 2 days, and differentiation medium for 1, 3, and 5 days. As
   depicted in Fig. [109]6E, the expression levels of circTGFβ2 escalated
   significantly during the differentiation phase, suggesting its pivotal
   role in muscle differentiation regulation. In addition, we investigated
   the subcellular localization of circTGFβ2 via RNA-FISH. The results
   indicated that circTGFβ2 predominantly resides in the cytoplasm
   (Fig. [110]6F), corroborating its hypothesized function as a miRNA
   sponge. A dual-luciferase reporter assay was conducted to confirm the
   interactions between circTGFβ2 and chi-miR-206/chi-miR-211. The
   findings indicated that chi-miR-206 and chi-miR-211 significantly
   decreased luciferase activity in comparison to the mutant group
   (Fig. [111]6G and H), thereby substantiating the binding relationship
   between chi-miR-206/chi-miR-211 and circTGFβ2.

Fig. 6.

   [112]Fig. 6
   [113]Open in a new tab

   The characterization and function of circTGFβ2 in MuSC differentiation.
   (A) The back-site junction of circTGFβ2 was identified using a
   divergent primer and sequenced by Sanger sequencing. (B) Agarose gel
   electrophoresis test for PCR product of circTGFβ2 using different
   templates. Random shows the cDNA template synthesized using random
   primers for RT-PCR, while Oligo(dT) denotes the cDNA synthesized using
   Oligo(dT) primers; gDNA represent genomic DNA; RNase R + indicates the
   samples were treated with exonuclease. (C) qRT-PCR analysis of the
   expression of circTGFβ2 and TGFβ2 in goat myoblasts treated with RNase
   R. (D) The abundance of circTGFβ2 in various tissues of goats on the
   third day after birth. (E) The expression level of circTGFβ2 during
   proliferation and differentiation of goat skeletal muscle satellite
   cells. (F) RNA FISH assay was performed to determine the subcellular
   localization of circTGFβ2 in goat skeletal muscle cells, blue indicates
   nuclei stained with DAPI; red indicates RNA probes recognizing
   circTGFβ2. Scale bar, 10 μm. (G) Luciferase reporter activity of
   circTGFβ2 in MuSCs co-transfected with chi-miR-206 mimics or mimics NC.
   (H) Luciferase reporter activity of circTGFβ2 in MuSCs co-transfected
   with chi-miR-211 mimics or mimics NC. (I) CircTGFβ2 was successfully
   overexpressed in MuSCs of goat. (J) CircTGFβ2 was successfully
   inhibited in MuSCs of goat. (K) The expression of TGFβ2 mRNA was
   detected by qRT-PCR after circTGFβ2 knockdown. (L) The expression of
   MyoD, MyoG, and MyHC was detected by qRT-PCR after transfection with
   pLV-circTGFβ2 and pLV-CiR vector. (M) The expression of MyoD, MyoG, and
   MyHC was detected by qRT-PCR after transfection with si-1, si-2, and
   si-3. *P < 0.05, **P < 0.01

The function of circTGFβ2 in MuSC differentiation

   To elucidate the function of circTGFβ2 in MuSC differentiation, we
   engineered an overexpression vector for circTGFβ2 (pLV-circTGFβ2) and
   designed three small interfering RNAs (siRNA-1, -2, -3) targeting the
   back-spliced junction of circTGFβ2. Quantitative reverse transcription
   PCR (qRT-PCR) analysis verified that the pLV-circTGFβ2 vector
   successfully facilitated the overexpression of circTGFβ2
   (Fig. [114]6I), Conversely, all three small interfering RNAs (siRNAs)
   markedly reduced circTGFβ2 expression in MuSCs (Fig. [115]6J).
   Importantly, the silencing of circTGFβ2 did not lead to any significant
   changes in TGFβ2 mRNA expression levels (Fig. [116]6K). Furthermore, we
   examined the effects of circTGFβ2 overexpression and knockdown on the
   differentiation of myoblasts. The expression levels of well-established
   myogenic differentiation markers, namely MyoD, MyoG, and MyHC, were
   quantified in goat MuSCs subsequent to transfection with pLV-circTGFβ2
   and siRNA constructs (si-1, si-2, si-3). The overexpression of
   circTGFβ2 resulted in a marked upregulation of the mRNA expression of
   MyoD, MyoG, and MyHC (Fig. [117]6L). In contrast, the knockdown of
   circTGFβ2 led to a significant downregulation of the mRNA expression of
   these markers (Fig. [118]6M). These findings suggest that circTGFβ2
   exerts a promotive effect on the differentiation of myoblasts in goats.

Discussion

   Myogenesis is a process including myoblast proliferation,
   differentiation and myotube formation and is controlled by a series of
   myogenic regulatory factors [[119]33, [120]34]. In recent years, an
   increasing number of non-coding RNAs capable of regulating myogenesis
   have been identified. As a unique form of non-coding RNA family,
   circRNAs are distinguished by their high abundance, high degree of
   stability, and tissue/spatiotemporal-specific expression, all of which
   play crucial roles in a variety of biological processes [[121]18,
   [122]35, [123]36]. Emerging evidence suggests that circRNAs play
   significant roles in the growth and development of skeletal muscle
   across multiple species, including pigs [[124]37, [125]38], cattles
   [[126]28, [127]39], and chickens [[128]40, [129]41]. Despite recent
   advancements, the expression profiles and biological functions of
   circRNAs in goat skeletal muscle stem cells remain largely unexplored.
   To elucidate the roles of circRNAs in myoblast differentiation, we
   conducted a comprehensive RNA sequencing analysis to profile circRNA
   expression in both proliferative and differentiated goat skeletal
   muscle satellite cells. We further predicted the functions of
   differentially expressed circRNAs through GO and KEGG pathway analyses,
   complemented by interaction network assessments.

   In this study, a total of 2,094 circRNAs were identified, with 84
   DEcircRNAs observed across the three stages. GO enrichment analysis
   indicated that the source genes of the DEcircRNAs were significantly
   enriched in processes related to the regulation of striated muscle
   contraction, skeletal muscle contraction, cell adhesion, muscle tissue
   development, and striated muscle hypertrophy, all of which are closely
   associated with muscle growth and development. In addition, KEGG
   pathway enrichment analysis indicated that the source genes of the
   DEcircRNAs were involved in Hippo, MAPK, TGF-beta and FoxO signaling
   pathways. Research has demonstrated that the Hippo signaling pathway
   plays a pivotal role in the regeneration and functional regulation of
   skeletal muscle [[130]42–[131]44]. Moreover, previous studies have
   established that the MAPK signaling pathway functions as the principal
   pathway during skeletal muscle differentiation [[132]45, [133]46].
   Additionally, the TGF-beta signaling pathway is acknowledged for its
   complexity and importance in the regulation of various cellular
   processes [[134]47]. It is well-characterized as a potent inhibitor of
   muscle cell differentiation due to its suppressive effect on the
   transcriptional activity of myogenic regulatory factors (MRFs)
   [[135]48]. Furthermore, the FoxO family of proteins, which are highly
   conserved transcription factors, play pivotal roles in regulating
   skeletal muscle plasticity. Specifically, FoxO1 and FoxO3 are essential
   for maintaining muscle energy homeostasis by modulating glycolytic and
   lipolytic fluxes as well as mitochondrial metabolism [[136]49].
   Collectively, these findings suggest that the corresponding circRNAs
   perform critical regulatory functions in the development of goat
   skeletal muscle. Additionally, the reliability of the circRNAs
   identified through RNA sequencing was corroborated by RT-PCR and
   qRT-PCR validation of six randomly selected DEcircRNAs.

   Previous studies have established that circRNAs can act as miRNA
   sponges by competitively binding to miRNAs, thereby influencing the
   activity of their target genes [[137]14, [138]50]. In order to
   elucidate the role of circRNAs in the differentiation of goat
   myoblasts, a circRNA-miRNA interaction network was constructed to
   predict potential interactions between circRNAs and miRNAs. Analysis of
   this network indicated that circRNAs are capable of interacting with
   multiple miRNAs. For instance, novel_circ_003559 was found to be
   associated with 18 miRNAs, whereas novel_circ_002145 was linked to 14
   miRNAs. While none of the circRNAs have been annotated in myoblasts,
   our observations indicate that several miRNAs, such as chi-miR-206,
   chi-miR-381, chi-miR-16b-5p, and chi-miR-15b-5p, have previously been
   shown to play functional roles in myogenesis [[139]30–[140]32]. These
   findings suggest that circRNAs may regulate gene expression related to
   myogenesis through interactions with miRNAs.

   Consequently, we selected circTGFβ2, a circRNA derived from the goat
   TGFβ2 gene, for validation experiments in MuSCs, due to its interaction
   with chi-miR-206, a well-established miRNA specific to skeletal muscle
   [[141]51–[142]53]. Before initiating the myoblast experiments, we
   employed RT-PCR, Sanger sequencing, and RNase R treatments to verify
   the presence of circTGFβ2. The results demonstrated that circTGFβ2 is
   abundantly expressed in skeletal muscle tissues and that its expression
   rapidly increases during the differentiation of MuSCs, suggesting a
   potential role in muscle differentiation. Furthermore, the subcellular
   localization of circRNAs is crucial in determining their molecular
   functions. Within the cytoplasm, circRNAs can function as competitive
   inhibitors of miRNAs by binding to them, thereby alleviating the
   miRNA-mediated repression of their target genes [[143]14, [144]17]. In
   the present study, circTGFβ2 was predominantly localized in the
   cytoplasm, suggesting a possible regulatory function through miRNA
   competition. We subsequently investigated the role of circTGFβ2 in the
   differentiation of goat myoblasts. The findings demonstrated that
   circTGFβ2 promotes myogenic differentiation in these myoblasts. Recent
   studies have elucidated the critical role of circRNAs in myogenesis.
   Specifically, circMEF2As [[145]54], circACTA1 [[146]55], circFgfr2
   [[147]56], and circPTPN4 [[148]57] have been confirmed to enhance
   myogenesis by positively regulating myoblast differentiation. In
   contrast, circMEF2D [[149]58], circHUWE1 [[150]59], circTAF8 [[151]60],
   and circCPE [[152]61] have been identified as inhibitors of
   differentiation. Collectively, these findings suggest that circTGFβ2
   promotes myogenic differentiation in goats, potentially via a miRNA
   sponge mechanism.

Conclusions

   In summary, this study characterized the expression profiles of
   circRNAs during the differentiation of myoblasts in goats. Enrichment
   analyses revealed that DEcircRNAs were significantly linked to a range
   of biological processes and pathways relevant to muscle development,
   such as muscle contraction, muscle tissue development, and the TGF-beta
   and MAPK signaling pathways. Moreover, the circRNA-miRNA interaction
   networks elucidated in this study provide a valuable resource of
   candidate circRNAs involved in myogenesis. Furthermore, we identified
   and validated a novel circRNA, circTGFβ2, which was observed to promote
   myogenic differentiation in goats. The findings of this study provide a
   foundation for future investigations into the mechanisms by which
   circRNA influence the proliferation and differentiation of goat
   skeletal muscle.

Materials and methods

Ethics statement

   The Animal Care and Use Committee of the College of Animal Science and
   Technology, Sichuan Agricultural University, Sichuan, China, approved
   all of the animal care, slaughter, and experimental procedures in
   accordance with the Regulations for the Administration of Affairs
   Concerning Experimental Animals (Ministry of Science and Technology,
   China) [Approval No. SAU20231023].

Sample preparation

   All experimental Chengdu Ma goats were managed and raised on a farm
   managed by the Chengdu Xilingxue Agricultural Development Co., Ltd.
   (Sichuan, China). Prior to sample collection, the goats were
   administered Jingsongling (batch number 870721, produced by Shandong
   Zibo Veterinary Medicine Factory, Shandong, China) via intramuscular
   injection in the hip region at a dosage of 2 mg/kg to achieve deep
   anesthesia. Once full anesthesia was confirmed, the goats were
   euthanized through arterial exsanguination, performed by trained
   slaughterhouse personnel (this method complies with euthanasia guide of
   the Chinese Society of Laboratory Animals, No. T/CALAS 31–2017). The
   primary skeletal muscle satellite cells (MuSCs) were isolated and
   cultured from longissimus dorsi muscle derived from a newborn goat
   (Chengdu Ma goat, female), as described previously [[153]62, [154]63].
   The MuSCs were seeded in 6-well plates (at a density of ~ 2 × 10^5
   cells/well) and cultured in Dulbecco’s Modified Eagle Medium (DMEM) in
   a 5% CO[2] incubator at 37℃, supplemented with 10% FBS (Gibco,
   Invitrogen, Carlsbad, CA, USA) and 2% penicillin/streptomycin
   (Solarbio, Beijing, China). When the confluence reached about 80–90%,
   the growth medium (GM) was replaced with differentiation medium (DM)
   containing DMEM, 2% horse serum (Gibco) and 2% penicillin/streptomycin
   in order to induce MuSC differentiation. The medium was replaced with
   fresh medium every 48 h. Proliferating MuSCs were labelled GM samples,
   while MuSCs differentiated for 1 and 5 days were labelled as DM1 and
   DM5 samples (three biological replicates at each time point). In
   addition, three female Chengdu Ma goats were sacrificed humanely on the
   third day after birth. Three skeletal muscles (longissimus dorsi
   muscle, psoas major muscle, and gastrocnemius muscle) and five internal
   organs (heart, liver, spleen, lung, and kidney) were collected. All of
   the samples were kept at -80℃ before RNA extraction.

RNA extraction, library construction, and sequencing

   Total RNAs was extracted using Trizol reagent kit (Invitrogen,
   Carlsbad, CA, USA) according to the manufacturer’s protocol. RNA
   quality was assessed on an Agilent 2100 Bioanalyzer (Agilent
   Technologies, PaloAlto, CA, USA) and checked using RNase free agarose
   gel electrophoresis. After extracted, total RNAs were treated with
   RNase R to degrade the linear RNAs, and purified using RNeasy MinElute
   Cleanup Kit (Qiagen, Venlo, The Netherlands). Next, strand-specific
   library was constructed using VAHTS Total RNA-seq (H/M/R) Library Prep
   Kit (Vazyme, Nanjing, China) for Illumina following the manufacture’s
   instructions. Briefly, ribosome RNAs were removed to retain circRNAs.
   The enriched circRNAs were fragmented into short fragments by using
   fragmentation buffer and reverse transcribed into cDNA with random
   primers. Second-strand cDNA were synthesized by DNA polymerase I, RNase
   H, dNTP (dUTP instead of dTTP) and buffer. Next, the cDNA fragments
   were purified with VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China),
   end repaired, poly(A) added, and ligated to Illumina sequencing
   adapters. Then UNG (Uracil-N-Glycosylase) was used to digest the
   second-strand cDNA. The digested products were purified with VAHTSTM
   DNA Clean Beads, PCR amplified, and sequenced using Illumina
   Novaseq60000 by Gene Denovo Biotechnology Co. (Guangzhou, China).

Quality control and circRNAs identification

   Raw reads were filtered to acquire clean reads by removing
   adapter-containing reads, poly-N-containing reads (over 10%), and
   low-quality reads (Q-value ≤ 20) using the fastp (version 0.18.0)
   [[155]64]. The clean reads were aligned to the goat reference genome
   (ARS1, GCF_001704415.1) using the TopHat2 (version 2.1.1) [[156]65].
   CircRNAs were then identified using find_circ [[157]66] and CIRI
   [[158]67]. The anchor alignments were then extended such that the
   complete read aligns and the breakpoints were flanked by GU/AG splice
   sites. A candidate circRNA was called if it was supported by at least
   two unique back spliced reads at least in one sample. The identified
   circRNAs were subjected to statistical analysis of type, chromosome
   distribution and length distribution. The expression level of each
   annotated circRNA was normalized by calculating the Reads Per Million
   mapped reads (RPM). The edgeR package (version 3.12.1) was applied to
   identify differentially expressed circRNAs (DEcircRNAs) in each
   comparison, with |log[2](Fold Change)| > 1 and P < 0.05.

GO Enrichment, KEGG Pathway, and CircRNA-miRNA network

   Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)
   pathway enrichment analyses were performed using the DAVID tool
   [[159]68, [160]69]. In the DAVID tool, the Benjamini method requests
   adjusted p-values using the linear step-up procedure of Benjamini and
   Hochberg, which was utilized to identify statistically significant GO
   terms and KEGG pathways [[161]70]. The genome of the goat (Capra
   hircus) served as the background parameter for this analysis. The
   hypergeometric test was employed to evaluate statistical significance,
   with a threshold significance level set at P < 0.05 for identifying
   significantly enriched GO terms and KEGG pathways. CircRNAs are
   proposed to regulate gene expression through competing endogenous RNAs
   acting as miRNA sponges [[162]29]. The target miRNAs, which could be
   bound to the circRNAs, were predicted using the miRanda [[163]71] and
   TargetScan [[164]72] software with a set of default parameters.
   Cytoscape [[165]73] was used for network visualization.

Validation of CircRNAs by qRT-PCR, sequencing, and RT-PCR

   The circRNA sequencing data was validated by randomly selecting six
   differentially expressed circRNAs to assess the reliability of the
   obtained results. Given the unique head-to-tail junction characteristic
   of circRNAs, their authenticity was confirmed using RT-PCR. The RT-PCR
   amplicons were analyzed via 2% agarose gel electrophoresis and
   subsequently sequenced using Sanger sequencing. To evaluate the
   sensitivity of circRNAs to RNase R, RT-PCR was conducted on RNA samples
   treated with and without RNase R. Additionally, qRT-PCR was employed to
   verify the expression levels of circRNAs, using GAPDH as the reference
   gene. Total RNA was extracted from tissues or cells using TRIzol
   Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s
   instructions. cDNA was synthesized using a RT reagent kit (Takara,
   Dalian, China) and qRT-PCR for circRNA and mRNA was performed on the
   CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) using
   SYBR Premix Ex Taq™ II (Takara, Dalian, China). Each experiment was
   conducted independently three times with three biological replicates.
   The 2^−ΔΔCt method [[166]74] was used to calculate relative expression
   levels. The primer sequences used in this study are listed in
   Supplementary Table [167]7.

Vector construction and siRNA synthesis

   To construct the novel_circ_000975 (designated as circTGFβ2)
   overexpression vector, exons 4, 5, 6, and 7 of the caprine TGFβ2 gene
   were amplified utilizing cDNA derived from MuSCs and subsequently
   inserted into a pLV-CiR vector (Geneseed, Guangzhou, China) at the
   EcoRI and BamHI restriction sites. Additionally, three small
   interfering RNAs (siRNAs) targeting circTGFβ2 were designed and
   synthesized by RiboBio (Guangzhou, China). The sequence of siRNAs are
   listed in Supplementary Table [168]7.

Cell culture and transfection

   The goat MuSCs were seeded in 6-well (~ 2 × 10^5 cells per well) plates
   containing growth medium (GM) (DMEM containing 10% FBS, Hyclone) and
   cultured at 37 °C in a 5% CO[2] atmosphere. When MuSCs achieved 80%
   confluence, GM was replaced with differentiation medium (DM) containing
   2% horse serum (Hyclone) to stimulate myoblast differentiation. The
   medium was replaced every 48 h. For transfection, MuSCs were seeded in
   6-well plates and allowed to reach 70–80% confluence. Transfection was
   performed using Lipofectamine 3000 (Invitrogen, USA) according to the
   manufacturer’s protocol, with a minimum of three replicates. After
   24 h, the transfected cells were collected for qRT-PCR assay. The final
   concentration of the pLV-circTGFβ2 vector and siRNAs used in
   transfections was 4 µg/mL and 50 nM, respectively.

RNA FISH

   For the synthesis of RNA-FISH probes and the subsequent RNA-FISH assay
   on goat MuSCs using the FISH assay kit (Gefan Biotechnology Co.,
   Shanghai, China), cells were seeded into 6-well plates and cultured
   until they reached 60 ~ 70% confluence. The cells were then fixed with
   4% paraformaldehyde. Following fixation, the cells were permeabilized
   with 0.3% Triton X-100 and incubated with a 20 mg/mL probe at 65 °C for
   48 h. The nuclei were counterstained with DAPI, and fluorescence images
   were captured using a fluorescence microscope (Olympus, Tokyo, Japan).

Luciferase activity assay

   To further elucidate the binding interaction between circTGFβ2 and
   chi-miR-206/chi-miR-211, both wild-type and mutant circTGFβ2 fragments
   were amplified and subcloned into the psiCHECK-2 vector. The circTGFβ2
   mutant (circTGFβ2-Mut) was generated by altering the binding site using
   the TaKaRa MutanBEST Kit. A list of primers utilized for plasmid
   construction is provided in Supplementary Table [169]7. In the
   luciferase reporter assays, MuSCs were cultured in 24-well plates and
   subsequently co-transfected with either circTGFβ2-WT or circTGFβ2-Mut
   constructs, along with chi-miR-206/chi-miR-211 mimics or a negative
   control mimic (NC mimics). Following a 48-hour post-transfection
   period, the activities of firefly and renilla luciferases were
   quantified using a dual luciferase reporter assay kit (Promega,
   Madison, WI, USA), in accordance with the manufacturer’s protocol.

Statistical analysis

   All data were expressed as mean ± standard error, with at least three
   biological replicates. The student’s t-test was used to analyze
   statistical significance for two-group comparisons, and one-way ANOVA
   with Tukey’s post-hoc test for multiple-group comparisons using SAS
   software version 9.2 (SAS, Cary, NC, USA). The significance level was
   set at *P < 0.05 and **P < 0.01.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [170]Supplementary Material 1^ (12.4KB, xlsx)
   [171]Supplementary Material 2^ (329.8KB, xlsx)
   [172]Supplementary Material 3^ (28.9KB, xlsx)
   [173]Supplementary Material 4^ (305.9KB, xlsx)
   [174]Supplementary Material 5^ (25.6KB, xlsx)
   [175]Supplementary Material 6^ (12.7KB, xlsx)
   [176]Supplementary Material 7^ (13.7KB, xlsx)
   [177]Supplementary Material 8^ (15.7KB, docx)
   [178]Supplementary Material 9^ (9.4MB, zip)

Acknowledgements