Front. Genet.Frontiers in GeneticsFront. Genet.1664-8021Frontiers Media S.A.10.3389/fgene.2020.587559GeneticsOriginal ResearchIdentification of Candidate Circular RNAs Underlying Intramuscular Fat Content in the DonkeyLiBojiang1FengChunyu1ZhuShiyu1ZhangJunpeng1IrwinDavid M.2ZhangXiaoying3WangZhe1ZhangShuyi1*1College of Animal Science and Veterinary Medicine, Shenyang Agricultural University, Shenyang, China2Department of Laboratory Medicine and Pathobiology, University of Toronto, Toronto, ON, Canada3Liaoning Province Engineering Center of Modern Agricultural Production Base, Shenyang, China
Edited by: Xiao Li, Northwest A&F University, China
Reviewed by: Zhuanjian Li, Henan Agricultural University, China; Renata Veroneze, Universidade Federal de Viçosa, Brazil
*Correspondence: Shuyi Zhang, szhang@syau.edu.cn
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Intramuscular fat (IMF) content is a crucial indicator of meat quality. Circular RNAs (circRNAs) are a large class of endogenous RNAs that are involved in many physiological processes. However, the expression and function of circRNA in IMF in the donkey remains unresolved. Here we performed an expression profiling of circRNAs in the donkey longissimus dorsi muscle and identified 12,727 candidate circRNAs. Among these, 70% were derived from the exons of protein genes. Furthermore, a total of 127 differentially expressed (DE) circRNAs were identified in high (H) and low (L) IMF content groups, including 63 upregulated and 64 downregulated circRNAs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the host genes of the DE circRNAs showed that the host genes were enriched in lipid metabolism related GO terms (e.g., fatty acid beta-oxidation using acyl-CoA dehydrogenase and MLL3/4 complex), and signaling pathways (e.g., TGF-beta and lysine degradation signaling pathway). Further analyses indicated that 127 DE circRNAs were predicted to potentially interact with miRNAs, leading to the construction of circRNA-miRNA regulatory network. Multiple circRNAs can potentially function as sponges of miRNAs that regulate the differentiation of adipocytes. Our results provide valuable expression profile information for circRNA in the donkey and new insight into the regulatory mechanisms of circRNAs in the regulation of IMF content.
circular RNAexpression profileIMFregulatory networkdonkeyDepartment of Education of Liaoning Province10.13039/501100007620Introduction
The donkey (Equus asinus) is an important livestock animal in many countries including China and Italy (Polidori et al., 2015; Zhang X. et al., 2019), and plays a crucial role in human agricultural society by providing various products (meat, milk, and leather), acting as a draft force and transportation (Xia et al., 2019). Donkey meat was previously obtained from animals that were slaughtered at the end of their working lives, leading to poor quality meat with bad sensorial and nutritional characteristics (Lorenzo et al., 2014). In recent years, with the increasing mechanization of the world, donkey meat production from young males has improved the quality of the meat and attracted more customers for its consumption (Lorenzo et al., 2014). Therefore, donkey has increasingly become an important meat producing livestock. Donkey meat has recently been recognized as a nutritious food for human consumption, as it contains high-quality protein, vitamins, and minerals (Polidori et al., 2015). For example, the protein content of donkey meat is 22.8%, and the potassium and phosphorus content are 343 and 212 mg per100 g, respectively (Polidori et al., 2008).
IMF, corresponds to the fat within muscles, and its amount is controlled by the number and size of intramuscular adipocytes (Li B. et al., 2018). A previous study indicated that the IMF content plays a key role in various quality traits of meat in many species (Hocquette et al., 2010). Increased levels of IMF content can positively influence sensory quality traits such as tenderness, juiciness, taste, and flavor (Bahelka et al., 2009). IMF content has a relatively high heritability in cattle (h2 = 0.51) (Nogi et al., 2011) and pig (h2 = 0.4–0.7) (Pena et al., 2016), which indicates that these animals can be selected for higher IMF content and bred to improve this trait in the next generation. However, there are currently no reports on the heritability of IMF in donkey. IMF content is complex quantitative trait, which is affected by multiple genetic components, environmental conditions, cellular signals and hormones (Cho et al., 2019). Therefore, identifying candidate genes and molecular markers could be useful in selection programs to improve IMF content in the donkey.
To date, groups that have studied donkey meat and carcass performance have focused on these traits (Polidori et al., 2008; De Palo et al., 2016), rather than the mechanisms that control them. Recently, many previous studies have shown that coding and non-coding RNAs (e.g., miRNA and lncRNA) regulate IMF formation in the pig, cattle and chicken (Li B. et al., 2018; Zhang et al., 2018, 2020; Zou et al., 2018). However, these mechanisms of coding or non-coding RNAs in the regulation of IMF deposition in the donkey is unknown. Circular RNAs (circRNAs) are a class of endogenous non-coding RNA that are processed from precursor mRNA (pre-mRNA) by back-splicing with a covalent linkage between the 3′ and 5′ ends (Dong et al., 2017). They are extensively distributed in mammalian cells and lack the typical 5′ caps and 3′ poly(A) tails of coding mRNA (Devaux et al., 2017). Recent studies have provided evidence that the expression levels of circRNAs are highly conserved among species (Legnini et al., 2017). Due to their non-linear structure, circRNAs have higher stability than linear RNAs and are involved in many different physiological processes (Rybak-Wolf et al., 2015; Sun et al., 2019). It has been shown that circRNAs play important roles in biological and physiological processes in livestock (Bahelka et al., 2009; Li H. et al., 2018; Cao et al., 2019). For example, 828 circRNAs were found to be significantly differently expressed between embryonic and adult bovine muscle tissues and that circLMO7 regulates myoblasts differentiation and survival through the sponging of miR-378a-3p (Wei et al., 2017). A recent study by Liu et al. (2018) demonstrated that a large number of circRNAs are significantly differently expressed during adipogenesis in subcutaneous fat in pigs (Liu et al., 2018). However, it is unknown whether circRNAs have any biological role during IMF formation in donkeys. Therefore, here we systematically investigated the expression profile of circRNAs in donkey longissimus dorsi muscles.
First, we performed expression profiling of circRNAs in the donkey longissimus dorsi muscles of donkey using RNA-seq. We then identified differentially expressed circRNAs between high (H) and low (L) IMF content. GO and KEGG analyses were conducted of the host genes of the differentially expressed circRNAs between H and L IMF content. Finally, we predicted that specific miRNAs were adsorbed by the differentially expressed circRNAs, allowing us to construct the circRNA-miRNA regulatory network in IMF.
Materials and MethodsEthics Statement
All animal procedures described in this study were conducted according to the animal husbandry guidelines of Shenyang Agricultural University. The studies in these animals were reviewed and approved by the Ethics Committee and Experimental Animal Committee of Shenyang Agriculture University.
Animals, Samples Collection, and Phenotypes Measurement
The animals used in this study were derived from 30 Liaoxi donkeys, and provided by Fuxin City Lv Xian Yuan Meat Food Co., Ltd. (Fuxin, China), which had been raised under standard conditions with ad libitum access to a mixture of cereal straw and grass hay, maize, bran, peas, minerals, vitamins, and water. At an age of about 15 months, all of the donkeys were slaughtered in the same abattoir (Fuxin, China). Longissimus dorsi muscle was collected and used for the measurement of IMF content and RNA extraction. Tissue samples for RNA extraction were immediately snap-frozen in liquid nitrogen, and then stored at −80 °C until use. IMF content of the longissimus dorsi muscle was measured using the Soxhlet extraction method as described previously (Li B. et al., 2018). Of the 30 samples tested, the three with the H IMF content and the three with the L IMF content were chosen for RNA extraction. The IMF content of the H and L groups were 6.34% (SEM = 0.47%) and 3.04% (SEM = 0.12%), respectively, which was significantly different (P ≤ 0.05), however, the body weights of the animals in these two groups were not different.
RNA Isolation, Library Construction, and RNA Sequencing (RNA-Seq)
Total RNA was isolated from each sample using TRIzol (Invitrogen, Carlsbad, CA, United States) according to manufacturer’s protocol. RNA concentration and purity were assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). The integrity of the RNA was assessed using an Agilent 2,100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, United States). Ribosomal RNA (rRNA) was removed using a Ribo-Zero Magnetic Gold Kit (Epicenter, Madison, WI, United States). The linear-stranded RNA was removed by RNase R (Epicenter, Madison, WI, United States). Sequencing libraries were generated using NEBNext® UltraTM Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, United States) following manufacturer’s recommendations. Sequencing of the libraries was performed on an Illumina Novaseq 6,000 system with PE150 sequencing mode by Novogene Co., Ltd. (Beijing, China).
CircRNA Sequencing Analysis
Clean reads were obtained by removing those composed of adapters, contain ploy-N, and low-quality sequences (containing more than 50% low-quality bases) from the raw data. All of the downstream analyses were based on the clean high-quality data. Clean reads were then aligned with the donkey reference genome (ASM303372v1) using the Burrows-Wheeler Aligner (BWA)-MEM (Li and Durbin, 2009). CIRI2 (Gao et al., 2018) software was applied to obtain back-spliced junction reads for circRNA prediction. The expression levels of individual circRNA were calculated as RPM (reads per million mapped reads). DESeq2 (Love et al., 2014) was used to identify differentially expressed (DE) circRNAs between the H and L groups. We defined circRNAs that had a fold change ≥ 2 or < 0.5 with a Benjamini-Hochberg method corrected p ≤ 0.05 between the two groups as significant differentially expressed circRNA.
GO and KEGG Pathway Enrichment Analysis
Gene Ontology (GO) enrichment analysis of the host genes of all differentially expressed circRNAs was performed by the clusterProfiler R package using default parameters (Yu et al., 2012). GO terms with a Benjamini-Hochberg method corrected p < 0.05 were considered significantly enriched. We used KOBAS software (Wu et al., 2006) with default parameters to test the statistical enrichment of all host genes of the differentially expressed circRNAs in KEGG pathways. Pathways with a Benjamini-Hochberg method corrected p < 0.05 were considered significantly enriched.
Prediction of CircRNA-miRNA Interactions
We used miRanda software (John et al., 2004) with “-sc 140 -en -10 -scale 4 -strict” to identify putative circRNA-miRNA interactions and Cytoscape software (Shannon et al., 2003) was used to construct the circRNA-miRNA networks.
Quantitative Reverse-Transcription PCR Analysis
Total RNA was extracted from longissimus dorsi muscle used in the RNA-seq using TRIzol (Invitrogen, Carlsbad, CA, United States), and then reverse-transcribed into complementary DNA (cDNA) using Primescript RT Master Kit (Takara, Dalian, China) according to the manufacturer’s instructions. According to the method of selecting DE circRNA in other studies (Zhang G. et al., 2019; Tian et al., 2020), we randomly select 5 DE circRNAs from the up-regulated and down-regulated DE circRNAs using sample function in R software, respectively. Primers for the differentially expressed circRNAs were designed using primer5 (Singh et al., 1998). All the primers used are listed in Supplementary Table S1. Quantitative Reverse-Transcription PCR (qRT-PCR) was performed with AceQ qPCR SYBR Green Master Mix (Vazyme, Nanjing, China) in a reaction volume of 20 μL. The cycling parameters were as follows: 95°C for 5 min, followed by 40 amplification cycles, each at 95°C for 10 s, then 60°C for 30 s. All reactions were performed in triplicate for each sample. Relative expression levels of the differentially expressed circRNAs were calculated by the 2–ΔΔCt method as reported previously (Wei et al., 2017; Cheng et al., 2019). The glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was used as a reference to normalize the relative expression of the circRNAs.
Statistical Analyses
Statistical analyses were conducted using SPSS 20.0 (SPSS Inc., Chicago, IL, United States). Statistical significance of the difference between the two groups was calculated using a Student’s t-test. A value of p < 0.05 was considered to represent a statistically significant difference.
ResultsCharacterization of CircRNAs in Longissimus Dorsi Muscles
To identify the circRNAs expression profile in donkey longissimus dorsi muscles, we performed RNA-seq and mined the data for circRNAs. In total, 663,324,766 raw reads (327,437,586 for H and 335,887,180 for L group) were generated from the six sequencing libraries (Supplementary Table S2). After removing low-quality, poly-N containing, and adapter-containing reads from the raw reads, 318,388,574 and 326,916,290 clean reads were obtained, for the H and L groups, respectively (Supplementary Table S2). An average of 82.49% (range: 78.58–85.87%) of the reads were mapped to the donkey genome (ASM303372v1) (Supplementary Table S3). A total of 12,727 circRNAs that were supported by at least two junction reads were identified from longissimus dorsi muscle and spliced from 3,251 genes (Supplementary Table S4). Approximately 70% of the circRNAs from both groups were derived from the exons of protein genes, while some of them were from intronic or intergenic regions (Figure 1A). The lengths of most of the identified circRNAs were less than 1,000 nucleotides (nt) (Figure 1B). Among the circRNA-producing genes, most host genes preferred to generate a single circRNA instead of multiple circRNA (Figure 1C). The majority of the identified circRNA species that originate from protein-coding genes contain two or three exons (Figure 1D).
Characteristics of the circRNAs identified in donkey muscle tissue. (A) Genomic origin of the identified circRNAs. (B) Length distribution of the identified circRNAs. (C) Distribution of the number of circRNA per gene. (D) The percentage of circRNAs was calculated based on the number of exons each circRNA in the exon-derived circRNAs.
Identification of DE CircRNAs Between H and L IMF Muscles
To identify candidate circRNAs affecting IMF content in donkey longissimus dorsi muscle, we calculated the expression level of each identified circRNAs and tested for differential expression of circRNAs between the H and L IMF content groups. Our results indicate that the density distribution of circRNA expression was not significantly different between the H and L IMF content groups (Figure 2A), however, 127 DE circRNAs were identified, including 63 upregulated and 64 downregulated circRNAs in the H IMF content group (Figures 2B,C). Detailed information on each DE circRNAs is provided in Supplementary Table S5. Among the candidates, novel_circ_0000323 was the most upregulated circRNA and novel_circ_0000319 was the most downregulated circRNA in the H IMF content group (Supplementary Table S5). Figure 2D shows a heatmap of the expression pattern of the DE circRNAs from the six samples, which indicates that the samples from the H IMF content group could be clearly distinguished from those from the L IMF content group.
Differential expression analysis of circRNAs between H and L IMF groups. (A) Density plot of the expression density distribution of circRNAs in each sample. (B) Volcano plot of all DE circRNAs between the H and L IMF groups. X-axis represents the value of log2 (H/L) and the Y-axis represents the value of –log10 (p-value). (C) Number of total, up-regulated and down-regulated circRNAs between the H and L IMF groups. (D) Clustered heatmap of the expression patterns of the DE circRNAs from three H and three L IMF samples.
Experimental Validation of the DE CircRNAs
To verify the reliability of the RNA-seq data, we use qRT-PCR to detect the expression levels of DE circRNAs. We tested the expression of randomly selected five upregulated circRNAs (novel_circ_0010172, novel_circ_0007969, novel_circ_0011073, novel_circ_0002126, and novel_circ_0010184) and downregulated circRNAs (novel_circ_0012311, novel_circ_0007411, novel_circ_0002621, novel_circ_0009905, and novel_circ_0002071). Divergent primers for each circRNA were designed to amplify the back-splice sequences (Figure 3A). PCR products of the divergent primers for each circRNA were confirmed by Sanger sequencing (Figure 3B and Supplementary Figure S1). Moreover, qRT-PCR results indicate that these circRNAs were significantly differentially expressed between H and L IMF groups (p < 0.05; Figure 3C), and the Pearson correlation coefficient of the log2(fold change) data between the qRT-PCR and RNA-Seq was 0.95 (p < 0.05; Figure 3D), suggesting that the DE circRNAs identified by RNA-Seq were reliable.
Validation of DE circRNAs by qRT-PCR. (A) Schematic of the divergent primer design for circRNAs. (B) Representative PCR products sequenced to validate the backsplice sequence of the circRNAs. Arrow represents the backsplice junction. (C) Expression level of DE circRNAs in H and L IMF samples were determined by qRT-PCR. *p < 0.05; **p < 0.01, two-tailed t-test. (D) Pearson correlation coefficient of the log2(fold change) data between qRT-PCR and RNA-Seq.
Enrichment Analysis of the Host Genes of DE CircRNAs
A previous study demonstrated that circRNA abundance is negatively correlated with their linear host gene mRNA and that there is a competition between pre-mRNA splicing and circRNA production (Ashwal-Fluss et al., 2014). To explore and analyze the potential biological function of these DE circRNAs, we performed GO and KEGG enrichment analysis of the host genes of the DE circRNAs. GO enrichment analysis indicated that 48, 37, and 36 GO terms were significantly enriched in biological process, cellular component, and molecular function, respectively (Supplementary Table S6). Among the most enriched GO terms, some of them were associated with lipid metabolism (e.g., fatty acid beta-oxidation using acyl-CoA dehydrogenase and MLL3/4 complex) (Figure 4A). Moreover, the KEGG enrichment analysis indicated that the host genes of the DE circRNAs were significantly enriched in 12 pathways, including related to lipid metabolism (e.g., TGF-beta and lysine degradation signaling pathway) (Figure 4B).
Gene Ontology (GO) and KEGG enrichment analysis of host genes of the DE circRNAs. (A) Most enriched GO terms in biological process, cellular component, and molecular function. X-axis represents GO terms and the Y-axis represents the value of -log10 (p-value). (B) Significantly enriched signaling pathways of the host genes of the DE circRNAs. The X-axis represents rich factor and the Y-axis represents pathway. Size and color of the bubble represent the number of host genes enriched in the pathway and enrichment significance, respectively.
Putative Functions of the DE CircRNAs as MiRNA Sponges
A previous study has shown that circRNA can act as miRNA sponges by acting as binding sites (Li X. et al., 2018). Some microRNAs (miRNAs) function as stimulators or inhibitors in adipocyte differentiation (Arner and Kulyte, 2015). To determine whether the DE cricRNA identified in this study can potentially function as miRNA sponges to regulate adipogenesis, we tested the ability of the DE circRNAs to bind to miRNAs. The circRNA-miRNA regulatory networks showed that 17,088 circRNA-miRNA interactions could be predicted based on 127 DE circRNAs and 690 miRNAs (Supplementary Table S7). Interestingly, we found that many of the upregulated circRNAs potentially function as sponges for many miRNA genes associated with adipogenesis such as miR-429, miR-224, miR-125a-5p, miR-223, miR-145, and miR-302a (Figure 5A). Additionally, many of the downregulated circRNAs also potentially bind miRNAs such as to miR-181a, miR-144, miR-199a-5p, and miR-127 (Figure 5B).
With the advancement of high-throughput sequencing technologies and bioinformatics, a large number of circRNAs have been identified in the human (Li S. et al., 2018), mouse (Legnini et al., 2017), chicken (Shen et al., 2019), pig (Hong et al., 2019; Li et al., 2020), cattle (Li H. et al., 2018), and sheep (Li X.Y. et al., 2019) genomes. Recently, the circRNA database CircAtlas1 has integrated over one million circRNAs from 6 species (human, macaca, mouse, rat, pig, and chicken) (Wu et al., 2020). However, the expression profiles and biological characteristics of circRNAs in the donkey has not been examined. Here we performed a comprehensive profile of circRNAs in the donkey using RNA-seq and identified 12,727 unique circRNAs in longissimus dorsi muscle samples. These results provide the foundation for establishing a donkey circRNA database and allowing future research on the potential role of circRNA in muscle development in the donkey. In this study, the major circRNAs were found to be derived from back-splicing of exons, which is consistent with previous studies in pig (Liang et al., 2017), chicken (Ouyang et al., 2017), and cattle (Wei et al., 2017) muscle. This result suggests that the circRNAs are highly similar in distribution pattern and evolutionarily conserved among different species. Moreover, our results show that the identified circRNAs contain multiple exons, with most containing two or three exons, which is similar to the circRNAs in many other species (Zhang et al., 2014; Wei et al., 2017), suggesting many host preferentially produce two or three circRNAs in donkey.
Previous studies have demonstrated that circRNAs are a large class of widespread and diverse endogenous RNAs in animals that play critical roles in many physiological and pathological processes, including fat (Jiang et al., 2020), muscle (Zhang P. et al., 2019), and cancer (Xie et al., 2020) development. For example, Jiang et al. (2020) characterized the expression profile of circRNA in fat tissue during developmental stages of the calf and adult cow and demonstrated that the circFUT10 regulates adipocyte cell proliferation and differentiation by sponging the miRNA let-7c (Jiang et al., 2020). Moreover, several studies have found that circRNAs are involved in muscle traits in the pig (Liang et al., 2017; Sun et al., 2017) and chicken (Chen et al., 2019). Recent studies have shown a potential role for circRNAs in IMF deposition in the yak based on the differentially expression of circRNAs between muscle and adipose tissue (Wang et al., 2020). In this study, we obtained the expression profiles of circRNA from muscle with differing IMF content, and identified 127 differentially expressed circRNAs between three H and L IMF samples. However, the expression level of these DE circRNAs in a large number of individuals are still needed to be confirmed. Furthermore, 10 DE circRNAs detected by qRT-PCR are significantly differentially expressed between H and L groups, indicating the DE circRNA results are reliable. These results suggest that these DE circRNAs play an important role in IMF deposition, and that circRNAs may affect livestock meat quality through the formation of IMF.
Since circRNAs are produced from precursor mRNA (pre-mRNA) by back-splicing, the production of circRNA can reduce the amount of linear transcripts by competing with linear splicing (Ashwal-Fluss et al., 2014; Liu et al., 2018). This suggests that circRNAs can act in an important role under physiological conditions to regulate the expression of their host genes. In this study, GO enrichment of host genes of DE circRNAs yielded some GO terms that are involved in lipid metabolism including MLL3/4 complex, fatty acid beta-oxidation using acyl-CoA dehydrogenase, fatty acid metabolic process, methyltransferase complex and histone deacetylase binding. Histone methyltransferases MLL3/4 play a crucial role in adipogenesis (Lee et al., 2020) and the deletion of MLL3 in mice results in significantly decreased amounts of white fat (Lee et al., 2008). For example, novel_circ_0009716, a down-regulated circRNA from the H IMF group, is derived from the histone mono-methyltransferase MLL4 (KMT2D), suggesting that it can regulate IMF formation through the MLL3/4 complex. Moreover, novel_circ_0007020 is derived from the MIER1 gene, which is enriched in the histone deacetylase binding GO term. A previous study indicated that MIER1 can bind to BAHD1 to form a hub for histone deacetylases and methyltransferases that are involved in energy metabolism (Lakisic et al., 2016). Therefore, novel_circ_0007020 might be involved in IMF formation by repressing the levels of its linear mRNA for MIER1. In addition, we found that host genes of the differently expressed circRNAs were significantly enriched in 12 pathways, including related to lipid metabolism (e.g., TGF-beta and lysine degradation signaling pathway). TGFβ and its downstream effector have been reported to play important roles in regulating glucose and energy homeostasis (Yadav et al., 2011). Previous study has shown that increased IMF deposition is due to lysine restriction in pig diets (Madeira et al., 2013). These evidences indicate that DE cricRNAs can regulate TGF-beta and lysine degradation signaling pathway by competing with host mRNA splicing, which ultimately affect IMF deposition. Together, all these results suggest that these differentially expressed circRNAs affect the physiological functions of their host genes leading to the regulation of IMF deposition.
It has been shown that circRNA can act as competing endogenous RNAs (ceRNAs) to sponge miRNAs and prevent them from binding and suppressing their target mRNAs (Chen, 2020). For example, circHIPK2 regulates astrocyte activation, via regulation of autophagy and endoplasmic reticulum (ER) stress, through targeting of MIR124-2HG and SIGMAR1 (Huang et al., 2017). We identified 17,088 circRNA-miRNA interactions in this study and these results suggest that DE circRNAs have a regulatory role in IMF deposition via the sequestration of miRNAs. MicroRNAs (miRNAs) are a class of small non-coding RNA that regulate target gene expression, and can stimulate or inhibit the differentiation of adipocytes in adipose tissue (Arner and Kulyte, 2015). Our current study revealed that many DE circRNAs may act as endogenous sponges of miRNAs including miR-429 (Peng et al., 2016), miR-224 (Zhang Y. et al., 2019), miR-125a-5p (Du et al., 2018), miR-223 (Li F. et al., 2019), miR-145 (Guo et al., 2012), miR-302a (Jeong et al., 2014), miR-181a (Zhang Z. et al., 2019), miR-144 (Shen et al., 2018), miR-199a-5p (Shi et al., 2014), and miR-127 (Gao et al., 2019) that are involved in the regulation of adipogenesis. miR-125a-5p is negatively involved with IMF content by targeting KLF13 and ELOVL6 (Du et al., 2018). We predicted that fifteen upregulated DE circRNA, such as novel_circ_0010184, can able to bind miR-125a-5p. Moreover, novel_circ_0010184 was upregulated in H IMF content compared to L IMF content using qRT-PCR in this study. Therefore, these results suggest that these circRNAs might be a sponge of miR-125a-5p and promote increases in IMF content. A downregulated DE circRNA, novel_circ_0012311, might function as sponge for miR-181a that facilitates porcine preadipocyte differentiation by targeting TGFBR1, indicating that novel_circ_0012311 might inhibit IMF preadipocyte differentiation through the miR-181a/TGFBR1 axis in donkey. Consequently, our results suggest that DE circRNAs have important roles in the regulation of IMF adipogenesis by acting as miRNA sponges to inhibit the miRNA targeting of genes.
Conclusion
We identified 12,727 circRNAs in the donkey longissimus dorsi muscle, thus, expanding our understanding of the complexity of the donkey transcriptome. The host genes were enriched in lipid metabolism related GO terms and signaling pathways, and DE circRNAs act as adipogenesis related miRNAs sponge, which provide insight into our understanding of the functions of circRNA in IMF content. Our results provide evidence that circRNAs have an important regulatory role in IMF content and this information might be useful for future research into circRNA and the regulation of IMF deposition.
Data Availability Statement
The data used in our study have been deposited in NCBI SRA (accession codes PRJNA647167).
Ethics Statement
The animal study was reviewed and approved by the Ethics Committee and Experimental Animal Committee of Shenyang Agriculture University. Written informed consent was obtained from the owners for the participation of their animals in this study.
Author Contributions
SZha conceived the study. BL, CF, SZhu, JZ, and XZ conducted the experiments. ZW analyzed the data. BL wrote the manuscript. DI revised the manuscript. All authors contributed to the article and approved the submitted version.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Funding. This work was supported by the grants from the Educational Department of Liaoning Province (Climbing Scholar), the Organization Department of Liaoning Provincial Committee of China (No. XLYC1907018), and the Shenyang Agricultural University Research Start-up Funding (No. 880418062).
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.587559/full#supplementary-material
ReferencesArnerP.KulyteA. (2015). MicroRNA regulatory networks in human adipose tissue and obesity.Nat. Rev. Endocrinol.11276–288. 10.1038/nrendo.2015.2525732520Ashwal-FlussR.MeyerM.PamudurtiN. R.IvanovA.BartokO.HananM. (2014). circRNA biogenesis competes with pre-mRNA splicing.Mol. Cell.5655–66. 10.1016/j.molcel.2014.08.01925242144BahelkaI.OravcovaM.PeskovicovaD.TomkaJ.HanusovaE.LahuckyR. (2009). Comparison of accuracy of intramuscular fat prediction in live pigs using five different ultrasound intensity levels.Animal31205–1211. 10.1017/S175173110900448022444851CaoZ. B.GaoD.XuT. T.ZhangL.TongX.ZhangD. D. (2019). Circular RNA profiling in the oocyte and cumulus cells reveals that circARMC4 is essential for porcine oocyte maturation.Aging Us118015–8034. 10.18632/aging.10231531562810ChenB.YuJ.GuoL.ByersM. S.WangZ.ChenX. (2019). Circular RNA circHIPK3 promotes the proliferation and differentiation of chicken myoblast cells by sponging miR-30a-3p.Cells8:177. 10.3390/cells802017730791438ChenL. L. (2020). The expanding regulatory mechanisms and cellular functions of circular RNAs.Nat. Rev. Mol. Cell. Biol.21475–490. 10.1038/s41580-020-0243-y32366901ChengZ.YuC.CuiS.WangH.JinH.WangC. (2019). circTP63 functions as a ceRNA to promote lung squamous cell carcinoma progression by upregulating FOXM1.Nat. Commun.10:3200. 10.1038/s41467-019-11162-431324812ChoI. C.ParkH. B.AhnJ. S.HanS. H.LeeJ. B.LimH. T. (2019). A functional regulatory variant of MYH3 influences muscle fiber-type composition and intramuscular fat content in pigs.PLoS Genet.15:e1008279. 10.1371/journal.pgen.100827931603892De PaloP.MaggiolinoA.MilellaP.CentoducatiN.PapaleoA.TateoA. (2016). Artificial suckling in Martina Franca donkey foals: effect on in vivo performances and carcass composition.Trop. Anim. Health Prod.48167–173. 10.1007/s11250-015-0940-226510946DevauxY.CreemersE. E.BoonR. A.WerfelS.ThumT.EngelhardtS. (2017). Circular RNAs in heart failure.Eur. J. Heart Fail.19701–709. 10.1002/ejhf.80128345158DongY. P.HeD.PengZ. Z.PengW.ShiW. W.WangJ. (2017). Circular RNAs in cancer: an emerging key player.J. Hematol. Oncol.10:2. 10.1186/s13045-016-0370-228049499DuJ.XuY.ZhangP.ZhaoX.GanM.LiQ. (2018). MicroRNA-125a-5p affects adipocytes proliferation, differentiation and fatty acid composition of porcine intramuscular fat.Int. J. Mol. Sci.19:501. 10.3390/ijms1902050129414921GaoY.WangY.ChenX.PengY.ChenF.HeY. (2019). MiR-127 attenuates adipogenesis by targeting MAPK4 and HOXC6 in porcine adipocytes.J. Cell. Physiol.23421838–21850. 10.1002/jcp.2866031206669GaoY.ZhangJ.ZhaoF. (2018). Circular RNA identification based on multiple seed matching.Brief. Bioinform.19803–810. 10.1093/bib/bbx01428334140GuoY.ChenY.ZhangY.ZhangY.ChenL.MoD. (2012). Up-regulated miR-145 expression inhibits porcine preadipocytes differentiation by targeting IRS1.Int. J. Biol. Sci.81408–1417. 10.7150/ijbs.459723197937HocquetteJ. F.GondretF.BaezaE.MedaleF.JurieC.PethickD. W. (2010). Intramuscular fat content in meat-producing animals: development, genetic and nutritional control, and identification of putative markers.Animal4303–319. 10.1017/S175173110999109122443885HongL.GuT.HeY.ZhouC.HuQ.WangX. (2019). Genome-wide analysis of circular RNAs mediated ceRNA regulation in porcine embryonic muscle development.Front. Cell. Dev. Biol.7:289. 10.3389/fcell.2019.0028931803743HuangR.ZhangY.HanB.BaiY.ZhouR.GanG. (2017). Circular RNA HIPK2 regulates astrocyte activation via cooperation of autophagy and ER stress by targeting MIR124-2HG.Autophagy131722–1741. 10.1080/15548627.2017.135697528786753JeongB. C.KangI. H.KohJ. T. (2014). MicroRNA-302a inhibits adipogenesis by suppressing peroxisome proliferator-activated receptor gamma expression.FEBS Lett.5883427–3434. 10.1016/j.febslet.2014.07.03525109777JiangR.LiH.YangJ.ShenX.SongC.YangZ. (2020). circRNA profiling reveals an abundant circFUT10 that promotes adipocyte proliferation and inhibits adipocyte differentiation via sponging let-7.Mol. Ther. Nucleic Acids20491–501. 10.1016/j.omtn.2020.03.01132305019JohnB.EnrightA. J.AravinA.TuschlT.SanderC.MarksD. S. (2004). Human MicroRNA targets.PLoS Biol.2:e363. 10.1371/journal.pbio.002036315502875LakisicG.LebretonA.PourpreR.WendlingO.LibertiniE.RadfordE. J. (2016). Role of the BAHD1 chromatin-repressive complex in placental development and regulation of steroid metabolism.PLoS Genet.12:e1005898. 10.1371/journal.pgen.100589826938916LeeJ.SahaP. K.YangQ. H.LeeS.ParkJ. Y.SuhY. (2008). Targeted inactivation of MLL3 histone H3-Lys-4 methyltransferase activity in the mouse reveals vital roles for MLL3 in adipogenesis.Proc. Natl. Acad. Sci. U.S.A.10519229–19234. 10.1073/pnas.081010010519047629LeeJ. E.ChoY. W.DengC. X.GeK. (2020). MLL3/MLL4-associated PAGR1 regulates adipogenesis by controlling induction of C/EBPbeta and C/EBPdelta.Mol. Cell. Biol.40:e00209-20. 10.1128/MCB.00209-2032601106LegniniI.Di TimoteoG.RossiF.MorlandoM.BrigantiF.SthandierO. (2017). Circ-ZNF609 is a circular RNA that can be translated and functions in myogenesis.Mol. Cell.6622–26. 10.1016/j.molcel.2017.02.01728344082LiB.WengQ.DongC.ZhangZ.LiR.LiuJ. (2018). A key gene, PLIN1, can affect porcine intramuscular fat content based on transcriptome analysis.Genes9:194. 10.3390/genes904019429617344LiB.YinD.LiP.ZhangZ.ZhangX.LiH. (2020). Profiling and functional analysis of circular RNAs in porcine fast and slow muscles.Front. Cell. Dev. Biol.8:322. 10.3389/fcell.2020.0032232528948LiF.LiD.ZhangM.SunJ.LiW.JiangR. (2019). miRNA-223 targets the GPAM gene and regulates the differentiation of intramuscular adipocytes.Gene685106–113. 10.1016/j.gene.2018.10.05430389563LiH.DurbinR. (2009). Fast and accurate short read alignment with burrows-wheeler transform.Bioinformatics251754–1760. 10.1093/bioinformatics/btp32419451168LiH.WeiX. F.YangJ. M.DongD.HaoD.HuangY. Z. (2018). circFGFR4 Promotes Differentiation of Myoblasts via Binding miR-107 to Relieve Its Inhibition of Wnt3a.Mol. Ther. Nucleic Acids11272–283. 10.1016/j.omtn.2018.02.01229858062LiS.MaY.TanY.MaX.ZhaoM.ChenB. (2018). Profiling and functional analysis of circular RNAs in acute promyelocytic leukemia and their dynamic regulation during all-trans retinoic acid treatment.Cell Death Dis.9:651. 10.1038/s41419-018-0699-229844435LiX. Y.LiC. Y.WeiJ. C.NiW.XuY. R.YaoR. (2019). Comprehensive expression profiling analysis of pituitary indicates that circRNA participates in the regulation of sheep estrus.Genes10:90. 10.3390/genes1002009030696117LiX.YangL.ChenL. L. (2018). The biogenesis, functions, and challenges of circular RNAs.Mol. Cell.71428–442. 10.1016/j.molcel.2018.06.03430057200LiangG.YangY.NiuG.TangZ.LiK. (2017). Genome-wide profiling of Sus scrofa circular RNAs across nine organs and three developmental stages.DNA Res.24523–535. 10.1093/dnares/dsx02228575165LiuX.LiuK. Q.ShanB. S.WeiS. J.LiD. F.HanH. Y. (2018). A genome-wide landscape of mRNAs, lncRNAs, and circRNAs during subcutaneous adipogenesis in pigs.J. Anim. Sci. Biotechnol.9:76. 10.1186/s40104-018-0292-730410752LorenzoJ. M.SarriesM. V.TateoA.PolidoriP.FrancoD.LanzaM. (2014). Carcass characteristics, meat quality and nutritional value of horsemeat: a review.Meat Sci.961478–1488. 10.1016/j.meatsci.2013.12.00624423453LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.Genome Biol.15:550. 10.1186/s13059-014-0550-825516281MadeiraM. S.CostaP.AlfaiaC. M.LopesP. A.BessaR. J.LemosJ. P. (2013). The increased intramuscular fat promoted by dietary lysine restriction in lean but not in fatty pig genotypes improves pork sensory attributes.J. Anim. Sci.913177–3187. 10.2527/jas.2012-542423572257NogiT.HondaT.MukaiF.OkagakiT.OyamaK. (2011). Heritabilities and genetic correlations of fatty acid compositions in longissimus muscle lipid with carcass traits in Japanese Black cattle.J. Anim. Sci.89615–621. 10.2527/jas.2009-230021036930OuyangH.ChenX.WangZ.YuJ.JiaX.LiZ. (2017). Circular RNAs are abundant and dynamically expressed during embryonic muscle development in chickens.DNA Res.2571–86. 10.1093/dnares/dsx03929036326PenaR. N.Ros-FreixedesR.TorM.EstanyJ. (2016). Genetic marker discovery in complex traits: a field example on fat content and composition in pigs.Int. J. Mol. Sci.17:2100. 10.3390/ijms1712210027983643PengY.ChenF. F.GeJ.ZhuJ. Y.ShiX. E.LiX. (2016). miR-429 inhibits differentiation and promotes proliferation in porcine preadipocytes.Int. J. Mol. Sci.17:2047. 10.3390/ijms1712204727941616PolidoriP.PucciarelliS.ArianiA.PolzonettiV.VincenzettiS. (2015). A comparison of the carcass and meat quality of Martina Franca donkey foals aged 8 or 12 months.Meat Sci.1066–10. 10.1016/j.meatsci.2015.03.01825863189PolidoriP.VincenzettiS.CavallucciC.BeghelliD. (2008). Quality of donkey meat and carcass characteristics.Meat Sci.801222–1224. 10.1016/j.meatsci.2008.05.02722063861Rybak-WolfA.StottmeisterC.GlazarP.JensM.PinoN.GiustiS. (2015). Circular RNAs in the mammalian brain are highly abundant, conserved, and dynamically expressed.Mol. Cell58870–885. 10.1016/j.molcel.2015.03.02725921068ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks.Genome Res.132498–2504. 10.1101/gr.123930314597658ShenL.LiQ.WangJ.ZhaoY.NiuL.BaiL. (2018). miR-144-3p Promotes Adipogenesis Through Releasing C/EBPalpha From Klf3 and CtBP2.Front. Genet.9:677. 10.3389/fgene.2018.0067730619490ShenX.LiuZ.CaoX.HeH.HanS.ChenY. (2019). Circular RNA profiling identified an abundant circular RNA circTMTC1 that inhibits chicken skeletal muscle satellite cell differentiation by sponging miR-128-3p.Int. J. Biol. Sci.152265–2281. 10.7150/ijbs.3641231592238ShiX. E.LiY. F.JiaL.JiH. L.SongZ. Y.ChengJ. (2014). MicroRNA-199a-5p affects porcine preadipocyte proliferation and differentiation.Int. J. Mol. Sci.158526–8538. 10.3390/ijms1505852624830555SinghV. K.MangalamA. K.DwivediS.NaikS. (1998). Primer premier: program for design of degenerate primers from a protein sequence.Biotechniques24318–319. 10.2144/98242pf029494736SunJ.XieM.HuangZ.LiH.ChenT.SunR. (2017). Integrated analysis of non-coding RNA and mRNA expression profiles of 2 pig breeds differing in muscle traits.J. Anim. Sci.951092–1103. 10.2527/jas.2016.086728380516SunL. F.ZhangB.ChenX. J.WangX. Y.ZhangB. W.JiY. Y. (2019). Circular RNAs in human and vertebrate neural retinas.RNA Biol.16821–829. 10.1080/15476286.2019.159103430874468TianJ.FuY.LiQ.XuY.XiX.ZhengY. (2020). Differential expression and bioinformatics analysis of CircRNA in PDGF-BB-induced vascular smooth muscle cells.Front. Genet.11:530. 10.3389/fgene.2020.0053032547599WangH.ZhongJ.ZhangC.ChaiZ.CaoH.WangJ. (2020). The whole-transcriptome landscape of muscle and adipose tissues reveals the ceRNA regulation network related to intramuscular fat deposition in yak.BMC Genomics21:347. 10.1186/s12864-020-6757-z32381004WeiX. F.LiH.YangJ. M.HaoD.DongD.HuangY. Z. (2017). Circular RNA profiling reveals an abundant circLMO7 that regulates myoblasts differentiation and survival by sponging miR-378a-3p.Cell Death Dis.8:e3153. 10.1038/cddis.2017.54129072698WuJ.MaoX.CaiT.LuoJ.WeiL. (2006). KOBAS server: a web-based platform for automated annotation and pathway identification.Nucleic Acids Res.34W720–W724. 10.1093/nar/gkl16716845106WuW.JiP.ZhaoF. (2020). CircAtlas: an integrated resource of one million highly accurate circular RNAs from 1070 vertebrate transcriptomes.Genome Biol.21:101. 10.1186/s13059-020-02018-y32345360XiaX.YuJ.ZhaoX.YaoY.ZengL.AhmedZ. (2019). Genetic diversity and maternal origin of Northeast African and South American donkey populations.Anim. Genet.50266–270. 10.1111/age.1277430854699XieM.YuT.JingX.MaL.FanY.YangF. (2020). Exosomal circSHKBP1 promotes gastric cancer progression via regulating the miR-582-3p/HUR/VEGF axis and suppressing HSP90 degradation.Mol. Cancer19:112. 10.1186/s12943-020-01208-332600329YadavH.QuijanoC.KamarajuA. K.GavrilovaO.MalekR.ChenW. (2011). Protection from obesity and diabetes by blockade of TGF-beta/Smad3 signaling.Cell. Metab.1467–79. 10.1016/j.cmet.2011.04.01321723505YuG.WangL. G.HanY.HeQ. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters.OMICS16284–287. 10.1089/omi.2011.011822455463ZhangG.DiaoS.ZhangT.ChenD.HeC.ZhangJ. (2019). Identification and characterization of circular RNAs during the sea buckthorn fruit development.RNA Biol.16354–361. 10.1080/15476286.2019.157416230681395ZhangH. M.XiaH. L.JiangH. R.MaoY. J.QuK. X.HuangB. Z. (2018). Longissimus dorsi muscle transcriptomic analysis of Yunling and Chinese simmental cattle differing in intramuscular fat content and fatty acid composition.Genome61549–558. 10.1139/gen-2017-016429883552ZhangM.HanY.ZhaiY.MaX.AnX.ZhangS. (2020). Integrative analysis of circRNAs, miRNAs, and mRNAs profiles to reveal ceRNAs networks in chicken intramuscular and abdominal adipogenesis.BMC Genomics21:594. 10.1186/s12864-020-07000-332847498ZhangP.ChaoZ.ZhangR.DingR.WangY.WuW. (2019). Circular RNA regulation of myogenesis.Cells8:885. 10.3390/cells808088531412632ZhangX. O.WangH. B.ZhangY.LuX.ChenL. L.YangL. (2014). Complementary sequence-mediated exon circularization.Cell159134–147. 10.1016/j.cell.2014.09.00125242744ZhangX.LiH.YuJ.ZhouX.JiC.WuS. (2019). Label-free based comparative proteomic analysis of whey proteins between different milk yields of Dezhou donkey.Biochem. Biophys. Res. Commun.508237–242. 10.1016/j.bbrc.2018.11.13030482389ZhangY.WangY.WangH.MaX.ZanL. (2019). MicroRNA-224 impairs adipogenic differentiation of bovine preadipocytes by targeting LPL.Mol. Cell. Probes4429–36. 10.1016/j.mcp.2019.01.00530703449ZhangZ.GaoY.XuM. Q.WangC. J.FuX. H.LiuJ. B. (2019). miR-181a regulate porcine preadipocyte differentiation by targeting TGFBR1.Gene68145–51. 10.1016/j.gene.2018.09.04630266501ZouC.LiL.ChengX.LiC.FuY.FangC. (2018). Identification and functional analysis of long intergenic non-coding RNAs underlying intramuscular fat content in pigs.Front. Genet.9:102. 10.3389/fgene.2018.0010229662503