BACKGROUND
feed efficiency  is one of the important traits directly related to profitability, productivity, and sustainability in the pork industry  <cit> .
while many non-genetic strategies have been developed to improve feed efficiency of pigs  <cit> , improving feed efficiency by genetic selection can be a sustainable alternative.
an increasingly popular measure of fe is residual feed intake , which is the difference between the actual and expected feed intake of an animal for production and maintenance  <cit> .
importantly, rfi is moderately heritable and responds well to genetic selection in pigs  <cit> .
over the last decade, significant efforts have been made to understand the molecular, genetic and physiological basis of rfi in pigs.
researchers have found many interesting differences across multiple tissues between pigs with divergent rfi phenotypes  <cit> .
for example, as compared to high rfi pigs , low rfi pigs  have altered feeding behaviors  <cit> , slightly lower growth rate  <cit> , less back fat  <cit> , lower protein turnover rate in the muscle  <cit> , altered mitochondrial protein profiles  <cit> , less mitochondrial reactive oxygen species  production  <cit> , and lower levels of leptin  <cit>  and juvenile igf- <dig> in circulating blood  <cit> .
to explore the genetic basis of rfi in pigs, several genome-wide association studies  have been conducted .
some chromosomal regions tagged by single nucleotide polymorphisms  have been found to be associated with rfi, but these associations were not consistent across studies and explained only small portions of the genetic variance for rfi .
the wide range of differences across multiple tissues between lines of pigs with divergent rfi and the lack of snps with major effects on rfi suggest that rfi is a highly polygenic, quantitative trait with multiple tissues contributing to its variation.
global gene expression profiling technologies have also been used to explore the molecular basis of rfi in pigs.
by profiling the transcriptomes of the adipose tissue of two lines of pigs divergently selected for rfi with gene expression microarrays, lkhagvadorj et al.
<cit>  found that genes involved in the lipid metabolic pathway were overrepresented among the differentially expressed genes  that had lower expression in low versus high rfi pigs, and genes involved in carbohydrate metabolism and response to stress were overrepresented among the degs that had higher expression in low versus high rfi pigs  <cit> .
they also found the leptin-related gene network to be different between the two lines  <cit> .
recently, vincent et al.
<cit>  and jing et al.
<cit>  profiled the transcriptome of the longissimus dorsi  muscle from pigs with divergent rfi by gene expression microarray and rna-seq, respectively.
using pigs from lines divergently selected for rfi, vincent et al.
<cit>  found genes involved in protein synthesis and glycolysis, and genes associated with mitochondrial energy/oxidative metabolism had higher and lower expression, respectively, in the low versus high rfi line.
using yorkshire barrows with extreme phenotypes for rfi, jing et al.
<cit>  found that genes involved in glycolysis had lower expression in the low versus high rfi group, while genes involved in muscle proliferation and differentiation had higher expression in the low versus high rfi group.
surprisingly, these two studies shared no degs and proposed opposite differences in glycolytic activities in the low versus high rfi pigs.
therefore, in consideration of the complexity of rfi and the inconsistency from study to study, the molecular mechanisms underlying rfi in pigs are still largely unclear.
to incorporate rfi into animal breeding programs for improving feed efficiency, feed intake, body weight gain and back fat depth must be recorded on individual pigs.
as it is very expensive and time-consuming to record feed intake on individual animals  <cit> , convenient and predictive biomarkers for rfi that can be measured at an early age are in demand.
in cattle, chen et al.
<cit>  successfully used  <dig> degs identified in the liver of angus bulls that were divergently selected for rfi to classify angus steers from the same divergent rfi lines, and al-husseini et al.
<cit>  developed a rfi predictor using  <dig> of these  <dig> degs and validated it in an independent angus population.
these biomarkers are, however, not very practical because invasive liver biopsies are needed.
the peripheral blood is an informative tissue not only because it carries a variety of cells directly involved in immunity and inflammation, but also because it interacts with every organ and tissue in the body via bioactive circulating factors, such as nutrients, metabolites, cytokines, hormones and exosomal cargoes, which are released from the same or different organs or tissues  <cit> .
these bioactive factors interact with blood cells and thus might modify the gene expression profiles of the blood cells dynamically.
molecular signatures in circulating blood, including gene expression profiles, have been shown to reflect the physiopathological status, growth stage and lifestyle of subjects .
due to its easy accessibility and informativeness, blood has become a popular sample  for disease diagnosis, prediction, prognosis, and biomarker discovery  <cit> .
interestingly, the concentration of igf- <dig> in serum at a young age has been shown to be different between animals with divergent rfi in poultry and livestock, including pigs  <cit> .
blood cell profiles at early growth stages have also been found to be different between livestock with divergent rfi phenotypes, including pigs  <cit>  and cattle  <cit> .
these results suggest that animals with divergent rfi phenotypes, measured later in life, have early physiological differences in circulating blood that may be reflected in blood gene expression profiles at these early stages.
in this study, our objective was to determine the blood transcriptomic differences between post-weaning pigs from two lines divergently selected for rfi, to explore potential molecular mechanisms underlying rfi in peripheral blood and to develop a list of candidate biomarkers for rfi prediction.
we hypothesized that post-weaning expression levels of some genes in whole blood were correlated with rfi phenotype measured during the grow-finish phase.
we identified  <dig> degs with q ≤  <dig>  and four co-expression modules that were differentially expressed between the low and high rfi groups.
a set of selected degs were validated by reverse transcription-quantitative pcr .
several interesting biological processes underlying degs and differential co-expression modules were suggested.
we also identified several candidate biomarkers for rfi.
methods
animals, blood sample collection and complete blood count  test
the experimental protocols for this study were approved by the institutional animal care and use committee  at iowa state university under permit number 11-1-4996-s. all pigs were from parity  <dig> of generation nine of the two lines divergently selected for residual feed intake: the low and high rfi lines  <cit> .
bunter et al.
found that the igf- <dig> concentration in blood measured between  <dig> to 42 days of age differed between the two lines and were genetically correlated with rfi measured later during the grow-finish phase  <cit> .
so blood samples were collected from the jugular vein into tempus™ blood rna tubes  for long-term storage at −80 °c from  <dig> post-weaning piglets of the two lines in this age range.
meanwhile, blood from those animals was also collected into edta tubes  and kept at 4 °c before cbc tests, which were performed on the bleeding day, if possible, or the next morning as described  <cit> .
differences in cbc profiles between the two lines have been published elsewhere  <cit> .
at  <dig>  ±  <dig>   days of age and  <dig>  ±  <dig>  kg of body weight ,  <dig> gilts and  <dig> barrows were randomly assigned to  <dig> mixed-line, mixed-sex finishing pens with electronic single-space feeders  for feed efficiency test, with  <dig> pens being randomly assigned to either of two diets: a high-fiber, low-energy diet  and a low-fiber, high-energy diet   <cit> .
the animals on feed efficiency test consisted of  <dig> barrows and  <dig> gilts from the low rfi line plus  <dig> barrows and  <dig> gilts from the high rfi line fed the hfd, and  <dig> barrows and  <dig> gilts from the low rfi line plus  <dig> barrows and  <dig> gilts from the high rfi line fed the lfd.
all pigs had ad libitum access to feed and water.
individual feed intake was real-time recorded, body weight was recorded biweekly, and backfat  depth above the 10th rib and loin muscle area were recorded at the end of the test, when the pigs were  <dig>  ±  <dig>  days of age and  <dig>  ±  <dig>  kg of body weight.
rfi of individual pigs were calculated as described  <cit>  with modifications and shown in fig.
<dig>  briefly, average daily feed intake  was estimated by fitting a quadratic polynomial regression model of the daily feed intake from the on-test day to the off-test day on the number of days on test for each pig; and individual average daily gain  was estimated as the slope from simple linear regression of bi-weekly bw on the number of days on test  <cit> .
a single trait animal model was used to analyze adfi with adjustments for fixed covariates of metabolic mid-body weight , adg, bf, deviation of the on-test weight from 50 kg, deviation of the off-test weight from 118 kg, and deviation of the on-test age from 90 days, and the random effect of pen.
the rfi value for each pig was estimated as the residual of the fitted model .
the growth performance and feed intake data of these animals have been published  <cit> .fig.
<dig> distribution of rfi values and sample selection for rna-seq and rt-qpcr assays.
post-weaning blood samples of  <dig> pigs  with extremely low rfi from gilts of the low rfi line  and  <dig> pigs  with extremely high rfi from gilts of the high rfi line  fed the high-fiber, low-energy diet  or the low-fiber, high-energy diet  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  <dig> 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  <dig> 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  <dig> animals originally selected for rna-seq.
the distribution of rfi values of barrows from each line by diet combination is also shown for reference
rna extraction and globin depletion
based on the rfi phenotype of individual pigs, the juvenile blood samples of  <dig> out of the  <dig> gilts were selected for rna-seq.
in detail, for the low rfi line, we selected eight samples from animals with extremely low rfi values on the lfd and eight samples from animals with extremely low rfi values on the hfd ; while for the high rfi line, we selected eight samples from animals with extremely high rfi values on the lfd and eight samples from animals with extremely high rfi values on the hfd  .
total rna extraction of the  <dig> samples was processed in four batches such that for each batch, two samples were randomly selected from each line by diet combination.
within each batch, the processing order of samples was randomized beforehand and followed in all subsequent procedures where blocking was necessary.
the total rna was extracted from blood samples using preserved blood rna purification kit i  by following the kit’s manual.
on-column dna digestion was performed as described using dnase i .
globin transcripts  were depleted by following an rnase h-mediated method  <cit> .
the quantity and integrity of the rna were determined by using nanodrop  <dig>  and bioanalyzer  <dig>  before and after globin depletion.
the efficiency of globin depletion of each sample was checked by conventional rt-qpcr with β-actin  and glyceraldehyde 3-phosphate dehydrogenase  as the internal reference genes.
the total rna of one selected blood sample from the low rfi line on the lfd was dropped from this study due to its low rna quality  was  <dig> , lower than our criterion of rin ≥  <dig>  before globin depletion).
the average rin for the  <dig> remaining samples before and after globin depletion were  <dig>  ±  <dig>  and  <dig>  ±  <dig> , respectively.
detailed information about the selected samples, including pedigree, rna quality, batch, and cbc is available in additional file 1: table s <dig>
rna-sequencing
library construction and sequencing were performed by the dna facility at iowa state university.
briefly, the rna-seq libraries were constructed using the illumina truseq rna sample preparation kit v <dig>  according to manufacturer’s instructions.
for each sample, poly -containing transcripts was enriched with oligo-dt-coated magnetic beads from  <dig>  to 2 μg of total rna with globin transcript depleted.
the enriched rna was fragmented by heat and reverse transcribed with hexamer random primers.
for each sample, adapters with unique barcodes were ligated to the end-polished cdna fragments.
the libraries were linearly amplified by pcr, size selected and quantitated.
the individual libraries were diluted to  <dig> nm and pooled in approximately equimolar amounts according to the processing batches mentioned above with eight libraries per pool, except for one pool of seven libraries.
one hundred base paired-end sequencing was run on an illumina hiseq <dig> platform with one pool per lane on a flow cell.
quality control, preprocessing and alignment of rna-seq reads
read quality was checked by using fastqc   <cit> .
adapters and low quality bases were trimmed by running trimmomatic   <cit>  in the paired-end mode with the following options, illuminaclip:adapters.fa:2:30:10:1:true leading: <dig> trailing: <dig> slidingwindow:4: <dig> leading: <dig> trailing: <dig> minlen: <dig>  such that the average base quality was not lower than  <dig> for every sliding window of  <dig> bases and the minimum length of kept reads was  <dig> bases.
for each set of raw paired-end reads, trimmomatic outputted a pair of files for the kept paired-end reads, a file for unpaired forward reads , and a file for unpaired backward reads .
trimmed paired-end and unpaired reads were separately aligned to the pig reference genome sscrofa  <dig>   using 2-pass rna-star , using the default settings  <cit> .
read counts per gene per library were summarized by using featurecounts   <cit> , with the resulting sam files for uniquely mapped, paired-end and unpaired reads as inputs, separately, and using the pig genome gtf file  as the genomic annotation reference file.
the default settings for other featurecounts options were used  <cit> .
a final table of read counts per gene of the  <dig> samples  was generated by summing up the individual count tables for paired-end and unpaired reads.
prior to differential expression analysis, hemoglobin genes  and genes with few reads  were removed from the count table to obtain a final count table with  <dig> genes.
this count table was used for the subsequent differential expression analysis and weighted gene co-expression network analysis  after further transformation and adjustment.
differential expression analysis
the statistical programming language r  was used for all statistical analyses, unless indicated otherwise.
differential expression analysis was carried out by using the r package “quasiseq”   <cit> .
for each of the  <dig> genes in the final count table, we used quasiseq to fit a full generalized linear model with a negative binomial response and a log link that included an upper-quartile normalization offset  <cit>  and the fixed effects of rfi group , diet , batch , and the linear covariates of rfi phenotype ,pre_conc , pre_rin , post_conc , post_rin , and concentrations of neutrophils, lymphocytes, monocytes, eosinophils and basophils.
note that, although blood was collected before pigs were fed the different diets, a diet effect was included in our initial full model because diet affected the component traits that were used to estimate rfi and thus the estimation of rfi.
because not all the variables included in the initial full model may be associated with transcript levels, we used a backward variable selection algorithm to identify the most relevant variables  <cit> .
the final model included rfi group, batch, pre_conc, post_rin, and the concentrations of neutrophils, lymphocytes, monocytes and basophils as independent variables.
the default settings for arguments in all function calls were used unless specified otherwise.
in the ql.fit function, the “method” argument was set to “optim”.
the reported p-values, q-values and log <dig>  associated with all tests of significance were calculated by using the qlspline method.
weighted gene co-expression network analysis
before co-expression analysis using the r package “wgcna”   <cit> , the expression levels for the  <dig> genes in the count table were adjusted for all independent variables in the final model used for differential expression analysis except rfi group.
briefly, log-counts per million  were calculated using the voom function of the bioconductor package “limma”  with the upper-quartile normalized counts as input  <cit> .
the lmfit function was used to fit a linear model with log-cpm per gene feature as the responsible variable, and rfi group, batch, and the linear covariates of pre_conc, post_rin, and concentrations of neutrophils, lymphocytes, monocytes and basophils as independent variables.
effects associated with relevant variables  were subtracted from the original log-cpm, to create adjusted transformed gene expression values.
the data matrix consisting of adjusted log-cpm per gene for the  <dig> samples  was used as the input for wgcna.
wgcna was performed by following tutorial i  <cit> , with slight modifications as needed.
in wgcna, all correlation coefficients between gene pairs were calculated by using pearson’s method.
a soft-thresholding power of seven was used by assuming the topology of the unsigned weighted gene co-expression network was scale-free.
the average linkage method was used for all clustering procedures.
the function cuttreedynamic was used for identification of modules.
only modules with a minimum of  <dig> genes were considered.
modules with eigengene correlations no less than  <dig>  were merged using the mergeclosemodules function with cutheight =  <dig> .
the eigengene of a module is the first principal component of the gene expression values of that module and can be considered as a representation of the expression profiles of genes in the module  <cit> .
we then fitted linear regression models with expression levels of the module eigengenes as the response variable and rfi group as the independent variable.
the estimated effect of rfi group and the associated p-value for the null hypothesis that the rfi group term was not useful in explaining expression of the eigengene of the module were used to quantify the strength and significance of the association between the eigengene of a module and rfi group.
hierarchical clustering, generation of heat map, and multi-dimension scaling analysis
the adjusted gene expression matrix was used for hierarchical clustering, heatmap generation, and multi-dimensional scaling  analyses.
spearman correlation coefficients for gene expression between samples were calculated and  <dig> minus this correlation coefficient was used as the distance between a pair of samples for both hierarchical clustering and heatmap construction.
the ward method was used in the function hclust for hierarchical clustering of the samples.
a heatmap was generated with the heat.map <dig> function to visualize the degs .
mds with the first two dimensions was used to visualize the relationships of samples with each other by using the function plotmds from the bioconductor package “limma”   <cit> .
for mds, the distance between each pair of samples was the euclidian distance between them based on the expression of all  <dig> genes.
gene ontology term and pathway enrichment analyses
for all go term and ingenuity pathway analysis -based pathway enrichment analyses, unadjusted p-values were reported  <cit> .
we used bioconductor package “topgo”   <cit>  to perform gene ontology  term analysis, including go biological process , go molecular function  and go cellular component .
go terms associated with each gene were downloaded from ensembl biomart .
the “classical” algorithm, which treats all go terms to be independent of each other, and fisher’s exact test were used to estimate significance of such enrichment using the function runtest.
for go term enrichment analysis of degs , we analyzed degs with higher and lower expression in the low versus high rfi group separately, while for go term enrichment analysis of wgcna modules, genes in a whole module were analyzed together.
the maximum subset of the  <dig> genes that was associated with at least one go-bp, go-mf or go-cc term, respectively, was used as the reference set  in the corresponding go term enrichment analyses.
significantly enriched go terms associated with more than  <dig> annotated genes in the pig genome annotation  were reported.
we performed other enrichment analyses of the degs  using ingenuity pathway analysis  and the integrated pathway analysis database  for systematic enrichment analysis  <cit> .
for ipa-based analysis,  <dig> of the  <dig> genes could be mapped to ipa identifiers via the gene symbols of pig genes and these genes were used as the reference set.
for networks and upstream regulator analysis, both direct and indirect relationships were considered.
the options, “all data sources”, “confidence” and “mutation”, were checked.
for species, “all mammals” was checked.
for tissues and cell lines, only data on tissues and primary cells were considered.
the cutoff for the log <dig>  was set to  <dig> and the q-value cutoff was set to  <dig> .
all degs were analyzed together, with  <dig> of  <dig> degs mapped to ipa identifiers.
for ipad-based analysis of degs, the ipad web server  <cit>  was used.
of the  <dig> degs ,  <dig> genes with human gene symbols were analyzed together.
the default reference set was used, as the ipad server does not allow the user to provide a reference set.
the raw p-values were corrected for multiple testing using the “bh” method  <cit>  for ipad-based analysis.
validation by rt-qpcr
we attempted to validate degs between the rfi groups with q ≤  <dig> , |log <dig> | ≥  <dig> and averaged fpkm   <cit>  within either rfi group ≥  <dig>  which resulted in  <dig> genes.
all pairs of primers corresponding to the  <dig> degs and six internal reference genes were designed and synthesized by the fluidigm corporation , such that the two primers of each pair were separated by exon-exon boundaries and could amplify all isoforms of the target gene if possible .
the efficiency of each primer pair was tested by conventional rt-qpcr on the dna engine opticon  <dig> system  by using the standard curve method  <cit>  and only primer pairs with amplification efficiency no less than  <dig>  were further considered.
as for specificity, only primer pairs that gave products with single peaks in melting curve analyses were used for the downstream rt-qpcr assays.
this resulted in  <dig> pairs of primers for degs and two for internal reference genes  passing the test .
the samples we used for rt-qpcr validation of the degs included  <dig> of the  <dig> rna samples that we had sequenced by rna-seq, and another  <dig> novel samples from gilts from the same contemporary group of pigs, with  <dig> samples from each line by diet combination.
the  <dig> novel blood samples were selected such that the rfi phenotypes of the corresponding animals were evenly distributed across the distribution of rfi phenotypes that was not covered by the rfi phenotypes of the  <dig> animals originally selected for rna-seq .
detailed information about the  <dig> samples used for validation of degs is in additional file 1: table s <dig>  rna was extracted from the novel samples as above in three batches, with two samples from each line by diet combination per batch.
the concentration and quality of  <dig> novel rna samples were determined as above.
total rna without undergoing globin depletion was used for cdna synthesis.
by following the fluidigm user guide for real-time pcr analysis  <cit> , real time-qpcr was done on a  <dig>  dynamic array chip  using the biomarker hd system .
data were analyzed with the fluidigm real-time pcr analysis software with the default settings, to obtain raw ct values.
the raw ct values were corrected for differences in the amount of input rna by using the geometric mean of the ct values of the two internal reference genes for the same rna samples to get – ∆ct = −   <cit> .
differential expression analysis was performed by fitting linear models with – ∆ct values as the response variables and rfi group, rna extraction batch, and the linear covariates of the concentrations of neutrophils, lymphocytes, monocytes and basophils as independent variables.
with the high rfi group as the reference, the estimated effects on – ∆ct for rfi group were defined as the –∆∆ct values.
the p-values associated with the effect of rfi group were adjusted to get q-values by using the “bh” method  <cit> .
if the q-value of the significance test of the rfi group was less than  <dig> , the corresponding gene was considered differentially expressed between the low and high rfi groups by rt-qpcr.
because the amplification efficiencies  of the primers were close to  <dig>  the fold changes of gene expression between the low and high rfi group were calculated as 2–∆∆ct  <cit> .
association analysis of gene expression with rfi phenotype
for association analysis of the expression levels of eigengenes of modules identified by wgcna with rfi phenotype, the eigengene expression levels were first analyzed with a linear model that included rfi group and diet as fixed effects and rfi phenotype as a covariate, along with the two-way and three-way interactions among these three factors.
we did not include rna processing batch, rin, and concentrations of the blood cell types in the linear model because the eigengene expression levels were calculated based on gene expression values that had already been adjusted for these effects.
gene expression determined by rt-qpcr was also used to identify genes associated with rfi phenotype.
the – ∆ct values for the  <dig> target genes in the  <dig> samples were calculated as above and analyzed with a linear model that included rfi group, rna extraction batch and diet as fixed effects, and rfi phenotype and the concentrations of neutrophils, lymphocytes, monocyte and basophils as covariates, along with the two-way and 3-way interactions among rfi group, diet and rfi phenotype.
however, none of the interaction terms were significant after correcting for multiple testing with the “bh” method , thus interaction terms were removed from the models.
the p-values associated with the regression coefficient on the rfi phenotype covariate were adjusted to get q-values by the “bh” method  <cit> .
if the q-value was less than  <dig> , the association between the gene and the rfi phenotype was considered significant.
RESULTS
differentially expressed genes between low and high rfi groups
mauch et al.
<cit>  showed that the low rfi line had significantly lower rfi than the high rfi line on the lfd , but the low rfi line only tended to have lower rfi than the high rfi line on the hfd .
thus, to maximize the contrast, we selected blood samples for rna-seq from gilts of the low rfi line with extremely low rfi values when fed the lfd or the hfd , and from gilts of the high rfi line with extremely high rfi values when fed the lfd or the hfd  .
the metadata, including rna concentration, rins before and after globin depletion, rfi phenotype and cbc test results are in additional file 1: table s <dig>  the distribution of the rfi values for all gilts of each line by diet combination is in fig.
<dig>  to increase the power to detect lowly expressed genes in the whole blood transcriptome, we depleted the most highly expressed transcripts, hemoglobin a and b , using the rnase h-mediated method  <cit> .
alignment of rna-seq reads to the reference genome showed that globin transcripts had been effectively reduced .
after globin depletion, on average, only  <dig>  ±  <dig>  and  <dig>  ±  <dig>  % of the trimmed reads mapped to the hba and hbb genes, respectively.
the numbers of raw reads, trimmed reads, and mapped reads for each sample are in additional file 3: table s <dig>  in summary,  <dig>  ±  <dig>  million pairs of 100-bp raw reads were sequenced per sample.
after trimming,  <dig>  ±  <dig>  % of raw paired-end reads were kept as paired-end reads, while  <dig>  ±  <dig>  % of raw forward reads  and  <dig>  ±  <dig>  % of raw reverse reads  were unpaired reads.
91 ± 1 % of trimmed fragments were mapped to the pig reference genome, with  <dig>  ±  <dig>  %  of trimmed fragments being uniquely mapped.
after removing genes with extremely low expression , and hba and hbb from the count table, we had expression data for  <dig> genes for downstream analyses.
mpetile et al.
<cit>  showed that pigs at  <dig> to 42 days of age from the low rfi line had lower concentrations of lymphocytes, monocytes and basophils than those from the high rfi line in their peripheral blood.
therefore, we decided to account for the concentration of different blood cell types in our analyses so that we could adjust gene expression for these concentration differences.
although the animals were randomly assigned to one of the two diets after the blood samples were collected, we started with a full model that included diet and rfi phenotype  and performed backward selection to establish a final model for the following reasons: 1) the samples for rna-seq were selected based on the rfi phenotype of the animals; 2) the rfi phenotype depended on both line and diet effects; 3) we were interested in identifying genes that might be associated with rfi phenotype.
however, after accounting for rfi group  and other variables identified via backward selection , we found the rfi phenotype, diet, pre_rin, post_conc and the concentration of eosinophils were not significantly associated with the expression levels of most genes .
with our final selected model, we found expression levels of  <dig>   <dig>   <dig>  and  <dig> genes to be significantly associated with the concentrations of neutrophils, lymphocytes, monocytes and basophils, respectively  .
however, we were most interested in degs whose expression differences were primarily dependent on line differences.
therefore, for all downstream analysis, we only considered degs between the low and high rfi groups, while accounting for the effects of relevant variables: batch, post_rin, and the concentrations of neutrophils, lymphocytes, basophils and monocytes.
the numbers of degs between the low and high rfi groups based on different stringencies in terms of q-value and fold change cutoffs are in table  <dig>  a full list of the  <dig> degs for q ≤  <dig>  is in additional file 4: table s <dig>  differential expression of the  <dig> degs with q ≤  <dig>  is shown in fig.
2a.
as the volcano plot  shows, the fold change between the low and high rfi groups was small for most genes.
for a q-value cutoff of  <dig> , only  <dig> degs had a fold change ≥  <dig> or a fold change ≤  <dig> .
however, the transcriptomes of the two groups of animals were collectively different, as the  <dig> samples were grouped into two clusters by rfi group based on hierarchical clustering and mds analyses .table  <dig> the number of differentially expressed genes identified at different q-value and fold change cutoffs
│log2│ >  <dig> a
afc, 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
fig.
<dig> differentially expressed genes and transcriptomic differences between the low and high rfi groups.
a heatmap showing  <dig> degs  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 , while animals with sample names in black were fed the low-fiber, high-energy diet .
the relative orders of genes and samples were determined by two-way hierarchical clustering based on the adjusted transformed gene expression of the  <dig> 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 │log2│ =  <dig>  while red horizontal dash line correspond to q-value of  <dig> .
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  <dig> rna-seq samples.
the samples were hierarchically clustered by using ward method with  <dig> minus spearman correlation as distance.
the spearman correlations between pairs of samples were calculated based on the adjusted transformed expression of the  <dig> genes
annotation of differentially expressed genes by go term and pathway enrichment analyses
to understand the biological differences between rfi groups based on the differential expression analyses, we carried out go term analysis of degs using a less stringent cutoff of q-value  such that we could detect a broader group of significant go terms.
these analyses were conducted separately for degs with higher and lower expression in the low versus high rfi group.
go-bp term analysis showed that genes involved in small molecule biosynthesis, carboxylic acid biosynthesis, organic acid biosynthesis, steroid biosynthesis, antigen processing and presentation of peptide antigen via mhc  class i, and organic hydroxy compound biosynthesis were overrepresented among degs with higher expression in the low versus high rfi group .
among the degs with lower expression in the low rfi group, genes involved in signal transduction, bone mineralization, regulation of phosphorylation, and phosphorylation were overrepresented .
interestingly, go-cc term analysis showed that genes functioning in the proteasome complex or the mitochondrion were enriched among the degs with higher expression in the low versus high rfi group .table  <dig> go terms overrepresented among differentially expressed genes
anumber of genes detected in the blood and associated with a given go term
bnumber of degs associated with a given go term
cexpected number of degs associated with a given go term
in addition, we conducted pathway enrichment analyses of the degs  by assuming pig genes have similar biological functions as their human orthologs.
detailed results from these ipa-based analyses are in additional file 6: table s <dig>  in summary, ipa-based analysis suggested genes functioning in several canonical pathways tended to be enriched among the degs, including p <dig> signaling , t cell receptor signaling , antigen presentation , il- <dig> signaling  and il- <dig> signaling .
but after correcting for multiple testing, none of these pathways were significantly enriched among the degs .
in addition, ipa-based analyses suggested that activities of genes involved in regulating the quantity of t lymphocytes , consumption of oxygen , insulin resistance , cell survival , glucose metabolism and transport , respiration of mitochondrion  and body mass index  were different between the two rfi groups.
ipa-based analysis also suggested several potential upstream regulators of the degs, including il <dig>  which had a higher expression in the low versus high rfi group .
the inferred regulatory network for il <dig> is shown in additional file 7: figure s <dig>
using the integrated pathway analysis database  for systematic enrichment analysis  <cit> , we found genes involved in the immune system   and metabolism   were enriched among the degs.
ipad-based disease-associated gene enrichment analysis suggested genes involved in several diseases were overrepresented among degs , including taste disorder, eating disorder, anorexia, hyperphagia, obesity, insulin resistance, mitochondrial diseases and lymphocytosis.
details are in additional file 8: table s <dig>  this is consistent with the phenotypic differences in feeding behavior  <cit> , body composition  <cit> , and growth rate  <cit>  that have been observed between the two rfi lines.
validation of degs
we used rt-qpcr to validate the differential expression of  <dig> degs that were selected based on average expression levels within rfi group , differential expression levels between rfi groups | ≥  <dig> and q ≤  <dig> ) and primer performance.
to test whether the expression differences we detected were due to differences between the lines, rather than due to the comparison of pigs of the low rfi group to pigs of the high rfi group, we selected another  <dig> samples from the same parity of the 9th generation of the two rfi lines.
these  <dig> samples were selected such that the rfi values of the corresponding animals were roughly evenly distributed across the ranges of rfi values for each line by diet combination, which were not covered by the rfi values of the  <dig> animals originally selected for rna-seq .
the rna from these new samples as well as the  <dig> samples used for rnaseq analyses were used in rt-qpcr analyses.
data from the original  <dig> samples and the new  <dig> samples were analyzed independently as well as jointly.
a total of  <dig> of the  <dig> degs were confirmed in the joint analysis of all  <dig> samples , while  <dig> of the  <dig> degs were confirmed when analyzing the original  <dig> samples .
however, only  <dig> of the  <dig> degs were validated in the analysis of the novel  <dig> samples .fig.
<dig> validation of degs by rt-qpcr.
<dig> of the  <dig> selected degs between low and high rfi groups were confirmed by rt-qpcr when using the  <dig> samples,  <dig> of which were used for rna-seq and another  <dig> of which were novel, as shown in this figure ;  <dig> of the  <dig> selected degs were confirmed by rt-qpcr when using the  <dig> samples that were used for rna-seq ; but only  <dig> of the  <dig> selected degs were validated by rt-qpcr when using the  <dig> novel samples  .
for comparison, the log <dig>  of these genes determined by rna-seq were also displayed.
genes were ordered based on their log <dig>  as determined by rt-qpcr for display.
error bars for rt-qpcr assays show the standard errors of mean log <dig> .
degs not confirmed by rt-qpcr are labeled in red.
genes without corresponding human orthologs are labeled with the last  <dig> digits of their ensembl gene ids, with common prefix “enssscg000000” omitted for simplicity.
for example, the id for “29500” should be enssscg00000029500
based on the rt-qpcr results of the  <dig> samples, we found the expression of  <dig> of the  <dig> degs  to be significantly associated with rfi phenotype .
they were lrp <dig> , enssscg <dig> , enssscg <dig>  pdl <dig>  and enssscg <dig> .
co-expressed gene modules associated with rfi
since conventional differential expression analysis considers each gene independently and suffers from loss of power due to correction for multiple testing, we performed weighted gene co-expression gene network analysis   <cit>  based on the adjusted gene expression matrix  to identify modules  differentially expressed between the rfi groups and associated with rfi phenotype.
we identified four modules, designated as c1-lightcyan, c2-darkturqoise, c3-skyblue <dig>  and c4-black, whose eigengene expression levels were significantly differentially expressed between the low and high rfi groups .
the eigengene of a module is a weighted average of the expression profiles of genes in the module, calculated as the first principal component score.
the expression levels of the eigengenes for the  <dig> samples are shown in fig.
4a.
among the genes in module c1-lightcyan, those involved in lipid metabolism-related biological processes, such as lipid metabolic process, lipid biosynthesis and steroid biosynthesis were overrepresented .
in module c2-darkturqoise, genes involved in biological processes related to bone mineralization, immunity and stress response and lipid metabolism were overrepresented .
we also found that degs identified by rna-seq  significantly overlapped with each of the four differentially expressed modules , which suggests that the modules identified by wgcna are not computational artifacts as they were enriched for the degs between the rfi lines.
in addition, we found that the expression levels of the eigengenes of modules c3-skyblue <dig> were significantly associated with rfi phenotype  .table  <dig> summary of wgcna modules differentially expressed between the low and high rfi groups
areg.
coef., regression coefficient estimated by regressing the expression level of the eigengene of a module on rfi groups, with the high rfi group as the reference
fig.
<dig> 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
discussion
small but significant blood transcriptomic differences between the low and high rfi groups
by transcriptomic analyses, we identified  <dig>   <dig> and  <dig> degs between the low and high rfi groups at q-value cutoffs of  <dig> ,  <dig>  and  <dig> , respectively.
we validated  <dig> of the  <dig> selected degs by rt-qpcr in a joint analysis of  <dig> samples:  <dig> of the  <dig> samples used for rna-seq plus  <dig> novel samples from the same population .
we validated  <dig> of the  <dig> selected degs by rt-qpcr when using  <dig> of the  <dig> samples that were used for rna-seq, but only  <dig> of the  <dig> selected degs were validated when using only the  <dig> novel samples .
our validation rate by rt-qpcr of degs detected by rna-seq was low compared to the high validation rate reported by wang et al.
<cit> , especially when based on the  <dig> novel samples.
there are several possible reasons for failure to validate some of the  <dig> degs that were detected by rna-seq by rt-qpcr in this study.
first, some of the selected degs could be false positives, especially as we selected degs with q ≤  <dig>  for validation, although we did not detect an obvious association between the validation rate and the q values for the degs.
second, the primers used for rt-qpcr maybe did not quantify all isoforms that were measured by rna-seq.
although the primers were designed to amplify all isoforms, these designs were based on poorly annotated gene models: most of the  <dig> genes only have one isoform in the ensembl pig genome annotation .
third, the  <dig> novel samples used for rt-qpcr were not the exact biological replicates of the  <dig> samples used for rna-seq given their different rfi phenotype.
in addition, as observed in other differential expression analyses between low and high rfi groups conducted in poultry, pigs and cattle , the magnitude of the differential expression in terms of fold change was generally small for most genes.
this could be because the experimental units for the low and high rfi groups were genetically heterogeneous within each group, which cancelled out some differences between extreme individuals of the two groups.
in addition, the experiments were conducted without strong external stimuli, thus the expression levels of genes were in their normal physiological range.
it was also possible that the differences in gene expression in blood at the post-weaning stage were not as dramatic as during the grow-finish phase when feed efficiency was tested.
however, it is possible that small changes in the transcript abundance of some genes involved in metabolism, transcription regulation, and signal transduction, could have significant effects on rfi phenotype.
as our differential analysis indicated, a significant portion of genes involved in several biosynthetic processes, signal transduction, and regulation of phosphorylation were differentially expressed between the low and high rfi groups.
lastly, we cannot exclude the possibility that there is a weak or indirect relationship between early blood gene expression and later rfi phenotype.
nevertheless, based on the adjusted global gene expression profiles of all samples, hierarchical clustering and mds analyses revealed that the post-weaning blood transcriptome was different between the two rfi groups.
potential relationship of the degs and differentially co-expressed gene modules in blood with rfi
although we did not attempt to validate most of the degs or the differentially co-expressed gene modules between the two rfi groups due to the generally small fold changes, we found some interesting biological processes and cellular components underlying them, which suggested their potential relationships with rfi phenotype.
go term enrichment analysis of degs suggested that genes involved in small molecule  biosynthesis, antigen processing and presentation of peptide antigen via mhc class i, and steroid biosynthesis were enriched among degs with higher expression in the low versus high rfi group.
degs involved in small molecular biosynthesis included fads <dig>  fads <dig>  and elovl <dig>  fads <dig> and fads <dig> are important genes in regulating the synthesis of polyunsaturated fatty acids, which have pleiotropic influences on health and diseases by functioning in several pathways, including metabolism and immunity  <cit> .
of note, fads <dig> also had a higher expression level in the liver in the low rfi versus high rfi group of nelore steers  <cit> .
elovl <dig> encodes one of the rate-limiting enzymes in elongation of very long chain fatty acids .
in mice, a proposed physiological function of elovl <dig> is to maintain lipid homeostasis by replenishing the intracellular pool of triacylglycerol  <cit> .
mat2b, is one of the genes associated with steroid biosynthesis but had lower expression in the low versus high rfi group.
it encodes a regulatory subunit of methionine adenosyltransferase, mat ii, which catalyzes the synthesis of s-adenosyl methionine .
sam is a key methyl donor in transmethylation reaction and polyamine biosynthesis and also functions as a cofactor in key metabolic pathways.
down-regulation of matiib expression causes a  <dig> to10-fold increase in intracellular sam levels  <cit> .
taken together, the expression difference of genes involved in small molecule biosynthesis and steroid biosynthesis might have significant effects on overall metabolism, immunity and beyond.
seven of the  <dig> genes associated with go term go: <dig>  had higher expression in the low rfi group.
these included sla- <dig> and tap <dig>  it has been shown that in post-weaning pigs resistant to escherichia coli f <dig>  sla- <dig> has higher expression level than in susceptible pigs  <cit> .
tap <dig> transports antigens from the cytoplasm to the endoplasmic reticulum for loading antigen peptide onto mhc class i molecules  <cit> .
in addition, we found  <dig> of the  <dig> degs associated with go term go: <dig>  to have higher expression in the low versus high rfi group.
the only exception was psmd <dig>  which had lower expression in the low versus high rfi group.
it has been shown that overexpression of psmd <dig> inhibits assembly and activity of 26s proteasome  <cit> , which is an important component in processing and presenting intracellular peptide antigen via mhc class i to cd8+ t cells.
antigen presenting cells  have both constitutive proteasomes and immunoproteasomes, both of which are important for antigen peptide processing  <cit> .
notably, an immunoproteasome subunit psmb <dig>  had the largest difference in expression level among the  <dig> degs that function in the proteasome complex.
given the increased expression levels of non-inhibitory proteasomal components in blood of the low rfi group, it is likely that proteasomal activity is higher in blood from the low versus high rfi group.
taken together, the expression differences of genes involved in antigen processing and peptide antigen presentation suggest that low rfi animals might have a more active system of antigen peptide processing and presentation, which might improve the robustness of their immune system of the low rfi animals.
consistent with this speculation, dunkelberger et al.
<cit>  recently showed that the low rfi line was less affected by experimental infection with porcine reproductive and respiratory syndrome virus , which mainly infects mature macrophages, derivatives of monocytes.
in view of the lower concentration of monocytes in the low rfi line, the putative higher activity of antigen processing and peptide antigen presentation might lead to an increased efficiency of innate immunity in the low rfi pigs.
it would be interesting to test the activity of the blood proteasome in more detail in these lines for a role in feed efficiency.
another interesting finding was that genes functioning in the mitochondrion were significantly overrepresented among degs with higher expression in the low versus high rfi group.
of the  <dig> blood genes associated with go-cc term mitochondrion ,  <dig> were degs with higher expression in the low versus high rfi group.
mitochondria play important roles in many biological processes, including fatty acid metabolism, amino acid metabolism, and steroidogenesis  <cit> .
importantly, mitochondrial function has been associated with feed efficiency in poultry, cattle and lambs  <cit> .
as well, it has been shown that the mitochondrial proteome profiles are different between these low and high rfi lines of pigs and the mitochondria from the liver and ld muscle of the low rfi pigs produced less ros .
interestingly, we found that four genes involved in detoxification of ros, gpx <dig>  sod <dig> and cat, had higher expression in the low versus high rfi group.
another interesting gene that functions in mitochondria is g0s <dig>  which is a noncompetitive inhibitor of adipose lipase, a rate-limiting lipase of triglyceride hydrolysis  <cit> .
g0s <dig> had higher expression in the low versus high rfi group.
it has been shown that g0s <dig> can inhibit ros production in endothelial cells  <cit> .
in addition, we found higher expression in blood from the low rfi group for genes that encode  <dig> mitochondrial ribosomal proteins , nine components of mitochondrial complex i , one component of complex ii , cytochrome c  and four components of complex iv , two components of the complex v , four key enzymes of the tca cycle , four enzymes involved in mitochondrial dna repair , and two enzymes involved in nucleotide metabolism , and three translocases of the inner/outer membranes .
taken together, it is reasonable to speculate that the low rfi group has more efficient mitochondria and can more efficiently handle oxidative stress.
among the  <dig> degs ,  <dig> genes were associated with go: <dig> , and genes involved in signal transduction were significantly enriched among degs with lower expression in the low versus high rfi group.
degs associated with signal transduction were mainly cytokines, receptors and kinases.
signal transduction-related genes with much lower expression in the low versus high rfi group included lrp <dig>  wnt10b, and fzd <dig>  lrp <dig> belongs to the low-density lipoprotein receptor family and plays a key role in lipoprotein endocytosis and as a co-receptor in wnt/β–catenin signaling.
lrp <dig> is also involved in regulating lipid homeostasis and body fat mass  <cit> , while lrp <dig>  wnt10b and fzd <dig> can function together in the wnt/β-catenin signaling pathway.
besides playing important roles in bone metabolism, the wnt/β-catenin signaling pathway also takes part in glycolysis and regulates mitochondrial physiology and insulin sensitivity and is thus linked to metabolic diseases  <cit> .
on the other hand, we found a few genes involved in signal transduction but with much higher expression in the low versus high rfi group, such as il <dig> and irs <dig>  il <dig> plays important roles as a pleiotropic cytokine in innate and adaptive immunity  <cit> .
irs <dig> is an important player in both the insulin signaling and the igf- <dig> signaling pathways and igf- <dig> is an important anabolic hormone, playing a key role in growth.
note that juveniles from the low rfi line have been shown to have lower serum concentration of igf- <dig> and the low rfi line grows only slightly slower than the high rfi line.
the higher expression of irs <dig> might partially compensate for the lower level of igf- <dig> in the low rfi line.
it is interesting to investigate the expression and activity of irs <dig> in other tissues of the two lines, such as the muscle.
the differential expression of many genes associated with signal transduction might be related to differences in metabolism and immune response between the rfi lines.
lastly, we found four co-expression modules to be differentially expressed between the rfi groups, and these modules shared a significant portion of their genes with the list of degs.
the eigengene expression of module c1-lightcyan was lower in the low rfi group than in the high rfi group.
the top highly connected genes in this module were mat2b, cpt1a, and lrp <dig>  these genes were also differentially expressed between the low and high rfi groups.
since we have already discussed mat2b and lrp <dig>  we focus on cpt1a herein.
cpt1a had lower expression in the low versus high rfi group and plays an important role in importing long chain fatty acids into mitochondria, by catalyzing the primary step of mitochondrial fatty acid oxidation  <cit> .
inhibition of hypothalamic cpt1a has been shown to decrease feed intake and glucose production in rats  <cit> .
thus it will be interesting to investigate its expression in other tissues including the hypothalamus.
in addition, go-bp term analysis showed that genes involved in lipid metabolism-related processes were enriched in module c1-lightcyan.
because low rfi pigs have less body fat and are leaner  <cit>  and that deposition of energy as fat costs more energy than as protein  <cit> , genes in this module might be relevant to differences in rfi.
potential predictive biomarkers for rfi
degs whose expression levels are associated with the phenotype of interest, as well as genes whose expression profiles are similar to the eigengenes of co-expression modules associated with the phenotype are often considered as good candidate biomarkers for the phenotype.
based on our association tests, the best candidate biomarkers included lrp <dig>  enssscg <dig>  enssscg <dig>  pdl <dig> and enssscg <dig>  lrp <dig> and enssscg <dig> were ones of the degs between the two rfi groups and ones of the highly connected, eigengene-like genes in modules c1-lightcyan and c2-darkturqoise, respectively .
notably, we only found a very few genes whose expression levels were significantly associated with rfi phenotype, which might be because the estimated rfi values were residuals which represented both random errors and true feed efficiency and thus were not accurate, or because there were not many genes whose post-weaning expression levels in blood were well associated with the rfi values measured during the grow-finish phase in blood, or because there may be little or no relationship of blood gene expression in blood and rfi phenotype regardless of age.
CONCLUSIONS
as far as we know, this is the first study to explore transcriptomic differences in blood between pig lines with divergent rfi by rna-seq.
we found that the blood transcriptome was clearly different between the low and high rfi groups, although only a small number of genes showed large fold changes of expression between the two groups.
the two rfi groups may be different in mitochondrial and proteasomal activities, small molecule biosynthetic process, and signal transduction.
these blood transcriptomic differences may be related to the difference in feed efficiency between these two groups, although the observed relationships of post-weaning blood gene expression with rfi phenotype measured during the grow-finish phase were not strong.
the top candidate biomarkers for predicting rfi included lrp <dig>  enssscg <dig>  enssscg <dig>  pdl <dig> and enssscg <dig>
availability of supporting data
the rna-seq data sets supporting the results of this article is available in arrayexpress under accession number: .
additional files
additional file 1: table s <dig>  meta data for rna-seq and rt-qpcr assays.
additional file 2: table s <dig>  primers used to validate differentially expressed genes by rt-qpcr.
additional file 3: table s <dig>  statistical summary of raw, trimmed, mapped and residual globin reads.
additional file 4: table s <dig>  differentially expressed genes identified by rna-seq.
additional file 5: figure s <dig>  mds plot showing relationships of rna-seq samples.
additional file 6: table s <dig>  summary of ingenuity ipa-based pathway analyses.
additional file 7: figure s <dig>  subnetwork consisting of differentially expressed genes regulated by the inferred upstream regulator il- <dig>  nodes were colored based on the log ratio of averaged expression levels of genes between groups, with up-regulated genes in red and down-regulated genes in green.
additional file 8: table s <dig>  summary of ipad-based pathway and disease enrichment analyses.
additional file 9: table s <dig>  validation of differentially expressed genes by rt-qpcr assays.
additional file 10: table s <dig>  association of gene expression with rfi phenotype.
additional file 11: table s <dig>  differentially expressed wgcna modules between the rfi groups.
additional file 12: table s <dig>  go-bp terms enriched in the wgcna modules.
additional file 13: table s <dig>  association of module eigengenes with rfi phenotype.
abbreviations
adfiaverage daily feed intake
adgaverage daily gain
apcantigen presenting cells
bhbenjamini–hochberg
bwbody weight
cbccomplete blood count
cpmcount per million mapped reads
degdifferentially expressed gene
fcfold change
fefeed efficiency
firefeed intake recording equipment
fpkmfragments per kilobase of exon per million fragments mapped
go-bpcc and mf, gene ontology-biological process, cellular component, and molecular function
gwasgenome-wide association study
hfdhigh-fiber, low-energy diet
hrfihigh rfi
ipaingenuity pathway analysis
ipadintegrated pathway analysis database
ldlongissimus dorsi
lfdlow-fiber, high-energy diet
lrfilow rfi
mbwmetabolic mid-body weight
mdsmulti-dimensional scaling
mhcmajor histocompatibility complex
rfiresidual feed intake
rinrna integrity number
rna-seqrna-sequencing
rosreactive oxygen species
snpsingle nucleotide polymorphism
wgcnaweighted gene co-expression network analysis
competing interests
the authors declare that they have no competing interests.
authors’ contributions
jcmd developed the lines of pigs.
hl, yn, dn, jcmd and ckt conceived and designed the experiments; hl performed the experiments; hl and yn performed analysis and interpreted the results; hl, yn, dn, jcmd and ckt drafted and revised the manuscript.
all authors read and approved the final manuscript.

 feed efficiency is one of the important traits directly related to profitability, productivity, and sustainability in the pork industry. while many non-genetic strategies have been developed to improve feed efficiency of pigs, improving feed efficiency by genetic selection can be a sustainable alternative. an increasingly popular measure of fe is residual feed intake, which is the difference between the actual and expected feed intake of an animal for production and maintenance. importantly, rfi is moderately heritable and responds well to genetic selection in pigs. over the last decade, significant efforts have been made to understand the molecular, genetic and physiological basis of rfi in pigs. researchers have found many interesting differences across multiple tissues between pigs with divergent rfi phenotypes. for example, as compared to high rfi pigs, low rfi pigs have altered feeding behaviors, slightly lower growth rate, less back fat, lower protein turnover rate in the muscle, altered mitochondrial protein profiles, less mitochondrial reactive oxygen species production, and lower levels of leptin and juvenile igf- in circulating blood. to explore the genetic basis of rfi in pigs, several genome-wide association studies have been conducted. some chromosomal regions tagged by single nucleotide polymorphisms have been found to be associated with rfi, but these associations were not consistent across studies and explained only small portions of the genetic variance for rfi. the wide range of differences across multiple tissues between lines of pigs with divergent rfi and the lack of snps with major effects on rfi suggest that rfi is a highly polygenic, quantitative trait with multiple tissues contributing to its variation. global gene expression profiling technologies have also been used to explore the molecular basis of rfi in pigs. by profiling the transcriptomes of the adipose tissue of two lines of pigs divergently selected for rfi with gene expression microarrays, lkhagvadorj et al. found that genes involved in the lipid metabolic pathway were overrepresented among the differentially expressed genes that had lower expression in low versus high rfi pigs, and genes involved in carbohydrate metabolism and response to stress were overrepresented among the degs that had higher expression in low versus high rfi pigs. they also found the leptin-related gene network to be different between the two lines. recently, vincent et al. and jing et al. profiled the transcriptome of the longissimus dorsi muscle from pigs with divergent rfi by gene expression microarray and rna-seq, respectively. using pigs from lines divergently selected for rfi, vincent et al. found genes involved in protein synthesis and glycolysis, and genes associated with mitochondrial energy/oxidative metabolism had higher and lower expression, respectively, in the low versus high rfi line. using yorkshire barrows with extreme phenotypes for rfi, jing et al. found that genes involved in glycolysis had lower expression in the low versus high rfi group, while genes involved in muscle proliferation and differentiation had higher expression in the low versus high rfi group. surprisingly, these two studies shared no degs and proposed opposite differences in glycolytic activities in the low versus high rfi pigs. therefore, in consideration of the complexity of rfi and the inconsistency from study to study, the molecular mechanisms underlying rfi in pigs are still largely unclear. to incorporate rfi into animal breeding programs for improving feed efficiency, feed intake, body weight gain and back fat depth must be recorded on individual pigs. as it is very expensive and time-consuming to record feed intake on individual animals, convenient and predictive biomarkers for rfi that can be measured at an early age are in demand. in cattle, chen et al. successfully used degs identified in the liver of angus bulls that were divergently selected for rfi to classify angus steers from the same divergent rfi lines, and al-husseini et al. developed a rfi predictor using of these degs and validated it in an independent angus population. these biomarkers are, however, not very practical because invasive liver biopsies are needed. the peripheral blood is an informative tissue not only because it carries a variety of cells directly involved in immunity and inflammation, but also because it interacts with every organ and tissue in the body via bioactive circulating factors, such as nutrients, metabolites, cytokines, hormones and exosomal cargoes, which are released from the same or different organs or tissues. these bioactive factors interact with blood cells and thus might modify the gene expression profiles of the blood cells dynamically. molecular signatures in circulating blood, including gene expression profiles, have been shown to reflect the physiopathological status, growth stage and lifestyle of subjects. due to its easy accessibility and informativeness, blood has become a popular sample for disease diagnosis, prediction, prognosis, and biomarker discovery. interestingly, the concentration of igf- in serum at a young age has been shown to be different between animals with divergent rfi in poultry and livestock, including pigs. blood cell profiles at early growth stages have also been found to be different between livestock with divergent rfi phenotypes, including pigs and cattle. these results suggest that animals with divergent rfi phenotypes, measured later in life, have early physiological differences in circulating blood that may be reflected in blood gene expression profiles at these early stages. in this study, our objective was to determine the blood transcriptomic differences between post-weaning pigs from two lines divergently selected for rfi, to explore potential molecular mechanisms underlying rfi in peripheral blood and to develop a list of candidate biomarkers for rfi prediction. we hypothesized that post-weaning expression levels of some genes in whole blood were correlated with rfi phenotype measured during the grow-finish phase. we identified degs with q ≤  and four co-expression modules that were differentially expressed between the low and high rfi groups. a set of selected degs were validated by reverse transcription-quantitative pcr. several interesting biological processes underlying degs and differential co-expression modules were suggested. we also identified several candidate biomarkers for rfi. methods animals, blood sample collection and complete blood count test the experimental protocols for this study were approved by the institutional animal care and use committee at iowa state university under permit number 11-1-4996-s. all pigs were from parity of generation nine of the two lines divergently selected for residual feed intake: the low and high rfi lines. bunter et al. found that the igf- concentration in blood measured between to 42 days of age differed between the two lines and were genetically correlated with rfi measured later during the grow-finish phase. so blood samples were collected from the jugular vein into tempus™ blood rna tubes for long-term storage at −80 °c from post-weaning piglets of the two lines in this age range. meanwhile, blood from those animals was also collected into edta tubes and kept at 4 °c before cbc tests, which were performed on the bleeding day, if possible, or the next morning as described. differences in cbc profiles between the two lines have been published elsewhere. at  ±  days of age and  ±   kg of body weight, gilts and barrows were randomly assigned to mixed-line, mixed-sex finishing pens with electronic single-space feeders for feed efficiency test, with pens being randomly assigned to either of two diets: a high-fiber, low-energy diet and a low-fiber, high-energy diet. the animals on feed efficiency test consisted of barrows and gilts from the low rfi line plus barrows and gilts from the high rfi line fed the hfd, and barrows and gilts from the low rfi line plus barrows and gilts from the high rfi line fed the lfd. all pigs had ad libitum access to feed and water. individual feed intake was real-time recorded, body weight was recorded biweekly, and backfat depth above the 10th rib and loin muscle area were recorded at the end of the test, when the pigs were  ±   days of age and  ±   kg of body weight. rfi of individual pigs were calculated as described with modifications and shown in fig.  briefly, average daily feed intake was estimated by fitting a quadratic polynomial regression model of the daily feed intake from the on-test day to the off-test day on the number of days on test for each pig; and individual average daily gain was estimated as the slope from simple linear regression of bi-weekly bw on the number of days on test. a single trait animal model was used to analyze adfi with adjustments for fixed covariates of metabolic mid-body weight, adg, bf, deviation of the on-test weight from 50 kg, deviation of the off-test weight from 118 kg, and deviation of the on-test age from 90 days, and the random effect of pen. the rfi value for each pig was estimated as the residual of the fitted model. the growth performance and feed intake data of these animals have been published.fig. distribution of rfi values and sample selection for rna-seq and rt-qpcr assays. post-weaning blood samples of pigs with extremely low rfi from gilts of the low rfi line and pigs with extremely high rfi from gilts of the high rfi line fed the high-fiber, low-energy diet or the low-fiber, high-energy diet 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 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 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 animals originally selected for rna-seq. the distribution of rfi values of barrows from each line by diet combination is also shown for reference rna extraction and globin depletion based on the rfi phenotype of individual pigs, the juvenile blood samples of out of the gilts were selected for rna-seq. in detail, for the low rfi line, we selected eight samples from animals with extremely low rfi values on the lfd and eight samples from animals with extremely low rfi values on the hfd; while for the high rfi line, we selected eight samples from animals with extremely high rfi values on the lfd and eight samples from animals with extremely high rfi values on the hfd. total rna extraction of the samples was processed in four batches such that for each batch, two samples were randomly selected from each line by diet combination. within each batch, the processing order of samples was randomized beforehand and followed in all subsequent procedures where blocking was necessary. the total rna was extracted from blood samples using preserved blood rna purification kit i by following the kit’s manual. on-column dna digestion was performed as described using dnase i. globin transcripts were depleted by following an rnase h-mediated method. the quantity and integrity of the rna were determined by using nanodrop and bioanalyzer before and after globin depletion. the efficiency of globin depletion of each sample was checked by conventional rt-qpcr with β-actin and glyceraldehyde 3-phosphate dehydrogenase as the internal reference genes. the total rna of one selected blood sample from the low rfi line on the lfd was dropped from this study due to its low rna quality was, lower than our criterion of rin ≥  before globin depletion). the average rin for the remaining samples before and after globin depletion were  ±  and  ± , respectively. detailed information about the selected samples, including pedigree, rna quality, batch, and cbc is available in additional file 1: table s rna-sequencing library construction and sequencing were performed by the dna facility at iowa state university. briefly, the rna-seq libraries were constructed using the illumina truseq rna sample preparation kit v according to manufacturer’s instructions. for each sample, poly -containing transcripts was enriched with oligo-dt-coated magnetic beads from to 2 μg of total rna with globin transcript depleted. the enriched rna was fragmented by heat and reverse transcribed with hexamer random primers. for each sample, adapters with unique barcodes were ligated to the end-polished cdna fragments. the libraries were linearly amplified by pcr, size selected and quantitated. the individual libraries were diluted to nm and pooled in approximately equimolar amounts according to the processing batches mentioned above with eight libraries per pool, except for one pool of seven libraries. one hundred base paired-end sequencing was run on an illumina hiseq platform with one pool per lane on a flow cell. quality control, preprocessing and alignment of rna-seq reads read quality was checked by using fastqc. adapters and low quality bases were trimmed by running trimmomatic in the paired-end mode with the following options, illuminaclip:adapters.fa:2:30:10:1:true leading: trailing: slidingwindow:4: leading: trailing: minlen: such that the average base quality was not lower than for every sliding window of bases and the minimum length of kept reads was bases. for each set of raw paired-end reads, trimmomatic outputted a pair of files for the kept paired-end reads, a file for unpaired forward reads, and a file for unpaired backward reads. trimmed paired-end and unpaired reads were separately aligned to the pig reference genome sscrofa using 2-pass rna-star, using the default settings. read counts per gene per library were summarized by using featurecounts, with the resulting sam files for uniquely mapped, paired-end and unpaired reads as inputs, separately, and using the pig genome gtf file as the genomic annotation reference file. the default settings for other featurecounts options were used. a final table of read counts per gene of the samples was generated by summing up the individual count tables for paired-end and unpaired reads. prior to differential expression analysis, hemoglobin genes and genes with few reads were removed from the count table to obtain a final count table with genes. this count table was used for the subsequent differential expression analysis and weighted gene co-expression network analysis after further transformation and adjustment. differential expression analysis the statistical programming language r was used for all statistical analyses, unless indicated otherwise. differential expression analysis was carried out by using the r package “quasiseq”. for each of the genes in the final count table, we used quasiseq to fit a full generalized linear model with a negative binomial response and a log link that included an upper-quartile normalization offset and the fixed effects of rfi group, diet, batch, and the linear covariates of rfi phenotype,pre_conc, pre_rin, post_conc, post_rin, and concentrations of neutrophils, lymphocytes, monocytes, eosinophils and basophils. note that, although blood was collected before pigs were fed the different diets, a diet effect was included in our initial full model because diet affected the component traits that were used to estimate rfi and thus the estimation of rfi. because not all the variables included in the initial full model may be associated with transcript levels, we used a backward variable selection algorithm to identify the most relevant variables. the final model included rfi group, batch, pre_conc, post_rin, and the concentrations of neutrophils, lymphocytes, monocytes and basophils as independent variables. the default settings for arguments in all function calls were used unless specified otherwise. in the ql.fit function, the “method” argument was set to “optim”. the reported p-values, q-values and log associated with all tests of significance were calculated by using the qlspline method. weighted gene co-expression network analysis before co-expression analysis using the r package “wgcna”, the expression levels for the genes in the count table were adjusted for all independent variables in the final model used for differential expression analysis except rfi group. briefly, log-counts per million were calculated using the voom function of the bioconductor package “limma” with the upper-quartile normalized counts as input. the lmfit function was used to fit a linear model with log-cpm per gene feature as the responsible variable, and rfi group, batch, and the linear covariates of pre_conc, post_rin, and concentrations of neutrophils, lymphocytes, monocytes and basophils as independent variables. effects associated with relevant variables were subtracted from the original log-cpm, to create adjusted transformed gene expression values. the data matrix consisting of adjusted log-cpm per gene for the samples was used as the input for wgcna. wgcna was performed by following tutorial i, with slight modifications as needed. in wgcna, all correlation coefficients between gene pairs were calculated by using pearson’s method. a soft-thresholding power of seven was used by assuming the topology of the unsigned weighted gene co-expression network was scale-free. the average linkage method was used for all clustering procedures. the function cuttreedynamic was used for identification of modules. only modules with a minimum of genes were considered. modules with eigengene correlations no less than were merged using the mergeclosemodules function with cutheight = . the eigengene of a module is the first principal component of the gene expression values of that module and can be considered as a representation of the expression profiles of genes in the module. we then fitted linear regression models with expression levels of the module eigengenes as the response variable and rfi group as the independent variable. the estimated effect of rfi group and the associated p-value for the null hypothesis that the rfi group term was not useful in explaining expression of the eigengene of the module were used to quantify the strength and significance of the association between the eigengene of a module and rfi group. hierarchical clustering, generation of heat map, and multi-dimension scaling analysis the adjusted gene expression matrix was used for hierarchical clustering, heatmap generation, and multi-dimensional scaling analyses. spearman correlation coefficients for gene expression between samples were calculated and minus this correlation coefficient was used as the distance between a pair of samples for both hierarchical clustering and heatmap construction. the ward method was used in the function hclust for hierarchical clustering of the samples. a heatmap was generated with the heat.map function to visualize the degs. mds with the first two dimensions was used to visualize the relationships of samples with each other by using the function plotmds from the bioconductor package “limma”. for mds, the distance between each pair of samples was the euclidian distance between them based on the expression of all genes. gene ontology term and pathway enrichment analyses for all go term and ingenuity pathway analysis -based pathway enrichment analyses, unadjusted p-values were reported. we used bioconductor package “topgo” to perform gene ontology term analysis, including go biological process, go molecular function and go cellular component. go terms associated with each gene were downloaded from ensembl biomart. the “classical” algorithm, which treats all go terms to be independent of each other, and fisher’s exact test were used to estimate significance of such enrichment using the function runtest. for go term enrichment analysis of degs, we analyzed degs with higher and lower expression in the low versus high rfi group separately, while for go term enrichment analysis of wgcna modules, genes in a whole module were analyzed together. the maximum subset of the genes that was associated with at least one go-bp, go-mf or go-cc term, respectively, was used as the reference set in the corresponding go term enrichment analyses. significantly enriched go terms associated with more than annotated genes in the pig genome annotation were reported. we performed other enrichment analyses of the degs using ingenuity pathway analysis and the integrated pathway analysis database for systematic enrichment analysis. for ipa-based analysis, of the genes could be mapped to ipa identifiers via the gene symbols of pig genes and these genes were used as the reference set. for networks and upstream regulator analysis, both direct and indirect relationships were considered. the options, “all data sources”, “confidence” and “mutation”, were checked. for species, “all mammals” was checked. for tissues and cell lines, only data on tissues and primary cells were considered. the cutoff for the log was set to and the q-value cutoff was set to. all degs were analyzed together, with of degs mapped to ipa identifiers. for ipad-based analysis of degs, the ipad web server was used. of the degs, genes with human gene symbols were analyzed together. the default reference set was used, as the ipad server does not allow the user to provide a reference set. the raw p-values were corrected for multiple testing using the “bh” method for ipad-based analysis. validation by rt-qpcr we attempted to validate degs between the rfi groups with q ≤ , |log | ≥  and averaged fpkm within either rfi group ≥  which resulted in genes. all pairs of primers corresponding to the degs and six internal reference genes were designed and synthesized by the fluidigm corporation, such that the two primers of each pair were separated by exon-exon boundaries and could amplify all isoforms of the target gene if possible. the efficiency of each primer pair was tested by conventional rt-qpcr on the dna engine opticon system by using the standard curve method and only primer pairs with amplification efficiency no less than were further considered. as for specificity, only primer pairs that gave products with single peaks in melting curve analyses were used for the downstream rt-qpcr assays. this resulted in pairs of primers for degs and two for internal reference genes passing the test. the samples we used for rt-qpcr validation of the degs included of the rna samples that we had sequenced by rna-seq, and another novel samples from gilts from the same contemporary group of pigs, with samples from each line by diet combination. the novel blood samples were selected such that the rfi phenotypes of the corresponding animals were evenly distributed across the distribution of rfi phenotypes that was not covered by the rfi phenotypes of the animals originally selected for rna-seq. detailed information about the samples used for validation of degs is in additional file 1: table s rna was extracted from the novel samples as above in three batches, with two samples from each line by diet combination per batch. the concentration and quality of novel rna samples were determined as above. total rna without undergoing globin depletion was used for cdna synthesis. by following the fluidigm user guide for real-time pcr analysis, real time-qpcr was done on a dynamic array chip using the biomarker hd system. data were analyzed with the fluidigm real-time pcr analysis software with the default settings, to obtain raw ct values. the raw ct values were corrected for differences in the amount of input rna by using the geometric mean of the ct values of the two internal reference genes for the same rna samples to get – ∆ct = −. differential expression analysis was performed by fitting linear models with – ∆ct values as the response variables and rfi group, rna extraction batch, and the linear covariates of the concentrations of neutrophils, lymphocytes, monocytes and basophils as independent variables. with the high rfi group as the reference, the estimated effects on – ∆ct for rfi group were defined as the –∆∆ct values. the p-values associated with the effect of rfi group were adjusted to get q-values by using the “bh” method. if the q-value of the significance test of the rfi group was less than, the corresponding gene was considered differentially expressed between the low and high rfi groups by rt-qpcr. because the amplification efficiencies of the primers were close to the fold changes of gene expression between the low and high rfi group were calculated as 2–∆∆ct. association analysis of gene expression with rfi phenotype for association analysis of the expression levels of eigengenes of modules identified by wgcna with rfi phenotype, the eigengene expression levels were first analyzed with a linear model that included rfi group and diet as fixed effects and rfi phenotype as a covariate, along with the two-way and three-way interactions among these three factors. we did not include rna processing batch, rin, and concentrations of the blood cell types in the linear model because the eigengene expression levels were calculated based on gene expression values that had already been adjusted for these effects. gene expression determined by rt-qpcr was also used to identify genes associated with rfi phenotype. the – ∆ct values for the target genes in the samples were calculated as above and analyzed with a linear model that included rfi group, rna extraction batch and diet as fixed effects, and rfi phenotype and the concentrations of neutrophils, lymphocytes, monocyte and basophils as covariates, along with the two-way and 3-way interactions among rfi group, diet and rfi phenotype. however, none of the interaction terms were significant after correcting for multiple testing with the “bh” method, thus interaction terms were removed from the models. the p-values associated with the regression coefficient on the rfi phenotype covariate were adjusted to get q-values by the “bh” method. if the q-value was less than, the association between the gene and the rfi phenotype was considered significant.  differentially expressed genes between low and high rfi groups mauch et al. showed that the low rfi line had significantly lower rfi than the high rfi line on the lfd, but the low rfi line only tended to have lower rfi than the high rfi line on the hfd. thus, to maximize the contrast, we selected blood samples for rna-seq from gilts of the low rfi line with extremely low rfi values when fed the lfd or the hfd, and from gilts of the high rfi line with extremely high rfi values when fed the lfd or the hfd. the metadata, including rna concentration, rins before and after globin depletion, rfi phenotype and cbc test results are in additional file 1: table s the distribution of the rfi values for all gilts of each line by diet combination is in fig.  to increase the power to detect lowly expressed genes in the whole blood transcriptome, we depleted the most highly expressed transcripts, hemoglobin a and b, using the rnase h-mediated method. alignment of rna-seq reads to the reference genome showed that globin transcripts had been effectively reduced. after globin depletion, on average, only  ±  and  ±   % of the trimmed reads mapped to the hba and hbb genes, respectively. the numbers of raw reads, trimmed reads, and mapped reads for each sample are in additional file 3: table s in summary,  ±  million pairs of 100-bp raw reads were sequenced per sample. after trimming,  ±   % of raw paired-end reads were kept as paired-end reads, while  ±   % of raw forward reads and  ±   % of raw reverse reads were unpaired reads. 91 ± 1 % of trimmed fragments were mapped to the pig reference genome, with  ±   % of trimmed fragments being uniquely mapped. after removing genes with extremely low expression, and hba and hbb from the count table, we had expression data for genes for downstream analyses. mpetile et al. showed that pigs at to 42 days of age from the low rfi line had lower concentrations of lymphocytes, monocytes and basophils than those from the high rfi line in their peripheral blood. therefore, we decided to account for the concentration of different blood cell types in our analyses so that we could adjust gene expression for these concentration differences. although the animals were randomly assigned to one of the two diets after the blood samples were collected, we started with a full model that included diet and rfi phenotype and performed backward selection to establish a final model for the following reasons: 1) the samples for rna-seq were selected based on the rfi phenotype of the animals; 2) the rfi phenotype depended on both line and diet effects; 3) we were interested in identifying genes that might be associated with rfi phenotype. however, after accounting for rfi group and other variables identified via backward selection, we found the rfi phenotype, diet, pre_rin, post_conc and the concentration of eosinophils were not significantly associated with the expression levels of most genes. with our final selected model, we found expression levels of and genes to be significantly associated with the concentrations of neutrophils, lymphocytes, monocytes and basophils, respectively. however, we were most interested in degs whose expression differences were primarily dependent on line differences. therefore, for all downstream analysis, we only considered degs between the low and high rfi groups, while accounting for the effects of relevant variables: batch, post_rin, and the concentrations of neutrophils, lymphocytes, basophils and monocytes. the numbers of degs between the low and high rfi groups based on different stringencies in terms of q-value and fold change cutoffs are in table  a full list of the degs for q ≤  is in additional file 4: table s differential expression of the degs with q ≤  is shown in fig. 2a. as the volcano plot shows, the fold change between the low and high rfi groups was small for most genes. for a q-value cutoff of, only degs had a fold change ≥  or a fold change ≤ . however, the transcriptomes of the two groups of animals were collectively different, as the samples were grouped into two clusters by rfi group based on hierarchical clustering and mds analyses.table the number of differentially expressed genes identified at different q-value and fold change cutoffs │log2│ >  a afc, 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 fig. differentially expressed genes and transcriptomic differences between the low and high rfi groups. a heatmap showing degs 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, while animals with sample names in black were fed the low-fiber, high-energy diet. the relative orders of genes and samples were determined by two-way hierarchical clustering based on the adjusted transformed gene expression of the 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 │log2│ =  while red horizontal dash line correspond to q-value of. 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 rna-seq samples. the samples were hierarchically clustered by using ward method with minus spearman correlation as distance. the spearman correlations between pairs of samples were calculated based on the adjusted transformed expression of the genes annotation of differentially expressed genes by go term and pathway enrichment analyses to understand the biological differences between rfi groups based on the differential expression analyses, we carried out go term analysis of degs using a less stringent cutoff of q-value such that we could detect a broader group of significant go terms. these analyses were conducted separately for degs with higher and lower expression in the low versus high rfi group. go-bp term analysis showed that genes involved in small molecule biosynthesis, carboxylic acid biosynthesis, organic acid biosynthesis, steroid biosynthesis, antigen processing and presentation of peptide antigen via mhc class i, and organic hydroxy compound biosynthesis were overrepresented among degs with higher expression in the low versus high rfi group. among the degs with lower expression in the low rfi group, genes involved in signal transduction, bone mineralization, regulation of phosphorylation, and phosphorylation were overrepresented. interestingly, go-cc term analysis showed that genes functioning in the proteasome complex or the mitochondrion were enriched among the degs with higher expression in the low versus high rfi group.table go terms overrepresented among differentially expressed genes anumber of genes detected in the blood and associated with a given go term bnumber of degs associated with a given go term cexpected number of degs associated with a given go term in addition, we conducted pathway enrichment analyses of the degs by assuming pig genes have similar biological functions as their human orthologs. detailed results from these ipa-based analyses are in additional file 6: table s in summary, ipa-based analysis suggested genes functioning in several canonical pathways tended to be enriched among the degs, including p signaling, t cell receptor signaling, antigen presentation, il- signaling and il- signaling. but after correcting for multiple testing, none of these pathways were significantly enriched among the degs. in addition, ipa-based analyses suggested that activities of genes involved in regulating the quantity of t lymphocytes, consumption of oxygen, insulin resistance, cell survival, glucose metabolism and transport, respiration of mitochondrion and body mass index were different between the two rfi groups. ipa-based analysis also suggested several potential upstream regulators of the degs, including il which had a higher expression in the low versus high rfi group. the inferred regulatory network for il is shown in additional file 7: figure s using the integrated pathway analysis database for systematic enrichment analysis, we found genes involved in the immune system and metabolism were enriched among the degs. ipad-based disease-associated gene enrichment analysis suggested genes involved in several diseases were overrepresented among degs, including taste disorder, eating disorder, anorexia, hyperphagia, obesity, insulin resistance, mitochondrial diseases and lymphocytosis. details are in additional file 8: table s this is consistent with the phenotypic differences in feeding behavior, body composition, and growth rate that have been observed between the two rfi lines. validation of degs we used rt-qpcr to validate the differential expression of degs that were selected based on average expression levels within rfi group, differential expression levels between rfi groups | ≥  and q ≤  ) and primer performance. to test whether the expression differences we detected were due to differences between the lines, rather than due to the comparison of pigs of the low rfi group to pigs of the high rfi group, we selected another samples from the same parity of the 9th generation of the two rfi lines. these samples were selected such that the rfi values of the corresponding animals were roughly evenly distributed across the ranges of rfi values for each line by diet combination, which were not covered by the rfi values of the animals originally selected for rna-seq. the rna from these new samples as well as the samples used for rnaseq analyses were used in rt-qpcr analyses. data from the original samples and the new samples were analyzed independently as well as jointly. a total of of the degs were confirmed in the joint analysis of all samples, while of the degs were confirmed when analyzing the original samples. however, only of the degs were validated in the analysis of the novel samples.fig. validation of degs by rt-qpcr. of the selected degs between low and high rfi groups were confirmed by rt-qpcr when using the samples, of which were used for rna-seq and another of which were novel, as shown in this figure; of the selected degs were confirmed by rt-qpcr when using the samples that were used for rna-seq; but only of the selected degs were validated by rt-qpcr when using the novel samples. for comparison, the log of these genes determined by rna-seq were also displayed. genes were ordered based on their log as determined by rt-qpcr for display. error bars for rt-qpcr assays show the standard errors of mean log. degs not confirmed by rt-qpcr are labeled in red. genes without corresponding human orthologs are labeled with the last digits of their ensembl gene ids, with common prefix “enssscg000000” omitted for simplicity. for example, the id for “29500” should be enssscg00000029500 based on the rt-qpcr results of the samples, we found the expression of of the degs to be significantly associated with rfi phenotype. they were lrp, enssscg, enssscg pdl and enssscg. co-expressed gene modules associated with rfi since conventional differential expression analysis considers each gene independently and suffers from loss of power due to correction for multiple testing, we performed weighted gene co-expression gene network analysis based on the adjusted gene expression matrix to identify modules differentially expressed between the rfi groups and associated with rfi phenotype. we identified four modules, designated as c1-lightcyan, c2-darkturqoise, c3-skyblue and c4-black, whose eigengene expression levels were significantly differentially expressed between the low and high rfi groups. the eigengene of a module is a weighted average of the expression profiles of genes in the module, calculated as the first principal component score. the expression levels of the eigengenes for the samples are shown in fig. 4a. among the genes in module c1-lightcyan, those involved in lipid metabolism-related biological processes, such as lipid metabolic process, lipid biosynthesis and steroid biosynthesis were overrepresented. in module c2-darkturqoise, genes involved in biological processes related to bone mineralization, immunity and stress response and lipid metabolism were overrepresented. we also found that degs identified by rna-seq significantly overlapped with each of the four differentially expressed modules, which suggests that the modules identified by wgcna are not computational artifacts as they were enriched for the degs between the rfi lines. in addition, we found that the expression levels of the eigengenes of modules c3-skyblue were significantly associated with rfi phenotype.table summary of wgcna modules differentially expressed between the low and high rfi groups areg. coef., regression coefficient estimated by regressing the expression level of the eigengene of a module on rfi groups, with the high rfi group as the reference fig. 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 discussion small but significant blood transcriptomic differences between the low and high rfi groups by transcriptomic analyses, we identified and degs between the low and high rfi groups at q-value cutoffs of, and, respectively. we validated of the selected degs by rt-qpcr in a joint analysis of samples: of the samples used for rna-seq plus novel samples from the same population. we validated of the selected degs by rt-qpcr when using of the samples that were used for rna-seq, but only of the selected degs were validated when using only the novel samples. our validation rate by rt-qpcr of degs detected by rna-seq was low compared to the high validation rate reported by wang et al., especially when based on the novel samples. there are several possible reasons for failure to validate some of the degs that were detected by rna-seq by rt-qpcr in this study. first, some of the selected degs could be false positives, especially as we selected degs with q ≤  for validation, although we did not detect an obvious association between the validation rate and the q values for the degs. second, the primers used for rt-qpcr maybe did not quantify all isoforms that were measured by rna-seq. although the primers were designed to amplify all isoforms, these designs were based on poorly annotated gene models: most of the genes only have one isoform in the ensembl pig genome annotation. third, the novel samples used for rt-qpcr were not the exact biological replicates of the samples used for rna-seq given their different rfi phenotype. in addition, as observed in other differential expression analyses between low and high rfi groups conducted in poultry, pigs and cattle, the magnitude of the differential expression in terms of fold change was generally small for most genes. this could be because the experimental units for the low and high rfi groups were genetically heterogeneous within each group, which cancelled out some differences between extreme individuals of the two groups. in addition, the experiments were conducted without strong external stimuli, thus the expression levels of genes were in their normal physiological range. it was also possible that the differences in gene expression in blood at the post-weaning stage were not as dramatic as during the grow-finish phase when feed efficiency was tested. however, it is possible that small changes in the transcript abundance of some genes involved in metabolism, transcription regulation, and signal transduction, could have significant effects on rfi phenotype. as our differential analysis indicated, a significant portion of genes involved in several biosynthetic processes, signal transduction, and regulation of phosphorylation were differentially expressed between the low and high rfi groups. lastly, we cannot exclude the possibility that there is a weak or indirect relationship between early blood gene expression and later rfi phenotype. nevertheless, based on the adjusted global gene expression profiles of all samples, hierarchical clustering and mds analyses revealed that the post-weaning blood transcriptome was different between the two rfi groups. potential relationship of the degs and differentially co-expressed gene modules in blood with rfi although we did not attempt to validate most of the degs or the differentially co-expressed gene modules between the two rfi groups due to the generally small fold changes, we found some interesting biological processes and cellular components underlying them, which suggested their potential relationships with rfi phenotype. go term enrichment analysis of degs suggested that genes involved in small molecule biosynthesis, antigen processing and presentation of peptide antigen via mhc class i, and steroid biosynthesis were enriched among degs with higher expression in the low versus high rfi group. degs involved in small molecular biosynthesis included fads fads and elovl fads and fads are important genes in regulating the synthesis of polyunsaturated fatty acids, which have pleiotropic influences on health and diseases by functioning in several pathways, including metabolism and immunity. of note, fads also had a higher expression level in the liver in the low rfi versus high rfi group of nelore steers. elovl encodes one of the rate-limiting enzymes in elongation of very long chain fatty acids. in mice, a proposed physiological function of elovl is to maintain lipid homeostasis by replenishing the intracellular pool of triacylglycerol. mat2b, is one of the genes associated with steroid biosynthesis but had lower expression in the low versus high rfi group. it encodes a regulatory subunit of methionine adenosyltransferase, mat ii, which catalyzes the synthesis of s-adenosyl methionine. sam is a key methyl donor in transmethylation reaction and polyamine biosynthesis and also functions as a cofactor in key metabolic pathways. down-regulation of matiib expression causes a to10-fold increase in intracellular sam levels. taken together, the expression difference of genes involved in small molecule biosynthesis and steroid biosynthesis might have significant effects on overall metabolism, immunity and beyond. seven of the genes associated with go term go: had higher expression in the low rfi group. these included sla- and tap it has been shown that in post-weaning pigs resistant to escherichia coli f sla- has higher expression level than in susceptible pigs. tap transports antigens from the cytoplasm to the endoplasmic reticulum for loading antigen peptide onto mhc class i molecules. in addition, we found of the degs associated with go term go: to have higher expression in the low versus high rfi group. the only exception was psmd which had lower expression in the low versus high rfi group. it has been shown that overexpression of psmd inhibits assembly and activity of 26s proteasome, which is an important component in processing and presenting intracellular peptide antigen via mhc class i to cd8+ t cells. antigen presenting cells have both constitutive proteasomes and immunoproteasomes, both of which are important for antigen peptide processing. notably, an immunoproteasome subunit psmb had the largest difference in expression level among the degs that function in the proteasome complex. given the increased expression levels of non-inhibitory proteasomal components in blood of the low rfi group, it is likely that proteasomal activity is higher in blood from the low versus high rfi group. taken together, the expression differences of genes involved in antigen processing and peptide antigen presentation suggest that low rfi animals might have a more active system of antigen peptide processing and presentation, which might improve the robustness of their immune system of the low rfi animals. consistent with this speculation, dunkelberger et al. recently showed that the low rfi line was less affected by experimental infection with porcine reproductive and respiratory syndrome virus, which mainly infects mature macrophages, derivatives of monocytes. in view of the lower concentration of monocytes in the low rfi line, the putative higher activity of antigen processing and peptide antigen presentation might lead to an increased efficiency of innate immunity in the low rfi pigs. it would be interesting to test the activity of the blood proteasome in more detail in these lines for a role in feed efficiency. another interesting finding was that genes functioning in the mitochondrion were significantly overrepresented among degs with higher expression in the low versus high rfi group. of the blood genes associated with go-cc term mitochondrion, were degs with higher expression in the low versus high rfi group. mitochondria play important roles in many biological processes, including fatty acid metabolism, amino acid metabolism, and steroidogenesis. importantly, mitochondrial function has been associated with feed efficiency in poultry, cattle and lambs. as well, it has been shown that the mitochondrial proteome profiles are different between these low and high rfi lines of pigs and the mitochondria from the liver and ld muscle of the low rfi pigs produced less ros. interestingly, we found that four genes involved in detoxification of ros, gpx sod and cat, had higher expression in the low versus high rfi group. another interesting gene that functions in mitochondria is g0s which is a noncompetitive inhibitor of adipose lipase, a rate-limiting lipase of triglyceride hydrolysis. g0s had higher expression in the low versus high rfi group. it has been shown that g0s can inhibit ros production in endothelial cells. in addition, we found higher expression in blood from the low rfi group for genes that encode mitochondrial ribosomal proteins, nine components of mitochondrial complex i, one component of complex ii, cytochrome c and four components of complex iv, two components of the complex v, four key enzymes of the tca cycle, four enzymes involved in mitochondrial dna repair, and two enzymes involved in nucleotide metabolism, and three translocases of the inner/outer membranes. taken together, it is reasonable to speculate that the low rfi group has more efficient mitochondria and can more efficiently handle oxidative stress. among the degs, genes were associated with go:, and genes involved in signal transduction were significantly enriched among degs with lower expression in the low versus high rfi group. degs associated with signal transduction were mainly cytokines, receptors and kinases. signal transduction-related genes with much lower expression in the low versus high rfi group included lrp wnt10b, and fzd lrp belongs to the low-density lipoprotein receptor family and plays a key role in lipoprotein endocytosis and as a co-receptor in wnt/β–catenin signaling. lrp is also involved in regulating lipid homeostasis and body fat mass, while lrp wnt10b and fzd can function together in the wnt/β-catenin signaling pathway. besides playing important roles in bone metabolism, the wnt/β-catenin signaling pathway also takes part in glycolysis and regulates mitochondrial physiology and insulin sensitivity and is thus linked to metabolic diseases. on the other hand, we found a few genes involved in signal transduction but with much higher expression in the low versus high rfi group, such as il and irs il plays important roles as a pleiotropic cytokine in innate and adaptive immunity. irs is an important player in both the insulin signaling and the igf- signaling pathways and igf- is an important anabolic hormone, playing a key role in growth. note that juveniles from the low rfi line have been shown to have lower serum concentration of igf- and the low rfi line grows only slightly slower than the high rfi line. the higher expression of irs might partially compensate for the lower level of igf- in the low rfi line. it is interesting to investigate the expression and activity of irs in other tissues of the two lines, such as the muscle. the differential expression of many genes associated with signal transduction might be related to differences in metabolism and immune response between the rfi lines. lastly, we found four co-expression modules to be differentially expressed between the rfi groups, and these modules shared a significant portion of their genes with the list of degs. the eigengene expression of module c1-lightcyan was lower in the low rfi group than in the high rfi group. the top highly connected genes in this module were mat2b, cpt1a, and lrp these genes were also differentially expressed between the low and high rfi groups. since we have already discussed mat2b and lrp we focus on cpt1a herein. cpt1a had lower expression in the low versus high rfi group and plays an important role in importing long chain fatty acids into mitochondria, by catalyzing the primary step of mitochondrial fatty acid oxidation. inhibition of hypothalamic cpt1a has been shown to decrease feed intake and glucose production in rats. thus it will be interesting to investigate its expression in other tissues including the hypothalamus. in addition, go-bp term analysis showed that genes involved in lipid metabolism-related processes were enriched in module c1-lightcyan. because low rfi pigs have less body fat and are leaner and that deposition of energy as fat costs more energy than as protein, genes in this module might be relevant to differences in rfi. potential predictive biomarkers for rfi degs whose expression levels are associated with the phenotype of interest, as well as genes whose expression profiles are similar to the eigengenes of co-expression modules associated with the phenotype are often considered as good candidate biomarkers for the phenotype. based on our association tests, the best candidate biomarkers included lrp enssscg enssscg pdl and enssscg lrp and enssscg were ones of the degs between the two rfi groups and ones of the highly connected, eigengene-like genes in modules c1-lightcyan and c2-darkturqoise, respectively. notably, we only found a very few genes whose expression levels were significantly associated with rfi phenotype, which might be because the estimated rfi values were residuals which represented both random errors and true feed efficiency and thus were not accurate, or because there were not many genes whose post-weaning expression levels in blood were well associated with the rfi values measured during the grow-finish phase in blood, or because there may be little or no relationship of blood gene expression in blood and rfi phenotype regardless of age.  as far as we know, this is the first study to explore transcriptomic differences in blood between pig lines with divergent rfi by rna-seq. we found that the blood transcriptome was clearly different between the low and high rfi groups, although only a small number of genes showed large fold changes of expression between the two groups. the two rfi groups may be different in mitochondrial and proteasomal activities, small molecule biosynthetic process, and signal transduction. these blood transcriptomic differences may be related to the difference in feed efficiency between these two groups, although the observed relationships of post-weaning blood gene expression with rfi phenotype measured during the grow-finish phase were not strong. the top candidate biomarkers for predicting rfi included lrp enssscg enssscg pdl and enssscg availability of supporting data the rna-seq data sets supporting the results of this article is available in arrayexpress under accession number:. additional files additional file 1: table s meta data for rna-seq and rt-qpcr assays. additional file 2: table s primers used to validate differentially expressed genes by rt-qpcr. additional file 3: table s statistical summary of raw, trimmed, mapped and residual globin reads. additional file 4: table s differentially expressed genes identified by rna-seq. additional file 5: figure s mds plot showing relationships of rna-seq samples. additional file 6: table s summary of ingenuity ipa-based pathway analyses. additional file 7: figure s subnetwork consisting of differentially expressed genes regulated by the inferred upstream regulator il- nodes were colored based on the log ratio of averaged expression levels of genes between groups, with up-regulated genes in red and down-regulated genes in green. additional file 8: table s summary of ipad-based pathway and disease enrichment analyses. additional file 9: table s validation of differentially expressed genes by rt-qpcr assays. additional file 10: table s association of gene expression with rfi phenotype. additional file 11: table s differentially expressed wgcna modules between the rfi groups. additional file 12: table s go-bp terms enriched in the wgcna modules. additional file 13: table s association of module eigengenes with rfi phenotype. abbreviations adfiaverage daily feed intake adgaverage daily gain apcantigen presenting cells bhbenjamini–hochberg bwbody weight cbccomplete blood count cpmcount per million mapped reads degdifferentially expressed gene fcfold change fefeed efficiency firefeed intake recording equipment fpkmfragments per kilobase of exon per million fragments mapped go-bpcc and mf, gene ontology-biological process, cellular component, and molecular function gwasgenome-wide association study hfdhigh-fiber, low-energy diet hrfihigh rfi ipaingenuity pathway analysis ipadintegrated pathway analysis database ldlongissimus dorsi lfdlow-fiber, high-energy diet lrfilow rfi mbwmetabolic mid-body weight mdsmulti-dimensional scaling mhcmajor histocompatibility complex rfiresidual feed intake rinrna integrity number rna-seqrna-sequencing rosreactive oxygen species snpsingle nucleotide polymorphism wgcnaweighted gene co-expression network analysis competing interests the authors declare that they have no competing interests. authors’ contributions jcmd developed the lines of pigs. hl, yn, dn, jcmd and ckt conceived and designed the experiments; hl performed the experiments; hl and yn performed analysis and interpreted the results; hl, yn, dn, jcmd and ckt drafted and revised the manuscript. all authors read and approved the final manuscript. 
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 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. 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. a weighted gene co-expression network analysis identified four co-expression modules that were differentially expressed between the low and high rfi groups. the post-weaning blood transcriptome was clearly different between the low and high rfi groups. we identified differentially expressed genes between the low and high rfi groups of animals by using rna sequencing technology. 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 to 42 day-old pigs with extremely low and high rfi 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. we validated of selected degs by reverse transcription-quantitative pcr in a joint analysis of of the samples already used for rna-seq plus novel samples from the same contemporary group of pigs. using an analysis of the novel samples alone, only nine of the selected degs were validated. genes involved in small molecule biosynthetic process, antigen processing and presentation of peptide antigen via major histocompatibility complex class i, and steroid biosynthetic process were overrepresented among degs that had higher expression in the low versus high rfi animals. currently, the molecular mechanisms underlying rfi are largely unknown. improving feed efficiency of pigs by genetic selection is of economic and environmental significance. 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. we also found five of the degs and one of the co-expression modules were significantly associated with the rfi phenotype of individual animals. 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. sus scrofafeed efficiencyresidual feed intakeblood transcriptomerna-sequsda national institute of food and agriculture agriculture and food research initiativenifa-afri-2011-68004-30336dekkers jack c. m. issue-copyright-statement© the author 2016 additionally, to incorporate rfi into animal breeding programs, feed intake must be recorded on individual pigs, which is costly and time-consuming. the degs significantly overlapped with genes associated with disease, including hyperphagia, eating disorders and mitochondrial diseases.
 improving feed efficiency of pigs by genetic selection is of economic and environmental significance. an increasingly accepted measure of feed efficiency is residual feed intake. 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 to 42 day-old pigs with extremely low and high rfi 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.  we identified differentially expressed genes between the low and high rfi groups of animals by using rna sequencing technology. we validated of selected degs by reverse transcription-quantitative pcr in a joint analysis of of the samples already used for rna-seq plus novel samples from the same contemporary group of pigs. using an analysis of the novel samples alone, only nine of the selected degs were validated. genes involved in small molecule biosynthetic process, antigen processing and presentation of peptide antigen via major histocompatibility complex 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. a weighted gene co-expression network analysis 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.  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. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords sus scrofafeed efficiencyresidual feed intakeblood transcriptomerna-sequsda national institute of food and agriculture agriculture and food research initiativenifa-afri-2011-68004-30336dekkers jack c. m. issue-copyright-statement© the author 2016 