Functional Genes in Relation to Residual Feed Intake in Murrah Buffalo Heifers ()
1. Introduction
Feed cost contributes to more than 80% of the total maintenance cost in large animals [1] and Residual feed intake (RFI) [1] is determined as the difference between the actual and predicted feed intake [2] , and also defined as the feed conversion efficiency (FCE) in cattle and buffalo [3] . Low RFI heifers end up with additional genetic gain and profitability by reducing 10% to 14% feed consumption, i.e., cutting short on input cost [4] on achieving a comparable growth rate as of high RFI counterparts. RFI has been used as a promising tool to identify feed efficient animals [5] , however, reports in buffaloes are scanty. Feed efficiency trait heritability is determined as 0.01 - 0.3 [1] and is improvable through genetic selection [6] . Variation in RFI has been reported ranging between (−0.20) and high (+0.18) in growing calves [3] and Murrah heifers [7] . Feed intake, nutrient partitioning, immune function and metabolic adaptation in different body tissues are the factors known for influencing feed utilization in animals. A reduction in dry matter intake without compromising the weight gain is expected to improve the overall productivity of animal. Recording of feed intake for an individual animal is arduous and biological markers have been reported (Sikka et al., 2020 [7] ; Baban et al., 2021 [8] ) in calves to enable marker-based selection of feed efficient calves. However, the identification of variants in this respect remains a challenge [8] for employing genomic selection for this trait. The present study is conducted with the aim of determining differentially expressed genes (DEGs) and comparing high and low feed utilizing heifers. Transcriptomic data was collected from low and high RFI heifers using whole blood. The study catalogued the functional genes that were differentially expressed between high and low RFI young buffalo heifers.
2. Material and Methods
2.1. Study Location, Animals and Sample Collection
Animal experiments were performed under approval and review by the Institutional Animal Ethics Committee (IAEC) at ICAR-Central Institute for Research on Buffaloes Hisar, Govt. of India. Heifers of 9 to 11 months age were subjected to control feeding in order to record individual dry matter intake (DMI) to determine residual feed intake (RFI) by estimating the difference between predicted and actual DMI for an individual animal. Dietary regimen of heifers included Jowar fodder and concentrate mixture. Nutrient value of feed was determined using standard methods. Chemical composition and nutrient digestibility of feed (g/kg) and nutritive value of diet were determined by reference methods.
Nutrient digestibility (Mean ± SE) of dry matter was 60.4 ± 2.43. Percentage of crude protein (CP), digestible crude protein (DCP) and total digestible nutrients (TDN) was determined as 14.2 ± 0.13, 12.5 ± 0.33 and 61.2 ± 1.77, respectively. Blood biochemical profiling revealed IGF1 and Glucose as significant biological attributes related to RFI efficiency in heifers [8] .
aValues represent hex plicate assays of each material; bIngredient composition of concentrate mixture: [maize grain, 175 g/kg; barley grain, 175 g/kg; wheat bran, 270 g/kg; mustard cake, 200 g/kg; cotton seed cake, 150 g/kg; mineral mixture, 20 g/kg; common salt, 10 g/kg].
Daily feed intake and fortnightly body weight were recorded for heifers during the period of 78 days feeding trial. Average DMI (Kg/h/d) and daily weight gain (ADG) were determined. Average BW was attained as 200 Kg in heifers after 78 days of feed trial initiating from 155 Kg. Mean ADG and DMI per animal was recorded as 577 g/h/d and 4.8 kg/day, respectively (Table 1). Variation in daily weight gain and DMI was recorded ranging from 382 to 897 g/day and 3.3 to 6.0 kg/d over the 42 heifers during the trial. Mean DMI was regressed with RFI at (R2 0.970) in the range of 0.359381 to −0.43778 at high and low extremes (Figure 1). DMI was recorded as 2.41 Kg lesser in efficient feed conversant group of heifers (low RFI) as compared with less efficient (high RFI). BW gain in these heifers was recorded as 614.84 g/h/d. Low RFI heifers gained 16.4% higher BW through lesser DMI by 13.53% (Figure 1).
A large variation in DMI per unit BW gain over was determined in heifers (n = 42) allocated to low and high RFI feed efficiency (Table 1). DGEs in respect of FCE were identified in heifers comparing transcriptomic data obtained from high and low extremes of RFI efficient heifers through high throughput (HTP) sequencing.
2.2. RNA Extraction, cDNA Library Preparation and Sequencing
Total RNA from blood cells collected through veni-puncture was extracted from
![]()
Table 1. Quality check statistics for RFI transcriptome in buffalo heifers.
![]()
Figure 1. Regression of actual with predicted dry matter intake (Kg/day) in Buffalo heifers (n = 42).
four heifers (two showing high RFI and two showing low RFI). 1 μg of intact RNA having integrity number (RIN) of 8.0 (Agilent 2100 Bio analyzer) was used to purify mRNA using oligo-dT beads (TruSeq RNA Sample Preparation Kit, Illumina) for each sample. cDNA library preparation and paired-end sequencing were done using denatured (90˚C) mRNA. The RNA fragments were reverse transcribed using random hexamers and Superscript II Reverse Transcriptase (Life Technologies). Second strand cDNA was synthesized on this first strand template, using RNaseH and DNA polymerase I. Clean cDNAs were obtained using Beckman Coulter Agencourt Ampure XP SPRI beads after end-repair and the addition of an “A” base. The clean cDNA molecules were ligated with Illumina adapters. Further, the cDNA library was amplified using PCR for the enrichment of the adapter-ligated fragments. The individual libraries were quantified using NanoDrop spectrophotometer (Thermo Scientific) and validated for quality by Bio analyzer (Agilent Technologies) for subsequent sequencing, using the Illumina HiSeq 2500 platform. Paired-end FASTQ files were subjected to standard quality control. High quality (HQ) filtered reads (Phred scores > 20) were selected using the NGSQC Tool Kit [9] .
2.3. Assembly, Annotation and Global Expression Profiling
RNA extracted from blood of the four heifers (Bubalus bubalis), covering extremely high and low RFI efficiency as estimated based on DMI/h/g BW gain as FCE1:FCE2: vs FCE3:FCE4 was subjected to standard paired end RNA-Seq library(s) preparation as per Illumina recommended protocol using 100 bp paired-end module. Approximately 40 - 60 million paired end sequence reads (Phred score > Q30) were obtained for each library and the same were subjected to stringent quality test using NGSQC tool kit [9] .
The improved transcripts resulting from the trans-assembly were subjected to the CD-HIT EST clustering pipeline. Redundant transcripts were removed and reference unigene transcriptome was generated from each library, all converging finally into a single non-redundant resultant transcriptome assembly by merging and clustering the transcripts from individual libraries. Resultant transcriptome assembly was aligned with the reference genome sequence, using Bowtie tool for quantification of expressed transcripts. Read count for each transcript was further normalized to determine value of RPKM. Transcripts with an RPKM of ≥1 were categorized as expressed transcript. Treatment specific transcripts were annotated with Gene Ontology and Pathways identified from expressed transcripts by DAVID Functional Annotation Tool [DAVID Bioinformatics Resources 6.8, NIAID/NIH [10] .
2.4. RNASeq Data Analysis
The Reference genome of Bubalus bubalis (https://www.ncbi.nlm.nih.gov/assembly/GCF_003121395.1) was used for read alignment and identification of transcripts coding regions using Kallisto followed by quantification and annotation [11] . Differentially expressed (DE) transcripts (fold-change ≥ 2.0, P-value < 0.05) were identified by DESEQ2 [12] which revealed underlying genetic make-up of efficient RFI, compared with high RFI as control.
2.5. Biological Pathways and Gene Ontology Analysis
The molecular functions, biological mechanisms and gene networks that emerged from the transcripts which were over-represented were determined. Dataset containing gene identifiers, corresponding expression and P-values were uploaded into IPA. Identification of de-regulated genes was done based on significant P-adj value (P < 0.05). The “focus” genes were overlaid onto a global molecular network using Cytoscape software [13] to lineate bio-molecular interaction networks into a unified conceptual framework by evaluating a score ranking the network obtained according to their degree of relevance to the eligible molecules in the dataset. The score is the result of the number of knowledge base considerations of the network eligible molecules into account out of the total number of network eligible molecules that emerged and the total number of possible molecules to be included in the network.
Ethics Approval: The protocols followed in the present animal study were reviewed and approved by the Institute Animal Ethics Committee (IAEC) of CIRB Hisar. India (Reg. No. 406/GO/RBI/L/01/CPCSEA).
3. Results and Discussion
3.1. Sequencing and Mapping Statistics Summary
An average of 51.9 million (varying between 36.9 - 72 million) paired end raw RNA Seq reads were obtained per sample using Illumina HiSeq platform. Average size of reads qualifying for all the QC criteria was 210 bp (Table 1). FASTQ sequences were filtered using NGSQC tool kit thus, an average of 52 million high quality reads per sample (~92.7%) were accounted for alignment and variant detection analysis. Around 88.7% of the HQ reads of created assembly were mapped on the human reference genome using KALLISTO pipeline, which confirms optimum alignment (Figure 2) of the two genomes.
3.2. Transcripts Expression in High and Low RFI Heifers
Per sample high quality reads obtained from all the heifers were comparable covering nearly 9% of the total reads and having average GC content of 47.22% (Table 1).
RNAseq dataset obtained from selected heifers were evaluated for DEGs including 49781 significant transcripts (Padj = 2.44995163703878E-54 to 0.0504821201656885). Percent alignment of reads with reference genome of Bubalus bubalis varied between 64.1% to 67.7% with respect to heifers (Figure 2) in the study. Differentially expressed (P < 0.05) transcripts were 850. Correlation matrix coefficients were determined as 0.821*** between replicates of high RFI, i.e. FCE1:FCE2 and 0.967*** between replicates in low RFI individuals i.e. FCE3: FCE4 (Table 2). It indicates high similarity in genetic make-up underlying the selected heifers as per the estimated RFI efficiency. Correlation coefficients while regressing the transcript data across high and low RFI category of heifers was significantly low (r2 0.851***) for FC1:FC3 and FC2:FC4 (r2 = 0.816***) while comparing within high or low RFI.
RFI of an individual heifer is influenced by digestibility, tissue metabolism and protein turnover [14] under the changing patterns in environmental stress, thermoregulation, heat increment, and body composition [15] while evaluating RFI efficiency having uniform feeding regimen.
![]()
Figure 2. Alignment percentage of HQ reads of Low RFI [FCE 1: FCE 2] and high RFI [FCE 3: FCE 4] Murrah heifers Referring Bubalus bubalis.
3.3. Differential Expression of Transcripts
Distinctive pattern of transcripts expression is depicted through Heatmaps (Figure 3(a) and Figure 3(b)) generated from high and low RFI efficient heifers. De-regulated genes (P-adjusted < 0.05) are presented as fold change value (>2 or <−2) expressing variation in RFI (Table 2 and Table 3). A total of 850 DEGs (p < 0.05) were detected based on the Cuffdiff analysis (FDR = 0.00174998654944857 to 58.5661012869287 (FDR < 0.05). 427 genes were down-regulated and 420 genes were up-regulated in the study, where top thirty genes in each category of deregulation are shown in Table 2 and Table 3.
3.3.1. Total 850
DEGs were identified from 49,781 transcripts, while comparing transcriptomic data obtained from high RFI, ranging between 0.26711 to 0.359381 and low RFI (−0.43778 to −0.20372) category of heifer. The top DEGs are mostly annotated under the Biological processes (>25.08%), Cellular component (21.67%) or Molecular function (20.64%) (Figure 5) of the GO terms (P 1.31569713353195E−06 to 0.0598850395935369) reported in Table 4.
Down and Up-Regulated DEGs in Relation to RFI
Top down-regulated genes (Padj < 0.05) out of 421 down-regulated genes related to RFI (low efficiency) are discussed as given in Table 2. Cell growth regulator (earlier reports annotated by GO terms related to cell cycle (negative regulation of cell growth, negative regulation of chromosome condensation, telomeric heterochromatin assembly, negative regulation of G0 to G1 transition), Ring finger Protein 213 and E3 ubiquitin-protein ligase gene related to adaptive immunity were emerged as least fold change (log2FoldChange = −11.535) and cell signaling (endocytosis) inducing Pleckstrin homology domain containing A2
![]()
Figure 3. Hierarchal clustering of deregulated transcripts in high RFI (a) and low RFI (b) Murrah Buffalo Heifers.
![]()
Table 2. Top (adj p < 0.05) down-regulated genes in relation to Residual Feed Intake (RFI) in Murrah Heifers.
gene showing maximum change (log2 FoldChange = −2.461) for down-regulation in efficient feed conversant(low RFI) heifers in comparison to less efficient heifers (high RFI) in this study.
Down-regulated Pleckstrin homology domain containing A2, an endocytsis & cell signalling [PLEKHA2]gene, Ring finger protein 213-like [ATPase, adaptive Immunity] [LOC512869], IFI27, Dipeptidyl peptidase 8 [blocking T-cell activation and immune, MHC class I heavy chain, Sterile alpha motif domain containing 9 [TNF-alpha signalling] and Lipopolysaccharide-binding protein [Gm –ve bacterial infections] genes fighting out the gram negative bacterial infection.
![]()
Table 3. Top (adjp < 0.05) Up-regulated genes in relation to Residual Feed Intake (RFI) in Murraj Heifers.
Down-regulation of Thyroid hormone receptor inter actor 12 [Innate Immune System gene associated with presenting Class I MHC mediated antigen and Innate Immune System; IFI44L gene regressing the Interferon stimulating antiviral proteins and CD274 gene, blocking T-cell activation building-up innate immunity and T cell activation for developing feed efficiency in heifers. (RNF144A) ring finger protein 144A promoting ubiquitination, down-regulated transcription regulator genes as (LOC407163) trappin 5 rate of transcription, 5OH methyl cytosine binding, (ES cell specific), zinc finger CCCH-type containing 11a genes having role in n mRNA export were down-regulated expressing the lower rate of chromosome segregation for replication due to down-regulating BUB3 gene and transforming growth factor beta (TGF-β) protein expression consequent upon lower expression of ribosomal protein S15a gene thus, hampering the cell growth
![]()
Table 4. Gene ontology terms and biological pathways identified from differential genes in high and low residual feed intake efficiency buffalo heifers.
in high RFI (less efficient) heifers. Ribosomal protein L17 protein is a cell growth inhibitor protein which may improve the muscle growth while down-regulated. Down-regulation of the cell signalling genes may re-model mitochondrial membrane in order to improve ROS through electron transport chain and other GTPases based energy balancing functions in cytoplasm. Down-regulation of SEPTIN 2 (GTPase) and Rabaptin (GTPase binding protein 1) genes may affect the cell signalling due to altered cell membrane and cell binding effector proteins. Down-regulation of regulation on GTPase signalling (G protein 14) may deplete energy in cytoplasm, thus affecting cell functions and growth in heifers. Down-regulation of LOC100295883 gene locus may implicate affected endometrial development and embryo implantation due to its role in lipid synthesis and Cytochrome P450 proteins (monooxygenase). Down-regulation of G protein-coupled membrane receptors (GPCR) genes may affect binding of receptors-transcription factors (tf) do alter the downstream regulation of nuclear mRNA transport, cell signalling thus, altering endocytosis hampering protein degradation and cell homeostasis in low efficient heifers. Pre-mRNA modification in terms of down-regulated miRNA biosynthesis might lower the rate of translation, and even transcription regressing the overall feed efficiency RNA binding protein (RBP) transcripts as down-regulated Zinc finger CCCH (ZC3H7A) gene in high RFI heifers. Down-regulated LOC508666 gene may hamper development and optimize the functioning of reproductive system in high RFI heifers. These findings corroborate the earlier reports [16] .
Low rate of oocyte cell division and mitogenesis may be the response to the same effect. Down-regulation of genes governing steroid synthesis, adaptive/innate immunity, endometrium development co-expressed with the low level of interferon induction against viral infection, immunologic response to gram-negative bacterial infections, MHC class I heavy chain and TNF signaling were profound in expression in 9 to 11 months age heifers. Down-regulation of receptors mediated desensitization of signaling molecules for T-cell activation and cytokine production maintain homeostasis of the immune response in growing heifers.
3.3.2. Up-Regulated Genes
Total of 428 up-regulated genes were identified (Padj < 8.81509190325442E−41 to 0.000226709239889707). Top up-regulated genes (Padj < 0.05) comparing high RFI (low efficiency) heifers with low RFI heifers as control are listed in Table 3 and discussed ahead.
Up regulation of RPL 26 protein, a 60 s ribosomal subunit gene indicates enhanced mRNA translation in muscle in efficiently growing animals as reported in Angus [17] .
Up-regulation of RBPJ gene in this study may influence its negative control on phagocyte oxidative burst due to repression in NADPH oxidase transcription in response to bacterial infection. Endogenous transmembrane of CLM-1 (LOC522174) gene remodelling to hold preferred energy deriving cycles in mitochondria. Up-regulation of Innate Immune System representing gene locus (LOC101903261) which regulates TNF-induced cellular inflammation response; T cell signaling and its activation (TEC); expression of cytokines IL-2, IL-3, IL-4, IL-5, granulocyte-macrophage colony-stimulating factor, Bubalus bubalis CMRF 35 like mol. 6, an immune-regulatory signaling entity, CD 46 mol. propagating T-cell proliferation/differentiation may compensate for body growth [18] . Ubiquitination Protein gene (USP25) and a stress responsive cell cycle regulator, peptidase like activity (USPLI) genes was up-regulated in efficient heifers.
Up-regulation of ATP2C1 gene (Table 3) known for ATP binding signal molecule as ligand for Golgi complex by promoting Ca transport and its secretion in cytoplasm, Up-regulation of cell growth enhancers for promoting transport of glucose, other sugars, bile salts, organic acids, metal ions and amine compounds, Solute carrier family 9 member A7 [SLC9 A7] protein gene involving endocytic pathways do establish homeostatic balance in cells. Up-regulation of Cardiac muscle contrac 2, ANKFY1 and Recombinant Signal binding protein for Ig Kappa J region [RBPJ], might improve energy metabolism in liver through enhancing Glutamate dehydrogenase 1 [GLUD1] activity and advocating Ammonia detoxification by deamination of Glutamate to 2-oxoglutarate (ATP generation), channeling dietary lipids and adipose triglycerides to the mitochondrial respiratory chain due to induced oxidation-reduction process and ETC, generating high levels of energy against the oxygen stress [18] was negatively correlated (P < 0.05) with RFI efficiency and growth. Enhanced muscle mitochondrial respiration is known to associate with high RFI Angus steers [19] . High serum SGPT in this study [7] is linked with energy (NADPH) generation through Pyruvate diversion to Citric acid cycle [20] to compensate for growth, expectedly, in less efficient (high RFI) heifers. More than six GO terms related to glucose metabolism and signaling (glycolytic process, gluconeogenesis, pyruvate metabolic process were identified in response to insulin signaling (bta 04910) and Insulin resistance (bta 04931) identified in this study (Table 4) as ATP binding, AMPK signaling (bta 004152), GTPase activation (GO: 0090630, GO: 0005096, positive regulation of MAP Kinase (GO: 0043406). Also, GO terms related to oxidative metabolism (oxidation-reduction process, tricarboxylic acid cycle, proton transmembrane transport, oxaloacetate metabolic process, 2-oxoglutarate metabolic process, mitochondrial ATP synthesis coupled proton transport, and ATP biosynthetic process (Intracellular signal transduction GO: 0035556, microtubule cytoskeleton, GO; 0015630, Regulation of actin cytoskeleton (bta04810, GO: 0031532, GO: 0005884)and related genes which were identified in this study do corroborate with other reports [21] confirming the up-regulation of energy conserving genes and reduction in mitochondrial Oxidative phosphorylation do support feed efficient (low-RFI) animals.
Induced DNA-binding activity of transcription factors, e.g., MGA and NCOR1-HDAC3 genes (Table 4) may promote protein translation. Circadian expression of the core clock gene VIARTNL/BMAL1 nuclear receptor co-repressor 1 (NCOR1) gene known for altering lipid metabolism in liver may optimize body weight and growth related functions in response to oxidative stress. Fork head box transcription factor (FOXO), a Growth factor having potential role in expression of adherent genes and stress regulating transcription factors was up-regulated in this study thus, regulating cell cycle (bta 04514, GO: 0008283), energy intake and metabolic rate through Adipocytokine signaling pathway (bta04920) to combat with stress. FOXO, another transcription factor regulates Cortisol: glucocorticoid receptor complexing (GO: 2000323) in cytosol, translocating the same to nucleus and modulating transcription of a large battery of genes as fatty acid oxidation in muscles altering cell phenotype. mTOR protein kinase—a growth factor that coordinates cell growth (bta 04150) involving TNF signaling pathway (bta04668), immunity functions regulation by stimulating transcription factor Nuclear factor-kappa B (NF-kappa B) gene which regulates inflammation and cell survival. ErbB receptor tyrosine kinases pathway (bta 04012, Table 4) with downstream Mitogen-activated protein kinase (MAPK) pathway and phosphatidylinositol-3-kinase (PI-3K) pathway (bta: 04151, Table 4), for cell proliferation and differentiation through coupling extracellular growth factor ligands with intracellular signaling.
Recombination signal binding protein for immunoglobulin kappa J region (RBPJ) gene (Table 4) which negatively regulates the phagocyte oxidative burst in response to bacterial infection by repressing transcription of NADPH oxidase may add to feed efficiency in efficient heifers. Tubulin beta 1 (TUBB1) gene indicates need for higher expression of platelets production and platelet release (Table 4). LOC10084s 8700 gene locus up-regulation might be responsive towards stress in response to the toxic elements in liver of efficient, low RFI heifers symbolizing n adaptation of cell metabolism to combat conditions of oxidative stress. LOC523126 MRP4, a functional Prostaglandin carrier molecule known as regulator of estrous cycle up-regulated may be indicative of oocytes differentiation due to enhanced expression of XPO6 gene in efficient heifers which might induce early puberty. High metabolic activity is correlated with higher DMI for partitioning of energy for puberty attainment in less efficient heifers (22). CD46 gene is known for governing spermatozoa: Oocyte fusion which supports its role in fertilization and conception in buffalo heifers, while up-regulated. Present study identified Dynein gene (Axonemal dynein light chain domain containing 1 protein, AXDND1) up-regulation (log2FoldChange = 4.814, 4.814) (Table 3 and Table 4) in heifers, playing a major role in energy production for sperm motility, [22] [23] thus having its suggestive role in reproduction.
3.4. Gene Enrichment Analysis
Significant GO terms enriched and the pathways identified (https://david.ncifcrf.gov/) in this study are listed in Table 4. Selected GO terms (Padj ≤ 0.05) were used as input to plot gene network (Figure 4) in respect of FCE trait. Clusters (Padj < 0.05) that emerged in Cytoscape showed Intracellular
![]()
Figure 4. Gene network related to residual feed intake (feed conversion efficiency) trait in buffalo heifers.
Transduction signals (GO: 0035556) as pronounced function in Biological Process (4.3%) followed by Cell migration and Cell proliferation (Figure 5).
Cytoplasm functions (GO: 0005737) were found more pronounced than Nucleoplasm (GO: 0005654) due to a higher number of genes (103) annotated with former than the latter (65 genes), respectively (Table 4), sharring 9.8% of Cellular component (Figure 5) term, however, transcripts of genes governing transcriptional regulation in nucleus were emerged and annotated by Nucleoplasm (GO: 0005654) significantly (Padj < 0.05) (Table 4). GO: terms annotated by high rate of energy dependent translational changes as ATP binding (GO: 0005524), Metal Ion binding (GO: 0046872), Zinc Ion binding (GO; 0008270, bta04151) (Table 4), corroborate with the energy related genes reported earlier as Rap 1 gene, [24] which harvests MAP kinase activity [25] to regulate body weight and feed intake in cattle [26] [27] . The Rap1 (bta 04015) and Ras signalling (bta 04014) pathways emerged as enriched (P < 0.05, FDR in range of 0.00174998654944857 to 58.5661012869287) pathways (Table 4) as regulator of fatty acid content in response to retinoic acid and glycerol transport [28] in higher feed efficiency heifers.
![]()
Figure 5. Top 10 Significant GO and pathway for feed conversion efficiency trait.
Up-regulation of regulatory function genes for ATP binding (bta 04152: AMPK signalling pathway (Table 4) in B cell activation in present study along with MHC class I protein complex (presenting intra-cellular degraded proteins to cytotoxic T cells); GTPase activity, cell proliferation and differentiation, proteins translocation, signal transduction through transmembrane receptors, followed by subsequent protein synthesis at ribosome, as reported earlier [29] and Platelet activation (bta04611, Table 4) [30] are expected physiological processes in efficient feed utilizing (low-RFI) heifers for achieving higher weight gain by reduced feed intake thus, lowering the input price as compared to rearing the less efficient counterparts.
3.5. System Biology Analysis
Biological analysis of transcripts (Padj < 0.05) obtained from system biology analysis was summarized in GO summary (GO: terms and pathways) (Table 4). The DEGs identified were annotated by 93 GO terms/pathways enriched, allocated to the category of post transcriptional regulation and chromatin organization in Nucleus, Cell energy balancing, Immunity, hormonal regulation functions and Cell signaling including the biological pathways as Rap1, Notch, Prolactin, Osteoclast differentiation, m TOR, Ras, PI3K-AKT, FoxO, TNF, ErbB and Sphingolipid. High RFI heifers were associated with high IGF1 (p < 0.05) and other blood parameters such as Cholesterol and Triglycerides (p < 0.001) as reported in buffalo heifers [7] and cattle [31] .
Feed utilization is an energy dependent function requiring availability of higher reducing power to maintain optimum productivity [7] . Metabolic profile differs in high and low RFI heifers. High level of IGF1 has a positive correlation with LDL, HDL fat and SGPT enzyme in this study [7] and might favor the deposition of fat by channeling free triglycerides (P < 0.05) toward muscle tissue to gain weight in less efficient animals. Lower blood level of non-esterified fatty acids and higher muscle growth has been reported in low RFI cattle [32] and pig [20] . Thus, IGF-1 is recognized as a potential biological marker for FE/RFI, in dairy heifers [7] and beef cattle [33] determining a negative correlation between RFI efficiency and circulatory level of IGF-1 [34] in heifers. However, the systemic IGF-I concentration is known to be influenced by environment also [7] [14] . In the contrary, a positive genetic relationship between plasma IGF-I level and RFI in sheep [35] has been reported. RFI DEGs may not be comparable with respect to age, breed and physiological status. However, feeding management and climate change may be the additive factors influencing the DEG patterns.
Differentially expressed genes in respect of RFI efficiency trait are listed the first time in Murrah dairy growing heifers. DG patterns were obtained comparing blood transcriptome [21] data obtained from low and high RFI heifers tried on controlled feeding trials. Feed efficiency is a complex functional expression resulting from the synergism between energy built-up (ATP synthesis) and growth reflected through basal metabolic rate, homeostatic control of body, immune response through lipid metabolism and hepatic inflammation respectively [18] .
A large variation (p < 0.05) in DMI and average daily weight gain (ADG) over young heifers maintained under common management was attributed to the difference in their respective feed utilization efficiency. Feed efficiency was estimated as the residual feed intake (RFI Kg/h/d) in this study. DMI (kg/h/d) was 13.14% lower in the low RFI heifers as compared to high RFI heifers. It hypothesized feed efficiency as a selection trait due to underlying variation expected in genetic make-up of these animals. Feed input cost can be significantly reduced by selective breeding of feed efficient heifers. Variation in RFI is pertinent to metabolism [36] of individual animal which is translated into low consistency in phenotypic growth and RFI.
Hierarchical clustering of differentially expressed transcripts obtained from high and low RFI heifers showed no common transcripts in categories of up and down regulated genes. It supports the selection proposal of the animals in high and low RFI subgroups, especially having common transcripts between the two groups being non-significant in enrichment analysis. The identified set of deregulated genes that were obtained from the analysis are considered as a fingerprint of differences in metabolic homeostasis. This leads to variation in growth rates that are governed by the underlying genetic potential of the heifers. Hence, these transcripts can be used to identify extreme FE phenotypes in buffaloes.
Immune response emerges as a potential body function underlying the biological variation in RFI in this study. Immunity is an energetically costly physiological process [37] as higher incidences of chronic inflammation are expected in high-RFI cattle, consequently making less energy available for growth. Down regulation of genes governing the immune functions indicates low immunity in the high-RFI heifers as reported earlier [19] , thus affecting the efficiency of feed utilization. It is consistent with the regular/increased level of immunity governing gene expression in energy efficient manner as reported in low-RFI cattle [38] [39] .
Efficient (low RFI) heifers showing down-regulated ROS function genes suggest the lower energy expenditure for metabolic function to achieve growth in feed efficient heifers. Genes related to lipid transport and energy production were up-regulated probably, in response to high ROS in the high-DMI group of less efficient animals. High-ADG is consistently increased in response to carbohydrate and lipid transport in favor of cell differentiation.
4. Conclusion
The present study highlights the genes underpinning RFI efficiency in buffalo growing heifers, which certainly is a less studied trait in Indian buffaloes. The DEGs identified from the study can be used as a resource of biomarkers for developing molecular signature markers for the selection of RFI traits. Underlying biological processes of FCE are complex and are influenced by diverse climates, feed properties consequent to different gut microbiota, and individual genetic predisposition. FCE determinant, i.e., RFI being a promising tool for selection has moderate genomic heritability in a range between 0.18 to 0.57 in cattle, but is independent of age, growth and body weight (BW) traits where, the feed intake, body weight, and weight gain do carry moderate heritability. Thus, the validation of DEs can be done by including animals of different developmental and physiological stages, with a view to examining their ability to maintain energy balance and sustain immunity. For more effective breeding plans, RFI DEGs and pathways identified in heifers will improve the knowledge of QTL links in favour of related functional traits such as fertility and growth predictions in heifers. Information on genomic regions harbouring the variant genes [40] can be utilized to analyze additional markers like CNV (Copy Number Variation) for developing genomic selection programmes.
Ethics Statements
The animal study was reviewed and approved by the Institute Animal Ethics Committee (IAEC) of CIRB Hisar. India (Reg. No. 406/GO/RBI/L/01/CPCSEA).
Supplementary Material and Data Access/Availability
RNA Sequence data has been submitted to the GenBank NCBI databases under SRA ID SUB5692555. SUB5733459 and the Bio Project ID: PRNA546485 Title: Locus tag prefixes: FHP27 (samn11959988) None (SAMN120577830).
Funding
The research work [ICAR-CIRB-IASRI Collaboration] was supported by the Network Project on Agricultural Bioinformatics and Computation Biology under the Centre for Agricultural Bioinformatics Scheme, ICAR-IASRI, Indian Council of Agricultural Research (ICAR), New Delhi, INDIA as/Project Code AGENIASRICIP201500900046.