Post-weaning blood transcriptomic differences between Yorkshire pigs divergently selected for residual feed intake

Abstract

Background: Improving feed efficiency (FE) of pigs by genetic selection is of economic and environmental significance. An increasingly accepted measure of feed efficiency is residual feed intake (RFI). Currently, the molecular mechanisms underlying RFI are largely unknown. Additionally, to incorporate RFI into animal breeding programs, feed intake must be recorded on individual pigs, which is costly and time-consuming. Thus, convenient and predictive biomarkers for RFI that can be measured at an early age are greatly desired. In this study, we aimed to explore whether differences exist in the global gene expression profiles of peripheral blood of 35 to 42 day-old pigs with extremely low (more efficient) and high RFI (less efficient) values from two lines that were divergently selected for RFI during the grow-finish phase, to use such information to explore the potential molecular basis of RFI differences, and to initiate development of predictive biomarkers for RFI.

Results: We identified 1972 differentially expressed genes (DEGs) (q ≤ 0.15) between the low (n = 15) and high (n = 16) RFI groups of animals by using RNA sequencing technology. We validated 24 of 37 selected DEGs by reverse transcription-quantitative PCR (RT-qPCR) in a joint analysis of 24 (12 per line) of the 31 samples already used for RNA-seq plus 24 (12 per line) novel samples from the same contemporary group of pigs. Using an analysis of the 24 novel samples alone, only nine of the 37 selected DEGs were validated. Genes involved in small molecule biosynthetic process, antigen processing and presentation of peptide antigen via major histocompatibility complex (MHC) class I, and steroid biosynthetic process were overrepresented among DEGs that had higher expression in the low versus high RFI animals. Genes known to function in the proteasome complex or mitochondrion were also significantly enriched among genes with higher expression in the low versus high RFI animals. Alternatively, genes involved in signal transduction, bone mineralization and regulation of phosphorylation were overrepresented among DEGs with lower expression in the low versus high RFI animals. The DEGs significantly overlapped with genes associated with disease, including hyperphagia, eating disorders and mitochondrial diseases (q < 1E-05). A weighted gene co-expression network analysis (WGCNA) identified four co-expression modules that were differentially expressed between the low and high RFI groups. Genes involved in lipid metabolism, regulation of bone mineralization, cellular immunity and response to stimulus were overrepresented within the two modules that were most significantly differentially expressed between the low and high RFI groups. We also found five of the DEGs and one of the co-expression modules were significantly associated with the RFI phenotype of individual animals (q < 0.05).

Conclusions: The post-weaning blood transcriptome was clearly different between the low and high RFI groups. The identified DEGs suggested potential differences in mitochondrial and proteasomal activities, small molecule biosynthetic process, and signal transduction between the two RFI groups and provided potential new insights into the molecular basis of RFI in pigs, although the observed relationship between the post-weaning blood gene expression and RFI phenotype measured during the grow-finish phase was not strong. DEGs and representative genes in co-expression modules that were associated with RFI phenotype provide a preliminary list for developing predictive biomarkers for RFI in pigs.

Fig. 1

Distribution of RFI values and sample selection for RNA-seq and RT-qPCR assays. Post-weaning blood samples of 16 pigs (8 per diet) with extremely low RFI from gilts of the low RFI line (LRFI) and 16 pigs (8 per diet) with extremely high RFI from gilts of the high RFI line (HRFI) fed the high-fiber, low-energy diet (HFD) or the low-fiber, high-energy diet (LFD) were selected for RNA-seq. One sample in the LRFI-LFD group was excluded from RNA sequencing because of low quality RNA, leaving a total of 31 samples for RNA-seq. The green and red dots represent individuals selected for RNA-seq. Twelve samples from each line by diet combination were selected for RT-qPCR validation of DEGs such that the corresponding RFI phenotypes were representative. The red and blue dots represent the samples selected for RT-qPCR assays. The 24 novel blood samples were selected such that the RFI values of the corresponding animals were roughly evenly distributed across the ranges of RFI not covered by the RFI phenotypes of the 24 animals originally selected for RNA-seq. The distribution of RFI values of barrows from each line by diet combination is also shown for reference

Fig. 2

Differentially expressed genes and transcriptomic differences between the low and high RFI groups. a Heatmap showing 454 DEGs (q ≤ 0.05) between low and high RFI groups identified by RNA-seq. Sample names are designated as RFI line followed by the pig identifier. LRFI, low RFI line; HRFI, high RFI line. Animals with sample names in blue were fed the high-fiber, low-energy diet (LFD), while animals with sample names in black were fed the low-fiber, high-energy diet (HFD). The relative orders of genes and samples were determined by two-way hierarchical clustering based on the adjusted transformed gene expression of the 454 DEGs. The adjusted gene expression was gene-wise standardized to get the z-score as displayed. b Volcano plot showing the magnitude and significance of differential expression of genes between low and high RFI groups. Black vertical dash lines correspond to │log₂(fold change)│ = 1, while red horizontal dash line correspond to q-value of 0.15. FC, fold change calculated as the ratio of mean gene expression in the low RFI group to mean gene expression in the high RFI group after accounting for the other relevant variables. c Hierarchical clustering showing relationship of the 31 RNA-seq samples. The samples were hierarchically clustered by using Ward method with 1 minus Spearman correlation as distance. The Spearman correlations between pairs of samples were calculated based on the adjusted transformed expression of the 12280 genes

Fig. 3

Validation of DEGs by RT-qPCR. 24 of the 37 selected DEGs between low and high RFI groups were confirmed by RT-qPCR when using the 48 samples, 24 of which were used for RNA-seq and another 24 of which were novel, as shown in this figure (q ≤ 0.15); 22 of the 37 selected DEGs were confirmed by RT-qPCR when using the 24 samples that were used for RNA-seq (q ≤ 0.15); but only 9 of the 37 selected DEGs were validated by RT-qPCR when using the 24 novel samples (q ≤ 0.15) (Addition file 9: Table S7). For comparison, the log₂ (fold change) of these genes determined by RNA-seq were also displayed. Genes were ordered based on their log₂ (fold change) as determined by RT-qPCR for display. Error bars for RT-qPCR assays show the standard errors of mean log₂ (fold change). DEGs not confirmed by RT-qPCR are labeled in red. Genes without corresponding human orthologs are labeled with the last 5 digits of their Ensembl gene IDs, with common prefix “ENSSSCG000000” omitted for simplicity. For example, the ID for “29500” should be ENSSSCG00000029500

Fig. 4

Co-expression modules differentially expressed between the low and high RFI group. a Distribution of eigengenes of modules highly associated with RFI groups. b Venn diagram showing overlapping between module genes and DEGs