Characterization of Circular RNAs in Chinese Buffalo (Bubalus bubalis) Adipose Tissue: A Focus on Circular RNAs Involved in Fat Deposition.

Circular RNAs (circRNAs) have been identified as a novel type of regulators involved in multiple biological processes. However, circRNAs with a potential function in fat deposition in buffalo are poorly understood. In this study, six RNA libraries of adipose tissue were constructed for three young and three adult Chinese buffaloes with paired-ends RNA sequencing using the Illumina HiSeq 3000 platform. A total of 5141 circRNAs were computationally identified. Among them, 252 circRNAs were differentially expressed (DE) between the young and adult buffaloes. Of these, 54 were upregulated and 198 were downregulated in the adult group. Eight DE circRNAs were further identified by quantitative real-time-PCR (qRT-PCR) and Sanger sequencing. Co-expression analysis revealed that 34 circRNAs demonstrated a strong correlation with fat deposition-associated genes (|r| > 0.980). Among these, expressional correlation between two circRNAs (19:45387150|45389986 and 21:6969877|69753491) and PR/SET domain 16 was further verified using qRT-PCR, and a strong correlation was revealed (1 > |r| > 0.8). These results strongly suggest that circRNAs 19:45387150|45389986 and 21:6969877|69753491 are potential regulators of buffalo fat deposition. In summary, this study characterized the circRNA profiles of adipose tissues at different stages for the first time and revealed two circRNAs strongly correlated with fat deposition-associated genes, which provided new candidate regulators for fat deposition in buffalo.

1. Introduction

Circular RNAs (circRNAs) are well known as a class of non-coding RNAs with a closed loop. More than 40 years ago, viral circRNAs were identified [1]. Subsequently, circRNAs were observed and identified in eukaryotes [2,3,4]. Due to the rapid advancement of RNA sequencing methodologies and bioinformatics, a great number of circRNAs have been discovered in many species in recent years. Originally, circRNAs were considered byproducts generated from an error in alternative splicing [5]. In recent years, circRNAs have been found to be involved in multiple biological processes [6,7,8]. In humans, many tissue-specific circRNAs have been identified [9,10,11], and some of them play roles in disease [12] or can be cancer biomarkers [13,14]. In mice, circRNA_010567 promotes myocardial fibrosis via miR-141 suppression by targeting TGF-β1 [15]. Similarly, circRNA_100782 regulates proliferation of pancreatic carcinoma [16]. Interestingly, expression of circRNAs can be induced by chemicals [17], implicating circRNAs in the response to external injury. In livestock animals, a great number of circRNAs have been identified by RNA sequencing technology, and some of them are associated with economically important traits. Based on RNA sequencing, 13,950 circRNAs have been identified in the preovulatory ovarian follicles of goats [18]. Among these, chi_circ_0008219 is a potential circRNA associated with goat reproductive traits [18]. In swine, 6,130 circRNAs have been identified in the muscle, and most of them are uniquely expressed in one breed [19], suggesting that circRNAs may be involved in maintaining characteristics of a breed. In cattle, 12,918 circRNAs have been identified in the muscle tissue [20]. The circLMO7 is abundant in muscle tissue and regulates myoblast differentiation and survival [20]. Characterizing the circRNA profiles of specific tissues or cells is a promising way to reveal functional circRNAs.

Beef (cattle meat) is popular all over the world for its high protein and vitamin content. Meat tenderness and juiciness are affected by the intramuscular fat (IMF) content. In recent years, cattle meat has been imported to meet consumer demand in China [21]. Although buffalo (Bubalus bubalis) is abundant in China, buffalo meat has not been widely accepted by customers because of the low IMF content [22]. Besides, the IMF level is too limited to be sampled for RNA sequencing analysis. Many studies demonstrate that back subcutaneous fat is significantly associated with IMF [23,24], which indicates that the character of backfat tissue may be used as an indicator of IMF. However, development of backfat tissue is poorly understood and requires further investigation in buffalo. Recently, circRNA expression profiles have been characterized in adipose tissues or adipose-associated cells of humans [25], mice [26], and swine [27]. These results provide vital information for further studies on fat deposition. However, to date, there has been no study characterizing the circRNA profiles in buffalo adipose tissue. Adipose tissue is formed at specific times to meet the requirement of an organism’s development. Generally, capacity of fat deposition in adult animals is higher than that of young animals. In the present study, circRNA profiles of subcutaneous adipose tissue of Chinese buffalo at two developmental stages (i.e., young and adult) were identified using paired-end sequencing using the Illumina HiSeq 3000 platform. We aimed to identify circRNAs with a significant association with fat deposition for further research on adipose tissue development in buffalo.

