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

Simple Summary Buffalo play a vital role in several southeastern and middle-eastern Asian countries and Africa. Fat deposition has received increased attention because of its importance to meat quality. However, information on the development of buffalo adipose tissue is poorly understood and requires further investigation. This study characterized circular RNA (circRNA) profiles of adipose tissues in different developmental stages in buffalo. Co-expression analysis and functional enrichment revealed a considerable number of circRNAs which may function in fat deposition. Further qRT-PCR analysis identified two circRNAs with a significant association with PR/SET domain 16, a fat metabolism gene. Our results identified candidate regulators of fat deposition in buffalo, which may be valuable for buffalo breeding. Abstract 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.


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.

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.

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 Animals 2019, 9, 403 3 of 12 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.

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.

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.

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.

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.

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).

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. Animals 2019, 9, x 5 of 13

Identification of DE circRNAs and Functional Enrichment
Of the 5,141 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 foldchange ≥ 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).

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)

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).

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.

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).

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 5,141 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 5,141 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

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 Animals 2019, 9, 403 9 of 12 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.
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.

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.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-2615/9/7/403/s1, Table S1: Eleven genes (mRNAs) associated with fat deposition showing DE between young and adult buffalo adipose tissues, Table S2: Details of quality control and mapping of RNA sequencing, Table S3: circRNAs identified in adipose tissue of Chinese buffalo based on RNA sequencing, Table S4: DE circRNAs between young and adult buffalo adipose tissue.