CircRNA Expression Profile during Yak Adipocyte Differentiation and Screen Potential circRNAs for Adipocyte Differentiation

The yak (Bos grunniens) is subjected to nutritional deficiency during the whole winter grazing season; deciphering the adipose metabolism and energy homeostasis under cold and nutrients stress conditions could be a novel way to understand the specific mechanism of energy metabolism. Circular RNAs (circRNAs) have elucidated that they play a key role in many biological events, but the regulatory function of adipose development remains mostly unknown. Therefore, the expression pattern of circRNAs were identified for the first time during yak adipocyte differentiation to gain insight into their potential functional involvement in bovine adipogenesis. We detected 7203 circRNA candidates, most of them contained at least two exons, and multiple circRNA isoforms could be generated from one parental gene. Analysis of differential expression circRNAs displayed that 136 circRNAs were differentially expressed at day 12 (Ad) after adipocyte differentiation, compared with the control at day 0 (Pread 0), while 7 circRNAs were detected on day 2. Sanger sequencing validated that six circRNAs had head-to-tail junction, and quantitative real-time PCR (qPCR) results revealed that the expression patterns of ten circRNAs were consistent with their expression levels from RNA-sequencing (RNA-seq) data. We further predicted the networks of circRNA-miRNA-gene based on miRNAs sponging by circRNAs, in which genes were participated in the adipocyte differentiation-related signaling pathways. After that, we constructed several adipocyte differentiation-related ceRNAs and revealed six circRNAs (novel_circ_0009127, novel_circ_0000628, novel_circ_0011513, novel_circ_0010775, novel_circ_0006981 and novel_circ_0001494) were related to adipogenesis. Furthermore, we analyzed the homology among yak, human and mouse circRNAs and found that 3536 yak circRNAs were homologous to human and mouse circRNAs. In conclusion, these findings provide a solid basis for the investigation of yak adipocyte differentiation-related circRNAs and serve as a great reference to study the energy metabolism of high-altitude animals.


Introduction
Yak is an iconic symbol of the Qinghai Tibetan Plateau [1]. Currently, there are more than 15 million yaks living on the Qinghai Tibetan Plateau, which represents about 90 percent of the world's yak population. Yaks are essential for Tibetans and other nomadic pastoralists living in high-altitude surroundings and furnish the primary resources such as meat, milk, transportation, dung

Ethics Statement
The handling of animals during the experiment was carried out in strict accordance with the Animal Ethics Procedures and Guidelines of the People's Republic of China, and all procedures were approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences of Chinese Academy of Agricultural Sciences (Permit No. SYXK-2014-0002).

Isolation of Pre-Adipocytes
Three healthy Datong yaks (3 days old) were offered by the Datong Yak Breeding Center (Datong County, Qinghai, China). Subcutaneous adipose tissue for cell isolation was harvested from three yaks at a local slaughterhouse and transported to the laboratory within sterile phosphate buffer saline (PBS) (HyClone, Thermo Fisher Scientific, Carlsbad, CA, USA) supplemented with 1% antibiotics (penicillin-streptomycin). Under a sterile environment, the potentially polluted epidermis, blood vessels and connective tissue were carefully removed. Subsequently, the remaining subcutaneous adipose tissue was washed with 1% antibiotic PBS several times. Then, samples were minced mechanically into approximately 1 mm 3 section under sterile environment. The sections were digested by collagenase Type I for about 60-90 min at 37 • C with constant agitation. The digested tissues were filtered by 40 µm nylon mesh screen, and the filtrates were centrifuged at 1400 × g for 5 min. Subsequently, pellets were incubated with erythrocyte lysis buffer (0.154 M NH 4 Cl, 10 mM KHCO 3, 0.1 mM EDTA) at room temperature for 10 min and filtered with 200 µm mesh screen and washed twice with serum-free medium. The adipocytes were collected by centrifugation at 1400 × g for 5 min, cell counting was done with a haemocytometer, and then resuspended in pre-adipocyte growth media (DMEM-F12 supplemented with 10% fetal bovine serum). Thereafter, pre-adipocytes were maintained in 5% CO 2 at 37 • C. After about 1 week, the cells were proliferated as pre-adipocyte naturally.