2. Materials and Methods

2.1. Animal Ethics

All animals were bred for commercial use rather than for experimental reasons, and they were slaughtered according to the food industry-approved halal food quality certified protocol by a Muslim cleric according to the law of Islam. Thus, no ethics approval was required by a specific committee.

2.2. Animals and Tissue Samples

For RNA sequencing, six Chinese buffaloes (swamp type) were sampled, including three young (6-month-old) and three adult (30 month old) individuals. In addition, another 49 buffaloes with variable months of age were also sampled for a subsequent validation assay. The buffaloes were raised at the Xinyang buffalo farm (Xinyang, Henan, China) with equivalent forage and feeding management conditions. Subcutaneous adipose was harvested after slaughter and frozen immediately in liquid nitrogen.

2.3. RNA Isolation and Sequencing

Total RNA was extracted using TRIzol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. The RNA quantity was verified with the NanoDrop 2000 (Nanodrop, Wilmington, DE, USA) and 1.5% agarose gel. The rRNA was removed by the Epicentre Ribo-zero rRNA Removal Kit (Epicentre, San Diego, CA, USA). The RNA was fragmented, and random primers were added to yield double-stranded cDNA. Then, the A-tail and ligating adapters were added to the cDNA, and PCR was performed to enrich the cDNA. The cDNA was precipitated with ethanol, and its quality was assessed using an Agilent Bioanalyzer 2100 system (Agilent, Santa Clara, CA, USA). Quantification of the sample was performed using the Quant-iT™ PicoGreen® dsDNA Assay Kit (Life, Southfield, MI, USA). Finally, the cDNA library was sequenced using the double terminal sequencing mode of the Illumina HiSeq 3000 platform. The RNA sequencing data were deposited in the NCBI Gene Expression Omnibus (GEO), and the accession number is GSE112744.

2.4. Quality Control, Transcriptome Assembly, and CircRNA Prediction

Trim Galore was used to remove adapter sequences and low-quality sequences [28]. The FastQC function in Trim Galore was used for quality control and calculating the proportion of Q20 and Q30 [28]. Reads with high quality were collected and used for analysis.

The cattle genome (UMD3.1) was used as a reference genome because the buffalo genome has not been successfully assembled yet. The cattle reference genome and the gene model annotation files were downloaded from the Ensemble database. High-quality reads for each sample were aligned to the reference genome by STAR [29]. For aligned sequence 6% mismatch was allowed; the maximum total length mapped to genome was 1,000,000 bp; minimum intron length was 20 bp; maximum intron length was 1,000,000 bp; other parameters were set to default values. For circRNA prediction, the CIRCexplorer2 was employed to identify the back-spliced junction of circRNAs [30]. The host genes of circRNAs, excluding the intergenic circRNAs and the circRNAs produced from transcripts without gene symbols, were statistically analyzed.

2.5. Differential Expression, Functional Enrichment, and Co-Expression Analysis

To calculate the expression levels of circRNAs, back-spliced reads per million mapped reads (BSRP) was used. To screen for differentially expressed (DE) circRNAs, DESeq2 software (Version 1.14.1) was used [31], and those with a fold-change ≥ 2 and a p-value ≤ 0.05 were considered as DE circRNAs. The host genes of DE circRNAs were used for functional enrichment by DAVID [32]. A p-value ≤ 0.05 was used as threshold to evaluate significant enrichment. Only items associated with lipid metabolism were retained.

