Comparative transcriptome analysis reveals that PCK1 is a potential gene affecting IMF deposition in buffalo

Background In China, although buffaloes are abundant, beef is mainly obtained from cattle, and this preference is mainly attributed to the low intramuscular fat (IMF) content of buffalo. Genetic factors are an important driver that affects IMF deposition. Results To reveal the intrinsic factors responsible for the low IMF content of buffalo, mRNA expression patterns in muscle and adipose tissue between buffalo and cattle were characterized by RNA sequencing analysis. The IMF content in Nanyang cattle was higher than that in Xinyang buffalo. A total of 1566 mRNAs expressed in adipose tissue showed differential expression between the longissimus dorsi muscles of buffalo and cattle. Functional annotation suggested a difference in the glycolysis/gluconeogenesis pathway between the two species. The results of RT-qPCR analysis and gain-of-function experiments confirmed the positive association between the IMF content and phosphoenolpyruvate carboxykinase 1 (PCK1) expression in buffalo. In both mouse C2C12 cells and cultured bovine myocytes, the activity of the PCK1 promoter in buffalo is lower than that in cattle. However, in mouse 3T3-L1 adipocytes and cultured bovine adipocytes, the activity of PCK1 in buffalo promoter is higher than that in cattle. Conclusions These results indicate the important role of PCK1 in buffalo IMF deposition and illustrate the differences between buffalo and cattle promoter activity that drive PCK1 expression. This research helps to establish a foundation for further studies investigating IMF deposition in buffalo.

next to cattle. Thus, buffalo hold promise for meat production. However, due to the low intramuscular fat (IMF) content, the quality of buffalo meat in terms of flavour and juiciness is poorer than that of cattle meat.
IMF is highly valued because it improves meat quality by enhancing taste, juiciness, and tenderness [7]. IMF deposition has attracted widespread attention for many years. Genetic, managerial, and nutritional factors can affect IMF deposition [8], and genetic factors present moderate or more pronounced effects [9,10]. For example, in Japanese Black cattle, notably high amounts of IMF (approximately 23%) are deposited in the longissimus dorsi muscle (LM), while in European breeds, the IMF content is less than 5% [11]. Comparisons of the gene expression pattern in LM between cattle breeds with high or low IMF content has revealed that many lipid metabolism-related genes are associated with IMF deposition, including adiponectin, C1Q and collagen domain containing (ADIPOQ) [12], peroxisome proliferator activated receptor gamma (PPARγ or PPARG) [13], thyroid hormone responsive [7], and fatty acid binding protein 4 (FABP4) [7]. Based on the high-density genotypes, a genome-wide association study suggested that PPARG coactivator 1 alpha, hepatocyte nuclear factor 4 gamma, and forkhead box P3 participated in IMF deposition in cattle [14]. Noncoding RNAs have also been reported to be involved in IMF deposition [15,16]. Although many studies have been performed on this subject, the molecular mechanism governing IMF deposition has not been elucidated.
In buffaloes, the low IMF content is mainly due to the long-term selection for draught power. Uncovering the genetic differences in the IMF content between buffalo and cattle may provide guiding significance for the improvement of buffalo meat quality. In this study, RNA sequencing was performed to compare mRNA expression patterns in muscle and adipose tissues between buffalo and cattle. A series of analyses revealed that the expression pattern of phosphoenolpyruvate carboxykinase 1 (PCK1), a gene with a significant association with IMF, differs between buffalo and cattle in adipose and muscle cells. Overexpression of PCK1 promotes lipid accumulation in buffalo intramuscular adipocytes. The results of this study suggest a positive role of PCK1 in buffalo IMF deposition, which can be used as a molecular marker for the improvement of buffalo meat quality through breeding.

IMF content of Xinyang buffalo and Nanyang cattle
To evaluate the effects of fattening in Xinyang buffalo and Nanyang cattle, the back fat thickness, LM area (LMA), and IMF content were measured at slaughter. As shown in Table 1, the LMA of cattle was 70.71 ± 11.19 cm 2 , which was larger than that of buffalo (59.61 ± 8.34 cm 2 ); the back fat thickness of cattle (4.03 ± 0.81 cm) was higher than that of buffalo (1.03 ± 0.24 cm). Correspondingly, the IMF content in cattle (2.43%) was higher than that in buffalo (0.51%), which could also be observed after slaughter (Additional file 1: Fig. S1).