Adipogenic Differentiation and Oil Red O Staining Identification
In order to induce adipogenic differentiation, cells reaching to confluence in growth media were induced with adipogenic agents composed of 3-isobuty-methylxanthine (MIX) (Sigma, St Louis, MO,USA), dexamethasone (Sigma, St Louis, MO, USA), rosiglitazone (Sigma, St Louis, MO, USA) and insulin (Sigma, St Louis, MO, USA) for 2 days. Then, cells were fed with the culture medium, which was changed between 2-3 days intervals. Generally, cells were washed twice with PBS and fixed with 4% formalin for 1 h. Cells were stained with a saturated Oil Red O solution for 30 min at room temperature. Subsequently, the cells were washed three times with distilled water, and images were taken with light microscopy.

High-throughput RNA Sequencing of circRNAs Relevant to Yak Adipocyte Differentiation
TRIzol reagent was used to extract total RNA from in vitro cultured yak pre-adipocytes and differentiated adipocytes at days 0, 2 and 12 from three biological replicates for each condition. The quality and quantity of total RNA were analyzed by Qubit™ RNA Assay Kit in Qubit™ 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and NanoPhotometer ® spectrophotometer (Implen, Westlake Village, CA, USA), respectively. For one sample, approximately 5 µg RNA was used for RNA sample preparations. Ribozero™ rRNA Removal Kit (Epicentre, USA) was used to remove ribosomal RNA from the total RNA. To digest the linear RNA 3 U of RNase R was used. According to the manufacturer's recommendations, NEBNext ® Ultra™ Directional RNA Library Prep Kit (NEB, Ipswich, MA, USA) was used to generate sequencing libraries for Illumina ® . The libraries were sequenced on an Illumina Hiseq 4000 platform (Novogene, Beijing, China). For quality control, Q20, Q30 and GC contents of clean data were calculated. Subsequently, find_circ and CIRI2 [10,13] were performed to detect and identify circRNAs. The Circos software [14] was used to construct the circos figure. Back-splice algorithm picked out the junctions of unmapped reads. These were considered the reference sequence for next analysis and their expression was detected by TPM (transcripts per million clean tags) [15]. TPM was used to normalize the raw counts. Normalized expression level = (readCount × 1,000,000)/libsize (libsize is the sum of circRNA read count). The data have been uploaded to the Sequence Read Archive (SRA) database. The valid accession number is PRJNA550036.

Annotating Host Liner Transcript and Differentially Expressed circRNAs
Linear transcripts were annotated based on the location of the chromosome where the circRNA sequence was overlapped by HISAT2 (version 2.0.4) [16]. CircRNA distribution in the genome can be explored by comparing circRNA with genetic elements. Based on the negative binomial distribution test, the DESeq R package (1.10.1) [17] was used to detect the differential expression analysis of circRNAs. The same circRNAs in two yak adipocyte differentiation groups were used at day 0 (Pread 0) group as a control. The resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. CircRNAs with an adjusted p-value < 0.05 were assigned as differentially expressed.

qRT-PCR Analysis for the Differential Expression of circRNAs
In order to verify the reliability of our analyzed data, ten differentially expressed circRNAs were selected. Total RNA was extracted from yak pre-adipocytes and differentiated adipocytes at days 0, Genes 2020, 11, 414 4 of 16 2, and 12 by using TRIzol reagent (Invitrogen, Carlsbad CA, USA) according to the manufacturer's protocol. The expression levels of circRNAs were tested by qRT-PCR. The qRT-PCR was performed with a total volume of 20 µL containing 10 µL 2 x SYBR Premix ExTaq II (Takara, China), 0.5 µL each primer (10 µM) and 1 µL diluted cDNA. PCR conditions were as follows: 95 • C for 30 s; 40 cycles of 94 • C for 15 s, primer specific Tm for 30 s, and 72 • C for 30 s. The 2 −∆∆Ct method was used to calculate the relative expression profiles of circRNAs [18] with β-actin as a reference gene. CircRNAs amplified primers are listed in Supplementary Table S1.