A number of DE genes were previously identified in adipose tissues of 6 and 30 month old buffaloes. Among them, 11 genes (peroxisome proliferator activated receptor gamma (PPARG), serum/glucocorticoid regulated kinase 1 (SGK1), cyclin dependent kinase 5 (CDK5), sterol regulatory element binding transcription factor 1 (SREBF1), thyroid hormone responsive (THRSP), FYN proto-oncogene (FYN), SMAD family member 3 (SMAD3), sterol regulatory element binding transcription factor 2 (SREBF2), secreted frizzled related protein 4 (SFRP4), PR/SET domain 16 (PRDM16), and bone morphogenetic protein 7 (BMP7)) involved in adipogenesis [33,34] were included (Table S1). To search for potential circRNAs involved in fat deposition, co-expression analysis was performed for these 11 DE genes and the DE circRNAs. The 1 > |r| > 0.8 means a strong correlation. To obtain a reliable result, only those mRNA–circRNA pairs with |r| > 0.980 were retained for the small sample size (n = 6) used in RNA sequencing.

2.6. Validation Assay

Quantitative real-time-PCR (qRT-PCR) and Sanger sequencing were performed for circRNA validation. Primers were designed by the Pick Primers function in NCBI (Table 1). Reverse transcription was performed using the PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). The qRT-PCR was performed to detect the expression levels of circRNAs and mRNAs using SYBR Green I (TaKaRa, Dalian, China) with two-step reactions, as recommended by the manufacturer’s protocol. The GAPDH was used as an internal control gene [35]. The cycle threshold (2−ΔΔCt) method was used to calculate relative expression levels of circRNA and mRNA. Relative expression levels of each tissue are presented as the means ± SEMs. Meanwhile, the PCR products were sequenced to validate the junctions of circRNAs.

Table 1

Primers used for the validation of differentially expressed circular RNAs (circRNAs) and circRNA–mRNA pairs.

CircRNA 1/mRNA	Forward primer (5′–3′)	Reverse primer (5′–3′)	Product Length (bp)
X:79761001|79786582	GCTAAGAACAGAAACACAAAATGCT	CCTCATCTCCAGGGTTTTCTTCA	225
16:34132574|34371163	TCGAACACTCTCTTCAGATGCAA	TGGGGTTTGGATTCTCTGCT	213
22:36903411|36917263	CAGACGGCTCTGGTGCAATA	GTCATCAGCCAGGCTACTCC	180
25:10461217|10480061	AAGGAGTGGGACCTAAAGCC	GGCATTGAAGCGGAAGAAGTC	183
10:44952496|44976492	AACACAGACCTACCGCATCC	GAAAGAGGGCGTAGGTGTCA	234
1:106999005|107007859	CCAGATGTGTGGAGATTGCC	TCCTGGTGGCTGGAAATACC	227
21:69696877|69753491	ATGCAAGTGGAGGTGAAGTGT	TAACTTCTTGCCGCCCATCT	117
19:45387150|45389986	TGGAAGAGGCTAGCAAACGA	GGGCTGAATCTGTCTCGCTG	152
PRDM16	TACAGGGTGGAGAAGCGGAA	GTACCTGCACGTGTATCGCT	264
GAPDH	CACTCACTCTTCTACCTT	GCCAAATTCATTGTCGTA	91

Indicates chromosomal locations in cattle.

3. Results

3.1. Overview of Sequencing and Quality Control

A total of six cDNA libraries of adipose tissues were constructed for RNA sequencing, including three young and three adult buffaloes. There were 561 million and 433 million clean reads obtained from young and adult buffaloes, respectively (Table S2). Since the buffalo genome assembly is not available, the cattle genome (UMD3.1) was used as the reference genome. An average of 80.96% of the clean reads were mapped to the reference genome (Table S2). Approximately 0.01% of the clean reads were identified as the junction reads for each sample (Figure 1A).

Figure 1

Profiling of circRNAs in young and adult Chinese buffalo adipose tissue. (A) Clean reads and junction reads obtained from young and adult buffalo adipose tissue. (B) circRNAs expressed in young and adult buffalo adipose tissue. (C) Overview of the types of total predicted circRNAs. (D) Splicing modes at the junctions. (E) Distribution of host genes producing different numbers of circRNAs. Different colors mean different numbers of circRNAs produced by the same host gene.

3.2. CircRNAs Expressed in Back Subcutaneous Adipose Tissue of Buffalo