Overview of mRNA expression profiles
To compare the gene expression pattern of IMF between buffalo and cattle, RNA sequencing was performed. The IMF in cattle and buffalo was too limited to be sampled for RNA sequencing analysis. Thus, both muscle and adipose tissues were used for analysis in this study. A total of 12 cDNA libraries were constructed, and 12 RNA sequencing datasets were obtained. After a highly stringent filtering pipeline and analyses, 18,713 and 18, 535 mRNAs were identified in buffaloes and cattle, respectively (Additional file 2: Table S1 and Additional file 3: Table S2). In buffaloes, 369 mRNAs were specifically expressed in muscle, and 17,906 coexpressed mRNAs were identified in both muscle and adipose tissues (Fig. 1a). In cattle, 97 mRNAs were specifically identified in muscle, and 18,220 coexpressed mRNAs were identified in both tissues (Fig. 1b).

Differential expression analysis
The analysis strategy used to identify intrinsic factors related to the low IMF content in buffaloes is presented in Fig. 2. mRNAs expressed in IMF should be involved in the coexpressed mRNAs of muscle and adipose tissues. In addition, IMF deposition can be affected by the muscle [18,19]. Therefore, a differential expression analysis was performed with the coexpressed mRNAs and muscle-specific mRNAs between buffalo and cattle (Fig.  2). In total, 1566 differentially expressed (DE) coexpressed mRNAs were ultimately identified, with 1143 upregulated mRNAs and 423 downregulated mRNAs in the muscle of cattle ( Fig. 3a and Additional file 4: Table  S3). Among these, adipose markers such as PCK1, fatty acid binding protein 3 (FABP3), FABP4, ADIPOQ, and perilipin 1 (PLIN1) were significantly upregulated in the  (Fig. 1d). In adipose tissue, however, significant differences were not found between buffalo and cattle (Fig. 1c). These results also indicated that the IMF level in cattle was higher than that in buffaloes, which was consistent with the slaughter traits presented in  coexpressed mRNAs and 70 DE muscle-specific mRNAs could distinguish buffalo muscle from cattle muscle (Fig.  3), which indicated the high quality of RNA sequencing.

RT-qPCR validation of differentially expressed mRNAs
To evaluate the quality of the results in the differential expression analysis, 6 DE coexpressed mRNAs and 2 DE muscle-specific mRNAs were randomly selected for validation via RT-qPCR. As shown in Fig. 4, the expression patterns of DE mRNAs in RT-qPCR were approximately similar to those of RNA sequences, although no significant difference was found for MYH6 and SLC16A6.
These results indicate the relatively high quality of the differential expression analysis; thus, the DE mRNAs can be used for the following analysis.

Functional enrichment
To characterize the putative functions of DE mRNAs, functional enrichment was performed for the 1566 DE coexpressed mRNAs and the 70 DE muscle-specific mRNAs. The results showed that the muscle-specific DE mRNAs were enriched in five Gene Ontology (GO) terms, namely, muscle filament sliding, muscle contraction, skeletal muscle contraction, regulation of striated muscle contraction, and sarcomere (Additional file 6: Table S5). Clearly, the terms were all associated with the development of muscle tissue. The aim of our functional enrichment was to find items associated with fat deposition. Thus, we did not investigate these results further. The coexpressed DE mRNAs were enriched in 41 GO terms and seven Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Additional file 7: Table S6). However, only glycolysis/gluconeogenesis (hsa00010), a pathway with a significant effect on fat deposition, was identified (Additional file 7: Table S6). Seventeen DE genes were included in the pathway. Among these genes, PCK1 has been demonstrated to be associated with fat deposition [20]. Overexpression of PCK1 in skeletal  UXT was used to normalize the expression level of genes, and their prime sequences were conserved between buffalo and cattle. Expression level data are presented as the mean ± S.E. (n = 6); * indicates p < 0.05; ** indicates p < 0.01 muscle can improve IMF levels in mice [17]. These results indicate that PCK1 may be a candidate gene that regulates IMF deposition in buffaloes.

