Identification and characterization of circular RNAs in Longissimus dorsi muscle tissue from two goat breeds using RNA-Seq

Circular RNAs (circRNAs) are a class of non-coding RNA that play crucial roles in the growth and development of skeletal muscle. However, little is known about the role of circRNAs in caprine skeletal muscle. In this study, the size of muscle fiber and the expression profiles of circRNAs were compared in Longissimus dorsi muscle of Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats with significant phenotypic differences in meat production performance, using hematoxylin and eosin staining and RNA-Seq, respectively. The size of muscle fiber in LC goats was larger than those in ZB goats (P < 0.05). A total of 10,875 circRNAs were identified and 214 of these were differentially expressed between the two caprine breeds. The parent genes of differentially expressed circRNAs were mainly enriched in connective tissue development, Rap1, cGMP-PKG, cAMP and Ras signaling pathway. In conclusion, circRNAs may play important roles in skeletal mass, meat production performance and meat quality traits in goats. The results provide an improved understanding of the functions of circRNAs in skeletal muscle development of goats.


Introduction
Skeletal muscle is the largest organ of animals and its growth and development directly determines skeletal muscle mass, meat production performance and meat quality in domestic animals. It is now found that many functional genes and noncoding RNAs are involved in the growth and development of skeletal muscle, so the identification of RNAs that regulate skeletal muscle activity can provide a chance to improve meat production. However, compared with our knowledge of skeletal muscle mRNAs (Hernández-Hernández et al. 2017;Taylor and Hughes 2017;Zammit 2017) and miRNAs (Tao et al. 2018;Cheng et al. 2020;Cai et al. 2021), the reports on the roles of circular RNA (circRNAs) in the skeletal muscle are very limited.
CircRNAs are a novel class of non-coding RNA and they are result from the covalently linage of 5′ and 3′ ends of linear RNA (Lasda and Parker 2014). Due to the closed loop structure, circRNAs are not susceptible to be affected by RNA exonucleases. It is, therefore, known that circR-NAs are more stable and evolutionarily conserved than the linear mRNAs (Jeck and Sharpless 2014). In recent years, the function of circRNAs is gradually uncovered. CircR-NAs mainly function as microRNA (miRNA) sponges. They Communicated by Joan Cerda.
1 3 thereby relieve the inhibition of target mRNAs by miRNAs, eventually resulting in an increase in expression level of target mRNAs (Hansen et al. 2013;Ouyang et al. 2018;Chen et al. 2019). In addition, exon-intron circRNAs predominantly localized in the nucleus can regulate the transcription of their parent genes in a cis-regulatory way (Li et al. 2015). Finally, other functions of circRNAs have also been reported, including regulating alternative splicing (Ashwal-Fluss et al. 2014), interacting with RNA-binding proteins (Hentze and Preiss 2014), and even being translated into protein .
To date, the studies of circRNAs in skeletal muscle tissues have mainly been focused on pigs (Hong et al. 2019;Shen et al. 2019a;Wang et al. 2019a;Li et al. 2020), chickens (Ouyang et al. 2018;Chen et al. 2019;Shen et al. 2019b) and cattle (Wei et al. 2017;Liu et al. 2020;Yan et al. 2020). These studies reported the expression profiles and characterization of circRNAs in skeletal muscle either from different breeds, or from different development stages, and also analyzed the function of circRNAs more by investigating the functions of their parent genes using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. For example, 13,377 circRNAs were identified in leg muscles of chickens between 11 embryo age, 16 embryo age and 1 day post hatch periods. Subsequently, one (cir-cRBFOX2s) of the differentially expressed circRNAs identified was found to promote the proliferation of chicken myoblast by binding with miR-206 (Ouyang et al. 2018). In addition, a total of 655 circRNAs were found to be differentially expressed in muscle tissues between Shandong black cattle and Luxi cattle with significant differences in the diameter, length and weight of muscle fibers, and their parent genes were mainly involved in muscle fiber development process, and MAPK and mTOR signaling pathways .
However, there are only two reports about the activity of circRNAs in muscle tissues in goats. Ling et al. (2020) identified 9090 circRNAs in skeletal muscle of Anhui white goats across seven development stages, and the parent genes of the differentially expressed circRNAs were involved in the regulation of myoblast differentiation, skeletal muscle maturation and hypertrophy. In addition, a circRNA CDR1as has been found to activate the differentiation of caprine skeletal muscle satellite cells by relieving the inhibition of IGF1R targeted by miR-7 .
Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats are famous local breeds in China and of economic importance for goat farmers. Compared with ZB goats, LC goats have higher carcass weight, content of intramuscular fat, and proportion of intramuscular collagen fibers, but poorer meat tenderness (P < 0.05) Wang et al. 2021). In this study, the circRNAs expression profiles of Longissimus dorsi muscle tissue were compared between LC and ZB goats using RNA-Seq. The GO enrichment and KEGG pathway were also analyzed for the parent genes of differentially expressed circRNAs between the two caprine breeds. The results will provide an improved understanding of the function of circRNAs in skeletal muscle growth and development processes in goats.