In total, 5141 circRNAs were computationally identified in the present study (Table S3). Among them, 501 circRNAs have been previously published for cattle [20], and the remainder were considered as novel bovine circRNAs. There were more circRNAs specifically expressed in the young group (1059) than that in the adult group (390), with 2,392 circRNAs being expressed in both developmental stages (Figure 1B). Based on the junction positions, the identified circRNAs were divided into four types: exon circRNA, exon–intron circRNA, intergenic circRNA, and intron circRNA. Most of the circRNAs identified were exon type (4,665, 90.74%), followed by exon–intron type (425, 8.27%), intergenic type (45, 0.88%), and intron type (6, 0.12%) (Figure 1C). Most of the junctions were flanked by the classical GT-AG splicing sites (98.40%, Figure 1D).

The circRNAs were produced from 2,187 host genes. A host gene could produce multiple circRNAs (Table S3, Figure 1E). About half of the host genes (1168, 53.40%) produced only a circRNA. Approximately 7% of host genes produced more than six circRNAs. The maximum number of circRNAs produced from a host gene was 32.

3.3. Identification of DE circRNAs and Functional Enrichment

Of the 5141 circRNAs, a total of 252 DE (p < 0.01) circRNAs were identified between the young and adult groups, including 54 upregulated (21.43%) and 198 downregulated (78.57%) circRNAs in adult buffaloes (Table S4). Among them, 20 circRNAs were downregulated with an absolute fold-change ≥ 4, and the maximum absolute fold-change was 6.18 (Table S4). However, all the upregulated circRNAs demonstrated a fold-change less than 4. Most of the DE circRNAs were exon type, and there was no intron circRNA, indicating that circRNAs produced from exons might play a more important role in the development of buffalo adipose tissue. Host genes of DE circRNAs were used for functional enrichment. The primary aim of this study was to identify circRNAs implicated in fat deposition. However, no term associated with lipid metabolism was found (data not shown).

3.4. Co-Expression Analysis Revealed circRNAs Associated with Fat Deposition

While the functional enrichment of the host genes did not detect any associations between specific circRNAs and fat deposition, the co-expression analysis revealed a total of 37 circRNA–mRNA pairs had an |r| > 0.980 (Table 2), including 10 pairs with a negative correlation and 27 pairs with a positive correlation. A total of 34 circRNAs and six mRNAs (CDK5, PPARG, PRDM16, SFRP4, SREBF2, and THRSP) were involved.

Table 2

The circRNA–mRNA pairs with potential effects on fat deposition.

Number	circRNA_ID 1	Symbol	r Value	Number	circRNA_ID 1	Symbol	r-Value
1	10:11563837|11593862	SREBF2	−0.9975	20	12:52715274|52733582	PRDM16	0.9834
2	15:83640019|83640860	THRSP	−0.9960	21	15:22884720|22887938	SREBF2	0.9841
3	12:89500403|89506932	PPARG	−0.9899	22	12:89492189|89525269	PRDM16	0.9841
4	2:37078234|37133108	THRSP	−0.9889	23	22:30718181|30723897	PRDM16	0.9843
5	12:23235488|23237718	THRSP	−0.9889	24	21:69696877|69753491	PRDM16	0.9852
6	12:2298522|2305539	PPARG	−0.9859	25	4:9762677|9766001	PRDM16	0.9854
7	2:37078234|37133108	CDK5	−0.9855	26	21:42001223|42024196	PRDM16	0.9855
8	12:23235488|23237718	CDK5	−0.9855	27	10:60619225|60625573	SFRP4	0.9861
9	15:83640019|83640860	CDK5	−0.9827	28	11:22513792|22530506	PRDM16	0.9865
10	9:84303407|84311225	THRSP	−0.9807	29	16:30873113|30903109	PRDM16	0.9867
11	X:79782632|79786582	PRDM16	0.9805	30	7:18298331|18317388	SFRP4	0.9903
12	17:40606888|40608971	PRDM16	0.9805	31	14:22698555|22703735	SFRP4	0.9906
13	14:28794206|28825387	PRDM16	0.9810	32	14:45838481|45839369	SFRP4	0.9928
14	11:100361051|100374772	PRDM16	0.9814	33	10:70853226|70859881	PRDM16	0.9933
15	X:79761001|79786582	PRDM16	0.9818	34	2:72283596|72308351	THRSP	0.9944
16	17:55314422|55316845	SFRP4	0.9828	35	2:37061458|37078299	PRDM16	0.9944
17	20:58640529|58642071	SFRP4	0.9828	36	5:10860435|10879727	THRSP	0.9948
18	9:15487323|15505719	PRDM16	0.9833	37	2:72283596|72308351	CDK5	0.9955
19	8:91848002|91864189	PRDM16	0.9834