Expression pattern of PCK1
According to the RNA sequencing analysis, PCK1 was upregulated in adipose tissue compared to muscle tissue in both buffalo and cattle (Fig. 5a). In addition, the expression of PCK1 in cattle muscle was higher than that in buffalo muscle (Fig. 5a). The mRNA expression level of PCK1 was then validated by RT-qPCR. The sequences of PCK1 primers were totally conserved between buffalo and cattle. Consistent with the results of RNA sequencing, PCK1 in adipose tissue was richer than that in muscle tissue in both buffalo and cattle (Fig. 5b, p < 0.01). In muscle, the expression of PCK1 was higher in cattle than that in buffalo (Fig. 5b, p < 0.01). In addition, compared to unfattened buffalo, the expression of PCK1 in muscle was upregulated in fattened buffalo (Fig. 5c, p < 0.05). These findings indicated a positive correlation between PCK1 and IMF content in buffalo.

Effects of PCK1 on the adipogenic differentiation of buffalo intramuscular preadipocytes
To detect the effect of PCK1 on the adipogenic differentiation of buffalo intramuscular adipocytes, the CDS region of buffalo PCK1 was inserted into the pHBAd vector and packaged into adenovirus for overexpression (Ad_PCK1), with EGFP as an internal indicator. Ad_EGFP was used as a control. EGFP was highly expressed on day 0 (24 h after overexpression) and day 6 of adipogenesis induction (Fig. 6a), which indicated a high level of overexpression. Indeed, the expression level of PCK1 in Ad_PCK1 cells was significantly higher than that in Ad_EGFP cells (Fig.  6b, p < 0.01). The adipogenic marker PPARG was upregulated with the overexpression of PCK1 on day 0 of adipogenesis induction (Fig. 6c, p < 0.05). On day 6, the fatty acid uptake gene (FABP4) was upregulated in Ad_PCK1 group (Fig. 6d, p < 0.01). No significant difference was detected for CCAAT enhancer binding protein alpha (C/ EBPα) (Fig. 6c), fatty acid translocase (FAT/CD36) (Fig.  6d), lipogenesis genes (glycerol-3-phosphate acyltransferase 4 (GPAT4) and diacylglycerol O-acyltransferase 1 (DGAT1)) ( Fig. 6g), and lipolysis genes (hormone-sensitive lipase (HSL), lipoprotein lipase (LPL), and peroxisome proliferator activated receptor alpha (PPARα)) ( Fig. 6h) between Ad_EGFP and Ad_PCK1 groups. Compared to the Ad_EGFP group, lipid accumulation was significantly enhanced in the Ad_PCK1 group on day 6 of adipogenesis induction ( Fig. 6e and f, p < 0.01).

Transcriptional activity of bovine PCK1 promoter
Although the expression pattern of PCK1 between muscle and adipose tissue in buffaloes was similar to UXT was used to normalize the expression level of PCK1, and their prime sequences were conserved between buffalo and cattle. Expression level data are presented as the mean ± S.E. (RNA sequencing, n = 3; RT-qPCR, n = 6); lower case indicates p < 0.05; upper case indicates p < 0.01 that in cattle, PCK1 in muscle was more abundant in cattle than in buffaloes, and the opposite finding was observed in adipose tissue ( Fig. 5a and 3b). These results suggested that PCK1 in muscle was more active in cattle; in contrast, PCK1 in adipose tissue was richer in buffaloes. Then, the transcriptional activity of the PCK1 promoter of buffalo and cattle was further analysed. The upstream putative promoter of PCK1 (Fig. 7a) was contained in the luciferase reporter construct and transfected into 3T3-L1 cells, C2C12 cells, primary bovine adipocytes, and primary bovine myoblasts. In 3T3-L1 cells and primary bovine adipocytes, the PCK1 promoter of buffalo demonstrated higher transcriptional activity than that of cattle ( Fig. 7b and d, p < 0.05). In C2C12 cells and primary bovine myoblasts, the PCK1 promoter of cattle presented higher activity than that of buffalo ( Fig. 7c and e, p < 0.05).