Ethics statement
All experimental procedures were approved by Faculty of Animal Science and Technology, Gansu Agricultural University, Lanzhou, China. The investigation for experimental animals was also coincide with the rules for animal care and use published by the Ministry of Science and Technology of the People's Republic of China (Approval number 2006-398).

Experimental animals and sampling
The same Longissimus dorsi muscle samples as described by Shen et al. (2021) were used in this study. Briefly, five 9-month-old LC rams and five 9-month-old ZB rams were selected from the Yongfeng Goat Breeding Company (Huan County, China) and then slaughtered. These goats were raised under the same feeding and management conditions. The carcass weight, muscle components and meat quality of these goats investigated in the study are presented in Table 1 Wang et al. 2021).
The Longissimus dorsi muscle samples from the area between 12 and 13th ribs on the left carcass were collected and then used for RNA isolation and hematoxylin and eosin staining. The samples for RNA isolation were frozen in liquid nitrogen, whereas the samples for hematoxylin and eosin staining was fixed with 4% paraformaldehyde.

Hematoxylin and eosin staining of Longissimus dorsi muscle
The Longissimus dorsi muscle samples fixed with 4% paraformaldehyde were treated using graded ethanol (75, 85, 95, and 100%) to remove moisture. Subsequently, the dehydrated specimens were embedded in paraffin and then cut into about 5 µm of thickness using Rotary cutting machine (Leica, Wetzlar, Germany). The paraffin sections were used for hematoxylin and eosin staining as suggested by Cao et al. (2014).
Micrographs (400×) of hematoxylin and eosin staining from three different fields of view for each sample were taken by BA200 Digital microscope (MOTIC, Xiamen, China). The diameter and cross-sectional area of muscle 1 3 fibers were then measured using Motic Images Advanced v3.2. The difference in these measurements between LC and ZB goats was analyzed using Student's t test in SPSS v24.0.

RNA samples preparation and sequencing
Total RNA was isolated from ten caprine Longissimus dorsi muscle samples using Trizol reagent (Invitrogen, Carlsbad, CA, United States). The integrity and concentration of these RNA samples were measured using Agilent 2100 Bioanalyzer (Agilent, CA, United States) and Nanodrop 2000 (Thermo Scientific, MA, United States), respectively. High-quality RNA samples were screened with a parameter of RNA integrity number > 7 being fitted as the threshold. The ribosomal RNA (rRNA) was removed from these high quality RNA samples using a Ribo-Zero Gold rRNA Removal Kit (Illumina, CA, United States). The remaining RNA was fragmented into 200-300 bp in length and then used for constructing cDNA libraries using a NEBNext Ultra RNA Library Prep Kit (New England Biolabs, MA, United States). The cDNA libraries were paired-end sequenced with a sequencing depth of 12 G clean data on a HiSeq™ 4000 sequencer (Illumina, CA, United States).