Indicates chromosomal locations in cattle.

3.5. Validation of DE circRNAs and circRNA–mRNA Pairs

To evaluate the reliability of the DE analysis results, eight DE circRNAs (Table 1) were selected randomly for validation by qRT-PCR. Fragments containing the junctions of circRNAs were clearly amplified (Figure 2A). The fragment sequences obtained using Sanger sequencing further confirmed the back-splice junctions (Figure 2B). The expression patterns of the eight circRNAs were consistent with those obtained from RNA sequencing (Figure 2C,D).

Figure 2

Validation of eight DE circRNA expression levels. (A) Electrophoretogram of the qRT-PCR products of the eight circRNAs obtained with divergent primers. M: marker 500. (B) Sanger sequencing of the eight circRNAs confirms the head-to-tail splicing. (C) Expression pattern of the eight circRNAs between young and adult Chinese buffalo adipose tissue based on RNA sequencing. (D) Validation of the expression pattern for the eight DE circRNAs by qRT-PCR. All of the eight circRNAs confirmed the RNA sequencing pattern. Data are presented as the means ± SEMs, n = 3.

In addition, the expressional correlation of two circRNAs (19:45387150|45389986 and 21:6969877|69753491) paired with the PRDM16 transcript were confirmed by qRT-PCR as well. First, a validation assay was performed using adipose tissues of three young and three adult buffaloes. In this part, r values for 19:45387150|45389986-PRDM16 and 21:6969877|69753491-PRDM16 were 0.963 and 0.893, respectively (Figure 3A,B, Table 3). Then, adipose tissue samples from another 49 buffaloes with variable months of age were further used for validation, and the r values for 19:45387150|45389986-PRDM16 and 21:6969877|69753491-PRDM16 were 0.889 and 0.893, respectively (Figure 3C,D, Table 3).

Figure 3

Validation of the two circRNAs with a strong correlation to PRDM16 by qRT-PCR. The y-axis shows the relative expression of circRNA and PRDM16. Pearson’s correlation coefficients between circRNA and PRDM16 based on their relative expression were calculated to yield an r-value. (A) Validation for 19:45387150|45389986-PRDM16 in young and adult buffalo, n = 6. (B) Validation for 21:6969877|69753491-PRDM16 in young and adult buffalo, n = 6. (C) Validation for 19:45387150|45389986-PRDM16 in buffalo with variable months of age, n = 49. (D) Validation for 21:6969877|69753491-PRDM16 in buffalo with variable months of age, n = 49.

Table 3

Two circRNA–mRNA pairs with potential effects on fat deposition.

CircRNA 1	CircRNA Type	Host Gene	mRNA	r-Value by RNA-Seq	r-Value by qRT-PCR
				Young and Adult(n = 6)	Young and Adult(n = 6)	Variable Months of Age(n = 49)
19:45387150|45389986	Exon	NMT1	PRDM16	0.968	0.963	0.889
21:69696877|69753491	Exon	MARK3	PRDM16	0.985	0.893	0.893

Indicates chromosomal locations in cattle.

4. Discussion

To identify circRNAs with a potential function in fat deposition in buffalo, circRNA expression profiles of adipose tissues were characterized by Ribo-Zero ribonucleic acid sequencing in Chinese buffalo. The RNA library construction by removing rRNAs and linear RNAs is more effective for circRNA identification [36]. In the present study, the RNA library was constructed and 5141 circRNAs were identified which was comparable with previous studies constructing libraries by removing rRNAs and linear RNAs [18,20]. Most of the circRNAs identified were exon type (Figure 1C), and ~98% of junctions reflected the classical GT-AG splicing mode (Figure 1D). In addition, 501 of the 5141 circRNAs were previously identified in the muscle tissue of Qinchuan cattle [20]. Our results further demonstrate that constructing an rRNA-depleted library without removing linear RNAs is practical for circRNA identification. Differential expression analysis revealed 252 DE circRNAs between the young and adult groups (Table S4). The validation results indicate a high reliability of the circRNA prediction.