Discussion
IMF deposition is a highly complex process that can be affected by genetic, managerial, and nutritional factors [8]. Among these factors, genetics is a significant driver. Japanese Black is the best beef breed and exhibits more than 20% IMF content at 24 months of age [11]. By comparison, other European breeds, such as Belgian Blue and German Angus, have only 0.6 and 4.4% IMF content, respectively [11]. IMF levels in Chinese breeds are similar to those of European breeds. In Qinchuan cattle, the IMF content ranges from~3% to~4% at 18 months of age [21]. Nanyang cattle at 24 months of age have an IMF content of (1.88 ± 0.64)% [22]. In the present study, the IMF content of Nanyang cattle at 30 months of age was 2.43% (Table 1). In contrast, Xinyang buffalo had only a 0.51% IMF content (Table 1). Accordingly, the back fat thickness and LMA in Nanyang cattle were higher than those in Xinyang buffalo. Xinyang buffalo and Nanyang cattle were raised with similar forage and feeding management conditions. Thus, the difference in IMF content in animals between the two groups may be attributed to genetic factors.
To identify candidate genes for IMF deposition, a large number of studies have been carried out in cattle [10]. Comparisons of gene expression between breeds with variable IMF contents have indicated that some adipogenic transcription factors are associated with IMF deposition in cattle [7,12,13]. Because the IMF content in the tissue was too limited to be sampled for transcriptome analysis or RT-qPCR assay, both muscle (containing very limited amounts of IMF) and adipose tissues were used to identify genes that explain the difference in the IMF content between buffalo and cattle according to the strategy shown in Fig. 2. Simply stated, the coexpressed mRNAs in muscle and adipose tissues were considered to be part of the mRNAs expressed in IMF tissue. Moreover, muscle-specific mRNAs also affect IMF deposition [18,19]. Both the coexpressed mRNAs and muscle-specific mRNAs that showed differential expression between buffalo and cattle muscle tissues were used for functional enrichment.
The number of mRNAs identified in buffalo muscle and adipose tissue is similar to that of cattle ( Fig. 1a and  b), which indicates a similar level of complexity between buffalo and cattle meat. However, there were 1566 coexpressed mRNAs (Additional file 4: Table S3) and 70 muscle-specific mRNAs (Additional file 5: Table S4) that showed differential expression between the two bovid species, suggesting a significant difference between buffalo and cattle IMF deposition. Of the 1566 DE mRNAs, 112 and 38 were also identified between Wagyu and Piedmontese [7] and between Japanese Black and Holstein Friesian breeds [23], respectively. However, only four of these DE genes were shared by the three studies, and they showed variable expression patterns among the three studies. These results indicate that the DE genes can be variable from different breed pairs and suggest a highly complex genetic regulation mechanism in IMF deposition.
Based on its specific location, IMF deposition can be affected by both the activity of adipocytes and musclespecific factors. Coculturing C2C12 and 3T3-L1 preadipocyte cells downregulates adipogenic marker gene expression [18]. Myostatin, which is secreted by myoblasts, decreases IMF deposition by inhibiting the differentiation of intramuscular preadipocytes [19]. However, a similar factor that is responsible for the difference in IMF deposition between buffalo and cattle was not identified in this study; all 5 identified GO items with Benjamini values (p) < 0.05 were associated with muscle development (Additional file 6: Table S5).
Lipid metabolism is notably different between buffalo and cattle, which is also reflected by the difference in milk fat content [24,25]. The percentage of milk fat in buffaloes is higher than that in cattle [24,26,27]. A previous study suggested that four genes (KCTD8, FOXO4, SSTR3, and LOC782855) related to carbohydrate metabolism and a gene (ESRRG) related to lipid metabolism may contribute to the differences between buffalo and cattle milk fat content [25]. In the present study, however, these five genes did not show differential expression between buffalo and cattle muscle tissues, and no item marked with "carbohydrate metabolism" was identified in our functional enrichment analysis. Instead, an energy metabolism pathway, the glycolysis/gluconeogenesis pathway (p < 0.01), was identified by the functional enrichment of the DE coexpressed mRNAs ( Fig. 2; Additional file 7: Table S6). Glycolysis is a catabolic pathway that reverses the carbon flow of anabolic pathways, gluconeogenesis and glyceroneogenesis [28]. Essentially, glyceroneogenesis is an abbreviated version of gluconeogenesis responsible for the de novo synthesis of glycerol-3-phosphate from the precursors of gluconeogenesis [29]. Glycerol-3-phosphate can esterify fatty acyl-CoAs into triglycerides (TGs). Thus, the carbon for TG comes from glycolysis. The hyperactivity of glyceroneogenesis enhances the storage of TG and results in fat deposition in muscle [30].
A total of 17 DE genes were involved in the glycolysis/ gluconeogenesis pathway. Among them, PCK1 drives a significant and irreversible reaction that is shared by glyceroneogenesis and gluconeogenesis [28]. In fact, PCK1 is considered an adipogenic marker, an obesity gene, and even an IMF deposition-related gene [17,31]. In mice, white adipose tissue-specific overexpression of PCK1 increases the rate of glyceroneogenesis in adipose tissue and results in obesity [32]. In contrast, knockdown of PCK1 expression in adipose tissue results in lipodystrophy in mice [20]. In skeletal muscle tissue, specific overexpression of PCK1 enhances IMF deposition in mice [17]. In pigs, the expression level of PCK1 demonstrates a positive association with the IMF content and fatty acid composition [31]. Functional mutation of the PCK1 gene results in a lower IMF content in pigs [33]. Interestingly, pigs with PCK1 specifically expressed in skeletal muscle show an enhanced IMF content [34]. In cattle, PCK1 is upregulated in Wagyu breed (high IMF) compared to Piedmontese breed (low IMF) [7]. In the present study, PCK1 was downregulated in buffalo muscle tissue ( Fig. 5a and b), which may be responsible for the lower IMF content in buffaloes (Table 1). Interestingly, overexpression of PCK1 promoted the expression of the adipogenic marker PPARG (Fig. 6c, p < 0.05) and the lipid accumulation in buffalo intramuscular adipocytes ( Fig. 6e and f, p < 0.01), suggesting a positive role of PCK1 in IMF deposition in buffalo. The expression level of fatty acid uptake gene FABP4 was upregulated by PCK1 on day 6 (Fig. 6d, p < 0.01). FABP4 is a fatty acid binding protein and promotes fatty acid transport [35]. As mentioned above, PCK1 drives glyceroneogenesis to produce glycerol-3-phosphate [28]. Both fatty acid and glycerol-3-phosphate can be utilized for TGs synthesis [28,29]. Therefore, overexpression of PCK1 enhances lipid accumulation by upregulating the expression of PPARG and FABP4 in buffalo intramuscular adipocytes, which is consistent with the effects of PCK1 in other species (20,35).
We hypothesized that there should be some intrinsic reasons for the lower expression of PCK1 in buffalo muscle than in cattle muscle. In terms of the PCK1 putative promoter sequence, differences can be found between buffalo and cattle, with 44 single-base, 2 double-base, and 1 tribasic base mutations (Additional file 8: Fig. S2). Accordingly, the activity of the PCK1 promoter demonstrated variable levels between buffalo and cattle (Fig. 5). Interestingly, the activity of the buffalo PCK1 promoter is higher than that of cattle in adipocytes, while the opposite occurs in myoblasts (Fig. 5). PCK1 is widely expressed in multiple tissues and has high levels in the liver and adipose tissue [36], which also suggests its important role in lipid metabolism. The molecular mechanism driving PCK1 expression is highly complex, and its expression in different tissues can be regulated in a variable manner according to its specific role in the tissue. The promoter of PCK1 has been widely studied, and some significant transcription factors have been identified, such as PPARG, sterol regulatory element binding protein 1, and CCAAT/enhancer binding protein beta [36]. Recently, a large number of noncoding RNAs have been demonstrated to participate in the regulation of coding genes via interaction with promoters [37,38]. In addition, methylation in the promoter sequence can influence the expression of a gene [39]. Therefore, the underlying molecular mechanism driving the expression of PCK1 in the promoter regions of buffalo and cattle, such as DNA methylation, protein-DNA interactions, and RNA-DNA interactions, needs further study.