Analysis of RNA-Seq data
The clean reads were obtained by removing the reads with quality scores < Q20, reads containing sequencing adapters and reads with > 10% unknown nucleotides from raw reads produced from sequencer, using fastp v0.18.0 . These clean reads were then mapped to Caprine Genome Assembly ARS1 (ftp:// ftp. ncbi. nlm. nih. gov/ genom es/ all/ GCF/ 001/ 704/ 415/ GCF_ 00170 4415.1_ ARS1) using HISAT2 v2.1.0 (Kim et al. 2015). For the reads that were unmapped against Caprine Genome Assembly ARS1, 20-mers from both ends were defined as anchor reads. The anchor reads were mapped to the reference genome again using bowtie2 v2.2.8 (Langmead and Salzberg, 2012) and circRNAs were identified using software Find_circ (Memczak et al. 2013). The identified circRNAs were characterized by counting their chromosomal distribution, type and length. The type of differentially expressed circRNAs was also counted. The expression level of each annotated circRNA was normalized by calculating the Reads Per Million mapped reads (RPM). CircRNAs with RPM > 0 were considered to be expressed. DESeq v2.0 (Love et al. 2014) was used to screen differentially expressed circRNAs in Longissimus dorsi muscle tissues between LC and ZB goats, with fold change > 2.0 and P value < 0.05. It is notable that the correction for multiple testing was not performed in this study as the calculation method of P value in DESeq2 is very strict and sufficient as the screening criterion for differentially expressed circRNAs. In addition, if the corrected P value was used, the number of differentially expressed circRNAs will sharply decrease, which was not conducive to analyze the function of the parent genes of differentially expressed circRNAs by GO and KEGG analysis.

Validation of the authenticity of circRNAs using reverse transcriptase-PCR and DNA sequencing
Based on the characteristic that circRNA has a unique head-to-tail junction, its authenticity was validated using reverse transcriptase-PCR (RT-PCR) and DNA sequencing. Briefly, 20 differentially expressed circRNAs between the two caprine breeds were randomly selected. The ten RNA samples that were used for RNA-Seq were also used to produce cDNA using an RT-PCR kit (Takara, Dalian, China). The cDNA products were then amplified using divergent primers (Table 2) in a 20-µL reaction, containing 0.8-µL of the cDNA, 0.5 U of Taq DNA polymerase (Takara, Dalian, China), 150 µM of each dNTP (Takara, Dalian, China), 2.5 mM Mg 2+ , 0.25 µM of each primer and 1× PCR buffer supplied with DNA polymerase enzyme. RT-PCR amplicons were checked in 1% agarose gels electrophoresis and then sequenced using Sanger sequencing. To confirm the authenticity of these circRNAs, the sequences from Sanger sequencing were aligned to the goat reference genome and RNA-Seq data to validate the location of the junction sites of these circRNAs.

Validation of the reliability of RNA-Seq results using RT-quantitative PCR
The same 20 differentially expressed circRNAs as those used for RT-PCR, were selected to confirm the reliability of RNA-Seq results using RT-quantitative PCR (RT-qPCR). The ten RNA samples that were used for RNA-Seq were also used to produce cDNA using an RT-PCR kit (Takara, Dalian, China). The RT-qPCR was performed in triplicate using a 2× ChamQ SYBR qPCR Master Mix (Vazyme, Nanjing, China) on an Applied Biosystems QuantStudio ® 6 Flex (Thermo Lifetech, MA, USA). Caprine GAPDH and β-actin were used as internal control for normalizing the expression levels of these circRNAs as suggested by Ling et al. (2020) and Wang et al. (2019a). A 2 −ΔΔCt method was finally used to calculate the relative expression levels.

Functional enrichment analysis of the parent genes of differentially expressed circRNAs and construction of circRNA-miRNA-mRNA network
The parent genes of circRNAs were annotated by aligning their sequence to the RefSeq non-redundant proteins (nr) database using BLASTX (E value < 10 -5 ). GO (Ashburner et al. 2000) analysis was used to analyze the main function of the parent genes of the differentially expressed circRNAs using DAVID (Dennis et al. 2003). The KOBAS 3.0 (Bu et al. 2021) was used to identify the pathway that included the parent genes of the differentially expressed circRNAs using KEGG database. The GO terms and KEGG pathways with P value < 0.05 were considered to be significantly enriched based on hypergeometric test. Meanwhile, five differentially expressed circRNAs were used to predict their miRNA binding sites using Mireap v0.2 , Miranda v3.3a (Turner 1985) and TargetScan v7.0 (Lewis et al. 2005). The three software were also used to predict the target genes of the miRNAs sponged by the five circRNAs and predicted results from the three kinds of software were compared. The Pearson's coefficients in expression levels between the five differentially expressed circRNAs and their corresponding mRNAs was calculated using SPSS v24.0. An interaction network among these circRNAs and their sponged miRNAs as well as corresponding mRNAs was constructed using StarBase v3.0 (Li et al. 2014) and then drawn using Cytoscape v3.5.1 (Shannon et al. 2003).

Comparison of the size of muscle fiber between LC and ZB goats
Muscle fibers are the main structural components of skeletal muscle tissues, accounting for about 70% of the total  TGC TTT CCG AAG TG  GCT TCA TTA AAG CCC ACC AA  circ_003757  CAG ATC CGA ATG CAA CTA  GCA ACA TCA ATG CTA CGC  circ_008447  TGC AGA CGA CGA TAA CTT GG  TCT TTA TCA CAC AGA ACT TGTTC  circ_002956  CAT GGT TCT TCT GGT TTT GGA  TAG GAA AGG AGT TCC CAG CA  circ_007137  CCA GGA AGT GAA GGA AGC AG  TGC CCA CTG TCC AAA CTT CT  circ_001875  CAC CTG CGA TCT CAA GCT CT  GCC CGA ATA ATG TCG TTG AT  circ_006718  AAG GGC CTA AGT GAC CGT TT  TGA CCC TGA ACC TGC TCA AT  circ_008092  TCA AGA CCC TAG AAG ATT TGCA  CAT GAA AAC GGA TGG TGG CA  circ_008770 CCC weight. The diameter and cross-sectional area are important morphology characteristics of muscle fibers and they play key roles in determining meat production performance of animals (Tuma et al. 1962;Larzul et al. 1997;Rehfeldt et al. 2000;Choi and Oh 2016). The comparison results of the diameter and cross-sectional area of Longissimus dorsi muscle fiber between LC and ZB goats are shown in Fig. 1. The diameter and cross-sectional area of muscle fibers from LC goats were 38.52 ± 2.20 μm and 1902.91 ± 156.92 μm 2 , respectively, which was larger than ZB goats with a diameter of 29.09 ± 3.81 μm and a cross-sectional area of 1117.72 ± 210.10 μm 2 (Fig. 1).

Longissimus dorsi muscle tissues
The raw reads obtained in the study from ten Longissimus dorsi muscle samples have been submitted to GenBank database (https:// www. ncbi. nlm. nih. gov/ sra) with accession numbers SRR13008213-SRR13008222. The RNA-Seq reads and mapped results to the goat reference genome have also been described in our previous study

Authentication of caprine circRNAs in the Longissimus dorsi muscle tissue
A total of 20 circRNAs were selected to validate the presence of their specific head-to-tail junctions using RT-PCR and Sanger sequencing. Agarose gel electrophoresis results showed that all 20 circRNAs were expressed and they produced a band with an expected size ( Fig. 2A). Sanger sequencing further affirmed the presence of head-to-tail splice junctions and the size for these caprine circRNAs, t-test and shown as mean ± SD (n = 5). *P < 0.05; **P < 0.01 which was consistent with those provided by RNA-seq (Fig. 2B).

Characterization of circRNAs identified in the caprine Longissimus dorsi muscle tissue
Of the six types of circRNAs, annot_exons were the most common sequences with an average proportion of 78.03%, followed by exon_intron, one_exon and antisense sequence types with an average proportion of 7.69%, 5.25%, and 3.97%, respectively. Intronic and intergenic sequences were the least common types with an average proportion of 3.14 and 1.91%, respectively (Fig. 3A). The circRNAs identified from Longissimus dorsi muscles were widely distributed across all the caprine chromosomes, with the exception of Y chromosome. The most circRNAs were distributed on chromosomes 2 and 1, while the least circRNAs were distributed on chromosome 27 (Fig. 3B). The length of most circRNAs was less than 1 kb (Fig. 3C).

Identification and validation of differentially expressed circRNAs
Of all the 10,875 circRNAs identified in caprine Longissimus dorsi muscle tissues, 214 circRNAs were found to be differentially expressed between LC and ZB goats, including 85 up-regulated circRNAs and 129 down-regulated circRNAs in Longissimus dorsi muscle of LC goats when compared to ZB goats (Supplementary File 1). The most up-regulated circRNA in LC goats was circ_008092 derived from signal transducer and activator of transcription 1 (STAT1) with a 19.4-fold increase in expression, while circ_003628 derived from myosin-4 (MYH4) was the most down-regulated cir-cRNA in LC goats with a 55.6-fold decrease in expression. Most notably, of 214 differentially expressed circRNAs identified, 25 circRNAs were only expressed in LC goats, while 39 circRNAs were only expressed in ZB goats. The type distribution of differentially expressed circRNAs was similar to the overall distribution of all circRNAs identified. Briefly, annot_exons were the most common sequences with an average proportion of 81.31%, followed by exon_intron, one_exon and intergenic sequences with an average proportion of 9.35%, 5.14%, and 2.34%, respectively. Antisense and intronic sequences were the least common types with an average proportion of 0.93% and 0.93%, respectively.
To validate the repeatability of RNA-seq data, 20 differentially expressed circRNAs that were used for the validation of the presence of circRNAs, were also subjected to RT-qPCR analysis. Because circ_006172 and circ_009217 were only expressed in LC and ZB goats, respectively, the log 2 fold-change for LC goats relative to ZB goats was infinity for these two circRNAs. Their relative expression levels are therefore not presented in Fig. 4. As shown in Fig. 4, the expression levels of circ_009262, circ_003757, circ_008447, circ_002956, circ_007137, circ_001875, circ_006718, circ_008092, circ_008770, circ_009387, and circ_002473 in LC goats were higher than those in ZB goats. In contrast, compared to ZB goats, the expression levels of circ_005348, circ_004394, circ_003976, circ_008047, circ_010303, circ_007151, and circ_008117 were lower in LC goats. The results suggest that RT-qPCR results for these circRNAs in the Longissimus dorsi muscle were in accordance with those obtained from RNA-Seq (Fig. 4), and that our RNA-Seq results are repeatable and reliable.

Functional analysis of circRNAs as miRNA sponges
For the 214 differentially expressed circRNAs identified, a total of 431 miRNA binding sites were predicted. For clearly presenting the interaction effect of circRNAs and miRNAs, 5 circRNAs with the type of annot_exons were further selected, including four up-regulated circRNAs (circ_002300, circ_006172, circ_008092, and circ_001875) and one down-regulated circRNA (circ_010839) in LC goats compared to ZB goats, given that annot_exons are mainly localized in the cytoplasm and they may, therefore, play biological function by acting miRNA sponges. There were 55 miRNA binding sites in total for these 5 circRNAs. Of these miRNAs, a total of eight sponged miRNAs were selected, including miR-34b-5p, miR-27b-3p, miR-34c-5p, miR-424-5p, miR-15b-5p, miR-148a-3p, miR-96 and miR-330-5p, based on the findings that miR-330-5p were related with adipose development (Shi et al. 2018), while other miRNAs have been reported to play important roles in the growth and development of skeletal muscle Ling et al. 2018;Hou et al. 2017;Connolly et al. 2018;Li et al. 2019b;Yin et al. 2020;Nguyen et al. 2020).In addition, the target genes of these miRNAs had positive correlations in expression with their corresponding circRNAs (r > 0.6, P < 0.05) ( Table 3 and Supplementary File 2). A total of 11 circRNA-miRNA-mRNA pairs were screened and then used to construct a circRNA-miRNA-mRNA network (Fig. 6).

Discussion
In this study, the diameter and cross-sectional area of Longissimus dorsi muscle fiber from LC goats were larger than those from ZB goats. It has been reported that the Fig. 5 Function annotation of the parent genes from which the differentially expressed circRNAs were derived between Liaoning cashmere (LC) goats and Ziwuling black (ZB) goats. A GO classification of the parent genes of the differentially expressed circRNAs between the two breeds. The X axis above represents the value of − Log10 P value, while X axis below indicates the number of the parent genes of differentially expressed circRNAs. B KEGG enrichment analysis for the parent genes of the differentially expressed circRNAs between the two breeds. The most significant GO terms and KEGG pathways with the lowest P values are shown 1 3 diameter or cross-sectional area reflects the size of muscle fiber, which directly determines skeletal muscle mass during postnatal period of animals (Glass 2005;Schiaffino et al. 2013). This may partly explain why in our findings LC goats had higher carcass weight than ZB goats. Choi and Oh (2016) also found that pigs with greater cross-sectional area of muscle fibers had higher carcass weight (P < 0.001). Besides carcass weight, cross-sectional area of muscle fiber has been reported to positively correlated with intramuscular fat content (r = 0.68) of pork (Larzul et al. 1997), and muscle fiber diameter was also positively correlated with shear force value (r = 0.63) and loin eye area (r = 0.56) of beef (Tuma et al. 1962). These studies further supported our observation that meat from LC goats with higher diameter and crosssectional area of Longissimus dorsi muscle fiber had higher muscle shear force value, intramuscular fat content and loin eye area than meat from ZB goats.
The number of circRNAs (10,875) identified in the study was higher than what was investigated by Ling et al. (2020), who described 9090 circRNAs in caprine Longissimus dorsi muscle tissues of Anhui white goats. In addition, 14,640 and 6988 circRNAs were found in muscle tissues from cattle ) and pigs , respectively. This likely revealed species-specific expression pattern of circRNAs. Our observation that most of circRNAs identified were the type of annot_exons, was in accordance with the  6 Interaction network of circRNA-miRNA-mRNA. The red triangles and inverted triangles represent up-regulated and down-regulated circRNAs in LC goats compared to ZB goats, respectively. The green squares represent predicted miRNAs, while the blue circles represent the target genes of the miRNAs 1 3 findings in muscle tissues in Anhui white goats (Ling et al. 2020), pigs (Liang et al. 2017), cattle (Wei et al. 2017), and chicken (Ouyang et al. 2018). Other types of circRNAs that have been identified in Longissimus dorsi muscle of Anhui white goats (Ling et al. 2020), chicken (Ouyang et al. 2018) and pigs , were also found in the study. Most of circRNAs identified in this study were less than 1 kb, which was also consistent with the length distribution of cir-cRNAs reported in skeletal muscle of cattle (Wei et al. 2017) and pigs (Cao et al. 2022). It is noteworthy that of multiple circRNAs produced by a single gene, there were only 1-2 circRNAs with a high expression level. This phenomenon was also observed in bovine muscle circRNAs (Wei et al. 2017). The uneven distribution of circRNAs in caprine chromosome observed in the study is perhaps unsurprising as caprine chromosomes 1 and 2 are the largest in size, while chromosome 27 is relatively small in the goat genome. Ling et al. (2020) found similar chromosome distribution of cir-cRNAs in caprine skeletal muscle tissues, and studies in cattle, chicken and pigs also confirmed that the numbers of circRNAs found was proportional to chromosome size (Wei et al. 2017;Shen et al. 2019b;Cao et al. 2022). It was noteworthy that the most highly expressed cir-cRNA in both LC and ZB goats was circ_001086 derived from LMO7. LMO7 is essential to skeletal muscle development as it maintains proper myoblast differentiation (Dedeic et al. 2011). Knockdown of LMO7 inhibited myogenesis by preventing myotube formation and decreasing the number of myoblast in chicken (Possidonio et al. 2016). In addition, circLMO7 produced from LMO7 promoted the proliferation of myoblast, but inhibited the differentiation and apoptosis of myoblasts by sponging miR-378a-3p in cattle (Wei et al. 2017). These suggest that circ_001086 play key roles for skeletal muscle development in both LC and ZB goats and it is worthy of further investigation.
Compared to ZB goats, the most up-regulated and downregulated circRNAs in LC goats were circ_008092 and circ_003628, respectively, which originated from STAT1 and MYH4, respectively. STAT1 has been reported to play a positive role in myoblast proliferation and hypertrophy of skeletal muscles (Sun et al. 2007;Begue et al. 2013). In addition, our circRNA-miRNA interaction network showed that circ_008092 would act as a miR-330-5p sponge (Fig. 6). The miR-330-5p negatively regulated ovine preadipocyte differentiation by targeting BCAT2 (Tao et al. 2018). It was therefore inferred that the up-regulated expression of circ_008092 in Longissimus dorsi muscle of LC goats may be involved in the observed higher carcass weight and intramuscular fat content. The circ_003628 was predicted to bind miR-208b and miR-1, which have been reported to be positively related to skeletal muscle cell proliferation and differentiation by targeting CDKN1A and BDNF, respectively (Varendi et al. 2014;Wang et al. 2019c). These findings suggest that the lower expression of circ_003628 in LC goats may contribute to its higher muscle mass by less sequestering with these important miRNAs.
Connective tissue development was the most significant GO term enriched for the parent genes of five differentially expressed circRNAs. The parent genes have been found to be associated with the synthesis of collagen fibers, which mainly form intramuscular connective tissue. For example, LOC102187872 encodes collagen alpha-1 (V) chain, and the protein SOX6 activated the expression of the gene encoding collagen type 2 by combing with a 48 bp enhancer (Lefebvre et al. 1998). Taken together, the differentially expression of these circRNAs may partly explain the difference in intramuscular collagen fiber content between LC and ZB goats.
Interestingly, in transcriptome profile of Longissimus dorsi muscle tissues between LC and ZB goats, LOC102187872 was found to be an up-regulated gene in LC goats . These findings suggest that circ_002339 had the opposite expression tendency with its parent gene LOC102187872. Genome wide analysis results showed that the expression level of 9.0% circRNAs identified in embryonic muscle tissue of pigs were negatively correlated with their parent genes (Hong et al. 2019). Unlike circ_002339, the expression tendency of circ_005286 and circ_007643 between LC and ZB goats was consistent with their parent genes LOC102191280 and LOC102190983, respectively, namely, these circRNAs and their parent genes were all down-regulated in Longissimus dorsi muscle from LC goats compared to ZB goats. Taken together these suggest that there is a complex relationship in expression between circRNA and its parent gene.
The skeletal muscle growth during postnatal period mainly depends on the hypertrophy of muscle fiber. The mTOR pathway was one of the most important factors promoting skeletal muscle hypertrophy (Glass 2005). It has been reported that the effect of mTOR pathway on muscle fiber depends on various GTPases (Kimball and Jefferson 2006;Bar et al. 2018). In the study, some differentially expressed circRNAs were enriched in GO terms related to GTPases, the circ_003976 and circ_007919 were one of them. FNIP1 is the parent gene of circ_003976, and involved in the accumulation of crucial muscle proteins, such as myosin heavy chain and troponins (Reyes et al. 2015). CYFIP1 producing circ_007919 promoted the remodeling of actin, which is one of the most important muscle proteins (De Rubeis et al. 2013). It was, therefore, inferred that differential expression of these circRNAs in Longissimus dorsi muscle tissues between LC and ZB goats may be involved in significant phenotypic differences in muscle mass and carcass weight originated from muscle fiber hypertrophy.
The cAMP signaling pathway, Ras signaling pathway and cGMP-PKG signaling pathway play important roles in growth and development of skeletal muscle and adipose tissue. The cAMP signaling and cGMP-PKG signaling pathways are crucial for skeletal muscle development (Berdeaux and Stewart 2012;Kuo and Ehrlich 2015) and also associated with adipocyte differentiation and lipolysis (Madsen and Kristiansen 2010;Nishikimi et al. 2009). Ras signaling was involved in inhibition of myoblast differentiation and skeletal myogenesis (Olson et al. 1987;Mitin et al. 2001), and regulation of adipocyte differentiation during brown adipogenesis (Murholm et al. 2010). As might be expected, the parent genes of nine differentially expressed circRNAs enriched in the pathways were also related with muscle hypertrophy and atrophy. For example, the proteins MAPK1 and AKT3 promoted the hypertrophy of postnatal skeletal muscle (Fadia and Adams 2004), and also participated in terminal differentiation and proliferation of myoblast (Li and Johnson 2006;Wei et al. 2013). The protein MET has been found to promote muscle hypertrophy by preventing apoptosis of myogenic progenitors (Ronzoni et al. 2017;Cassano et al. 2008). The significant role of MEF2A has well been established in controlling embryonic myogenesis, adult skeletal muscle growth, hypertrophy and regeneration (Schiaffino et al. 2018). The knockout of NFKB1 inhibited the unloading-induced muscle atrophy by increasing cross-sectional areas of muscle fiber (Hunter and Kandarian 2004). Meanwhile, the regulation effects of parent genes MAPK1, AKT3 and MEF2A on adipogenesis have also been described (Zhao et al. 2011;Ding et al. 2017;Wu et al. 2018). These indicate that the parent genes of these differentially expressed circRNAs detected contributed to the differences in carcass weight and content of intramuscular fat between LC and ZB goats.
It was notable that most differentially expressed circR-NAs screened in the study were the type of annot_exons (81.31%), suggesting that they may mainly be localized in the cytoplasm and can function as miRNA sponges to positively regulate the expression levels of the target genes. For the sponge mechanism of circRNA-miRNA-mRNA, there is a positive correlation in expression levels between cir-cRNA and corresponding mRNA (Quan et al. 2021). The positive correlation relationships were also observed in the 11 circRNA-miRNA-mRNA pairs constructed in the study, indicating possible sponge mechanism existed. In this context, the roles of circRNAs in various cell activities can be reflected by the functions of their binding miRNAs. In the study, some predicted binding miRNAs have previously been reported to be associated with skeletal muscle development and intramuscular fat deposition. For example, the miR-424-5p and miR-15b-5p would target with circ_001875 ( Fig. 6), which was up-regulated in Longissimus dorsi muscle of LC goats with higher carcass weight. Previous studies found that miR-424-5p inhibited human skeletal muscle mass by reducing protein synthesis (Connolly et al. 2018), and miR-15b-5p suppress myoblast proliferation and differentiation by targeting IRS1 . It was therefore inferred that the higher expression level of circ_001875 may contribute to higher muscle mass of LC goats by sequestering the negative effect of miR-424-5p and miR-15b-5p on skeletal development. Besides these miRNAs described above, miR-27b-3p that would be sponged by circ_002300 and circ_006172, also played important roles in proliferation and differentiation of myoblast by targeting PAX3 in goats (Ling et al. 2018). These results indicate that these differentially expressed circRNAs identified in the study may play key miRNAs sponge roles in regulating the differences in meat production performance between LC and ZB goats.
In this study, a total of 214 circRNAs were differentially expressed between the two goat breeds with significant phenotype difference in meat production performance. The parent genes and sponged miRNAs of differentially expressed circRNAs were associated with muscle development and intramuscular fat deposition. This study provides an improved understanding of the roles of circRNAs in skeletal muscle development of goats. Moreover, this study lays the foundation for further research into the function of individual circRNAs in the development of muscle and adipose tissues.