For genes, a specific expression profile indicates a specific role and function. The circRNAs were expressed more abundantly in young buffalo than in adult buffalo (Figure 1B). Meanwhile, most of the DE circRNAs (198/252) were downregulated in adult buffalo (Table S4). According to the mammalian development pattern, a young animal is less efficient at fat accumulation than an old one. Thus, the richer circRNA set in young buffalo indicates that circRNAs may mainly act as negative regulators of fat deposition in Chinese buffalo.

The primary aim of this study was to reveal circRNAs associated with fat deposition. Previous studies have suggested that the host genes of circRNAs can reveal important information on the function of circRNAs [30]. In this study, functional enrichment analysis for the host genes was performed. However, no term associated with fat deposition was identified (p < 0.05). In fact, circRNAs participate in biological processes with multiple modes of action [20,37,38], not limited to regulation of their host genes. Therefore, functional enrichment based on the host genes may not be enough. Besides, DAVID is a tool that systematically maps a large number of interesting genes to the associated biological annotation, and statistically highlights the most enriched biological annotation out of thousands of linked terms and contents [39]. However, some genes without strong neighbors will be omitted from the analysis. The limited annotation information of bovines may affect the result as well. Instead, co-expression analysis was used and 34 circRNAs with putative association with fat deposition were screened. Though without systematic and integrative analysis as is available in DAVID, co-expression analysis can directly capture a pair of transcripts with high correlation and should make the study of regulatory mechanisms of circRNAs more difficult.

Initially, 10 of the 37 circRNA–mRNA pairs were selected for validation by qRT-PCR. However, specific primers were available for only three pairs, and high correlations were validated for only two pairs, namely, 19:45387150|45389986-PRDM16 and 21:6969877|69753491-PRDM16. Previous studies have indicated that circRNAs can regulate the expression of their host genes [7,37,40]. The circRNAs 19:45387150|45389986 and 21:6969877|69753491 are produced by N-myristoyltransferase 1 (NMT1) and microtubule affinity regulating kinase 3 (MARK3), respectively (Table 3). N-myristoyltransferase is a key cellular enzyme that carries out lipid modification by facilitating the attachment of myristate (a 14-carbon saturated fatty acid) to the N-terminal glycine of several protein molecules and increases protein lipophilicity [41]. MARK3 is a vital regulator involved in energy metabolism. The Mark3 mice are protected against high-fat diet-induced obesity and display attenuated weight gain [42]. These lines of evidence further suggest that circRNAs 19:45387150|45389986 and 21:6969877|69753491 are likely to be involved in lipid metabolism.

In addition, PRDM16 is identified as a highly expressed gene in brown fat tissue and drives a brown fat-selective transcriptional program [43,44]. Brown adipose tissue can dissipate chemical energy in the form of heat (non-shivering thermogenesis) when encountering cold exposure and excessive feeding, which is considered to exert negative effects on feeding efficiency and inhibit fat deposition in animals [45,46]. In our study, PRDM16 is richer in young buffaloes than in adult individuals (Table S1, Figure 4). These results indicate that brown adipose tissue is more active in young individuals, which is consistent with the rule of fat development in buffaloes. Both 9:45387150|45389986 and 21:6969877|69753491 demonstrated a positive correlation with PRDM16, which indicates that they may inhibit fat deposition in buffalo. However, further functional identification is needed.

Figure 4

Expression of PRDM16 in young and adult buffalo adipose tissue. (A) Expression pattern analyzed by RNA sequencing. (B) Expression pattern analyzed by qRT-PCR. Data are presented as the means ± SEMs, n = 3. ** p < 0.01, compared with the young group.

5. Conclusions

This study characterized the circRNA expression profiles in adipose tissues of young and adult Chinese buffaloes for the first time. We primarily focused on circRNAs with a potential function in fat deposition. Two circRNAs showing a strong correlation with PRDM16 are likely to inhibit the deposition of white adipose tissue in buffalo. These results provide new candidate factors for further studies of the molecular regulatory mechanism of fat deposition in buffalo.