Abstract Background It is unclear whether improving feed efficiency by selection for low residual feed intake (RFI) compromises pigs’ immunocompetence. Here, we aimed at investigating whether pig lines divergently selected for RFI had different inflammatory responses to lipopolysaccharide (LPS) exposure, regarding to clinical presentations and transcriptomic changes in peripheral blood cells. Results LPS injection induced acute systemic inflammation in both the low-RFI and high-RFI line (n = 8 per line). At 4 h post injection (hpi), the low-RFI line had a significantly lower (p = 0.0075) mean rectal temperature compared to the high-RFI line. However, no significant differences in complete blood count or levels of several plasma cytokines were detected between the two lines. Profiling blood transcriptomes at 0, 2, 6, and 24 hpi by RNA-sequencing revealed that LPS induced dramatic transcriptional changes, with 6296 genes differentially expressed at at least one time point post injection relative to baseline in at least one line (n = 4 per line) (|log[2](fold change)| ≥ log[2](1.2); q < 0.05). Furthermore, applying the same cutoffs, we detected 334 genes differentially expressed between the two lines at at least one time point, including 33 genes differentially expressed between the two lines at baseline. But no significant line-by-time interaction effects were detected. Genes involved in protein translation, defense response, immune response, and signaling were enriched in different co-expression clusters of genes responsive to LPS stimulation. The two lines were largely similar in their peripheral blood transcriptomic responses to LPS stimulation at the pathway level, although the low-RFI line had a slightly lower level of inflammatory response than the high-RFI line from 2 to 6 hpi and a slightly higher level of inflammatory response than the high-RFI line at 24 hpi. Conclusions The pig lines divergently selected for RFI had a largely similar response to LPS stimulation. However, the low-RFI line had a relatively lower-level, but longer-lasting, inflammatory response compared to the high-RFI line. Our results suggest selection for feed efficient pigs does not significantly compromise a pig’s acute systemic inflammatory response to LPS, although slight differences in intensity and duration may occur. Keywords: Sus scrofa, Lipopolysaccharide, Systemic inflammation, Residual feed intake, RNA-seq Background Feed efficiency in pigs is a trait of economic, environmental and societal importance. One increasingly accepted measure of feed efficiency is residual feed intake (RFI), which is defined as the difference between an individual animal’s observed and expected feed intake for growth and maintenance [[39]1]. Thus pigs with a low RFI are more feed efficient than those with a high RFI. Pilot studies of divergent selection for RFI in pigs showed that RFI responds well to genetic selection [[40]2–[41]4]. Compared to high-RFI pigs, pigs selected for low RFI have reduced feed intake, but similar rate of growth [[42]2–[43]4]. This difference occurs likely because the low-RFI pigs are more efficient in allocating resources for production and maintenance [[44]5]. The immune response is a nutrient- and energy-demanding biological process and directly relates to pig health and performance [[45]6, [46]7]. Thus, one interesting question is whether improving feed efficiency by selection for low RFI affects the animal’s ability to respond to immune challenges. Based on resource allocation theory [[47]5], selection for high feed efficiency is expected to compromise the animal’s capacity to handle immune stimulation, such as the response that occurs during infectious diseases [[48]8]. This has been confirmed in studies on chickens and beef cattle, where selection for increased feed efficiency indeed negatively affected their immune system [[49]9]. Several experiments have investigated the potential side effects of selection for divergent RFI phenotypes on the immune response in pigs. First, a study of healthy pigs from the divergently selected RFI lines at Iowa State University (ISU) [[50]2–[51]4], from which representatives were used in the current study, showed that the low-RFI line had lower numbers of monocytes, lymphocytes, and basophils, but a higher hemoglobin concentration and red blood cell volume compared to the high-RFI line [[52]10]. Second, based on results from an experimental infection with the porcine reproductive and respiratory syndrome virus (PRRSV) in pigs from the ISU RFI lines [[53]2–[54]4], Dunkelberger et al. [[55]11] reported that pigs from the low-RFI line had a lower viral RNA load in the blood, a faster humoral immune response to PRRSV, and were less affected in terms of reduced growth rate than pigs from the high-RFI line. Third, in a parallel divergent selection experiment conducted at the French National Institute for Agricultural Research (INRA), the low-RFI line had a lower basal expression of many genes involved in immune and inflammatory response than the high-RFI line [[56]12, [57]13]. Fourth, to test the immune response in the divergently selected RFI lines developed by INRA, piglets from both lines were challenged with the Complete Freund’s Adjuvant (CFA) to induce a non-infectious pneumonia [[58]14–[59]17]. This work showed that both RFI lines handled the inflammatory challenge similarly, but did so by adopting different metabolic strategies [[60]15]. Interestingly, the protein abundance of inflammatory cytokines was lower in the low-RFI line in multiple tissues involved in the immune response 1 week after CFA injection [[61]16]. Lastly, Vigros et al. [[62]18] examined the expression profiles of a set of target genes related to intestinal inflammation in pigs with extremely divergent RFI phenotypes, both before and after an ex vivo LPS exposure of ileal and colonic tissue explants. No differentially expressed genes were detected in the un-stimulated explants. However, the mRNA levels of several proinflammatory cytokines (IL8, IL1, IL6, TNFα, IFNγ) and SOCS3 were lower in the low-RFI than the high-RFI explants following LPS challenge [[63]18]. These authors proposed that low-RFI pigs may adopt an energy saving mechanism during intestinal responses to an immune challenge [[64]18]. Taken together, although no significant negative side effects of selection for increased feed efficiency based on reduced RFI on the immune response in pigs have been detected, little is known about the effect of selection for RFI on the global transcriptomic profiles during the course of acute systemic inflammatory response in pigs. As a major component of the outer membrane of most gram-negative bacteria [[65]19], LPS has been widely used in vertebrates as an inflammatory immunostimulant. In vertebrates, LPS induces inflammatory response mainly via the TLR4-dependent NFκB signaling pathway [[66]20, [67]21], although a TLR4-independent host response to LPS has also been identified [[68]22]. Thousands of genes, including many pro-inflammatory and anti-inflammatory cytokines and chemokines, have been shown to be involved in the LPS-induced inflammatory response in multiple vertebrates, including pigs [[69]23–[70]28]. LPS stimulation can cause many physiological and behavioral changes, including elevated body temperature, dramatic hemodynamic, increased cytokine levels, reduced feed intake, and altered metabolism [[71]29–[72]32]. The objective of this study was to determine whether divergent selection for RFI significantly affects the pig’s systemic inflammatory/immune response to LPS. We induced acute systemic inflammatory response by intramuscular injection of LPS in two lines that were divergently selected for RFI [[73]2, [74]33], and then measured changes in body temperature, complete blood count (CBC), plasma cytokine levels, and peripheral blood transcriptome over time in these lines by using RNA-seq. We detected a few minor differences between the two RFI lines’ systemic inflammatory responses triggered by LPS, although acute responses were substantial in both lines. Our work suggests that the peripheral blood gene expression of low-RFI pigs is only slightly different from that of high-RFI pigs in response to an acute inflammatory challenge. Methods Animals, experimental design, and sample collection The pigs for this study were exclusively from Generation 8 of the Yorkshire pig lines divergently selected for RFI at ISU [[75]2, [76]33], which are developed and owned by Dr. Jack CM Dekkers, one of the co-authors. The experiment was performed in two independent replicates of a 2-by-2 factorial design with repeated measures. The sample size, LPS dosage, and injection route of LPS were determined by referring to our previous study [[77]34]. Figure [78]1a shows the experimental design for one replicate. A total of 28 gilts with an initial body weight (BW) of 63 ± 4 kg from the low-RFI and high-RFI lines (n = 14 per line) were randomly selected and utilized for the two replicates. Pigs were housed in randomly assigned individual metabolism crates, had ad libitum access to water, and were fed a corn-soy-based diet twice daily (8,00 am and 5:00 pm), with feed restriction (1.5 kg/day), as previously described [[79]35]. After a 9-day adaptation period, pigs within each line were randomly assigned to either a control (n = 6, three pigs per line) or LPS treatment (n = 8, four pigs per line) group. Pigs in the treatment group were then challenged with LPS using an established method [[80]34] via intramuscular injection of 30 μg/kg BW of LPS from E. coli O55:B5 (Sigma-Aldrich, St. Louis, MO, USA) dissolved in a endotoxin-free, sterile saline solution at baseline, 0 h post injection (hpi). Pigs in the control group were injected with an equivalent volume of endotoxin-free, sterile saline solution at the equivalent time. The rectal temperature of each pig was measured at 0, 2, 4, 6, and 24 hpi. At 0, 2, 6, and 24 hpi, blood samples were collected from the jugular vein into Tempus™ Blood RNA tubes (Life Technologies, Grand Island, NY, USA) for long-term storage at − 80 °C, into EDTA tubes (BD, Franklin Lakes, NJ, USA) for complete blood count tests and cytokine assays. Injection, temperature measurement, and blood collection followed the same order which was the predefined order of metabolism crates in the pen rooms. The sampling time points were determined by referring to a previous study where time series response of humans to LPS were investigated [[81]24]. Figure [82]1b show the number of blood samples collected at each time points from each replicate. At 168 hpi, all pigs were anesthetized by intraperitoneal injection of sodium pentobarbital (100 mg/kg BW), and exsanguinated by cutting their carotid and jugular vessels. Fig. 1. [83]Fig. 1 [84]Open in a new tab Experimental design and blood sample collection. a The experiment was performed in two replicates. Shown here is one of the two replicates. Fourteen gilts with a similar initial body weight (BW) from the low-RFI and high-RFI lines (n = 7 per line) were randomly selected and used for each replicate. Pigs were housed in individual metabolism crates, and had ad libitum access to water, but were restricted to feed intake. After a 9-day adaptation period, pigs within a line were randomly assigned to either a control (n = 6, three pigs per line) or LPS treatment (n = 8, four pigs per line) group. Pigs in the treatment group were then challenged with LPS via intramuscular injection of 30 μg/kg BW of LPS from E. coli O55:B5 dissolved in an endotoxin-free, sterile saline solution at baseline (0 hpi). Pigs in the control group were injected with an equivalent volume of endotoxin-free, sterile saline solution at the equivalent time. The rectal temperature of each pig was measured at 0, 2, 4, 6, and 24 hpi. At 0, 2, 6, and 24 hpi, blood samples were collected from the jugular vein into Tempus™ Blood RNA tubes for long-term storage at − 80 °C, into EDTA tubes for CBC tests and cytokine assays. For more details, see Materials and Methods. b Shown are the numbers of animals with blood samples collected from the two replicates at 0, 2, 6, and 24 hpi and the types of assays performed on different samples. Only blood samples collected from the LPS treated group from Replicate 2 were used for RT-qPCR and RNA-seq. Blood samples collected from all animals from both replicates were used for CBC tests and cytokine assays. The images of pigs were created by one of the co-authors, Anoosh Rakhshandeh and agreed to be published here. HRFI, high-RFI line; LRFI, low-RFI line RNA preparation Total RNA was extracted from the Tempus tube samples from Replicate 2 of the treatment groups for both pig lines by using preserved blood RNA purification kit I (Norgen Biotek Corp, Thorold, Ontario, CA) per the manufacturer’s instructions. On-column DNA digestion was performed using DNase I (Qiagen, Valencia, CA, USA). Globin transcripts (HBB and HBA) were depleted by following an RNase H-mediated method [[85]36]. The quantity and integrity of the RNA were monitored by using Nanodrop 2000 (Thermo Scientific, Waltham, MA, USA) and Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) before and after globin depletion. The efficiency of globin depletion of each sample was checked by conventional RT-qPCR with ACTB and GAPDH as the internal controls. Globin depletion efficiencies for all RNA samples were above 85%. Metadata, including RNA integrity numbers (RINs) and concentration of RNA post globin depletion, CBC, and sequencing batches are available in Additional file [86]1: Table S1. Complete blood count (CBC) tests and plasma cytokine assays CBC tests were performed for all except five clotted samples by the Pathology Lab, College of Veterinary Medicine at ISU as described [[87]10]. Eight cytokines (IFNα, IFNγ, IL1β, IL4, IL6, IL8, IL10 and TNFα) in the all 112 plasma samples were assayed by using Aushon SearchLight Arrays for pig cytokines (Aushon BioSystems, Billerica, MA, USA) per the manufacturer’s instructions. Verification of transcriptional response to LPS stimulation by RT-qPCR Primers for 47 candidate genes, which were porcine orthologs of human and murine genes responsive to LPS stimulation or important for Gram-negative sepsis control and resolution, and three house-keeping genes (ACTB, RPL32 and GAPDH), were designed and synthesized by Fluidigm Corporation (Fluidigm Corporation, San Francisco, CA, USA) such that two primers of each pair were separated by exon-exon boundaries and could amplify all known isoforms of the target gene, if possible. The specificities of the primer pairs were tested by conventional RT-qPCR using the DNA Engine Opticon 2 System (BioRad, Hercules, CA, USA) and only primer pairs that gave single peaks in melting curve analyses were kept. Additionally, the qPCR was performed without the melt curve analysis step and the amplicons were visualized on a 2% agarose gel for doublets, significant primer dimers, and confirmation of the amplicon size. Primer pairs producing amplicons of unexpected sizes were removed from the study. In total, this quality control scheme resulted in 36 and two primer pairs for genes of interest and internal controls (RPL32 and GAPDH) meeting the requirements, respectively (Additional file [88]2: Table S2). The 32 RNA samples from Replicate 2 of the treatment group were used for Fluidigm RT-qPCR without globin depletion. By following the Fluidigm user guide for Real-Time PCR analysis [[89]37], real time-qPCR was done on a 48.48 dynamic array chip (Fisher Scientific, Pittsburgh, PA, USA), along with reactions for assessing primer amplification efficiency, using the Biomarker HD system (Fluidigm Corporation, San Francisco, CA, USA). Data were analyzed with the Fluidigm Real-Time PCR analysis software using the default settings, to obtain raw C[t] (cycle of threshold) values. Since the expression levels of internal controls were not very stable for individual pigs during the time course, RT-qPCR data were analyzed by using the R package MCMC.qpcr (version 1.2.3) [[90]38]. With this method, internal reference genes are not mandatory but can be incorporated as Bayesian priors or as trackers of global effects when template abundances correlate with experimental conditions [[91]38]. Briefly, C[t] values were converted into the copy number of templates by incorporating the amplification efficiencies of the primers and then analyzed using a generalized linear mixed model, which assumed the copy number of the transcripts of a given gene follows a lognormal-Poisson distribution. In the generalized linear mixed model, effects for line, time, line-by-time interaction, and plate used for Fluidigm PCR were treated as fixed effects, while individual animal was included as a random effect. Bayesian MCMC prior distributions for fixed effects and the random effect were derived by using expression data of GAPDH and RPL32 during the 24-h time course. The p-values associated with the effects of line, time and line-by-time interaction were adjusted by using the Benjamini-Hochberg (BH) method [[92]39]. RNA-sequencing The 32 RNA samples used for the RT-qPCR assays above were also used for RNA-seq after depleting globin transcripts as described above. Library construction and sequencing were performed by the Beijing Genomics Institute (BGI, Hongkong, CN). Briefly, the RNA-seq libraries were constructed using the Illumina TruSeq RNA Sample Preparation Kit v2 (Illumina, San Diego, CA, USA) per the manufacturer’s instructions. Individual libraries were diluted to 2 nM and pooled in approximately equimolar amounts, with 8 libraries per pool. Paired-end sequencing (2 × 50 bases) was run on an Illumina Hiseq 2000 platform with one pool per lane. Quality control, preprocessing and alignment of RNA-seq reads Read quality was checked and filtered by BGI using their custom scripts. For a pair of reads, the whole pair was removed if either read met the following criteria: (1) either read had more than 50% of their bases aligned to the adapter sequences; (2) either read contained more than 10% of ‘N’ bases; (3) either read had more than 40% of bases with PHRED+ 64 quality scores lower than 20. The kept reads were aligned to the pig reference genome Sscrofa 11.1 (version 90, Ensembl) using 2-pass rna-STAR (version 2.5.3a) with the ENCODE standard option settings plus two explicit option settings: --outFilterMismatchNoverReadLmax 0.1 --outFilterIntronMotifs RemoveNoncanonical [[93]40, [94]41]. The resulting BAM files were further processed by using MMR to assign multi-mapper reads to their most likely locations such that the variances of local basewise coverage were minimized [[95]42]. Read counts per gene per library were summarized by using featureCounts (version 1.5.3) [[96]43] with explicit settings -d 30 -M, with the pig genome GTF file (version 90, Ensembl) as the genomic annotation reference file. Prior to differential expression analysis, hemoglobin genes (HBA and HBB) and genes with few reads were removed from the count table, such that only genes with count per million (cpm) mapped reads greater than one in at least four samples were kept. This analysis resulted in a final count table for 12,703 genes. This count table was used for subsequent differential expression analysis and clustering analysis after further transformation and adjustment (see below). Differential expression analysis of RNA-seq data Although in recent years, a few tools have been developed for time-course RNA-seq differential gene expression analysis [[97]44, [98]45], there is no generally accepted, applicable method to analyze RNAseq data from a short time-series experiment, with less than five time points and a small sample size as this study, in which within-individual measures are generally correlated. Therefore we were unable to take into account expected within-animal correlation in our differential expression analysis. We used the R/Bioconductor package DESeq2 (1.20.0) [[99]46] for differential expression analysis for two reasons: (1) DESeq2 adopts empirical Bayes shrinkage estimation for dispersions and fold changes, which improves stability and interpretability of estimates; (2) DESeq2 allows statistical tests of differential expression with a specified minimum effect size, which avoids the issue that post hoc filtering of differentially expressed genes (DEGs) based on a fold change threshold results in a false discovery rate (FDR) that is not easy to interpret [[100]46]. Additionally, known nuisance variables, such as RNA preparation batch, RIN, BW, and blood cell composition, could not well account for the treatment-unrelated variations of the complicated blood transcriptome. This is likely because several of them, such as the concentration of rare subtypes in peripheral blood (such as basophils, eosinophils, and monocytes), were not accurately measured. And more importantly, single cell RNA-seq data in human [[101]47] and pig peripheral blood mononuclear cells (PBMCs) (Crystal L. Loving, Haibo Liu et al., unpublished) revealed that the lymphocytes are composed of transcriptionally very heterogeneous subtypes of cells. Surrogate variable analysis has been shown to be a powerful method to detect and adjust for hidden variations in high throughput gene expression data [[102]48, [103]49]. Therefore, we decided to use surrogate variables estimated by the svaseq function of the R/Bioconductor package sva (v3.28.0) [[104]50] to account for the hidden nuisance variations in our RNA-seq data. Six surrogate variables were estimated by using a full model that included terms for an intercept, line, time and line-by-time interaction, and a reduced model that included only an intercept term. Differential expression analyses were conducted using DESeq2 as mentioned above. Briefly, a generalized linear model was fitted for each gene in the count table, with a negative binomial response and a log link that included a DEseq2 normalization offset and the effects of line, time and line-by-time interaction, and the six surrogate variables as estimated above. The nbinomWaldTest function was used to estimate and test the significance of regression coefficients with the following explicit parameter settings: betaPrior = FALSE, maxit = 5000, useOptim = TRUE, useT = FALSE, useQR = TRUE. Differentially expressed genes between conditions were identified by testing the significance of relevant contrasts and using the results function with the following explicit parameters: alpha = 0.01, lfcThreshold = log2(1.2), altHypothesis = “greaterAbs”, that is, testing whether the absolute values of the log[2] fold changes between conditions were greater than log[2](1.2). Thus, the estimates of the fold changes were shrunken by performing empirical Bayes shrinkage and their significances were tested by specifying a minimum effect size, which can improve the stability and interpretability of the estimates [[105]46]. Multiple testing correction was performed by using the BH method [[106]39]. Genes with absolute values of the log[2] fold change greater than log[2](1.2) and adjusted p values less than 0.05 were considered to be DEGs. Statistical analysis of body temperature, CBC and cytokine profile data Cytokine levels below the lower limit of detection were replaced by one half of the smallest positive values in the cytokine dataset. The CBC data and the imputed cytokine data were natural log-transformed for further analyses. Rectal temperature (all 140 data points), the transformed CBC data (107 data points) and cytokine data (all 112 data points) were analyzed by using the SAS PROC MIXED procedure (SAS Institute Inc., Cary, NC) for repeated measures analysis. The models used for analyzing CBC data and rectal temperature data contained RFI line (Line), sampling time (Time), treatment (Treatment), two-way and three-way interactions of the three former factors, body weight (BW) and age (AgeOfChallenge) at LPS injection, and pen by rooms (PenRoom) as fixed effects. For cytokine data analysis, plate used in the assays (Plate) was also included as a fixed effect, along with all the independent variables used for rectal temperature data analysis. Six commonly used residual correlation structures for repeated measures data analysis, compound symmetry (CS), heterogeneous compound symmetry (CSH), autoregressive (AR), heterogeneous autoregressive (ARH), spatial power (SP), and unstructured (UN) correlation structures, were considered. For each response variable, the correlation structure giving the smallest Akaike information criterion (AIC) was selected for the final analysis. The goodness of fit of the models was assessed as descried [[107]51]. Type III tests were performed for fixed effects and contrasts were adjusted using the Kenward-Roger method [[108]52]. Clustering of gene expression profiles The filtered RNA-seq count data were transformed to log (cpm) by using the voom function of the limma package (v3.36.2) [[109]53, [110]54]. The transformed gene expression levels were then adjusted for the nuisance variables, i.e., the six surrogate variables used in the model for RNA-seq differential expression analysis. The adjusted log (cpm) for all individuals of the two lines were altogether used as input for the software STEM (v1.3.8), which is a tool specifically developed for short time-series expression data mining [[111]55]. For parameter settings, see Supplementary Methods (Additional file [112]3). To test how likely the observed profiles were created at random, permutation tests were performed (For details, see Additional file [113]3). Expression profiles with more genes assigned than expected based on a null distribution derived from permutation, where gene expression values at 2, 6 and 24 hpi were randomly permuted within each animal for 50 iterations and used for STEM analysis (q < 0.05), were considered as significant profiles. Significant expression profiles were clustered together if the correlation coefficients of two mean profiles were no less than 0.6. GO term overrepresentation analysis and visualization were performed for genes in each cluster of profiles with a q value cutoff of 0.05 by using built-in functionality in STEM. Gene ontology term overrepresentation analysis (GOA) Gene ontology (GO) annotation for pig genes was downloaded from Ensembl BioMart (Release 90). GO terms associated with less than 10 or more than 500 genes were excluded. For the RNA-seq data, 12,703 Ensembl Gene IDs with detectable expression in the blood samples were used as the background references. Hypergeometric tests of overrepresentation of GO