Conclusions
This study evaluated the IMF content of fattened buffalo and cattle and compared the mRNA expression patterns in muscle and adipose tissue between the two bovid species. This study demonstrates that (1) the IMF content of fattened cattle at 30 months of age is higher than that of fattened buffalo, as well as other relative traits, including back fat thickness and LMA; (2) a large number of mRNAs have differential expression between buffalo and cattle muscles; (3) the glycolysis/gluconeogenesis pathway may regulate IMF deposition in buffalo; (4) PCK1 shows a positive association with the IMF content in buffalo; and (5) PCK1 promoter activity of buffalo is higher than that of cattle in adipocytes, which is in contrast to that of myoblasts.

Methods
The IMF content in buffalo is lower than that in cattle, and a genetic factor should be responsible for such a difference. The aim of this research is to identify the genetic factors responsible for the low IMF content in buffalo. Because the IMF content in tissue is too limited to be sampled for detection, both muscle and adipose tissues were used for high-throughput RNA sequencing and further validation.

Animals, sample collection, and fattening traits
Xinyang buffaloes (n = 12) and Nanyang cattle (n = 12) were raised in the breeding farm of Xinyang Buffalo (Guangshan, Henan Province, P.R. China) and the breeding centre of Nanyang cattle (Nanyang, Henan Province, P.R. China), respectively, with equivalent forage and feeding management conditions. They were weaned at 3 months of age, castrated at 6 months of age, started to fatten at 18 months of age, and slaughtered at 30 months of age. Details of the feeding conditions were described in a previous study [40]. In this study, only six buffaloes and six cattle were fattened, and the others were not fattened. For the unfattened animals, concentrate was replaced by roughage from 18 to 30 months of age.
The LM (muscle) and back subcutaneous fat (adipose) tissues were sampled immediately after slaughter. For determination of IMF, fresh muscle was kept in an icebox and taken back to the laboratory immediately. For RNA extraction, tissue was frozen immediately in liquid nitrogen and stored at − 80°C.
The back fat thickness, LMA, and IMF content of fattened animals were measured. The 12th-13th rib back fat thickness and LMA were measured by an Ultrasonic Backfat Tester (HS-1600 V; Honda Electronics, Toyohashi, Japan) at slaughter. The IMF content was measured by a near-infrared meat quality analyser (Series 3000 Food Analyser, NIR, Australia).

