Abstract
Systems biology approaches are used as strategy to uncover tissue-specific perturbations and regulatory genes related to complex phenotypes. We applied this approach to study feed efficiency (FE) in beef cattle, an important trait both economically and environmentally. Poly-A selected RNA of five tissues (adrenal gland, hypothalamus, liver, skeletal muscle and pituitary) of eighteen young bulls, selected for high and low FE, were sequenced (Illumina HiSeq 2500, 100 bp, pared-end). From the 17,354 expressed genes considering all tissues, 1,335 were prioritized by five selection categories (differentially expressed, harboring SNPs associated with FE, tissue-specific, secreted in plasma and key regulators) and used for network construction. NR2F6 and TGFB1 were identified and validated by motif discovery as key regulators of hepatic inflammatory response and muscle tissue development, respectively, two biological processes demonstrated to be associated with FE. Moreover, we indicated potential biomarkers of FE, which are related to hormonal control of metabolism and sexual maturity. By using robust methodologies and validation strategies, we confirmed the main biological processes related to FE in Bos indicus and indicated candidate genes as regulators or biomarkers of superior animals.
Introduction
Since the domestication of the first species, animal selection aims to meet human needs and their changes over time. The current main selection goals in livestock production are increase of productivity, reduction of the environmental impact and reduction of competition for grains for human nutrition (). Thus, feed efficiency (FE) has become a relevant trait of study, as animals considered of high feed efficiency are those presenting reduced feed intake and lower production of methane and manure without compromising animal’s weight gain (). However, the incorporation of FE as selection criteria in animal breeding programs is costly and time consuming. Daily feed intake and weight gain for a large number of animals need to be recorded for at least 70 days to obtain accurate estimates of FE ().
In the past years, several studies have been carried out with the aim to identify molecular markers associated with FE to enable a faster and cost-effective identification of superior animals (; ; ; ). However, for each population, different biological processes seem to be identified (; ; ; ). Probably, that is because FE is a multifactorial trait and many different biological mechanisms seems to be involved in its regulation (; ). It has been demonstrated that high FE animals present increased mitochondrial function (; ), less oxygen consumption () and delayed puberty (; ; ). On the other hand, low FE animals have increased physical activity, ingestion frequency and stress level (; ; ; ), increased leptin and cholesterol levels (; ; ; ), higher subcutaneous and visceral fat (; ; ), higher energy wastage as heat (; , ) and more hepatic lesions associated with inflammatory response (; ).
In the context of such a complex trait, we perform a multiple-tissue transcriptomic analysis of high (HFE) and low (LFE) feed efficient Nellore cattle across tissues related to endocrine control of hunger/satiety, hydric and energy homeostasis, stress and immune response, physical and sexual activity, as is the case of hypothalamus-pituitary-adrenal axis and organs as liver and skeletal muscle. Using gene co-expression across tissues and conditions, we derived a regulatory network revealing NR2F6 and TGFB1 signaling as key regulators of hepatic inflammatory response and muscle tissue development, respectively. Next, we applied advanced motif discovery methods which (i) validate that co-expressed genes are enriched for NR2F6 and TGFB1 signaling effector molecule SMAD3 binding sites in their 10 kb upstream regions and (ii) predict direct transcription factor (TF) – Target gene (TG) interactions at the sequence level. These binding interactions were experimentally validated with public TF ChIP-seq from ENCODE (; ). Regulatory activity in the tissues of interest was also confirmed by performing an enrichment analysis on open chromatin tracks and histone chromatin marks across cell types and tissues in the human and cow genome. Moreover, we propose a hormonal control of differences in metabolism and sexual maturity between HFE and LFE animals, indicating potential biomarkers for further validation such as adrenomedullin, FSH, oxytocin, somatostatin and TSH.
Results
Multi-Tissue Transcriptomic Data Reveal Differences Between High and Low Feed Efficient Animals
Feed efficiency is a complex trait characterized by multiple distinct biological processes including metabolism, ingestion, digestion, physical activity and thermoregulation (; ). To study FE at transcriptional level we performed RNAseq of five tissues (i.e., adrenal gland, hypothalamus, liver, muscle and pituitary) from nine male bovines of high feed efficiency [HFE, characterized by low residual feed intake (RFI) ()] and nine of low FE (LFE, characterized by high RFI). In total, we analyzed 18 samples of liver, hypothalamus and pituitary; 17 of muscle and 15 of adrenal gland, yielding 13 million reads per sample on average (Supplementary Table 1). Gene expression was estimated for 24,616 genes present in the reference genome (UMD 3.1) and after quality control (refer to Section “Materials and Methods”), 17,354 genes were identified as being expressed in at least one of the five tissues analyzed.
Differential expression (DE) analysis between HFE and LFE animals resulted in 471 DE genes across tissues (P < 0.001, Supplementary Image 1), namely, 111 in adrenal gland, 125 in hypothalamus, 91 in liver, 104 in muscle and 98 in pituitary (Supplementary Tables 2A–E). Although no significant functional enrichment was found for the 281 genes up-regulated in HFE group, the 248 genes down-regulated presented a significant enrichment of GO terms such as response to hormone (Padj = 5.43 × 10-6), regulation of hormone levels (Padj = 3.48 × 10-6), cell communication (Padj = 3.18 × 10-4), regulation of signaling receptor activity (Padj = 3.20 × 10-4), hormone metabolic process (Padj = 5.86 × 10-4), response to corticosteroid (Padj = 6.28 × 10-4), regulation of secretion (Padj = 7.2 × 10-4), response to lipopolysaccharide (Padj = 7.9 × 10-4) and regulation of cell proliferation (Padj = 1.86 × 10-3). Refer to Supplementary Image 2 to see all enriched terms.
Overlap Between Gene Selection Criteria Prioritizes Genes Associated With Feed Efficiency
The genetic architecture behind complex traits involves a large variety of genes with coordinated expression patterns, which can be represented by gene regulatory networks as a blueprint to study their relationships and to identify central regulatory genes (). Therefore, it is important to select relevant genes and gene families according to the phenotype of interest to be used for network analysis. We defined five categories of genes (see Section “Materials and Methods” for further information) for inclusion in co-expression analysis: (1) differentially expressed (DE), (2) genes harboring SNPs previously associated with FE (harboring SNP), (3) tissue specific (TS), (4) genes coding proteins secreted in plasma by any of the five tissues analyzed (secreted) and (5) key regulators.
As reported before, we have identified 471 DE genes between HFE and LFE animals (Figure 1A and Supplementary Table 3A). In addition, 267 genes were selected for harboring SNPs previously associated with FE, as not only differences in expression levels can influence the phenotype but also polymorphism in the DNA sequence that can alter the translated protein behavior (Supplementary Table 3B). Moreover, 396 were selected for being tissue specific (refer to Section “Materials and Methods” for definition); 22 in adrenal gland, 32 in hypothalamus, 215 in liver, 118 in muscle and 9 in pituitary (Supplementary Table 3C). A total of 244 genes coding proteins secreted in plasma were selected because of their potential as biomarkers of FE (Supplementary Table 3D). From those, 135 had liver as the tissue of maximum expression and were functionally enriched for GO terms such as complement activation (Padj = 1.82 × 10-19), regulation of acute inflammatory response (Padj = 1.89 × 10-14), innate immune response (Padj = 9.71 × 10-12), negative regulation of endopeptidase activity (Padj = 2.35 × 10-10), platelet degranulation (Padj = 1.08 × 10-10), regulation of coagulation (Padj = 3.39 × 10-9), triglyceride homeostasis (Padj = 1.23 × 10-6), cholesterol efflux (Padj = 1.03 × 10-5) (Supplementary Image 3). Finally, from 1570 potential regulators in publicly available Animal TFdb, 78 were identified as key regulators of the genes selected by all the other categories, i.e., 78 genes presented a coordinated expression level with many of the genes in the network reflecting a tight control of expression patterns across tissues (Supplementary Table 3E).
FIGURE 1
Considering all the inclusion criteria, 1,335 genes were selected to be included in co-expression network analysis (Figure 1B and Supplementary Table 4), some of them selected in more than one category (Figure 1C). Regarding DE genes, six of them were also reported before as harboring SNPs associated with the phenotype (LUZP2, MAOB, SFRS5, SLC24A2, SOCS3 and WIF1) (; ; ) and 13 of them were key regulators (HOPX, PITX1, CRYM, PLCD1, ND6, CYTB, ND1, MT-ND4L, ND5, ATP8, ND4, ENSBTAG00000046711 and ENSBTAG00000048135). Many of the genes that are both DE and regulators are involved in respiratory chain (ND6, CYTB, ND1, MT-ND4L, ND5, ATP8 and ND4) and were all up-regulated in HFE group.
Considering both DE and secreted genes, 18 were identified (NOV, SPP1, CTGF, OXT, PTX3, VGF, CCL21, COL1A2, PGF, SOD3, SERPINE1, PRL, PON1, SST, JCHAIN, PCOLCE, IGFBP6 and SCG2). In addition, four genes were DE, secreted and tissue specific, two from liver (CXCL3 and IGFBP1) and two from pituitary (NPY and CYP17A1). Genes RARRES2 and PENK (proenkephalin) were DE, secreted and had been previously reported as harboring SNP associated with FE (RARRES2:AnimalQTLdb Release 35 – QTL:20671, rs133399845; PENK: - rs136198266, rs134428213, rs137492938, rs132881564). Other DE genes worthy to highlight, due to their well-known role in metabolic processes, are AMH (anti-mullerian hormone), TSHB (thyroid stimulating hormone beta), FGF21 (Fibroblast growth factor 21) and FST (follistatin), up-regulated in HFE group, and PMCH (pro-melanin concentrating hormone), ADM (adrenomedullin) and FSHB (follicle stimulating hormone beta), up-regulated in LFE group.
Co-expression Network Reveals Regulatory Genes and Biological Processes Related to Feed Efficiency
The co-expression network (Figure 2) was composed of 1,317 significant genes and 91,932 connections, with a mean of 70 connections per gene (considering only genes with significant expression correlation ≥|0.90|). Most of the connections (51%) involved a DE gene and 23% of those were between two DE genes. Tissue specific (TS) genes were involved in 49% of the connections with 119 connections per gene on average, which was higher than the overall network mean and reflects the close relationship among genes involved in tissue specific functions. Key regulators were the least represented category in the network (only 78 genes) but accounted for 11% of the connections in the network with the highest value of mean connections per gene, 131 connections, which is in accordance with their regulatory role. Regarding the connections within tissues, when we ranked all the genes in the network by the number of connections and looked at the top 50 genes, 29 were from liver, 15 were from muscle and 3, 2 and 1 were from pituitary, adrenal gland and hypothalamus, respectively. These results indicate very well-coordinated expression patterns in liver and muscle that could be a reflection of the number of TS genes in those tissues and the presence of central regulatory genes coordinating the expression of many other genes.
FIGURE 2
In the network (Figure 2), genes were grouped together by tissue which was mostly driven by TS genes. As mentioned before, most of the secreted protein-coding genes were located in the liver. Most of the key regulators were located peripherally in relation to the clusters which could be reflecting their regulatory nature independent of tissue specificity. Despite that, some regulators draw attention because of their high number of connections.
The top five most connected regulators were EPC1, NR2F6, MED21, ENSBTAG00000031687 and CTBP1, varying from 317 to 284 connections. They were all first neighbors of each other and were connected mainly to genes with higher expression in liver and essentially enriched for acute inflammatory response (Padj = 4.5 × 10-13, Supplementary Image 4). The next most connected regulator is TGFB1 with 217 connections. It is mainly connected to genes from muscle that are primarily enriched for muscle organ development (Padj = 6.87 × 10-5) and striated muscle contraction (Padj = 1.39 × 10-5, Supplementary Image 5). Besides indicating main regulator genes, the gene co-expression networks approach can be useful to access the role of specific genes. For instance, gene FGF21, a hormone up-regulated in liver of HFE animals, is directly connected to genes enriched for plasma lipoprotein particle remodeling, regulation of lipoprotein oxidation and cholesterol efflux (Padj = 5.64 × 10-3, Supplementary Image 6). Indeed, according to the literature, this gene is associated with decrease in body weight, blood triglycerides and LDL-cholesterol ().
Motif Discovery Confirms NR2F6 as a Key Regulator of Liver Transcriptional Changes Between High and Low Feed Efficiency
By means of the power-law theory, co-expression networks present many nodes with few connections and few central nodes with many connections (), being the last ones indicated as central regulatory genes responsible for the transcriptional changes between the divergent phenotypes analyzed. In our study, the most connected regulators were indicated, together with their target genes, i.e., their first neighbors in the network. Those genes are a mixture of direct and indirect regulator targets. In order to validate the regulatory role of the most connected regulators in the network and identify their core direct targets, we performed motif discovery in their co-expressed target genes. It is noteworthy that motif discovery should confirm the presence of DNA motifs of a TF in the regulatory regions of co-expressed genes. From the top five most connected regulators from our previous co-expression analysis, only NR2F6 has the ability to bind DNA. In contrast, the other four regulators act mainly as cofactors (corepressor, i.e., CTBP1; coactivator, i.e., MED21; or histones modifier, i.e., EPC1), that is co-binding through protein–protein interactions.
The analysis of 313 co-expressed genes with NR2F6 (Figure 3A) yield the Nuclear Factor motif HNF4-NR2F2 (transfac_pro-M01031) as the second motif most enriched out of 9732 PWMs (position weight matrices) with a Normalized Enrichment Score (NES) of 7.98 (Figure 3B). In addition, a total of 19 motifs associated with HNF4-NR2F2 were enriched in the dataset, associating HNF4-NR2F2 to 168 direct target genes (Figure 3C). Due to motif redundancy or highly similarity between a plethora of TFs, these motifs can be associated with multiple TFs from HNF4 (direct) to several nuclear factors such as NR2F6 (motif similarity score FDR 1.414 × 10-5). However, our co-expression analysis strongly indicates that NR2F6 is the key TF, since it was the TF with the highest number of nodes in the co-expression network (Figure 3C) and neither HNF4 nor NR2F2 were prioritized by any selection category to be included in the network.
FIGURE 3
Each of the NR2F6 inferred direct target genes contain one or more predicted enhancers, i.e., regions with high-scoring motif binding sites for NR2F6 or TFs with highly similar motifs. To validate the binding of these genomic regions by NR2F6 or TFs with highly similar motifs to NR2F6, we performed a region enrichment analysis of our predicted NR2F6 binding sequences against public TF ChiP-seq bound regions in human cell lines from the ENCODE consortium (1394 TF binding site tracks, ; ). This analysis confirms the experimental binding of TFs with similar binding as NR2F6 in HepG2 cells. In particular, HNF4A on human HepG2 (ENCFF001UGH, GSM803460, NES = 9.57), HNF4G (ENCFF001UGI, GSM803404, NES = 7.83), RXRA (ENCFF001UHJ, GSM803404, NES = 6.85), and NR2F2 (ENCFF001UGV, GSM1010810, NES = 4.45,) as the most enriched tracks (Supplementary Data Sheet 1). Recent NR2F6 ChIP-seq data on HepG2 (ENCODE experiment ENCSR518WPL, GSE96210) also confirms an enrichment for NR2F6, indicating predicted NR2F6 binding regions are experimentally bound by NR2F6 in hepatocyte cell lines (Figure 3D).
Next, to validate that the NR2F6 binding in those regions is functional in liver we performed an enrichment analysis for open-chromatin (tracks = 655) and histone modifications (tracks = 2450) related to active regulatory elements (Supplementary Table 5). This analysis yielded DNA-seq on human hepatocytes (ENCFF001SOV, GSM816663, NES = 4.10), and H3K29ac and H3K4me3 in adult liver (; GSM621630, GSM537709, respectively) as the most enriched tracks, respectively, strongly indicating that not only predicted target enhancers are bound by NR2F6 in Hepatocyte cell lines, but these regulatory regions are functionally active in hepatocytes and human liver (Figure 3D).
Regarding the cow genome, a recent open-chromatin study () has mapped active promoters and enhancers by H3K4me3 and H3K27ac ChIP-seq in cow liver resulting in 13,796 promoter and 45,786 enhancers. We performed an enrichment analysis of predicted NR2F6 enhancers converted to cow coordinates (n = 779, Supplementary Table 6, Array Express Accession number E-MTAB-2633) resulting in 446 regions being identified as functional regulatory regions in cow liver. This number is significantly higher compared to the only 43 regions expected to overlap by random (1000 permutation tests) (Figure 3E).
Finally, in addition to NR2F6 motif, HNF1A motif was found as a potential co-regulator in liver, in particular swissregulon-HNF1A.p2 with a NES = 10.17 and in total 20 enriched motifs and 170 direct targets were associated to HNF1A (Figure 3B). HNF1 is a master regulator of liver gene expression (), thus making its finding justified.
Motif Discovery Validates TGFB1 Signaling Through SMAD3/MYOD1 Binding as Drivers of Transcriptional Differences in Muscle of Divergent Feed Efficient Cattle
The analysis of the 217 genes co-expressed with TGFB1 (Figure 4A) showed that most target genes motifs were enriched for master regulators of muscle differentiation, namely, MEF2 (NES = 10.42), a MADS box Transcription factor with 148 target genes, and MYOD1 (NES = 5.09), a bHLH transcription factor (CANNTG) with 136 direct target genes (Figure 4B and Supplementary Data Sheet 2). To evaluate the precision of our predicted MYOD1 (bHLH) target genes, we assessed how many of these TF-TG relationships had been previously experimentally reported. Based on MYOD1 ChIP-seq binding in mouse myotubules, 86 genes had already been associated with MYOD1 resulting in a 63% success rate (hypergeometric test 1.72 × 10-22). SMAD3, the effector molecule of TGFB1 signaling, is known to recruit MYOD1 to drive transcriptional changes during muscle differentiation (). Thus, we evaluated whether predicted MYOD1 target genes were enriched for known SMAD3 target genes resulting in 21 out of 135 MYOD1 predicted target genes presenting SMAD3 ChIP-seq binding in myotubes, thus indicating that there is a statistically significant association between MYOD1 target genes and SMAD3 target genes in myotubes (hypergeometric test 1.98 × 10-6) (Figure 4C) (). By contrast, no significant association was found between predicted MYOD1 target genes in this study and SMAD3 target genes in other cell lines, such as pro-B and ES cell (hypergeometric test 0.056 and 0.076, respectively) (). That is in agreement with the fact that the effect of TGFB1 signaling driven by SMAD3 DNA binding is tissue-specific (). Our analysis predicted 621 potential MYOD1 binding sites, of which 114 (18%, Supplementary Tables 7, 8) and 152 (24.5%, Supplementary Table 9) present a MYOD1 ChIP-seq signal in mouse C2C12 myotubes cells () and in primary myotubes (), respectively.
FIGURE 4
Finally, we evaluate whether predicted MYOD1 binding regions were regulatory regions active in muscle cells across different species, namely human, mouse and cow. To tackle this issue we performed an enrichment analysis across 2113 open-chromatin ENCODE tracks (; ). This analysis resulted in a clear enrichment of our predicted MYOD1 binding regions with H3K27ac (NES = 15.98) and H3K9ac (NES = 8.78) regions in the skeletal muscle (Figure 4D). Both chromatin marks are associated with active transcription, H3K27ac related to active enhancers and H3K9ac related to active gene transcription (), thus validating most of our enhancer predictions that MYOD1 in human is active in the skeletal muscles. In cow, we assessed the overlap of predicted MYOD1 enhancers and promoter regions in cow muscle experimentally detected with H3K4me3 (). This resulted in 275 regions out of 653 (42%) overlap when only 11 regions are expected to overlap by random 1000 permutation test) (Figure 4E, Supplementary Table 10).
Differential Co-expression
Although the general co-expression network provides important insights about regulatory genes and their behavior, by creating specific networks for HFE and LFE and comparing the connectivity of the genes in each one, we can identify genes that change their behavior depending on the situation, moving from highly connected to lowly connected and vice-versa. We were able to identify 87 differentially connected genes between HFE and LFE (P < 0.05); 63 mainly expressed in liver, 19 in muscle and 3, 1 and 1 in hypothalamus, adrenal gland and pituitary, respectively (Supplementary Table 11). Those genes were enriched for terms such as regulation of blood coagulation (Padj = 3.14 × 10-10), fibrinolysis (Padj = 7.71 × 10-7), platelet degranulation (Padj = 7.49 × 10-6), regulation of peptidase activity (Padj = 6.16 × 10-4), antimicrobial humoral response (Padj = 2.49 × 10-3), acute inflammatory response (Padj = 2.18 × 10-4) and induction of bacterial agglutination (Padj = 3.58 × 10-2) (Supplementary Image 7). It is important to highlight that 20 of the differentially connected genes were also differentially expressed (Table 1) and three of them, i.e., SST, JCHAIN and IGFBP1, were secreted in plasma as well, which make them very promising potential biomarkers.
Table 1
| Gene name | Number of connections | Category∗ | Tissue of maximum expression | Tissue of differential expression | |
|---|---|---|---|---|---|
| Low feed efficiency | High feed efficiency | ||||
| SST | 0 | 45 | DE, SEC | Hypothalamus | Hypothalamus |
| SNORA73 | 41 | 108 | DE | Liver | Liver |
| ENSBTAG00000047700 | 56 | 111 | DE | Liver | Liver |
| ENSBTAG00000047121 | 62 | 111 | DE | Liver | Liver |
| ENSBTAG00000047816 | 53 | 96 | DE | Liver | Liver |
| ENSBTAG00000039928 | 50 | 89 | DE | Liver | Liver |
| ANXA13 | 115 | 63 | DE | Liver | Liver |
| FST | 113 | 56 | DE | Liver | Liver |
| PBLD | 115 | 55 | DE | Liver | Liver |
| ENSBTAG00000021368 | 95 | 0 | DE | Liver | Liver |
| JCHAIN | 52 | 113 | DE, SEC | Liver | Liver |
| IGFBP1 | 55 | 0 | DE, TS, SEC | Liver | Liver |
| SBK2 | 0 | 70 | DE | Muscle | Muscle |
| ACTC1 | 54 | 0 | DE | Muscle | Muscle |
| MYH1 | 0 | 47 | DE, TS | Muscle | Muscle |
| HR | 119 | 50 | DE | Pituitary | Muscle |
| TAGLN | 83 | 31 | DE | Adrenal | Muscle, Pituitary |
| SFRP2 | 41 | 91 | DE | Hypothalamus | Pituitary |
| FN1 | 119 | 69 | DE | Liver | Pituitary |
| CAV1 | 98 | 50 | DE | Muscle | Pituitary |
Differentially connected and differentially expressed genes between high and low feed efficiency.
∗Differentially expressed genes between high and low feed efficiency (DE), tissue specific genes (TS) and genes encoding proteins secreted in plasma (SEC).
Discussion
Feed efficiency is a complex trait, regulated by several biological processes. Thus, the identification of genomic regions associated with this phenotype, as well as regulators genes and biomarkers to select superior animals and to direct management decisions, is still a great challenge. In this work, multi-tissue transcriptomic data of high and low feed efficient Nellore bulls were analyzed through robust co-expression network methodologies in order to uncover some of the biology that governs these traits and put forward candidate genes to be the focus of further research. In this sense, the validation of target genes of main transcription factors (key regulators) in the network by motif search proves the efficacy of the methodology for network construction and prioritizes some transcription factors as central regulators (; ; ). Moreover, the addition of a category of genes coding proteins secreted in plasma in the co-expression analysis highlights the genes with potential to be explored as biomarkers of feed efficiency. We were able to identify genes related to main biological processes associated with feed efficiency and indicate key regulators.
Firstly, it is important to state that the 98 animals used to select the HFE and LFE groups in this study have been previously analyzed with regard to several phenotypic and molecular measures (; ; ). It was observed that HFE and LFE groups had similar body weight gain, carcass yield and loin eye area but LFE animals had higher feed intake, greater fat deposition, higher serum cholesterol levels, as well as hepatic inflammatory response, indicated by transcriptome analysis of liver biopsy and proved by the higher number of periportal mononuclear infiltrate (histopathology) and increased serum gamma-glutamyl-transferase (GGT, a biomarker of liver injury) in this group (). In the present study, the simultaneous analysis of five distinct tissues revealed the importance of hepatic tissue. Liver presented the most connected genes in the network, the largest number of differentially connected genes and the largest number of secreted genes, which, although can be explained by its biological function, are enriched mostly for terms related to lipid homeostasis and inflammatory response. Moreover, the top five most connected regulators in the network are co-expressed mainly with genes highly expressed in liver and also enriched for inflammatory response.
The relationship between FE and genes or pathways related to immune response and lipid metabolism is becoming more evident, as recent studies also reported in beef cattle (; ; ; ; ) and pigs (; ). In our previous work (), we proposed that increased liver lesions associated with higher inflammatory response in the liver of LFE animals could be due to increased lipogenesis and/or higher bacterial infection in the liver. While further evidence is needed to test these hypotheses, the enrichment of terms such as induction of bacterial agglutination and response to lipopolysaccharide makes bacterial infection a strong possibility. Indeed, pigs with low FE were reported to have a increased risk of intestinal inflammation, higher neutrophil infiltration biomarkers and increased serum endotoxin (lipopolysaccharide and other bacterial products) which could be related to increased bacterial infection or to decreased capacity to neutralize endotoxins (). The authors hypothesized that differences in bacterial population could partially explain the increase in circulating endotoxins, which could also be true for cattle given that differences in intestinal and ruminal bacterial population between high and low FE animals have already been reported (, ). Furthermore, the literature reports lipopolysaccharides (LPS) may cause up-regulation of adrenomedullin (ADM) hormone (), an up-regulated gene in LFE individuals as showed here. It was also demonstrated in rats that intravenous infusion of LPS caused up-regulation of ADM in ileum, liver, lung, aorta, skeletal muscle and blood vessels () whereas in our study, ADM presented differential expression in muscle, but not in liver.
Against pathogen invasion, a tightly regulated adaptive immune response must be triggered in order to allow T lymphocytes to produce cytokines or chemokines and B cells to differentiate and produce antibodies (). This regulation is known to be strongly influenced by the expression level and transcriptional activity of several nuclear receptors, including the NR2F-family, which consists of three orphan receptors: NR2F1, NR2F2 and NR2F6 (). Those receptors present highly conserved DNA and ligand binding domains among each other and across species (), and all three are expressed in adaptive and immune cells (). In our study, NR2F6 appeared as the second most connected regulator gene in the network while the other family members, although present in our expression data, were not selected by any of our inclusion criteria, thus indicating they might not be so relevant in our conditions. Indeed, NR2F6 appears to be a critical regulatory factor in the adaptive immune system by directly repressing the transcription of key cytokine genes in T effector cells (; ). The role of NR2F6 as a key regulator of inflammatory response in our network was validated at gene level by the identification of the binding motif HNF4-NR2F2 (transfac_pro-M01031) as one of the most enriched in NR2F6 target genes, due to the high similarity between NR2F2 and NR2F6 binding sites. Furthermore, using open chromatin data publicly available, we provided experimental evidence of the binding of TFs with highly similar binding motifs as NR2F6 in hepatocyte cells in humans and in cattle, thus, indicating that predicted target enhancers are functional in this tissue.
Another regulator prioritized in our analysis is TGFB1, the sixth most connected gene in the co-expression network, and a potential driver of transcriptional changes between high and low FE cattle in muscle. This gene has been previously described as a master regulator of FE in beef cattle, using genomics and metabolomics data (). Moreover, our motif discovery analysis showed that TGFB1 co-expressed genes are mostly enriched for binding sites of master regulators of muscle differentiation such as MEF2 and MYOD. Indeed, public available data show many of TGFB1 target genes were associated with MYOD (). As it is known, signaling pathways are an effective mechanism for cells to respond to environmental cues by regulating gene expression. TGFB1 signaling triggers the phosphorylation of SMAD2/3 transcription factors, which co-bind with cell-type master regulators at the nuclear level allowing/triggering/leading to cell-type specific transcriptional changes (; ). In skeletal muscle cells, myoblasts and myotubes, SMAD3 co-binds with MYOD1 (). The overlap between MYOD1 and SMAD3 target genes demonstrate the significant association between both genes in skeletal muscle, in agreement with the tissue-specific TGFB1 signaling response (). The overlap percentage between our predicted binding sites and MYOD1 Chip-seq data (18 and 24.5%) confirms previous analyses in mice where they reported only 20% of experimental validated distal enhancers in mouse myotubes with a bHLH (MyoD1) binding were actually bound by MYOD1 ChIP-seq data (). Thus, suggesting that additional transcription-factors and/or histone modifications have a key role in MYOD1 binding. The SMAD3/MYOD1 co-bound regions for known target genes are also captured, such as the promoter regions of ACTA1 and ANKRD1, both genes involved in skeletal muscle differentiation. We also demonstrated predicted MYOD1 binding regions are enriched for muscle regulatory regions across species (human, mouse and cow).
Altogether, we showed that co-expressed genes with TGFB1 are enriched for SMAD3/MYOD1 binding sites, which we validated at the gene and enhancer level by proving not only MYOD1 and SMAD3 binding, but also their accessibility, in human, mouse and cow. In pigs, increased feed efficiency is associated with stimulation of muscle growth by TGFB1 signaling pathway (). Finally, although not directly co-expressed with TGFB1, oxytocin (OXT) was DE in muscle and despite the lack of knowledge about its role in this tissue, previous work in cattle showed a massive increase of OXT expression in the muscle of bovines chronically exposed to anabolic steroids (). It is not known yet if oxytocin alone has an anabolic activity, but in a context where muscle growth seems to be associated with high FE animals, this hormone should be the focus of further investigation.
From the 13 regulator genes that are DE between groups, six are involved in respiratory chain and are up-regulated in HFE group. Genes ND1, ND4, ND4L, ND5, ND6 and also ND2 (which is DE but not identified as a regulator) are core subunits of the mitochondrial membrane respiratory chain Complex I (CI) which functions in the transfer of electrons from NADH to the respiratory chain, while ATP8 is part of Complex V and produces ATP from ADP in the presence of the proton gradient across the membrane. Interestingly, greater quantities of mitochondrial CI protein were associated with high FE cattle by whereas found higher CI-CII and CI-CIII concentration ratios for the same group. Other studies demonstrated that HFE animals consume less oxygen () and present lower plasma CO2 concentrations, which suggests a decreased oxidation process (). In general, the literature suggests mitochondrial ADP has greater control of oxidative phosphorylation in high FE individuals () and their increased mitochondrial function may contribute to feed efficiency (). In pigs, differences in mitochondrial function were reported when analyzing muscle (), blood () and adipose tissue transcriptomes (). Differences in metabolic rate associated with FE have long been discussed () and here the hypothesis is corroborated by the up-regulation of TSHB in HFE animals, which stimulates production of T3 and T4 in thyroid, thus increasing metabolism. Metabolism is inhibited by SST, a down-regulated hormone in this group which was also found to be differentially connected between HFE and LFE.
Examining the DE genes, many hormones can be identified. Hormones are signaling proteins that are transported by the circulatory system to target distant organs in order to regulate physiology. Regarding the relationship between FE and other production traits of economic importance, FSHB, responsible for spermatozoa production by activating Sertoli cells in the testicles (), is up-regulated in LFE group and is inhibited by follistatin (FST), a gene found to be down-regulated in the same group. Moreover, in rats, it has been demonstrated that FSH secretion is stimulated by somatostatin expression, which is up-regulated in LFE animals (). In this scenario, one could argue that selection for high FE delay reproduction traits, something that could be related to the lower fat deposition in this group, as previously observed (; ; ). Indeed, differences in body composition and in intermediary metabolism can impact on reproductive traits () and it has been observed before that feed efficient bulls present features of delayed sexual maturity, i.e., decreased progressive motility of the sperm and higher abundance of tail abnormalities (; ). Moreover, high FE heifers presented lower fat deposition and later sexual maturity which results in calving later in the calving season than their low FE counterparts (; ). LFE animals also exhibit down-regulation of AMH and the decrease of this hormone in serum is an excellent marker of Sertoli cells pubertal development ().
Concerning the differences in lipid metabolism in divergent FE phenotypes, FGF21, a hormone up-regulated in liver of HFE animals, is associated in humans to decrease in body weight, blood triglycerides and LDL-cholesterol, with improvement in insulin sensitivity (). It is an hepatokine released to the bloodstream and an important regulator of lipid and glucose metabolism (). When we performed an enrichment analysis of co-expressed genes with FGF21, we indeed found terms related to plasma lipoprotein particle remodeling, regulation of lipoprotein oxidation and cholesterol efflux mostly due to FGF21 co-expression with the apolipoproteins APOA4, APOC3 and APOM. In the same context, pro-melanin-concentrating hormone (PMCH) encodes three neuropeptides: neuropeptide-glycine-glutamic acid, neuropeptide-glutamic acid-isoleucine and melanin-concentrating hormone (MCH), the last one being the most extensively studied (). MCH up-regulation has been related to obesity and insulin resistance, as well as increased appetite and reduced metabolism in murine models (; ). The PMCH gene is up-regulated in LFE animals and harbors SNPs found to be associated with higher carcass fat levels and marbling score (; ).
In this work, we were able to identify several biological processes known to be related to feed efficiency, which together with the validation of the main transcription factors of the network, demonstrate the quality of the data and the robustness of the analyses, giving us the confidence to identify candidate genes as regulators or biomarkers of superior animals for this trait. The regulatory genes NR2F6 and TGFB1 play central roles in liver and muscle, respectively, by regulating genes related to inflammatory response and muscle development and growth, two main biological mechanisms associated to feed efficiency. Likewise, hormones and other proteins secreted in plasma as oxytocin, adrenomedulin, TSH, somatostatin, follistatin and AMH are interesting molecules to be explored as potential biomarkers of feed efficiency.
Materials and Methods
Phenotypic Data and Biological Sample Collection
All animal protocols were approved by the Institutional Animal Care and Use Committee of Faculty of Food Engineering and Animal Sciences, University of São Paulo (FZEA-USP – protocol number 14.1.636.74.1). All procedures to collect phenotypes and biological samples were carried out at FZEA-USP, Pirassununga, State of São Paulo, Brazil. Ninety-eight Nellore bulls (16 to 20 months old and 376 ± 29 kg BW) were evaluated in a feeding trial comprised of 21 days of adaptation to feedlot diet and place and a 70-day period of data collection. Total mixed ration was offered ad libitum and daily dry matter intake (DMI) was individually measured. Animals were weighed at the beginning, at the end and every 2 weeks during the experimental period. Feed efficiency was estimated by RFI which is the residual of the linear regression that estimates DMI based on average daily gain and mid-test metabolic body weight (). 40 animals selected either as high feed efficiency (HFE) or low feed efficiency (LFE) groups were slaughtered on 2 days with a 6-day interval. Adrenal gland (longitudinal section), hypothalamus, liver (lateral portion of the left lobe), skeletal muscle (medial portion of Longissimus lumborum, close to 12th rib) and pituitary samples were collected from each animal, rapidly frozen in liquid nitrogen and stored at –80°C. Further information about management and phenotypic measures of the animals used in this study can be found in .
RNAseq Data Generation
Samples of nine animals from each feed efficiency group (high and low) were selected for RNAseq using RFI measure. For hypothalamus and pituitary, the nitrogen frozen tissue was macerated with crucibles and pistils to ensure all portions of the tissue were represented, and stored in aliquots at –80°C. Then, RNA was extracted using AllPrep DNA/RNA/Protein Mini kit (QIAGEN, Crawley, United Kingdom). For liver, muscle and adrenal gland, a cut was made in the frozen tissue and the RNA was extracted using RNeasy Mini Kit (QIAGEN, Crawley, United Kingdom). RNA quality and quantity were assessed using automated capillary gel electrophoresis on a Bioanalyzer 2100 with RNA 6000 Nano Labchips according to the manufacturer’s instructions (Agilent Technologies Ireland, Dublin, Ireland). Samples that presented an RNA integrity number (RIN) of less than 8.0 were discarded.
RNA libraries were constructed using the TruSeqTM Stranded mRNA LT Sample Prep Protocol and sequenced on Illumina HiSeq 2500 equipment in a HiSeq Flow Cell v4 using the HiSeq SBS Kit V4 (2×100 pb). Liver, pituitary and hypothalamus were sequenced on the same run, each one in a different lane. Muscle and adrenal gland were sequenced in a second run, in different lanes.
Gene Expression Estimation
The quality of the sequencing was evaluated using the software FastQC Version 31. Sequence alignment against the bovine reference genome (UMD3.1) was performed using STAR Version 2.2.1 (), according to the standard parameters and including the annotation file (Ensembl release 89) and secondary alignments, duplicated reads and reads failing vendor quality checks were removed using Samtools Version 1.9 (). Then, HTseq Version 0.6.0 () was used to generate gene read counts and expression values were estimated by fragments per kilobase of gene per million mapped reads (FPKM). Genes with average value lower than 0.2 FPKM across all samples and tissues were discarded.
Gene expression normalization was performed using the following mixed effect model ():
where, the log2-transformed FPKM value for i-th library (86 levels), j-th gene (17,354 levels), k-th tissue (5 levels), l-th RFI phenotype (2 levels), corresponding to Yijkl, was modeled as a function of the fixed effect of library (Li) and the random effects of gene (Gj), gene by tissue (GTjk) and gene by RFI phenotype (GPjl). Random residual (eijkl) was assumed to be independent and identically distributed. Variance component estimates and solutions to the model were obtained using VCE6 (). Normalized mean expression (NME) values for each gene were defined as the linear combination of the solutions for random effects.
The mixed model used to normalize the expression data explained 96% of the variation in gene expression, of which the largest proportion (0.30) was due to tissue-specificity. Contrariwise, differences between HFE and LFE represented no variation (0.27 × 10-11). For that reason, normalized mean expression (NME) was only used to identify tissue specific genes and the raw FPKM values were used for differential expression and co-expression analysis.
Gene Selection for Network Construction
In order to select a set of relevant genes for network analysis, we defined five categories based on the following inclusion criteria:
Differential Expression (DE)
The mean expression value of each gene, for each group (HFE and LFE) and each tissue was calculated and then the expression of LFE group was subtracted from the expression in HFE group. Next, genes were ranked according to their mean expression in all samples for each tissue and divided into five bins. Genes were considered differentially expressed when the difference between the expression in HFE and LFE groups were greater than 3.1 or smaller than – 3.1 standard deviation from the mean in each bin, corresponding to a t-test P < 0.001 ().
Harboring SNPs
Genes harboring SNPs associated with feed efficiency, mainly indicated by GWAS, were identified using the PubMed database2 and the AnimalQTL database – Release 353 and only bovine data were considered regardless of breed.
Tissue Specific (TS)
A gene was considered as tissue specific when the average NME in that tissue was greater than one standard deviation from the mean of all genes and the average NME in all the other four tissues was smaller than zero.
Secreted
The human secretome database4 (; ) was used to select genes encoding proteins secreted in plasma by any of the analyzed tissues (adrenal gland, hypothalamus, liver, muscle and pituitary).
Key Regulators
In order to identify key regulatory genes to be included in the co-expression network, a list of genes were obtained from the Animal Transcription Factor Database 2.05 () and it was compared to a set of potential target genes in each tissue, composed of the categories: TS, DE, harboring SNPs and secreted. The analysis was based on regulatory impact factor metrics (), which comprises a set of two metrics designed to assign scores to regulator genes consistently differentially co-expressed with target genes and to those with the most altered ability to predict the abundance of target genes. Those scores deviating ± 1.96 standard deviation from the mean (corresponding to P < 0.05) were considered significant. Genes presenting mean expression value less than the mean of all genes expressed were not considered in this analysis.
Some of the genes selected by the categories above were represented by more than one Ensembl ID. Those duplications were removed for further analysis, keeping only the expression value of the most meaningful Ensembl ID. Additionally, genes with mean expression across the samples equal to zero were also removed from further analysis.
Co-expression Network Analysis
For gene network inference, genes selected using the five categories described previously were used as nodes and significant connections (edges) between them were identified using the Partial Correlation and Information Theory (PCIT) algorithm (), considering all animals and all tissues. PCIT determinates the significance of the correlation between two nodes after accounting for all the other nodes in the network. Connections between gene nodes were accepted when the partial correlation was greater than two standard deviations from the mean (P < 0.01). The output of PCIT was visualized on Cytoscape Version 3.6.1 ().
Network Validation Through Transcription Factor Biding Motifs Analysis
Using the regulatory impact factor metric (RIF) we prioritized key regulator genes from gene expression data and predicted target genes based on co-expression network. In order to assess whether those target genes were enriched for motifs associated with the top most connected regulators in the network with a DNA binding domain (transcription factors – TF), we performed motif discovery analysis in the set of co-expressed target genes (first neighbors of the TF) using the i-cistarget method () and i-Regulon v1.3, a Cytoscape plug-in (). These tools use humans (hg19) as the reference species, therefore only genes with human orthology are assessed. Then, to validate the binding of the identified genomic regions by the TFs, we performed a region enrichment analysis across experimentally available TF bound regions from ChiP-seq in cell lines from the ENCODE consortium (1,394 TF binding site tracks, ; ). Briefly, the tools evaluate whether there is an over-representation of motifs in the set of co-expressed genes and across evolution. We examined 10 kb upstream of the gene transcription start site and their conservation in 7 vertebrate species, including cow. Thus, the tools provide over-represented motifs across evolution, allowing us to predict regulatory interactions TF to target gene in cow. In our analysis, we performed motif discovery using i-Regulon v1.3 () and i-cistarget database 3 (), that is using their 9713 motif collection. Both methods result in highly similar enrichments. Whereas i-Regulon is a user-friendly method to deliver a regulatory network, i-cistarget also yields the genomic position of the TF binding in the genome. Both i-Regulon and i-cistarget can be used to validate the TF binding on predicted genomic regions in the human genome. The tool contains a collection of TF ChIP-seq data in cell lines mostly from the ENCODE consortium (1,394 TF binding tracks), 2003 Histone modifications from the ENCODE consortium and Epigenomics roadmap and 908 Histone modification and open-chromatin. The tool allows to perform an enrichment of the different human tracks at the region level.
Finally, we converted identified enhancer regions into cow coordinates and searched for regions of open-chromatin using data from publicly available studies in cow tissues. Namely, cow liver promoters and enhancer from (Array Express Accession number E-MTAB-2633) and skeletal muscle cow promoters from (GSE61936). For MYOD1 and SMAD3 binding in myotubes and pro-B cells, data from was used (GEO: GSE21621); and for MYOD1 binding in primary myotubes, data from (GEO: GSE20059) was used.
Differential Connectivity
In order to explore differentially connected genes between HFE and LFE, two networks were created, one for each condition, using the same methodology described before. Then, the number of connections of each gene in each condition was computed and scaled so that connectivity varied from 0 to 1, making it possible to compare the same gene in the two networks. The connectivity in LFE group was subtracted from the connectivity in HFE group and results deviating ± 1.96 standard deviation from the mean were considered significant (P < 0.05).
Functional Enrichment
Functional enrichment analysis was performed on the online platform GOrilla (Gene Ontology enRIchment anaLysis and visuaLizAtion tool6), using all genes that passed FPKM filter as background, hypergeometric test and multiple test correction (FDR – false discovery rate). The human database was used to take advantage of a more comprehensive knowledgebase regarding gene functions. GO terms were considered significant when Padj < 0.05. For genes in co-expression networks, visualized using Cytoscape (), the functional enrichment was performed with BiNGO plug-in () using the same background genes and statistical test.
Statements
Data availability statement
Datasets supporting the results of this article are public available in the European Nucleotide Archive (ENA) as part of FAANG consortium under de study ID PRJEB27337 and can be accessed following the link https://www.ebi.ac.uk/ena/data/view/PRJEB27337. Moreover, any additional information as well as the scripts used to perform the analysis are available upon request, please email pamela.alexandre@usp.br. This manuscript has been published as a preprint in bioRxiv, doi: https://doi.org/10.1101/360396.
Author contributions
PA performed formal analysis, investigation and visualization of all presented data and wrote the original draft. HF was the overall project leader responsible for conceptualization, project administration and supervision of PA. AR performed bioinformatics analysis and investigation and supervised PA. MN-S performed bioinformatics analysis and investigation. JF was responsible for funding acquisition and project administration. LP-N provided computing resources and conceptualization. All authors contributed to reviewing and editing of this manuscript’s final version.
Funding
This study and PA scholarships were funded by São Paulo Research Foundation (FAPESP) – Proc. 2014/07566-2, 2014/02493-7, 2015/22276-3 and 2017/14707-0. MN-S was funded by the CSIRO Science Excellence Research Office.
Acknowledgments
The authors thank Mrs. Elisângela C. M. Oliveira and Ms. Arina L. Rochetti for all technical support.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2019.00230/full#supplementary-material
Footnotes
1.^http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
2.^www.ncbi.nlm.nih.gov/pubmed/
3.^www.animalgenome.org/cgi-bin/QTLdb/index
4.^www.proteinatlas.org/humanproteome/secretome
References
1
AertsS.QuanX. J.ClaeysA.SanchezM. N.TateP.YanJ.et al (2010). Robust target gene discovery through transcriptome perturbations and genome-wide enhancer predictions in drosophila uncovers a regulatory basis for sensory specification.PLoS Biol.8:e1000435. 10.1371/journal.pbio.1000435.g001
2
AlexandreP. A.KogelmanL. J. A.SantanaM. H. A.PassarelliD.PulzL. H.Fantinato-NetoP.et al (2015). Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle.BMC Genomics16:1073. 10.1186/s12864-015-2292-8
3
AndersS.PylP. T.HuberW. (2014). HTSeq - A Python framework to work with high-throughput sequencing data.Bioinformatics31166–169. 10.1093/bioinformatics/btu638
4
ArcherJ. A.ArthurP. F.HerdR. M.ParnellP. F.PitchfordW. S. (1997). Optimum postweaning test for measurement of growth rate, feed intake, and feed efficiency in British breed cattle.J. Anim. Sci.752024–2032. 10.2527/1997.7582024x
5
ArcherJ. A.RichardsonE. C.HerdR. M.ArthurP. F. (1999). Potential for selection to improve efficiency of feed use in beef cattle: a review.Aust. J. Agric. Res.50147–161. 10.1071/A98075
6
BlumR.VethanthamV.BowmanC.RudnickiM.DynlachtB. D. (2012). Genome-wide identification of enhancers in skeletal muscle: the role of MyoD1.Genes Dev.262763–2779. 10.1101/gad.200113.112
7
BolormaaS.HayesB. J.SavinK.HawkenR.BarendseW.ArthurP. F.et al (2011). Genome-wide association studies for feedlot and growth traits in cattle.J. Anim. Sci.891684–1697. 10.2527/jas.2010-3079
8
CafeL. M.RobinsonD. L.FergusonD. M.GeesinkG. H.GreenwoodP. L. (2011). Temperament and hypothalamic-pituitary-adrenal axis function are related and combine to affect growth, efficiency, carcass, and meat quality traits in Brahman steers.Domest. Anim. Endocrinol.40230–240. 10.1016/j.domaniend.2011.01.005
9
CaoP.HanaiJ.-I.TanksaleP.ImamuraS.SukhatmeV. P.LeckerS. H. (2009). Statin-induced muscle damage and atrogin-1 induction is the result of a geranylgeranylation defect.FASEB J.232844–2854. 10.1096/fj.08-128843
10
CaoY.YaoZ.SarkarD.LawrenceM.SanchezG. J.ParkerM. H.et al (2010). Genome-wide MyoD binding in skeletal muscle cells: a potential for broad cellular reprogramming.Dev. Cell.18662–674. 10.1016/j.devcel.2010.02.014
11
ChavesA. S.NascimentoM. L.TullioR. R.RosaA. N.AlencarM. M.LannaD. P. (2015). Relationship of efficiency indices with performance, heart rate, oxygen consumption, blood parameters, and estimated heat production in Nellore steers.J. Anim. Sci.935036–5046. 10.2527/jas.2015-9066
12
ChenL.MaoF.CrewsD. H.VinskyM.LiC. (2014). Phenotypic and genetic relationships of feeding behavior with feed intake, growth performance, feed efficiency, and carcass merit traits in Angus and Charolais steers.J. Anim. Sci.92974–983. 10.2527/jas.2013-6926
13
CheungB. M.DengH. B. (2014). Fibroblast growth factor 21: a promising therapeutic target in obesity-related diseases.Exp. Rev. Cardiovasc. Ther.12659–666. 10.1586/14779072.2014.904745
14
ConnorE. E.KahlS.ElsasserT. H.ParkerJ. S.LiR. W.Van TassellC. P.et al (2010). Enhanced mitochondrial complex gene function and reduced liver size may mediate improved feed efficiency of beef cattle during compensatory growth.Funct. Integr. Genomics1039–51. 10.1007/s10142-009-0138-7
15
DavisM. P.BrooksM. A.KerleyM. S. (2016). Relationship between residual feed intake and lymphocyte mitochondrial complex protein concentration and ratio in crossbred steers.J. Anim. Sci.941587–1591. 10.2527/jas2015-9843
16
De JagerN.HudsonN. J.ReverterA.WangY.NagarajS. H.CafeL. M.et al (2011). Chronic exposure to anabolic steroids induces the muscle expression of oxytocin and a more than fiftyfold increase in circulating oxytocin in cattle.Physiol. Genomics43467–478. 10.1152/physiolgenomics.00226.2010
17
de la FuenteA. (2010). From “differential expression” to “differential networking” - identification of dysfunctional regulatory networks in diseases.Trends Genet.26326–333. 10.1016/j.tig.2010.05.001
18
DobinA.DavisC. A.SchlesingerF.DrenkowJ.ZaleskiC.JhaS.et al (2013). STAR: ultrafast universal RNA-seq aligner.Bioinformatics2915–21. 10.1093/bioinformatics/bts635
19
Encode Project Consortium (2012). An integrated encyclopedia of DNA elements in the human genome.Nature48957–74. 10.1038/nature11247
20
FontouraA. B. P.MontanholiY. R.Diel de AmorimM.FosterR. A.ChenierT.MillerS. P. (2016). Associations between feed efficiency, sexual maturity and fertility-related measures in young beef bulls.Animal1096–105. 10.1017/S1751731115001925
21
FooteA. P.TaitR. G.KeislerD. H.HalesK. E.FreetlyH. C. (2016). Leptin concentrations in finishing beef steers and heifers and their association with dry matter intake, average daily gain, feed efficiency, and body composition.Domest. Anim. Endocrinol.55136–141. 10.1016/j.domaniend.2015.12.007
22
FranciscoC. L.ResendeF. D.BenattiJ. M. B.CastilhosA. M.CookeR. F.JorgeA. M. (2015). Impacts of temperament on Nellore cattle: physiological responses, feedlot performance, and carcass characteristics.J. Anim. Sci.93:5419. 10.2527/jas.2015-9411
23
GerberP. J.SteinfeldH.HendersonB.MottetA.OpioC.DijkmanJ.et al (2013). Tackling Climate Change Through Livestock - A Global Assessment of Emissions and Mitigation Opportunities.Rome: Food & Agriculture Organization.
24
GiraltM.Gavaldà-NavarroA.VillarroyaF. (2015). Fibroblast growth factor-21, energy balance and obesity.Mol. Cell. Endocrinol.41866–73. 10.1016/j.mce.2015.09.018
25
GomesR. C.SainzR. D.SilvaS. L.CésarM. C.BoninM. N.LemeP. R. (2012). Feedlot performance, feed efficiency reranking, carcass traits, body composition, energy requirements, meat quality and calpain system activity in Nellore steers with low and high residual feed intake.Livest. Sci.150265–273. 10.1016/j.livsci.2012.09.012
26
GonanoC.MontanholiY.SchenkelF.SmithB.CantJ.MillerS. (2014). The relationship between feed efficiency and the circadian profile of blood plasma analytes measured in beef heifers at different physiological stages.Animal81684–1698. 10.1017/S1751731114001463
27
GondretF.VincentA.Houée-BigotM.SiegelA.LagarrigueS.CauseurD.et al (2017). A transcriptome multi-tissue analysis identifies biological pathways and genes associated with variations in feed efficiency of growing pigs.BMC Genomics18:244. 10.1186/s12864-017-3639-0
28
GroeneveldE.KovačM.MielenzN. (2010). VCE User’s Guide and Reference Manual, Version 6.0.Neustadt: Institute of Farm Animal Genetics.
29
HayesB. J.LewinH. A.GoddardM. E. (2013). The future of livestock breeding: genomic selection for efficiency, reduced emissions intensity, and adaptation.Trends Genet.29206–214. 10.1016/j.tig.2012.11.009
30
HelgesonS. C.SchmutzS. M. (2008). Genetic variation in the pro-melanin-concentrating hormone gene affects carcass traits in Bos taurus cattle.Anim. Genet.39310–315. 10.1111/j.1365-2052.2008.01717.x
31
HerdR. M.ArthurP. F. (2009). Physiological basis for residual feed intake.J. Anim. Sci.87E64–E71. 10.2527/jas.2008-1345
32
HerdR. M.OddyV. H.RichardsonE. C. (2004). Biological basis for variation in residual feed intake in beef cattle. 1. Review of potential mechanisms.Aust. J. Exp. Agric.44423–430. 10.1071/ea02220
33
Hermann-KleiterN.BaierG. (2014). Orphan nuclear receptor NR2F6 acts as an essential gatekeeper of Th17 CD4+ T cell effector functions.Cell Commun. Signal.121–12. 10.1186/1478-811X-12-38
34
Hermann-KleiterN.GruberT.Lutz-NicoladoniC.ThuilleN.FresserF.LabiV.et al (2008). The nuclear orphan receptor NR2F6 suppresses lymphocyte activation and T Helper 17-dependent autoimmunity.Immunity29205–216. 10.1016/j.immuni.2008.06.008
35
HerrmannC.Van De SandeB.PotierD.AertsS. (2012). i-cisTarget: an integrative genomics method for the prediction of regulatory features and cis-regulatory modules.Nucleic Acids Res.40:e114. 10.1093/nar/gks543
36
ItoM.GomoriA.IshiharaA.OdaZ.MashikoS.MatsushitaH.et al (2003). Characterization of MCH-mediated obesity in mice.Am. J. Physiol. Endocrinol. Metab.284E940–E945. 10.1152/ajpendo.00529.2002
37
JankyR.VerfaillieA.ImrichováH.van de SandeB.StandaertL.ChristiaensV.et al (2014). iRegulon: from a gene list to a gene regulatory network using large motif and track collections.PLoS Comput. Biol.10:e1003731. 10.1371/journal.pcbi.1003731
38
JingL.HouY.WuH.MiaoY.LiX.CaoJ.et al (2015). Transcriptome analysis of mRNA and miRNA in skeletal muscle indicates an important network for differential residual feed intake in pigs.Sci. Rep.5:11953. 10.1038/srep11953
39
KarisaB.MooreS.PlastowG. (2014). Analysis of biological networks and biological pathways associated with residual feed intake in beef cattle.Anim. Sci. J.85374–387. 10.1111/asj.12159
40
KellyA. K.McGeeM.CrewsD. H.FaheyA. G.WylieA. R.KennyD. A. (2010). Effect of divergence in residual feed intake on feeding behavior, blood metabolic variables, and body composition traits in growing beef heifers.J. Anim. Sci.88109–123. 10.2527/jas.2009-2196
41
KitaokaM.TakanoK.KojimaI.OgataE. (1989). A stimulatory effect of somatostatin: enhancement of activin A-mediated FSH secretion in rat pituitary cells.Biochem. Biophys. Res. Commun.162958–962. 10.1016/0006-291X(89)90766-3
42
KlepschV.Hermann-KleiterN.BaierG. (2016). Beyond CTLA-4 and PD-1: orphan nuclear receptor NR2F6 as T cell signaling switch and emerging target in cancer immunotherapy.Immunol. Lett.17831–36. 10.1016/j.imlet.2016.03.007
43
KochR. M.SwigerL. A.ChambersD.GregoryK. E. (1963). Efficiency of feed use in beef cattle.J. Anim. Sci.22486–494. 10.2527/jas1963.222486x
44
LancasterP. A.CarstensG. E.MichalJ. J.BrennanK. M.JohnsonK. A.DavisM. E. (2014). Relationships between residual feed intake and hepatic mitochondrial function in growing beef cattle.J. Anim. Sci.923134–3141. 10.2527/jas.2013-7409
45
LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al (2009). The sequence alignment/map format and SAMtools.Bioinformatics252078–2079. 10.1093/bioinformatics/btp352
46
LiuD.BlackB. L.DerynckR. (2001). TGF-beta inhibits muscle differentiation through functional repression of myogenic transcription factors by Smad3.Genes Dev.152950–2966. 10.1101/gad.925901
47
LiuH.NguyenY. T.NettletonD.DekkersJ. C. M.TuggleC. K. (2016). Post-weaning blood transcriptomic differences between Yorkshire pigs divergently selected for residual feed intake.BMC Genomics17:73. 10.1186/s12864-016-2395-x
48
LouveauI.VincentA.TacherS.GilbertH.GondretF. (2016). Increased expressions of genes and proteins involved in mitochondrial oxidation and antioxidant pathway in adipose tissue of pigs selected for a low residual feed intake.J. Anim. Sci.945042–5054. 10.2527/jas2016-0619
49
LudwigD. S.TritosN. A.MastaitisJ. W.KulkarniR.KokkotouE.ElmquistJ.et al (2001). Melanin-concentrating hormone overexpression in transgenic mice leads to obesity and insulin resistance.J. Clin. Invest.107379–386. 10.1172/JCI10660
50
MaderC. J.MontanholiY. R.WangY. J.MillerS. P.MandellI. B.McBrideB. W.et al (2009). Relationships among measures of growth performance and efficiency with carcass traits, visceral organ mass, and pancreatic digestive enzymes in feedlot cattle.J. Anim. Sci.871548–1557. 10.2527/jas.2008-0914
51
MaereS.HeymansK.KuiperM. (2005). BiNGO: a cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks.Bioinformatics213448–3449. 10.1093/bioinformatics/bti551
52
ManiV.HarrisA. J.KeatingA. F.WeberT. E.DekkersJ. C. M.GablerN. K. (2013). Intestinal integrity, endotoxin transport and detoxification in pigs divergently selected for residual feed intake.J. Anim. Sci.912141–2150. 10.2527/jas.2012-6053
53
MontanholiY. R.FontouraA. B. P.Diel de AmorimM.FosterR. A.ChenierT.MillerS. P. (2016). Seminal plasma protein concentrations vary with feed efficiency and fertility-related measures in young beef bulls.Reprod. Biol.16147–156. 10.1016/j.repbio.2016.04.002
54
MontanholiY. R.SwansonK. C.PalmeR.SchenkelF. S.McBrideB. W.LuD.et al (2010). Assessing feed efficiency in beef steers through feeding behavior, infrared thermography and glucocorticoids.Animal4692–701. 10.1017/S1751731109991522
55
MontanholiY. R.SwansonK. C.SchenkelF. S.McBrideB. W.CaldwellT. R.MillerS. P. (2009). On the determination of residual feed intake and associations of infrared thermography with efficiency and ultrasound traits in beef bulls.Livest. Sci.12522–30. 10.1016/j.livsci.2009.02.022
56
MotaL. F. M.BonaféC. M.AlexandreP. A.SantanaM. H.NovaisF. J.ToriyamaE.et al (2017). Circulating leptin and its muscle gene expression in Nellore cattle with divergent feed efficiency.J. Anim. Sci. Biotechnol.8:71. 10.1186/s40104-017-0203-3
57
MukiibiR.VinskyM.KeoghK. A.FitzsimmonsC.StothardP.WatersS. M.et al (2018). Transcriptome analyses reveal reduced hepatic lipid synthesis and accumulation in more feed efficient beef cattle.Sci. Rep.8:7303. 10.1038/s41598-018-25605-3
58
MullenA. C.OrlandoD. A.NewmanJ. J.LovénJ.KumarR. M.BilodeauS.et al (2011). Master transcription factors determine cell-type-specific responses to TGF-β signaling.Cell147565–576. 10.1016/j.cell.2011.08.050
59
MyerP. R.SmithT. P. L.WellsJ. E.KuehnL. A.FreetlyH. C. (2015). Rumen microbiome from steers differing in feed efficiency.PLoS One10:e0129174. 10.1371/journal.pone.0129174
60
MyerP. R.WellsJ. E.SmithT. P. L.KuehnL. A.FreetlyH. C. (2016). Microbial community profiles of the jejunum from steers differing in feed efficiency.J. Anim. Sci.94:327. 10.2527/jas.2015-9839
61
Naval-SańchezM.PotierD.HaagenL.SánchezM.MunckS.Van De SandeB.et al (2013). Comparative motif discovery combined with comparative transcriptomics yields accurate targetome and enhancer predictions.Genome Res.2374–88. 10.1101/gr.140426.112
62
NkrumahJ. D.KeislerD. H.CrewsD. H.BasarabJ. A.WangZ.LiC.et al (2007). Genetic and phenotypic relationships of serum leptin concentration with performance, efficiency of gain, and carcass merit of feedlot cattle.J. Anim. Sci.852147–2155. 10.2527/jas.2006-764
63
NovaisF. J.PiresP. R. L.AlexandreP. A.DrommsR. A.IglesiasA. H.FerrazJ. B. S.et al (2019). Identification of a metabolomic signature associated with feed efficiency in beef cattle.BMC Genomics20:8. 10.1186/s12864-018-5406-2
64
OliveiraP. S. N.CesarA. S. M.NascimentoM. L.ChavesA. S.TiziotoP. C.TullioR. R.et al (2014). Identification of genomic regions associated with feed efficiency in Nelore cattle.BMC Genet.15:100. 10.1186/s12863-014-0100-0
65
ParadisF.YueS.GrantJ. R.StothardP.BasarabJ. A.FitzsimmonsC. (2015). Transcriptomic analysis by RNA sequencing reveals that hepatic interferon-induced genes may be associated with feed efficiency in beef heifers.J. Anim. Sci.93:3331. 10.2527/jas.2015-8975
66
PereiraF. A.TsaiM. J.TsaiS. Y. (2000). COUP-TF orphan nuclear receptors in development and differentiation.Cell Mol. Life Sci.571388–1398. 10.1007/PL00000624
67
PotierD.DavieK.HulselmansG.NavalSanchezM.HaagenL.Huynh-ThuV. A.et al (2014). Mapping gene regulatory networks in drosophila eye development by large-scale transcriptome perturbations and motif inference.Cell Rep.92290–2303. 10.1016/j.celrep.2014.11.038
68
Ramayo-CaldasY.BallesterM.SánchezJ. P.González-RodríguezO.RevillaM.ReyerH.et al (2018). Integrative approach using liver and duodenum RNA-Seq data identifies candidate genes and pathways associated with feed efficiency in pigs.Sci. Rep.8:558. 10.1038/s41598-017-19072-5
69
RamosM. H.KerleyM. S. (2013). Mitochondrial complex I protein differs among residual feed intake phenotype in beef cattle.J. Anim. Sci.913299–3304. 10.2527/jas.2012-5589
70
RandelR. D.WelshT. H. (2013). Joint Alpharma-Beef species symposium: interactions of feed efficiency with beef heifer reproductive development.J. Anim. Sci.911323–1328. 10.2527/jas.2012-5679
71
ReverterA.BarrisW.McWilliamS.ByrneK. A.WangY. H.TanS. H.et al (2005). Validation of alternative methods of data normalization in gene co-expression studies.Bioinformatics211112–1120. 10.1093/bioinformatics/bti124
72
ReverterA.ChanE. K. F. (2008). Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks.Bioinformatics242491–2497. 10.1093/bioinformatics/btn482
73
ReverterA.HudsonN. J.NagarajS. H.Pérez-EncisoM.DalrympleB. P. (2010). Regulatory impact factors: unraveling the transcriptional regulation of complex traits from expression data.Bioinformatics26896–904. 10.1093/bioinformatics/btq051
74
ReyR.Lordereau-RichardI.CarelJ. C.BarbetP.CateR. L.RogerM.et al (1993). Anti-müllerian hormone and testosterone serum levels are inversely during normal and precocious pubertal development.J. Clin. Endocrinol. Metab.771220–1226. 10.1210/jcem.77.5.8077315
75
Roadmap Epigenomics ConsortiumKundajeA.MeulemanW.ErnstJ.BilenkyM.YenA.et al (2015). Integrative analysis of 111 reference human epigenomes.Nature518317–329. 10.1038/nature14248
76
RolfM. M.TaylorJ. F.SchnabelR. D.McKayS. D.McClureM. C.NorthcuttS. L.et al (2011). Genome-wide association analysis for feed efficiency in Angus cattle.Anim. Genet.43367–374. 10.1111/j.1365-2052.2011.02273.x
77
SaatchiM.BeeverJ. E.DeckerJ. E.FaulknerD. B.FreetlyH. C.HansenS. L.et al (2014). QTLs associated with dry matter intake, metabolic mid-test weight, growth and feed efficiency have little overlap across 4 beef cattle studies.BMC Genomics15:1004. 10.1186/1471-2164-15-1004
78
SantanaM. H. A.RossiP.AlmeidaR.CuccoD. C. (2012). Feed efficiency and its correlations with carcass traits measured by ultrasound in Nellore bulls.Livest. Sci.145252–257. 10.1016/j.livsci.2012.02.012
79
SantanaM. H. A.UtsunomiyaY. T.NevesH. H. R.GomesR. C.GarciaJ. F.FukumasuH.et al (2014). Genome-wide association analysis of feed intake and residual feed intake in Nellore cattle.BMC Genet.15:21. 10.1186/1471-2156-15-21
80
SchmiererB.HillC. S. (2007). TGFβ-SMAD signal transduction: molecular specificity and functional flexibility.Nat. Rev. Mol. Cell Biol.8970–982. 10.1038/nrm2297
81
SeaburyC. M.OldeschulteD. L.SaatchiM.BeeverJ. E.DeckerJ. E.HalleyY. A.et al (2017). Genome-wide association study for feed efficiency and growth traits in U.S. beef cattle.BMC Genomics18:386. 10.1186/s12864-017-3754-y
82
ShafferK. S.TurkP.WagnerW. R.FeltonE. E. D. (2011). Residual feed intake, body composition, and fertility in yearling beef heifers.J. Anim. Sci.891028–1034. 10.2527/jas.2010-3322
83
ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD.et al (2003). Cytoscape: a software Environment for integrated models of biomolecular interaction networks.Genome Res.132498–2504. 10.1101/gr.1239303
84
ShinJ. H.LiR. W.GaoY.BaldwinR. VILiC. J. (2012). Genome-wide ChIP-seq mapping and analysis reveal butyrate-induced acetylation of H3K9 and H3K27 correlated with transcription activity in bovine cells.Funct. Integr. Genomics12119–130. 10.1007/s10142-012-0263-6
85
ShindoT.KuriharaH.KuriharaY.MoritaH.YazakiY. (1998). Upregulation of endothelin-1 and adrenomedullin gene expression in the mouse endotoxin shock model.J. Cardiovasc. Pharmacol.31(Suppl. 1)S541–S544. 10.1097/00005344-199800001-00156
86
ShojiH.MinaminoN.KangawaK.MatsuoH. (1995). Endotoxin markedly elevates plasma concentration and gene transcription of adrenomedullin in rat.Biochem. Biophys. Res. Commun.215531–537. 10.1006/bbrc.1995.2497
87
SloanC. A.ChanE. T.DavidsonJ. M.MalladiV. S.StrattanJ. S.HitzB. C.et al (2016). ENCODE data at the ENCODE portal.Nucleic Acids Res.44D726–D732. 10.1093/nar/gkv1160
88
SwamiM. (2009). Networking complex traits.Nat. Rev. Genet.10:2566. 10.1038/ng.332
89
TroncheF.YanivM. (1992). HNF1, a homeoprotein member of the hepatic transcription regulatory network.BioEssays14579–587. 10.1002/bies.950140902
90
UhlénM.BjörlingE.AgatonC.SzigyartoC. A.-K.AminiB.AndersenE.et al (2005). A human protein atlas for normal and cancer tissues based on antibody proteomics.Mol. Cell. Proteomics41920–1932. 10.1074/mcp.M500279-MCP200
91
UhlenM.FagerbergL.HallstromB. M.LindskogC.OksvoldP.MardinogluA.et al (2015). Tissue-based map of the human proteome.Science3471260419–1260419. 10.1126/science.1260419
92
VillarD.BerthelotC.AldridgeS.RaynerT. F.LukkM.PignatelliM.et al (2015). Enhancer evolution across 20 mammalian species.Cell160554–566. 10.1016/j.cell.2015.01.006
93
VincentA.LouveauI.GondretF.TréfeuC.GilbertH.LefaucheurL. (2015). Divergent selection for residual feed intake affects the transcriptomic and proteomic profiles of pig skeletal muscle.J. Anim. Sci.932745–2758. 10.2527/jas.2015-8928
94
WalkerW. H.ChengJ. (2005). FSH and testosterone signaling in Sertoli cells.Reproduction13015–28. 10.1530/rep.1.00358
95
WalterL. J.GaschC. A.MceversT. J.HutchesonJ. P.DefoorP.MarquessF. L. S.et al (2014). Association of pro-melanin concentrating hormone genotype with beef carcass quality and yield.J. Anim. Sci.92325–331. 10.2527/jas2013-6931
96
WeberK. L.WellyB. T.Van EenennaamA. L.YoungA. E.Port-NetoL. R.ReverterA.et al (2016). Identification of Gene networks for residual feed intake in Angus cattle using genomic prediction and RNA-seq.PLoS One11:e0152274. 10.1371/journal.pone.0152274
97
WidmannP.ReverterA.WeikardR.SuhreK.HammonH. M.AlbrechtE.et al (2015). Systems biology analysis merging phenotype, metabolomic and genomic data identifies Non-SMC condensin I complex, subunit G (NCAPG) and cellular maintenance processes as major contributors to genetic variability in bovine feed efficiency.PLoS One10:e0124574. 10.1371/journal.pone.0124574
98
ZarekC. M.Lindholm-PerryA. K.KuehnL. A.FreetlyH. C. (2017). Differential expression of genes related to gain and intake in the liver of beef cattle.BMC Res. Notes10:1. 10.1186/s13104-016-2345-3
99
ZhangH. M.LiuT.LiuC. J.SongS.ZhangX.LiuW.et al (2015). AnimalTFDB 2.0: a resource for expression, prediction and functional study of animal transcription factors.Nucleic Acids Res.43D76–D81. 10.1093/nar/gku887
100
ZhaoC.CarrilloJ. A.TianF.ZanL.UpdikeS. M.ZhaoK.et al (2015). Genome-wide H3K4me3 analysis in angus cattle with divergent tenderness.PLoS One10:e0115358. 10.1371/journal.pone.0115358
Summary
Keywords
feed efficiency, residual feed intake, Nellore (Zebu), Bos indicus, inflammation, muscle development, motif discovery, regulatory gene network
Citation
Alexandre PA, Naval-Sanchez M, Porto-Neto LR, Ferraz JBS, Reverter A and Fukumasu H (2019) Systems Biology Reveals NR2F6 and TGFB1 as Key Regulators of Feed Efficiency in Beef Cattle. Front. Genet. 10:230. doi: 10.3389/fgene.2019.00230
Received
09 August 2018
Accepted
04 March 2019
Published
22 March 2019
Volume
10 - 2019
Edited by
David E. MacHugh, University College Dublin, Ireland
Reviewed by
Elisabetta Giuffra, INRA Centre Jouy-en-Josas, France; Carolina Neves Correia, University College Dublin, Ireland
Updates
Copyright
© 2019 Alexandre, Naval-Sanchez, Porto-Neto, Ferraz, Reverter and Fukumasu.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Heidge Fukumasu, fukumasu@usp.br
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.