Analysis of Enrichment and microRNA Target Site
GO seq R package was used for the Gene Ontology (GO, http://www.geneontology.org/) enrichment analysis for the host genes that showed differential expression of circRNAs [19]. If p-value ≤ 0.05, GO terms were regarded as significant enrichment. Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/) database is a main resource for understanding high-level functions and utilities of biological systems. KOBAS software [20] was used to examine the statistical enrichment of genes or circRNA host genes that showed differential expression in KEGG pathways. MiRanda (http://www.microrna.org/) was used to identify the microRNA target site that located in exons of circRNA loci.

Co-expression Network Analysis of circRNA-miRNA-Gene during Yak Adipocyte Differentiation
The network of circRNA-miRNA-gene was established according to miRNA target site in exons of circRNA loci analysis. The putative correlations between circRNA, miRNA, and gene were ranked by miRnada and based on the hypergeometric distribution's p-value. The networks were constructed by Cytoscape software (wersion 3.5.1, BiNGO plug-in) [21] in which circle nodes represented miRNAs and rectangle nodes represented genes and triangle nodes represented circRNAs.

Homology Analysis of circRNAs among Yak, Human and Mouse
The analysis of homology among yak, human and mouse was performed by BLAST (e-value < 1.0 × e −5 ) [22]. The circRNAs data of human and mouse were downloaded from circBase (http://www.circbase.org).

High-throughput RNA Sequencing of Yak Pre-adipocytes and Adipocytes
After 12 days of induction with adipogenic agents, lipid droplets visibility in adipocytes increased significantly, as stained with Oil Red O ( Figure 1A-C). Adipocyte-specific marker genes PPARγ, C/EBPα, and FABP4, showed markedly higher expression on day 12 as compared to earlier stages (day 0, day 2) ( Figure 1D), indicating that our model system met the functional criteria for yak adipocyte differentiation.
In order to comprehensively identify the potential function of circRNAs during adipogenesis, high-throughput sequencing was performed to investigate the profiles of circRNAs by following a previously described pipeline [10]. The total numbers of reads obtained from nine samples are summarized in Supplementary Table S2. In total, 7203 circRNAs were identified during yak adipocyte differentiation at days 0, 2, and 12, of which 4139 circRNAs were detected at day 0, 4755 at day 2, 5505 at day 12, and 2737 at all three-time points (Figure 2A). The full length of adipocyte differentiation-related circRNAs ranged from 150 nt to over 15,000 nt. Among 7203 circRNAs, about 65.2% had a length of less than 15,000 nt, 8.61% had a length of less than 1000 nt, 29.64 % were between 1000 nt and 5000 nt in length, and 21.77% were between 5000 nt and 10,000 nt in length ( Figure 2B). The expression abundance of circRNAs was measured based on TPM [15] during adipocyte differentiation, and the results indicated no abnormal expression in all the three groups ( Figure 2C). This was highly consistent with the nine tested samples, as shown in Figure 2D.

High-throughput RNA Sequencing of Yak Pre-adipocytes and Adipocytes
After 12 days of induction with adipogenic agents, lipid droplets visibility in adipocytes increased significantly, as stained with Oil Red O ( Figure 1A-C). Adipocyte-specific marker genes PPARγ, C/EBPα, and FABP4, showed markedly higher expression on day 12 as compared to earlier stages (day 0, day 2) (Figure 1D), indicating that our model system met the functional criteria for yak adipocyte differentiation. In order to comprehensively identify the potential function of circRNAs during adipogenesis, high-throughput sequencing was performed to investigate the profiles of circRNAs by following a previously described pipeline [10]. The total numbers of reads obtained from nine samples are summarized in Supplementary Table S2. In total, 7203 circRNAs were identified during yak adipocyte differentiation at days 0, 2, and 12, of which 4139 circRNAs were detected at day 0, 4755 at day 2, 5505 at day 12, and 2737 at all three-time points (Figure 2A). The full length of adipocyte differentiation-related circRNAs ranged from 150 nt to over 15,000 nt. Among 7203 circRNAs, about 65.2% had a length of less than 15,000 nt, 8.61% had a length of less than 1000 nt, 29.64 % were between 1000 nt and 5000 nt in length, and 21.77% were between 5000 nt and 10,000 nt in length ( Figure 2B). The expression abundance of circRNAs was measured based on TPM [15] during adipocyte differentiation, and the results indicated no abnormal expression in all the three groups ( Figure 2C). This was highly consistent with the nine tested samples, as shown in Figure 2D.

Annotation of Host Linear Transcripts during Yak Adipocyte Differentiation
CircRNAs biogenesis occurs through back splicing and then the canonical splicing of linear RNAs, and they are commonly formed by the alternative splicing of pre-mRNA. We annotated linear transcripts of circRNAs from the correlative genes which revealed a comprehensive landscape of the relationship between circRNAs biogenesis and their parental linear transcripts. Furthermore, we investigated the distribution of circRNAs in the genome, based on the overlapping sequences of circRNAs on chromosomes (Supplementary Table S3). The patterns of circRNAs distribution were

Annotation of Host Linear Transcripts during Yak Adipocyte Differentiation
CircRNAs biogenesis occurs through back splicing and then the canonical splicing of linear RNAs, and they are commonly formed by the alternative splicing of pre-mRNA. We annotated linear transcripts of circRNAs from the correlative genes which revealed a comprehensive landscape of the relationship between circRNAs biogenesis and their parental linear transcripts. Furthermore, we investigated the distribution of circRNAs in the genome, based on the overlapping sequences of circRNAs on chromosomes (Supplementary Table S3). The patterns of circRNAs distribution were characterized in the genome ( Figure 3A). The mean transcript length of protein-coding genes is shown in Figure 3B. It revealed that the genomic loci of circRNAs were widely distributed across chromosomes ( Figure 3C). In addition, we found that one parental gene could generate multiple circRNA isoforms, in this study, 7203 circRNAs were formed from only 2909 host genes ( Figure 3D).

Expression Analysis of circRNAs during Yak Adipocyte Differentiation
Comparison of the expression of circRNAs between two stages (day 0 and 12) indicated that 136 circRNAs were differentially expressed at 12 days after adipocyte differentiation, of which 92 circRNAs were upregulated and 44 were downregulated ( Figure 4B), while at 2 days after adipocyte differentiation, 7 circRNAs were found to be differentially expressed, with 6 upregulated and 1 downregulated ( Figure 4A). This indicated that more circRNAs were upregulated at day 12, compared to day 2 after adipocyte differentiation.

Expression Analysis of circRNAs during Yak Adipocyte Differentiation
Comparison of the expression of circRNAs between two stages (day 0 and 12) indicated that 136 circRNAs were differentially expressed at 12 days after adipocyte differentiation, of which 92 circRNAs were upregulated and 44 were downregulated ( Figure 4B), while at 2 days after adipocyte differentiation, 7 circRNAs were found to be differentially expressed, with 6 upregulated and 1 downregulated ( Figure 4A). This indicated that more circRNAs were upregulated at day 12, compared to day 2 after adipocyte differentiation.

Quantitative Real-Time PCR Verification
We selected ten circRNAs and used their sequence-specific qPCR primers ( Figure 5A) to detect differential expression of circRNAs and amplified the junction regions. The circRNAs junctions were verified by Sanger sequencing of the amplified PCR products ( Figure 5B). The qPCR expression trends of circRNAs were similar to RNA sequencing data ( Figure 5C), which indicated that the sequencing data were reliable.

Quantitative Real-Time PCR Verification
We selected ten circRNAs and used their sequence-specific qPCR primers ( Figure 5A) to detect differential expression of circRNAs and amplified the junction regions. The circRNAs junctions were verified by Sanger sequencing of the amplified PCR products ( Figure 5B). The qPCR expression trends of circRNAs were similar to RNA sequencing data ( Figure 5C), which indicated that the sequencing data were reliable.

Quantitative Real-Time PCR Verification
We selected ten circRNAs and used their sequence-specific qPCR primers ( Figure 5A) to detect differential expression of circRNAs and amplified the junction regions. The circRNAs junctions were verified by Sanger sequencing of the amplified PCR products ( Figure 5B). The qPCR expression trends of circRNAs were similar to RNA sequencing data ( Figure 5C), which indicated that the sequencing data were reliable.

GO and KEGG Pathway Analysis the Host Genes of circRNAs
In order to explore the biological functions of differentially expressed circRNAs, we performed GO analysis for the host genes of circRNAs. GO analysis classified host linear transcripts in biological process (BP), cellular components (CC) and molecular function (MF) after 12 days of yak adipocyte differentiation, as shown in Figure 6A. The prediction terms (p-value ≤ 0.05) for processes were screened and ranked according to their p-value. Our findings showed that genes involved in every GO category were connected with the RNA polyadenylation, metabolic process, peptidase activity, and nucleic binding, which suggested that several circRNAs played a main role in the basic process of adipocyte differentiation. As shown in Figure 6B, the KEGG pathways were enriched in the adipocytokine signaling pathway, forkhead box O (FoxO) signaling pathway, extracellular matrix (ECM)-receptor interactions, focal adhesion, lysine degradation, inflammatory mediator regulation of transient receptor potential (TRP) channels, thyroid hormone signaling pathway, and regulation of actin cytoskeleton, these were all concerned with adipocyte differentiation and proliferation, which revealed that circRNAs might play a pivotal role in adipocyte differentiation. After 2 days of yak adipocyte differentiation, GO and KEGG results were as shown in Supplementary Figure S1A-B, respectively. The most enriched GO terms were serine-type endopeptidase inhibitor activity, endopeptidase regulator activity, as well as endopeptidase inhibitor activity. The most enriched KEGG pathways were acute myeloid leukemia, transcriptional misregulation, and pathways in cancer.
Genes 2020, 11, x FOR PEER REVIEW 9 of 18 circRNAs. (C) Change in circRNA levels between day 0, day 2, and day 12 groups. Day 2/day 0 and day 12/day 0 ratios for 10 circRNAs were based on the RNA-seq data and the expression of differentially expressed circRNAs as determined by qRT-PCR (three biological replicates, each done in triplicate).

GO and KEGG Pathway Analysis the Host Genes of circRNAs
In order to explore the biological functions of differentially expressed circRNAs, we performed GO analysis for the host genes of circRNAs. GO analysis classified host linear transcripts in biological process (BP), cellular components (CC) and molecular function (MF) after 12 days of yak adipocyte differentiation, as shown in Figure 6A. The prediction terms (p-value ≤ 0.05) for processes were screened and ranked according to their p-value. Our findings showed that genes involved in every GO category were connected with the RNA polyadenylation, metabolic process, peptidase activity, and nucleic binding, which suggested that several circRNAs played a main role in the basic process of adipocyte differentiation. As shown in Figure 6B, the KEGG pathways were enriched in the adipocytokine signaling pathway, forkhead box O (FoxO) signaling pathway, extracellular matrix (ECM)-receptor interactions, focal adhesion, lysine degradation, inflammatory mediator regulation of transient receptor potential (TRP) channels, thyroid hormone signaling pathway, and regulation of actin cytoskeleton, these were all concerned with adipocyte differentiation and proliferation, which revealed that circRNAs might play a pivotal role in adipocyte differentiation. After 2 days of yak adipocyte differentiation, GO and KEGG results were as shown in Supplementary Figure S1A-B, respectively. The most enriched GO terms were serine-type endopeptidase inhibitor activity, endopeptidase regulator activity, as well as endopeptidase inhibitor activity. The most enriched KEGG pathways were acute myeloid leukemia, transcriptional misregulation, and pathways in cancer.

Identification of CircRNA-MiRNA Axis/Pairs
The previous studies on functions of circRNAs were mainly focused on whether circRNAs could work as miRNA sponges to modulate the gene expression or not [10,23,24]. In this study, we performed miRanda to reveal that 131,884 interactions existed in the numerous miRNAs and 7203 circRNAs (Supplementary Table S4). Interestingly, we found that some well-known miRNAs were strongly associated with the differentiation of adipocytes. We considered that they would be promising candidates for subsequent research. CircRNAs (novel_circ_0009127, novel_circ_0000628,

Construction of the circRNA-miRNA-Gene Network
In order to elucidate the competing endogenous RNA network, the targets of differentially expressed circRNAs and downstream-regulated genes were predicted by miRanda, which was formed on the basis of circRNA-miRNA-gene connectivity (Supplementary Table S5). Further, to explore the bio-function of circRNAs involved in adipocyte differentiation, we used Cytoscape software (version 3.5.1, BiNGO plug-in) to construct a competing endogenous network of circRNA-miRNA-gene, which was based on the underlying effect of circRNAs (novel_circ_0006981, novel_circ_0009127, novel_circ_0000628, novel_circ_0011513, novel_circ_0010775 and novel_circ_0001494), and their putative miRNA targets and downstream-regulated mRNAs ( Figure 7A; Supplementary Figures  S2-S6). The potentially significant circRNAs contain numerous binding sites for the miRNAs associated with adipogenic differentiation. Besides, the function of putative mRNA targets in the networks were annotated by KEGG pathway, as indicated in Figure 7B and Supplementary Figures  S2-S6, which showed the KEGG pathways were enriched in phosphatidylinositol signaling system, nuclear factor kappa B (NF-kappa B) signaling pathway, Hedgehog signaling pathway, ECM-receptor interactions, steroid biosynthesis, pancreatic secretion, glycerolipid metabolism, cytokine-cytokine receptor interaction, Janus kinase (Jak)-signal transducer and activator of transcription (STAT) signaling pathway, Type II diabetes mellitus, FoxO signaling pathway and VEGF signaling pathway, these were all concerned with adipocyte differentiation and proliferation. This information provided a significant theoretical basis and reference for us to investigate the underlying mechanism of circRNA in the differentiation of yak adipocyte. Taken together, circRNAs may play an extensive endogenous regulatory role in the differentiation of yak adipocytes by interacting with multiple miRNAs. As a marked potential ceRNA, Novel circ 0006981 which controls adipocyte differentiation will be explored in further studies.
proliferation. This information provided a significant theoretical basis and reference for us to investigate the underlying mechanism of circRNA in the differentiation of yak adipocyte. Taken together, circRNAs may play an extensive endogenous regulatory role in the differentiation of yak adipocytes by interacting with multiple miRNAs. As a marked potential ceRNA, Novel circ 0006981 which controls adipocyte differentiation will be explored in further studies.

Homology Analysis of Yak circRNAs
The circ-RNAs of yak, compared with human and mouse, showed that the yak circRNAs shared 6636 (92.13%) and 3582 (49.73%) homology with human and mouse circRNAs, respectively. Most of the yak circRNAs displayed high homology with these two species, especially with a human. In total, 3536 yak circRNAs showed homology with human and mouse circRNAs (Figure 8).

Homology Analysis of Yak circRNAs
The circ-RNAs of yak, compared with human and mouse, showed that the yak circRNAs shared 6636 (92.13%) and 3582 (49.73%) homology with human and mouse circRNAs, respectively. Most of the yak circRNAs displayed high homology with these two species, especially with a human. In total, 3536 yak circRNAs showed homology with human and mouse circRNAs (Figure 8).

Discussion
The proliferation and differentiation of animal cells are mediated by a variety of regulatory factors, including growth factors, transcription factors, regulatory proteins and some hormone receptors, and so on. Previous studies on molecular mechanisms mainly focused on DNA, mRNA, and miRNA levels. Recently, the non-coding RNA of circRNA has changed its status from a rare curiosity to the central regulatory role in cell metabolism [10,24,25]. So far, transcriptome sequencing

Discussion
The proliferation and differentiation of animal cells are mediated by a variety of regulatory factors, including growth factors, transcription factors, regulatory proteins and some hormone receptors, and so on. Previous studies on molecular mechanisms mainly focused on DNA, mRNA, and miRNA levels. Recently, the non-coding RNA of circRNA has changed its status from a rare curiosity to the central regulatory role in cell metabolism [10,24,25]. So far, transcriptome sequencing as a preferred biotechnique permits us to identify the circRNAs in diverse species [26][27][28][29][30], some of the circRNAs proved to play a vital role in animal growth and development. Recently,circRNA expression profiles have been characterized in adipose tissues of swine [31], mice [32], humans [33,34], buffalo [35] and cattle [28]. Nevertheless, the circRNAs in yak remain unknown. In order to fully understand the circRNAs related to bovine adipogenesis, in this study, we used RNA-seq to analyze the non-coding regions of primary cultured yak pre-adipocytes and adipocytes and successfully identified 7203 circRNAs among numerous circRNAs.
The continued discoveries of functional circRNAs showed that circRNAs as a transcriptional product had an essential role in various tissues and cell types of animals [8,31,36]. A total of 136 and 14 differentially expressed circRNAs were identified at day 2 and 12 after yak adipocyte differentiation and pre-adipocyte differentiation at day 0, respectively. These circRNAs might have potential biological function in yak adipocytes differentiation. Adipocytes proliferation and differentiation involve numerous differentially expressed genes and non-coding RNAs [37][38][39][40]. For instance, Li et al. [11] reported that lncRNA ADNCR suppressed cattle adipocyte differentiation by inhibiting miR-204 as a competitive endogenous RNA, and Chen et al. [41] elucidated that Lnc-U90926 attenuated 3T3-L1 adipocytes differentiation by suppressing the transactivation of PPARγ2 or PPARγ. Recently, Li et al. [31] reported that 275 circRNAs were differentially expressed in subcutaneous adipose tissues of Laiwu and Large White pigs and Sun et al. [42] identified 4080 circRNAs were significantly expressed in human visceral preadipocytes and adipocyte. Additionally, Huang et al. [35], identified a total of 5141 circRNAs in buffalo, among them 252 circRNAs were differentially expressed between the adult and young buffaloes and concluded that circRNAs strongly correlated with fat deposition associated genes. These findings suggested that non-coding RNAs possibly participated in the development of animal adipose tissue and adipocyte differentiation through post-transcriptional regulation. Thus, we hypothesized that circRNAs might play a key role in post-transcriptional regulation during yak adipocyte differentiation.
Previous findings revealed that several biological pathways participate in adipocyte differentiation, including the signaling pathway of AMP-activated protein kinase (AMPK), Wnt, Hedgehog, insulin-like growth factor (IGF), adipocytokine, FoxO, mammalian target of rapamycin (mTOR), as well as ECM-receptor interactions [57][58][59][60][61][62]. However, there is limited knowledge available related to the function of circRNAs in yak adipocytes and adipose tissue. Our research clearly shows that circRNAs played a vital role in regulating the yak adipocytes differentiation by KEGG enriched pathways. Mainly, the phosphatidylinositol signaling system, NF-kappa B signaling pathway, Hedgehog signaling pathway, ECM-receptor interactions, steroid biosynthesis, pancreatic secretion, glycerolipid metabolism, cytokine-cytokine receptor interaction, Jak-STAT signaling pathway, and VEGF signaling pathway were significantly enriched.
Homologous analysis of circRNAs is very important to reveal their common characteristics. The sequence of circRNAs shared high identity among different species [63,64]. In our research, BLAST was employed to analyze the homology of circRNAs. The results showed that 3582 and 6636 yak circRNAs were homologous to mouse and human, respectively. Our findings are in agreement with Guo et al., who reported that the mouse circRNAs were homologous to human circRNAs [65], and Sun et al. revealed that human hsa_circ_0094183 and hsa_circ_0116913 were homologous with mice MMU_CIRCpedia_216382 and MMU_CIRCpedia_14213 [42]. These results implied that circRNAs might have similar functions in eukaryotes.

Conclusions
In our paper, the comprehensive landscape of circRNA expression profile was identified during yak adipocyte differentiation. For the first time, a circRNA-miRNA-gene interaction network was established for plateau animals. Furthermore, several core networks were established as a base for miRNA target site in exons of circRNA loci analysis. Meanwhile, we explored several interesting circRNAs (novel_circ_0009127, novel_circ_0000628, novel_circ_0011513, novel_circ_0010775, novel_circ_0006981 and novel_circ_0001494), which functioned in NF-kappa B signaling pathway, Hedgehog signaling pathway, ECM-receptor interactions, steroid biosynthesis, glycerolipid metabolism, cytokine-cytokine receptor interaction, Jak-STAT signaling pathway, Type II diabetes mellitus, FoxO signaling pathway and VEGF signaling pathway, these all were closely linked to adipocyte differentiation and proliferation. Undoubtedly, our findings will facilitate the research of the regulatory mechanism for circRNA in bovine adipogenesis, which may offer further insights into the adipocyte differentiation process and valuable resources for understanding the circRNA biology in plateau animals.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/4/414/s1, Figure S1: Annotations and enrichment of differentially expressed circRNAs between day 2 and day 0, Figure S2: CircRNA-miRNA-gene network with potentially effective novel_circ_0009127 (A) and KEGG pathway analysis (B), Figure S3: CircRNA-miRNA-gene network with potentially effective novel_circ_0000628 (A) and KEGG pathway analysis (B), Figure S4: CircRNA-miRNA-gene network with potentially effective novel_circ_0011513 (A) and KEGG pathway analysis (B), Figure S5: CircRNA-miRNA-gene network with potentially effective novel_circ_0010775 (A) and KEGG pathway analysis (B), Figure S6: CircRNA-miRNA-gene network with potentially effective novel_circ_0001494 (A) and KEGG pathway analysis (B), Table S1: Sequence specific primers were used in this study. F: forward, R: reverse. β-actin was used as endogenous control genes for miRNA and mRNA., Table S2: CircRNA-seq quality control statistics., Table S3: Distribution of circRNAs in the genome based on the overlapping sequences of circRNAs on chromosomes, Table S4: Interaction relationships between circRNAs and various miRNAs were found, Table S5: The circRNA-miRNA-gene interaction network between circRNA and its target gene.

Conflicts of Interest:
The authors have no conflict of interest.