RNA isolation and sequencing
Total RNA was extracted using TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. The quality and concentration of total RNA were determined by a NanoDrop 2000 (Nanodrop, Wilmington, DE) and 1% agarose gel. Total RNA from muscle and adipose tissues (three buffaloes and three cattle) of high quality was sent to Genergy Biotechnology Co., Ltd. (Shanghai, China) for library construction and RNA sequencing. Briefly, approximately 2 μg of total RNA per sample was used for library construction according to the protocol of the TruSeq® RNA LT Sample Prep Kit v2 (Illumina, USA). The cDNA library was sequenced using the double terminal sequencing mode of the Illumina HiSeq 3000 platform. A total of four groups that included the LM and the back fat tissues of buffalo and cattle were used for RNA sequencing. Three samples were obtained for each group.

Bioinformatics analysis of sequencing data
Low-quality reads and those containing adapters were removed to obtain clean reads using Trim Galore (Version 0.4.2, Babraham Institute, Cambridge, UK) [41]. High-quality data were used for the subsequent analysis. Since annotation information of the buffalo genome was not available, the cattle genome (UMD3.1) was used in the present study. Cattle genome and gene model annotation files were downloaded from the Ensemble database (http://www.ensembl.org/index.html). Reads were aligned to the reference genome using STAR [42]. Mapped reads were assembled using Cufflinks v2.2.1 [43]. These assemblies were merged to the reference gene model annotation by Cuffmerge [44]. Fragments per kilobase of transcript per million mapped reads (FPKM) was used as an index to calculate the expression level of each transcript in every library using the Cuffdiff program. The expression level was presented as log2(FPKM+ 1). The DE mRNAs were identified by DESeq2, with an absolute value of log2(fold change) ≥ 1 and an FDR ≤ 0.05 [44]. DE mRNAs were presented in heatmaps by the R package [45]. Functional enrichment for DE mRNAs was performed using DAVID bioinformatics resources (http://david.abcc.ncifcrf.gov/) [46]. The human database was used. GO terms (including molecular function (MF), biological process (BP), and cellular components (CC)) and KEGG pathways with a Benjamini value (p) < 0.05 were considered to be significantly enriched.

RT-qPCR
Primers were designed by Primer-BLAST (http://www. ncbi.nlm.nih.gov/tools/primer-blast/) ( Table 2). The ubiquitously expressed prefoldin-like chaperone (UXT) gene was used to normalize the expression level of the candidate gene [47,48]. Total RNA was transcribed into cDNA using the PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). qPCR was performed using SYBR Green I (TaKaRa, Dalian, China) with twostep reactions according to the manufacturer's recommended protocol. The efficiency of amplification was detected for each primer pair to ensure a reliable result. The cycle threshold (2 -ΔΔCt ) method was used to calculate the relative expression level of candidate genes. The RT-qPCR analysis was performed with six biological repeats, and each sample had three technological repeats.

Adenovirus packaging
The overexpression adenovirus packaging was performed by the Hanbio Biotechnology Co.; Ltd. (Shanghai, China) as previously described [40]. Briefly, the CDS of buffalo PCK1 was synthesized and ligated to the AdMax system. EGFP was included to indicate the transduction efficiency. Ad_EGFP was used as a negative control.

Isolation of buffalo intramuscular preadipocytes and cell culture
Buffalo intramuscular preadipocytes were isolated as described in a previous study [49]. Briefly, fresh LM tissue was sampled and the visible connective tissue was removed. Muscle tissue was finely minced and digested with 0.2% collagenase type II (Gibco, Grand Island, NY) in a shaking water bath at 37°C until complete digestion. The digested sample was filtered by 80 μm nylon mesh filters and washed three times with Dulbecco's modified Eagle's medium (DMEM) (Gibco, Grand Island, NY). Cells were then cultured with complete culture medium (DMEM with 10% foetal bovine serum and 1% penicillin-streptomycin (Gibco, Grand Island, NY)) under conditions of 37°C with 5% CO 2 for 1 h. The complete culture medium was changed to remove the nonadherent cells. The adherent cells (intramuscular preadipocytes) were continuously cultured with complete culture medium. When reaching 90% confluence, the cells were digested and collected for further culture.

Adenoviral transduction, adipogenic differentiation, oil red O staining, and quantification
The intramuscular preadipocytes were planted in twelvewell plates. Adenoviral transduction was conducted when the cells reached 80% confluence. Two days after transduction, the cells were induced to adipogenic differentiation. Six days later, Oil Red O staining was performed. The Oil Red O was eluted by isopropanol to quantify the lipid accumulation in cells. In addition, cells induced to differentiation for 6 days were collected for RT-qPCR analysis. Details of adenoviral transduction, adipogenic differentiation, Oil Red O staining, and quantification can be found in a previous study [40]. All these experiments were performed with six biological repeats.

Vector construction
The pGL3-Basic system (Promega, Madison, USA) was used for the promoter activity analysis of PCK1 as previously described [50]. The putative upstream promoter regions were cloned from the genomic DNA of cattle and buffalo ( Fig. 7a and Additional file 8: Fig. S2). Fragments were inserted into the KpnI and XhoI restriction sites of the pGL3-Basic vector to obtain recombinant plasmids (pGL3-cattle-PCK1 and pGL3-buffalo-PCK1). The pLR-TK vector was used as an internal control.
Cell culture, transfection, and luciferase reporter assays 3T3-L1 cells (ATCC, Shanghai, China), C2C12 cells (ATCC, Shanghai, China), primary bovine adipocytes, and primary bovine myoblasts were used to detect the activity of putative promoters. Isolation and culture of primary bovine adipocytes and primary bovine myoblasts were performed as previously described [40,51]. Cells were cultured with complete culture medium in 5% CO 2 at 37°C. Cells were plated in a 48-well plate at a density of 1 × 10 6 cells/well 24 h prior to transfection. Cells were cotransfected with 400 ng recombinant plasmid DNA and 16 ng pLR-TK plasmid using Lipofectamine® 3000 . The ratio of the firefly luciferase value to Renilla luciferase value was used to evaluate the activity of the PCK1 promoter fragment. Each promoter fragment plasmid was analysed with six repeats, and the assay was repeated three times.

Statistical analysis
Data were analysed by SPSS 19 software. Significant difference analysis was performed with one-way ANOVA. A p-value < 0.05 was considered to indicate significant differences. The results are presented as the means ± standard error (S.E.) using the OriginPro 8.5 program.