# Small networks of expressed genes in the whole blood and relationships to profiles in circulating metabolites provide insights in inter-individual variability of feed efficiency in growing pigs

Camille Juigné, Emmanuelle Becker, Florence Gondret

<figure><img src="image_1.png" /></figure>

## To cite this version:

Camille Juigné, Emmanuelle Becker, Florence Gondret. Small networks of expressed genes in the whole blood and relationships to profiles in circulating metabolites provide insights in inter-individual variability of feed efficiency in growing pigs. 2023. hal-04112110

HAL Id: hal-04112110

https://hal.science/hal-04112110

Preprint submitted on 31 May 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Small networks of expressed genes in the whole blood and
relationships to profiles in circulating metabolites provide insights
in inter-individual variability of feed efficiency in growing pigs
(preprint)

Camille Juigné<fnref n="1" />,<fnref n="2" />, Emmanuelle Becker<fnref n="2" />, and Florence Gondret<fnref n="1" />

**Abstract**

**Background:** Feed efficiency is a research priority to support a sustainable meat production. It is recognized as a complex trait that integrates multiple biological pathways orchestrated in and by various tissues. This study aims to determine networks between biological entities to explain inter-individual variation of feed efficiency in growing pigs.

**Results:** The feed conversion ratio (FCR), a measure of feed efficiency, and its two component traits, average daily gain and average daily feed intake, were obtained from 47 growing pigs from a divergent selection for residual feed intake and fed high-starch or high-fat high-fiber diets during 58 days. Datasets of transcriptomics (60 k porcine microarray) in the whole blood and metabolomics (1H-NMR analysis and target gas chromatography) in plasma were available for all pigs at the end of the trial. A weighted gene co-expression network was built from the transcriptomics dataset, resulting in 33 modules of co-expressed molecular probes. The eigengenes of eight of these modules were significantly ($P \le 0.05$) or tended to be ($0.05 < P \le 0.10$) correlated to FCR. Great homogeneity in the enriched biological pathways was observed in these modules, suggesting co-expressed and co-regulated constitutive genes. They were mainly enriched in genes participating to immune and defense-related processes. They were also generally associated with growth rate and percentage of lean mass. In the network, only one module composed of genes participating to the response to substances, was significantly associated with daily feed intake and body adiposity. The profiles in circulating metabolites and fatty acids of plasma were summarized by weighted linear combinations using a dimensionality reduction method. Close association was notably found between a module composed of co-expressed genes participating to T cell receptor signaling and cell development process in the whole blood and related to FCR, and the circulating concentrations of omega-3 fatty acids in plasma.

**Conclusion:** These systemic approaches have highlighted networks of entities driving key biological processes involved in the phenotypic difference in feed efficiency between animals. Connecting transcriptomics and metabolic levels together had some additional benefits.

**Introduction**

In the context of various geopolitical tensions and societal questions on the agri-agro-food systems, feed efficiency (FE) is a research priority to support food security and a sustainable meat production. Indeed, better FE is associated with a reduced amount of feeds needed for production and lower environmental wastes and emissions. Feed efficiency is measured on farms by the feed conversion ratio (FCR), an index calculated as the ratio of feed intake to body weight (BW) gain during a test period. Residual feed intake (RFI) has been also proposed as a refined measure of FE for genetic selection; it is calculated as the difference between observed feed intake and predicted feed intake from production and maintenance needs, which allows RFI to

<footnote n="1">¹ PEGASE, INRAE, Institut Agro, F-35590 Saint Gilles, France</footnote>

<footnote n="2">² Univ Rennes, Inria, CNRS, IRISA - UMR 6074, F-35000 Rennes, France</footnote>
reflect digestive and metabolic variabilities [1, 2]. A number of studies have been performed in the past years to depict and understand the biological bases of FE. Based on the comparison of animals with low or high FCR or RFI, they all concluded to the complex nature of FE since this trait integrates multiple biological pathways orchestrated in and by various tissues [2, 3, 4]. Among tissues, peripheral blood is a convenient and relatively easy sampling source of biological information that can highlight the variations in tissues metabolism and physiology to understand complex phenotypes [5] and with potential outcomes for applications in diagnostics and selection [6]. In pigs, about one thousand genes were found differentially expressed in the whole blood between two lines of pigs divergently selected for RFI [7, 8]. Gene set enrichment analysis on the whole blood transcriptome in beef cattle has also allowed identifying biological pathways associated with a divergent selection for low or high RFI [5]. Moreover, we recently showed that gene expression profiling in the whole blood is suitable to identify a few number of molecular candidate biomarkers for FCR in growing pigs, when gene expression levels were analyzed by machine learning algorithms based on classification or regression trees [9]. Likewise, metabolomic studies have shown that circulating concentrations of metabolites in the blood can be related to economically-important traits including FE [10]. However, all these approaches did not address the inter-individual variation in FE traits and did not attend to depict the interactions among the biological entities at a given level or different levels of omics organization. Therefore, a systemic approach considering the multiple relationships between molecules can add new insights in the architecture of complex traits such as FE.

Among various systems biology approaches, the weighted gene co-expression network analysis (WGCNA) has been proposed as a suitable method for defining interactions between transcripts of genes, by grouping them in modules of pairwise correlations to reveal the higher-order organization of the transcriptome [11]. Based on RNA-sequencing data in the liver, two co-expression networks were identified as associated with high or low FE in dairy cows [12]. Another approach based on linear models allowed to combine gene expression data and high throughput metabolomics data in skeletal muscle [13] ; in this study, pairs of metabolite-transcript associated with sphingolipid catabolism, multicellular organismal process, and purine metabolic processes were associated with differences in FE between two pig breeds and between two groups of pigs of low or high FE values.

The aim of the present study was to depict the biological bases underlying inter-individual variability in FE and other related traits in growing pigs by identifying small networks of interconnected gene transcripts in the whole blood and their relationships with global profiles of circulating metabolites. Eight molecular modules mainly composed of interconnected genes involved in immune and defense-related processes, were related to variability between pigs in FCR. Other important biological processes represented in these gene networks were the response to organic substance, ribosome biogenesis and translation, and cell development and learning, respectively. One module of inter-connected expressed genes related to immune process was associated with circulating concentrations of omega3 fatty acids in the plasma, thus connecting transcripts to metabolites in the determinism of variability in FE.

Material and methods

**Ethics** We reused transcriptomics and metabolomics datasets acquired in the whole blood from purebred French Large White pigs produced in a divergent selection experiment for RFI [7, 14], and that have been previously analyzed separately and for the line-associated effects only. The animal phenotypes have been published by [15]. These data were complemented by data on circulating concentrations of fatty acids (FA) in the plasma of the same pigs to specify lipid-related processes, that have not been previously published. In the original publications, the care and use of pigs were performed in compliance with the European Union legislation (directive 2010/63/EU). The protocol was approved by an Ethics Committee in Animal Experiment (Comité Rennais d'Ethique en matière d'Expérimentation Animale, CREEA N°007, agreement N°07-2012). All animals were reared and killed in compliance with the national regulations and according to procedures approved by the French veterinary Services. All methods were reported in accordance with ARRIVE guidelines. In the present study, reusing published data to perform new analyses for different animal traits perfectly fits with the 3R (Replacement, Reduction and Refinement) principles.
**Origin of phenotypic data** Full description of the experimental design that provided the original datasets is referenced by [15]. Briefly, data were obtained from a total of 48 growing pigs of two lines in the 7th generation of divergent selection for RFI, and fed diets formulated at isocaloric and isoproteic bases but differing in energy source and nutrients (lipids and fibers vs. starch), were tested from $74d \pm 0.3d$ of age to $132.5d \pm 0.5d$ (SEM) of age. From this publication, we considered body weight at slaughter (BW in kg), age at slaughter (in days), ADFI (in g/day), ADG (in g/day), FCR calculated as the ratio between the average daily feed intake ADFI and average daily gain ADG during the feeding trial, and the percentage (relative to carcass weight) of the dorsal subcutaneous adipose tissue (%backfat) and of the loin muscle cut (%loin) as surrogates of body composition.

**Transcriptomic dataset** The transcriptomic data were retrieved from NCBI's Gene Expression Omnibus (GEO) Subserie accession number GSE70838 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70838). These data have been obtained from the whole blood collected in pigs at $132.5d \pm 0.5d$ (SEM) of age, by using a custom porcine microarray ($8 \times 60K$, GPL16524, Agilent Technologies France, Massy, France) containing 60,306 porcine probes. Full description of the methods to produce the raw microarray data can be found in [7]. After quality filtration based on four criteria (background intensity value, diameter, saturation and uniformity of the spot), the original dataset contained 26,322 annotated probes expressed in the whole blood. There were approximately 2.2 replicates per unique gene in the transcriptomic dataset (min: 1; max: 33).

**Metabolomic dataset** The metabolomic dataset consisted in high resolution 1H-NMR spectra generated in plasma of the 48 pigs, and was retrieved from Jegou et al. [14]. A total of 94 buckets, and consequently of 94 variables, was considered for dedicated analyses in the present study. Buckets are individual metabolites or groups of two metabolites. They were notably assigned to different amino acids, creatine, lactate and glucose.

When using the 1H-NMR approach, lipids were all grouped as a single spectrum. In the present study, to investigate more deeply the role of the lipid family in FE variability, the circulating concentrations of fatty acids (FA) in plasma were newly determined in the 48 pigs. For that, lipids were extracted from plasma as previously described for tissues [15], and methylated. Dedicated analyses were then performed on a gas chromatograph (Nelson Analytical, Manchester, NH) equipped with a fused-silica capillary column ($30m \times 0.25mm$ internal diameter), with a base-deactivated silica stationary phase (a $0.25-\mu m$ film thickness) filled with a stationary phase (80% biscyanopropyl and 20% cyanopropylphenyl) and using margaric acid (C17:0) as the internal standard. The furnace temperature was $180^{\circ}C$, and injector and detector temperatures were $240^{\circ}C$. Retention times and peak areas were determined. Peaks were identified by comparison with the retention times of standard FA methyl esters. Individual FAs were then quantified as percentages of the sum of FA identified in each sample. To facilitate the biological interpretation, FA have been then grouped into families of saturated FA with 14C or less (sC14:0), polyunsaturated FA of the omega-6 family (ssn-6c), and polyunsaturated FA of the omega-3 family (ssn-3). All other identified FA were kept as these. This led to a total of 14 variables representing circulating FA concentrations for the 48 pigs (Supplementary Table 1).

**Construction of the weighted gene co-expression networks**

Starting from a matrix whose individuals are pigs and features are probes expression levels (48 pigs $\times$ 26,322 probes), we performed a hierarchical clustering to identify outlier individuals as recommended [16]. One pig was detected as an outlier and further removed from the dataset due to aberrant values. Then, we quantified the number of probes that were significantly linked to the animal phenotypic traits of interest (linear models with a *P*-value < 0.05 as cut-off). The results of the linear regressions showed that age at slaughter and line were the factors affecting the most the expression levels of molecular probes (19% of the probes significantly affected by age at slaughter and 14% by the line effect) (Supplementary Table 2). For next steps of analysis, we used the residuals of these linear regressions for the age effect, but preserved the intra- and inter-line variability for FE.
Corrected probe expression levels were then filtered. Only probes *i* whose log fold-change $FC_i$ was greater than 1 were selected, with

$$log_2(FC_i) = max(log_2(expression_i)) - min(log_2(expression_i)) \quad (1)$$

Indeed, we considered that the variant probes formed a reduced dataset in which connected probes can be more easily found when reducing noise and spurious associations. A smaller dataset is also associated with a less computationally demanding analysis for analyzing gene networks. This reduced dataset included 16,190 probes for the 47 pigs.

To calculate the co-expression network, we used the Weighted Gene Co-expression Network Analysis (WGCNA) step-by-step method [11], performed with R 4.2.2. The first step consisted in calculating a measure of co-expression similarity $s_{i,j}$ between each pair of probes to highlight the pairs of probes whose expression varies in a similar way. The adjacency matrix $A = [a_{i,j}]$ was then constructed by raising the co-expression similarity measure $s_{i,j}$ to the power $\beta$, using the signed hybrid method (only positive correlations were kept) [equation 2].

$$signed \ hybrid \ a_{ij} = \begin{cases} cor(x_i, x_j)^\beta & cor(x_i, x_j) > 0 \\ 0 & cor(x_i, x_j) < 0 \end{cases} \quad (2)$$

where $a_{ij}$ is the element $(i, j)$ of the adjacency matrix $A$, $\beta$ is the soft-thresholding, $x_i$ is the level expression of the ith probe, and $cor(x_i, x_j)$ is the Pearson correlation between expression profiles of the i and j probes.

$\beta$ is a non-dichotomic soft thresholding that allows to evaluate connection between probes without losing the continuous character of the co-expression. Low correlations are better masked with high $\beta$ values. We set $\beta$ to 6 according to the criterion of the “approximate scale free topology” [17].

Considering probes as nodes, a weighted co-expression network can be then deduced from the adjacency matrix, by adding edges of weight $a_{ij}$ between pairs of probes whose weight is strictly positive.

## Detection of modules of co-expressed probes and their relationships with animal phenotypic traits

To detect modules in this network of co-expressed probes, a proximity measure was calculated using the Topological Overlap Measure [18] and a hierarchical clustering was performed. This results in a dendrogram, in which clusters can be detected using the dynamic TreeCut algorithm [11] and a sensitivity threshold set to 2. To keep modules of highly co-expressed genes, we set the minimum module size to 25 probes. Each module was then summarized by its first principal component called the module eigengene (ME), that is a mathematical solution to condense the expression profile of the probes in the module.

A heatmap of correlations between ME and phenotypic traits was then produced. The heatmap can be examined to find the most significant associations. In this study, we considered $P \le 0.05$ as significant association and $0.05 < P \le 0.10$ as a tendency.

We also calculated the gene significance (GS) as (the absolute value of) the correlation between the probe and the phenotypic trait. For each module, we defined the quantitative measure of module membership (MM) as the correlation of ME and the gene expression profile. Using both the GS and MM measures, we can identify genes that have a high significance for FE traits (central players) and high module memberships in the modules. Finally, for a subset of modules, we provide a graphical representation of the sub-networks considering the annotated probes only. From the adjacency matrix, pairs of probes with correlation coefficients greater than the 95th percentile in a given module were selected, and represented using Cytoscape with nodes corresponding to probes and edges corresponding to the adjacency matrix values. For probes with low correlation with other probes, edges were represented in a different color.

## Biological functional enrichment in modules of co-expressed probes

For modules that were significantly related to FCR based on the heatmap examination between module eigen-genes (ME) and animal traits, we matched their constitutive molecular probes to the corresponding unique
genes (official gene symbol), using the annotation provided by the manufacturer of the expression microarray. When applicable, the gene ontology (GO) terms for biological processes were then automatically searched in each module, using the Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics knowledgebase v2022q4 released (http://david.abcc.ncifcrf.gov/). The GOBP terms FAT were selected to filter the broadest terms. The results were downloaded using the "Functional annotation clustering" option of the DAVID tool. Only clusters of terms with an enrichment score (measured by the geometric mean of the EASE score of all enriched annotations terms) greater than 1.2 were considered. Within each cluster, the top GO term was listed together with its own enrichment score and the associated modified Fisher exact P-value.

## Establishing profiles of circulating metabolites and evaluating connections between metabolic and transcriptomic levels

The second and third datasets considered in this study encompassed the circulating concentrations of metabolites and of FA, respectively. To reduce the dimension and facilitate correlation analyses with gene expression networks, these datasets were first summarized by Principal Component Analysis (PCA) using the R packages FactoMineR and factoextra [19, 20]. This would avoid bias when considering hundreds metabolic variables with only few modules of highly correlated genes. A PCA was used to summarize the profile in metabolites identified after 1H-NMR analysis (94 variables) and another PCA was used for circulating concentration of FA (14 variables). The first components of PCA were called Dimi\_Metab for the metabolomic 1H-NMR table and Dimi\_FA for the FA table, respectively, with *i* = 1 to 5. These dimensions were linear combinations of the original variables.

To connect information at the two omics levels, the dimensions of each PCA were then correlated to the eigengenes of the WGCNA modules (ME) by using Pearson correlations. The correlations were represented by heatmaps to facilitate the description.

## Results

Data obtained in a total of 47 growing pigs with inter-individual differences in FCR (i.e., the measure of FE on farms) due to genetic selection for RFI and to the diet received during a test period of 58 days, were considered in the present study. In addition to FCR, the average daily gain (ADG) and average daily feed intake (ADFI) (i.e., the two components of FCR), and body composition estimated by percentage (relative to carcass weight) of backfat (%backfat) and of the loin cut (%loin) were also obtained (Supplementary Table 1).

## Definition of gene co-expression network in the whole blood of pigs

A network was built with the WGCNA package from the expression levels of 16,190 molecular probes expressed in the whole blood, where nodes correspond to the expression profile of the molecular probes, and edges are determined by the pairwise correlations between probes expression (the adjacency matrix is available at https://data-access.cesgo.org/index.php/s/YPz0J2ItxIEuN5M). This network was then analysed to find modules defined as groups of co-expressed probes that may represent the molecular architecture behind the animal phenotypic traits.

After excluding the grey module which is used to hold all the probes that do not clearly belong to any other modules, we show that the co-expression network was composed of 33 modules composed of 27 to 3,829 molecular probes (Supplementary Table 3). Annotations were used when applicable to identify the corresponding genes. The distribution in the number of probes and their corresponding unique genes per module is shown in Figure 1.

The module eigengene (ME) was then calculated as the representative of expression profiles of the genes in the module. The module membership (MM) was calculated by correlating ME and the expression profile of the probes within each module. By definition, the closer to 1 or -1 is MM, the higher is the gene connected
to ME. The medians of the MM values indicated a satisfactory clustering from the whole network of the molecular probes expressed in the whole blood of the 47 pigs. The values are available in Supplementary Table 3.

## Relationships between modules of co-expressed genes and animal phenotypic traits

The heatmap of the correlation coefficients between the module eigengene (ME) of each WGCNA module and the animal phenotypic traits is presented in Figure 1. For six modules, the ME was significantly correlated ($P-value \le 0.05$) with Feed Conversion Ratio (FCR): they were the modules violet (probes: 32, unique genes: 16, $P-value = 8e-04$), darkred (probes: 82, unique genes: 53, $P-value = 2e-04$), royalblue (probes: 111, unique genes: 81, $P-value = 0.04$), lightcyan (probes: 180, unique genes: 114, $P-value = 0.008$), white (probes: 59, unique genes: 28, $P-value = 0.04$), and darkorange (probes: 60, unique genes: 21, $P-value = 0.05$). For two other modules, ME tended ($P-value < 0.1$) to be correlated with FCR: they were the modules darkolivegreen (probes: 27, unique genes: 13, $P-value = 0.08$) and steelblue (probes: 45, unique genes: 23, $P-value = 0.09$).

For the modules violet, darkred, royalblue, lightcyan, and white, the correlations between ME and ADG were also significant, and for the module darkorange, there was a trend for correlation between ME and ADG. As expected, the signs of correlation between ME and ADG and between ME and FCR (which is the ratio between ADFI and ADG during the test period), were opposite. However, for the modules darkolivegreen and steelblue, the correlation coefficients between ME and ADG did not reach statistical significance.

Four of the eight modules associated with FCR also displayed significant correlations with %loin (with opposite signs of correlation): the modules violet, darkred, lightcyan and white. For these modules, there was no significant correlation between ME and %backfat.

Finally, none of the eight modules associated with FCR were significantly related to ADFI. In the network, only the darkgreen module (probes: 82, unique genes: 43, $P-value = 8e-04$) was highly correlated with ADFI; it was also significantly related with %loin ($P-value = 0.004$) and %backfat ($P-value = 2e-04$). The saddlebrown module (probes: 46, unique genes: 12) tended to be correlated with ADFI ($P-value = 0.1$) and with %backfat ($P-value = 0.1$). Finally, the eigengenes of the modules green (probes: 1, 113, unique genes: 412, $P-value = 0.07$), brown (probes: 1, 441, unique genes: 773, $P-value = 0.02$) and skyblue (probes: 46, unique genes: 45, $P-value = 0.05$) displayed significant associations with %loin, but without significant correlations with any other animal phenotypic traits.

## Close-vicinity of the different modules of co-expressed genes

To evaluate the connectivity between modules in the network, a hierarchical clustering was performed between the eigengenes (ME) of the 34 modules. The resulting dendrogram is shown in Figure 2, in which the eight modules that were significantly associated or tending to be associated with FCR are enlightened. Among these eight modules, two clusters were identified. The first cluster associated the lightcyan (114 genes), steelblue (23 genes) and darkolivegreen (13 genes) modules. The second cluster associated the darkred (53 genes), violet (16 genes) and royalblue (81 genes) modules. The white (28 genes) and the darkorange (21 genes) modules were isolated in the dendrogram. In addition, the green (411 genes) and brown (772 genes) modules, that were significantly associated with %loin, were clustered together (Figure 2).

## Functional enrichment of the modules in biological processes

For the eight modules identified as significantly correlated or tending to be correlated with FCR, an enrichment analysis was performed to find the main biological processes shared by the co-expressed gene transcripts within each module (Table 1). Annotations of the probes were first retrieved, and the corresponding gene name was associated to each probe when applicable. The DAVID tool was used on the gene list uploaded for each module.
<table><caption>Table 1: Functional enrichment in biological pathways for molecular modules related to FCR and body composition.</caption><thead><tr><th>Module</th><th>probes</th><th>unique genes</th><th>match ID DAVID</th><th>GO/terms</th></tr></thead><tbody><tr><td>darkolivegreen</td><td>27</td><td>13</td><td>9</td><td>GO:0006954 inflammatory response E = 15.6 P &lt; 0.05</td></tr><tr><td>darkorange</td><td>60</td><td>21</td><td>15</td><td>GO:0010033 response to organic substance E = 4.4 P &lt; 0.001, GO:0002376 immune system process E = 2.9 P &lt; 0.05</td></tr><tr><td>darkred</td><td>82</td><td>53</td><td>44</td><td>GO:0042273 ribosomal large subunit biogenesis E = 18.8 P &lt; 0.01, GO:0002181 cytoplasmic translation E = 16.5 P &lt; 0.01</td></tr><tr><td>lightcyan</td><td>180</td><td>114</td><td>96</td><td>GO:0007166 cell surface receptor signaling pathway E = 2.3 P &lt; 0.001, GO:0032502 developmental process E = 1.5 P &lt; 0.001, GO:0006955 immune response E = 2.74 P &lt; 0.001, GO:0006909 phagocytosis E = 4.5 P &lt; 0.01, GO:0042113 B cell activation E = 9.8 P &lt; 0.001, GO:0070887 cellular response to chemical stimulus E = 1.7 P &lt; 0.01, GO:0045088 regulation of innate immune response E = 3.4 P &lt; 0.05, GO:0050727 regulation of inflammatory response E = 3.1 P &lt; 0.05</td></tr><tr><td>royalblue</td><td>111</td><td>81</td><td>71</td><td>GO:0050852 T cell receptor signaling pathway E = 21.1 P &lt; 0.001, GO:0030155 regulation of cell adhesion E = 5.6 P &lt; 0.001, GO:0048869 cellular developmental process E = 2.0 P &lt; 0.001</td></tr><tr><td>steelblue</td><td>45</td><td>23</td><td>10</td><td>GO:0008104 protein localization E = 3.9 P &lt; 0.05</td></tr><tr><td>violet</td><td>32</td><td>16</td><td>12</td><td>GO:0045087 innate immune response E = 9.4 P &lt; 0.001, GO:0016567 protein ubiquitination E = 7.5 P &lt; 0.05, GO:0006952 defense response E = 6.0 P &lt; 0.001</td></tr><tr><td>white</td><td>59</td><td>28</td><td>23</td><td>GO:0072359 circulatory system development E = 4.4 P &lt; 0.01, GO:0007612 learning E = 16.9 P &lt; 0.01, GO:0030036 actin cytoskeleton organization E = 5.0 P &lt; 0.05</td></tr></tbody></table>

The DAVID tool was used to identify the top enriched pathways across the list of unique annotated genes within each module. The gene ontology (GO) terms for biological processes are indicated together with enrichment score (E) of the process and Fisher P value.

The *lightcyan* which was the biggest module correlated to FCR (180 probes corresponding to 114 unique genes), displayed a large number of different clusters of biological processes. On the opposite, for the other seven modules considered, there was/were only one to three clusters (E > 1) identified among the GO terms in each module (Table 1). This indicates a good consistency of the biological processes shared by the intra-connected genes within each module. The modules *violet*, *royalblue*, *darkorange*, *lightcyan* and *darkolivegreen* showed a predominance of immune, inflammatory and defense-related pathways across their constitutive genes. The *darkorange* module also included genes involved in the response to organic substance. The *darkred* module was rather oriented towards ribosome biogenesis and the process regulating translation. The *white* module was related to circulatory cell development and to learning.

The biological processes identified in the other modules related to feed intake (ADFI) or body composition (%loin and %backfat), but not to FCR, were described in Table 2.

The module *saddlebrown* that tended to be correlated to ADFI and %backfat was composed of 16 unique genes participating to the response to stimulus. The *skyblue* module that was significantly correlated to muscle mass (%loin), showed a predominance of genes related to protein metabolism among its 44 unique constitutive genes (Table 2). The *green* module was a large module of 411 unique genes related to various processes such as epigenetics processes (chromatin organization, histone methylation), to cellular responses to compounds (nitrogen, organic cyclic) for stress and defense, and to circadian regulation. The *brown* was also a big module of 772 unique genes that regulated cell differentiation, protein and lipid processes, and
<table><caption>Table 2: Functional enrichment of molecular modules correlated with feed efficiency or body composition.</caption><thead><tr><th>Module</th><th>probes</th><th>unique genes</th><th>match ID DAVID</th><th>GO terms</th></tr></thead><tbody><tr><td>brown</td><td>1,441</td><td>772</td><td>623</td><td>GO:0009966 regulation of signal transduction E = 1.4 P &lt; 0.001,<br>GO:0050790 regulation of catalytic activity E = 1.4 P &lt; 0.001,<br>GO:0007169 transmembrane receptor protein tyrosine kinase signaling pathway E = 1.9 P &lt; 0.001,<br>GO:0036211 protein modification process E = 1.3 P &lt; 0.001,<br>GO:0042886 amide transport E = 1.4 P &lt; 0.001,<br>GO:0007167 enzyme linked receptor protein signaling pathway E = 1.7 P &lt; 0.001,<br>GO:0031400 negative regulation of protein modification process E = 1.9 P &lt; 0.001,<br>GO:0042692 muscle cell differentiation E = 1.8 P &lt; 0.01,<br>GO:0071363 cellular response to growth factor stimulus E = 1.9 P &lt; 0.001,<br>GO:0031331 positive regulation of cellular catabolic process E = 1.7 P &lt; 0.05,<br>GO:0008654 phospholipid biosynthetic process E = 2.5 P &lt; 0.01,<br>GO:1903320 regulation of protein modification by small protein conjugates or removal E = 2.2 P &lt; 0.01,<br>GO:0045595 regulation of cell differentiation E = 1.4 P &lt; 0.01,<br>GO:0034284 response to monosaccharide E = 2.1 P &lt; 0.01</td></tr><tr><td>darkgreen</td><td>82</td><td>42</td><td>33</td><td>GO:0009636 response to toxic substance E = 9.9 P &lt; 0.01,<br>GO:0032496 response to lipopolysaccharide E = 7.9 P &lt; 0.01,<br>GO:0051128 regulation of cellular component organization E = 2.2 P &lt; 0.05,<br>GO:0042063 gliogenesis E = 8.9 P &lt; 0.01</td></tr><tr><td>green</td><td>1,113</td><td>411</td><td>356</td><td>GO:0006325 chromatin organization E = 2.0 P &lt; 0.001,<br>GO:0009891 positive regulation of biosynthetic process E = 1.7 P &lt; 0.001,<br>GO:0034968 histone lysine methylation E = 5.3 P &lt; 0.001,<br>GO:0033555 multicellular organismal response to stress E = 4.8 P &lt; 0.001,<br>GO:1901699 cellular response to nitrogen compound E = 2.1 P &lt; 0.001,<br>GO:0032922 circadian regulation of gene expression E = 5.7 P &lt; 0.01,<br>GO:0006622 protein targeting to lysosome E = 9.37 P &lt; 0.001,<br>GO:0016050 vesicle organization E = 2.4 P &lt; 0.001,<br>GO:0071407 cellular response to organic cyclic compound E = 2.0 P &lt; 0.01,<br>GO:0016071 mRNA metabolic process E = 2.0 P &lt; 0.001,<br>GO:0009611 response to wounding E = 2.1 P &lt; 0.01</td></tr><tr><td>saddlebrown</td><td>46</td><td>16</td><td>13</td><td>GO:0006357 regulation of transcription E = 3.72 P &lt; 0.01,<br>GO:0009628 response to abiotic stimulus E = 9.6 P &lt; 0.001</td></tr><tr><td>skyblue</td><td>46</td><td>44</td><td>36</td><td>GO:0006807 nitrogen compound metabolic process E = 1.80 P &lt; 0.001,<br>GO:0045184 establishment of protein localization E = 2.7 P &lt; 0.01,<br>GO:0034660 ncRNA metabolic process E = 4.8 P &lt; 0.05</td></tr></tbody></table>

The DAVID tool was used to identify the top enriched pathways across the list of unique annotated genes within each module. The gene ontology (GO) terms for biological processes are indicated together with enrichment score (E) of the process and Fisher P value.
<table><caption>Table 3: Top genes in the molecular modules related to feed conversion ratio (FCR).</caption><thead><tr><th>darkolivegreen</th><th>darkorange</th><th>darkred</th><th>lightcyan</th><th>royalblue</th><th>steelblue</th><th>violet</th><th>white</th></tr></thead><tbody><tr><td>PTGER3</td><td>BMP6</td><td>SLCO2B1</td><td>ITGAD</td><td>POFUT1</td><td>UQCC2</td><td>SLPI</td><td>IGDCC3</td></tr><tr><td>TUBB6</td><td>PADI2</td><td>E4</td><td>DDC</td><td>LPAR3</td><td></td><td>P2RY1</td><td>TMEM14C</td></tr><tr><td>SLA-DRB1</td><td>DNAJB9</td><td>CCR3</td><td>HTRA1</td><td>CCR7</td><td></td><td>MUC4</td><td>KIAA0247</td></tr><tr><td></td><td>WBSCR27</td><td>SLC46A2</td><td>EBF1</td><td>PTTG1</td><td></td><td>MED8</td><td>HRH1</td></tr><tr><td></td><td>TAOK3</td><td>EIF1B</td><td>CBR3</td><td>NIPSNAP3B</td><td></td><td>HSPA1B</td><td>HTR7</td></tr><tr><td></td><td>SERHL2</td><td>ZFAND6</td><td>MYBL1</td><td>EPHB6</td><td></td><td>TR10D</td><td>AFF1</td></tr><tr><td></td><td></td><td>CD300C</td><td>C4BPB</td><td>SLC4A11</td><td></td><td>HSP70.2</td><td>ZUFSP</td></tr><tr><td></td><td></td><td>FAM102A</td><td>KIAA0556</td><td>STRN</td><td></td><td></td><td>LDB3</td></tr><tr><td></td><td></td><td>RPL14</td><td>MS4A1</td><td>ZCCHC10</td><td></td><td></td><td>PROX1</td></tr><tr><td></td><td></td><td>WWP1</td><td>SLA-DOA</td><td>FMNL3</td><td></td><td></td><td>ACVR2B</td></tr><tr><td></td><td></td><td>KRTCAP3</td><td>HMG20A</td><td>SKAP1</td><td></td><td></td><td>DGKA</td></tr><tr><td></td><td></td><td>POLB</td><td>CYAC3</td><td>NPY</td><td></td><td></td><td>LUM</td></tr><tr><td></td><td></td><td>PIGL</td><td>RALGPS2</td><td>PPP1R26</td><td></td><td></td><td>TBC1D19</td></tr><tr><td></td><td></td><td>TMEM52B</td><td>ELL3</td><td>IZUMO4</td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>PLA2G12A</td><td>PIKFYVE</td><td>NIPSNAP3A</td><td></td><td></td><td></td></tr></tbody></table>

The unique genes corresponding to the annotated probes were listed in each module according to their GS.FCR value. In the table, only genes with a value greater than 0.3 for GS.FCR are indicated. The GS.FCR is the correlation between Gene Significance (GS) of the module eigengene and FCR.

signaling pathways.

## Hierarchy of expressed genes in the modules related to feed efficiency traits

To determine which expressed genes accounted the most in the correlations between the module eigengene (ME) and FCR, we calculated the Gene Significance (GS) in each module, and expressed the GS as a function of FCR. The top genes are listed in Table 3 and the all data are provided in Supplementary Table 3.

For the white module, IGDCC3, TMEM14C, HRH1, HTR7 and AFF1 were notably found as top genes. For the violet module, SLPI, P2RY1, MUCA, MED8, HSPA1B and HSP70.2 were the most important genes triggering the correlation of the module with FCR. For the royalblue module, POFUT1, LPAR3, CCR7, PTTG1, STRN, NPY and PPAR26 were pointed as important in the correlation with FCR. For the darkred module, EIF1B, RPL14 and KRTCAP3 were among the top 15 genes. For the lightcyan module, ITGAD, DDG, HTRA1, EBF1 and CBR3 were notably listed. A graphical representation of the importance of the genes in the modules is also provided. For the royalblue (Figure 3), this shows that the probes that were highly correlated to other probes were listed in the top list of probes based on their module membership (MM).

## Metabolic profiles in the whole blood

The second part of the present study addressed the metabolic level by considering several variables obtained in the circulating blood of the 47 pigs, after non-target analysis (1H-NMR spectra) of plasma or specific analysis of the lipidic fraction (gas chromatography-derived information) of plasma.

The relationships between the 94 variables obtained from the 1H-NMR spectra were summarized by the first five dimensions of a principal component analysis (PCA). We showed that the two first dimensions represented 41.4% and 26.8% of the total variability, respectively. The Figure 4 shows the corresponding correlation circle. The first dimension of this PCA mainly opposed lactate to the majority of the identified amino acids (lysine, tyrosine, valine, phenylalanine, methionine, leucine, isoleucine, glutamine-glutamate) and to high density lipoproteins (HDL). The second dimension mainly opposed glucose, eventually combined
<table><caption>Table 4: Relationships between transcriptomic and metabolomic levels.</caption><thead><tr><th colspan="6">Dimensions of the principal component analysis of the metabolic level in blood</th></tr><tr><th>Modules</th><th>Dim1 met</th><th>Dim2 met</th><th>Dim3 met</th><th>Dim4 met</th><th>Dim5 met</th></tr></thead><tbody><tr><td>Significantly correlated (Pvalue ≤ 0.05)</td><td>1</td><td>2</td><td>8</td><td>0</td><td>11</td></tr><tr><td>Trend (Pvalue ≤ 0.1)</td><td>3</td><td>5</td><td>6</td><td>7</td><td>6</td></tr></tbody><thead><tr><th colspan="6">Dimensions of the principal component analysis of the fatty acids in blood</th></tr><tr><th>Modules</th><th>Dim1 FA</th><th>Dim2 FA</th><th>Dim3 FA</th><th>Dim4 FA</th><th>Dim5 FA</th></tr></thead><tbody><tr><td>Significantly correlated (Pvalue ≤ 0.05)</td><td>2</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>Trend (Pvalue ≤ 0.1)</td><td>4</td><td>0</td><td>3</td><td>3</td><td>0</td></tr></tbody></table>

Modules of co-expressed genes were identified from microarray data in the whole blood by using weighted gene correlation network analysis (WGCNA). Circulating biochemical molecules in the plasma were analysed by 1H-NMR (metabolites; met) or target gas chromatography for the lipid fraction subset (FA) and the data were summarized by weighted linear regressions (dim) using principal component analysis (PCA). Correlations were calculated between the module eigengenes (ME) and the first five dimensions of each PCA. The table indicates the number of significant correlations or trends between modules of coexpressed genes and PCA dimensions.

with other molecules, to circulating lipoproteins (very low density lipoproteins VLDL, low density lipoprotein LDL, and lipids) and to threonine. The third, fourth and fifth dimensions represented respectively 7.9%, 6.0% and 3.8% of the total variability. The third dimension mainly opposed circulating concentrations of glutamine (Gln), glutamate (Glu) and proline (Pro) on one hand, and beta-hydroxybutyrate on the other hand. The fifth dimension opposed circulating concentrations of betaine (bet) and trimethylamine N-oxide (TMAO; i.e., a metabolite produced by the liver and associated with microbiote metabolism), to VLDL and inositol concentrations. Correlation circles for the third, fourth and fifth are available in Supplementary Figures 1 and 2.

Considering specifically the lipid fraction of the plasma, we analyzed the fatty acids (FA) composition by target methodology. From the 30 FA (10 to 22 carbons) that can be analyzed, some of them were present in negligible concentrations in the plasma or even cannot be detected from the background in some pigs (e.g.; C10:0; C12:0; C18:4 n-3). Therefore, parts of the FA were grouped in biologically relevant families (saturated FA with 14 carbons or less, n-6 FA family; n-3 FA family). This led to a total of 14 variables representing single FA or groups of FA. They were then represented by a second PCA to summarize the profiles in circulating FA among the 47 pigs. The Figure 5 shows the corresponding correlation circle. The first dimension represented 42.5% of the total variability and opposed omega-6 (n-6) family of FA and to a lesser extent omega-3 (n-3) FA, to saturated family of FA. The second dimension represented 13.5% of the total variability and opposed C15:0 to C20:0 FA. The third, fourth and fifth dimension represented respectively 8.0%, 7.5% and 7.1% of the total variability. The third dimension opposed C20:1 to C20:2 FA on one hand and C20:0 FA on the other hand, whereas the fourth dimension opposed the sum of n-3 FA to C20:1 FA. Correlation circles for the third, fourth and fifth are available in Supplementary Figures 3 and 4.

## Connecting the two omics levels

We calculated the correlations between the eigengenes of the molecular modules (ME) and the profiles in circulating metabolites and fatty acids represented by the different dimensions of the two PCA. This allows connecting the transcriptome (via the WGCNA modules) and the metabolome (via the principal components).

The numbers of modules for which ME was correlated with at least one of the five dimensions of each PCA are presented in Table 4. A heatmap representing the correlation coefficients calculated between module eigengene (ME) of all modules identified from microarray data and the first five dimensions of the PCA, is available in Supplementary Figures 5 and 6.
There were few associations between molecular modules (33 modules) and metabolic and lipid profiles (as summarized by the first dimensions of the two different PCA). A summary is presented in Figure 6 considering only the list of modules of interconnected genes that have been found to be associated with the phenotypic traits of interest in subsection **Definition of gene co-expression network in the whole blood of pigs**. There were no significant correlations between the eigengenes of the modules associated with FCR and the profiles in circulating metabolites. Only the **darkorange** module tended to be associated with the second dimension (dim2 met) which opposed circulating concentrations of glucose to circulating concentrations of LDL and VLDL lipoproteins. In addition, the **darkgreen** module which was significantly associated with ADFI, was highly correlated with dim3, which summarized the circulating concentrations in some amino acids (Gln, Glu and Pro) and hydroxybutyrate. The **saddlebrown** module which tended to be associated with ADFI and %backfat, was significantly correlated to the second dimension (dim2 met), and also tended to be associated with the fifth dimension (dim5 met) which opposed betaine and TMAO circulating concentrations to inositol concentration. Interestingly, more correlations were generally observed between dimi met and modules related to %loin: the **darkgreen** module was highly correlated with dim3 met, and the **brown** module was correlated with dim3, and to a lesser extent with dim1 met, dim4 met and dim5 met.

Regarding the PCA summarizing the plasma concentrations of FAs, we observed that the module eigen-
gene (ME) of the **royalblue** module, that was negatively correlated with **FCR** and positively correlated with **ADG**, was highly positively correlated to the first dimension (dim1_FA). The **darkorange** module that was also negatively correlated with **FCR**, tended to be correlated with the fourth dimension (dim4_FA). The ME of the **green** also tended to be positively correlated with the first dimension (dim1_FA), whereas ME of the **darkolivegreen** tended to be negatively correlated.

Discussion

**Analyzing inter-individual variability in feed efficiency**

In animal and plant breeding, there has been an increasing interest in intermediate omics traits such as metabolomics and transcriptomics that mediate the effect of genetics on the phenotype of interest [21]. This study confirms that analyzing transcriptome in the whole blood and metabolome in plasma of growing pigs enables to depict the biological molecular pathways involved in various phenotypic traits related to FE, namely feed conversion ratio (FCR) (i.e., the on-farm measure of feed efficiency), and average daily gain (ADG) (i.e., one of the component of FCR describing growth rate during the test period), and to a lesser extent, body composition (%loin, %backfat). This study is a step ahead for the understanding of the relationships between entities that can act in the inter-individual variability in these traits. Indeed, previous studies have addressed differentially-expressed (DE) genes [7, 8], gene set enriched pathways [22] or metabolic signatures [14, 5, 23] between the lowest and the highest FE animals for pigs, cattle or sheep. They have thus compared a gene or a molecule with itself in different conditions such as the response to selection for RFI (a measure of net feed efficiency) or extreme groups based on RFI or FCR. This does not enlighten or explain the interactions between entities (architecture of the trait) and the behavior of regulatory genes in a complex system [24]. Moreover, the aforementioned studies often analyzed only a single type of omics data. Even when they included serum biochemistry in addition to transcriptomics [25] to illustrate consequences of variations in gene expression profiles, they did not intent to depict the correlations between the two levels of life organization (i.e., gene expression profiles and metabolites).

Networks approaches are rather based on the assumption that the effect of the change in the expression level of one entity can be propagated through the interactions on other entities to orchestrate complex phenotypes [26]. The present study identified a total of 33 networks of co-expressed genes in the whole blood across 47 pigs. The network was built with a low threshold set for the minimum number of co-expressed entities in the modules, with the assumption that this may facilitate the identification of co-expressed but also co-regulated genes. The fact that only one to three clusters of enriched biological processes were found in the majority of the molecular modules built in this study argues for the homogeneity in biological processes shared by the co-expressed genes participating to each module. Furthermore, modules in close-vicinity in the dendrogram (hierarchical clustering), such as the **darkolivegreen**, **lightcyan** and **steelblue** modules or
the royalblue, darkred and violet modules, respectively, did not share identical GO terms, which argues for keeping these modules separated rather than merging them together. Altogether, this procedure may then facilitate the identification of co-regulated entities in the different modules.

**Enriched pathways in co-expressed genes modules** From the 33 modules, we found six modules that were significantly associated with FCR and two modules that tended to be related to FCR. In the whole blood of young pigs, [8] previously identified four co-expression modules (minimum of 30 genes per module as threshold, leading to 89 to 786 genes per modules) that were DE between low and high RFI groups, and indicated that these DE genes overlapped with each of the four differentially expressed modules; altogether, they found only one module that was significantly correlated to the RFI phenotype. In addition, [24] identified 34 modules of co-expressed genes from RNA-seq analysis in the liver of low vs high RFI cattle (using a threshold of 30 genes as the minimum in each module), out of which four modules showed significant correlations to RFI. In the present study, the majority of the modules related to FCR were also significantly related to ADG and to %loin, whereas none of them were associated with feed intake (ADFI). This suggests that the main molecular mechanisms in the whole blood explaining the inter-individual variations in FE are involved in the determinism of lean growth potential rather than acting in the regulation of feeding behaviour. The pigs used in this study originated from a divergent selection for RFI and fed different diets during the test period. Because ADG is an independent variable in the regression that estimates predicted feed intake, RFI and ADG have no correlation. The situation is however quite different for FCR, since ADG is part of the ratio in the calculation. Although Gilbert and colleagues [2] indicated low responses (although statistically significant) to RFI selection on lean meat content across generations of pigs, other studies have consistently found higher proportion of lean pieces in the most feed efficient pigs as induced by genetic selection [27, 28] or by management strategies [29]. Skeletal muscle is the largest organ in the body and plays important roles in the utilisation and storage of a large proportion of the energy from feed, which could explain why molecular modules related to FCR were also partly associated with %loin.

The whole blood reflects the concerted actions of the different tissues. In the present study, five of the aforementioned eight modules related to FCR, the violet, lightcyan, darkorange, royalblue and darkolivegreen, were significantly enriched in immune and defense-related processes. In the whole blood, [8] also pointed out immunity and stress response as biological pathways involved in the difference between low and high efficient pigs; however in their study, the module that correlated to RFI phenotype rather included genes participating to cell adhesion, apoptotic process and immunoglobulin production. Other studies on DE genes in pigs have also reported differences in expression levels of genes involved in defense pathways in different tissues between divergent lines for RFI [28]. In cattle, [25] also found an enrichment of the transcriptomic networks in the inflammatory response, regulation of monocyte differentiation, proliferation and differentiation of T lymphocytes, but also apoptosis in the liver from animals with low or high RFI. Although the contribution of the different tissues to the whole blood transcriptome is not known, this suggests that immune-related pathways are probably shared mechanisms orchestrated by tissues that determine inter-individual differences in FE. Defense mechanisms can trigger nutrients to basal metabolism rather than to production performance [30]. Among defense mechanisms, we showed that T cell signaling (royalblue) was negatively related to FCR (i.e., the higher was expression levels of genes involved in T cell signaling, the higher was pig FE), whereas B cell activation (lightcyan) and inflammation (darkolivegreen) were positively associated to FCR (i.e., the higher was inflammation process, the lower was animal FE). Among the modules related to FCR, three modules were rather composed of co-expressed genes that participate to translation (darkred), protein localization (steelblue) or system development (white). In the whole blood of cattle, a set of genes associated with the metabolism of proteins was the most enriched pathway identified as differentially inhibited or activated in high-RFI compared to low-RFI beefs [22]. In our study, these pathways could account for the regulation of lean growth rate, because correlations with ADG and %loin were also identified. Finding nitrogen metabolic process (skyblue), protein modification and muscle cell differentiation (brown) as enriched processes in modules related to %loin further support the assumption that whole blood can encompass molecular mechanisms involved in muscle development and metabolism. When looking at the hierarchy of genes in the molecular modules underlying FCR, we observed that some of
the genes contributing the most to the association between transcriptomics level and animal trait (GS.FCR) have been previously identified as top-ranked genes with fold changes in expression greater than |2| between low and high (RFI) pigs [7]: SLPI (violet module), EIF1B (darkred module) and HTRA1 (ligthcyan. This validates the procedure used in the present study for ranking particular genes in the network architecture associated with feed efficiency.

Finally, two modules, the **darkgreen** and **saddlebrown**, were or tended to be related to both ADFI and **%backfat**, a surrogate of body adiposity. The **darkgreen** module encompassed genes involved in the response to toxic substances. In accordance, there is generally a marked reduction in voluntary feed intake in disease-challenged pigs [31].

**Relationships between transcriptomic and metabolic levels to the definition of feed efficiency or related traits**
The present study finally addressed the interactions between two levels of organization in the circulating entities, i.e., transcriptome and metabolome. Indeed, combined phenotype-metabolome-genome analysis by inferring gene networks based on partial correlation and information theory approaches has been valuable to look on pathways driving FE variation in bovine [32]. In our study, variations at the metabolome level were first summarized by linearly transforming the data into few new coordinates that explained most of the total variance, by using principal component analysis (PCA). This was an exploratory analysis for the hundred variables collected at the metabolic level, that was close-similar to the approach used at the transcriptomic level by calculating the eigengenes (ME) of the WGNA modules in the molecular network of expressed probes. Then, the PCA coordinates were correlated to the eigengenes of the 33 modules. This approach would allow to determine whether pigs with the same profiles in circulating fatty acids (FA) or metabolites may have certain groups of co-expressed and co-regulated genes. Altogether, only few significant correlations were identified between the two organization levels. Nevertheless, an important finding is that the eigengene of the **royalblue** module was highly correlated to the first dimension of the PCA that opposed saturated FA to polyunsaturated FA (omega-6 FA, and to a lesser extent, omega-3 FA families). This means that higher expression levels of genes involved in T signaling, cell adhesion and cell developmental process, were associated with higher proportions in saturated FA and lower proportions in PUFAs in the whole blood, respectively, and this also accounted for a better FE (lower FCR) and higher ADG in the pigs. In addition, the eigengene of the **darkorange** module, which was composed of co-expressed genes also related to immune process and was negatively correlated to FCR, tended to be correlated with the circulating concentration of omega-3 FA family. In pigs, Liu et al. have found that GO terms for lipid metabolism-related processes were enriched in a module of genes that were DE in the whole blood between low and high RFI pigs [8]. Moreover, muscle of high-FE pigs exhibited lower proportion of saturated FA and an enhanced proportion of polyunsaturated FA [27]. Other evidence of altered lipid metabolism between high and low feed efficient groups were found in co-expression analysis in the liver of steers [25]. It has been shown that omega-3 (n-3) PUFA can suppress T cell antigen presentation, activation, proliferation and cytokine expression [33]. In addition, high fat western diet promotes inflammation and modifies immunity [34]. These elements may contribute to the correlation found between the networks of genes circulating concentrations of polyunsaturated FA.

In addition, we suggest a significant association between the regulation of feed intake (ADFI) and the cir-culating concentration of beta-hydroxybutyrate. This was expected since concentration in this ketone body is rising up after fasting and in situation of energy deficit. Circulating concentrations in specific amino acids (AA) such as Pro, Glu and Gln were varying in an opposite synergy with beta-hydroxybutyrate. Because sensing of AA can likely act on the hypothalamic control of food intake [35], this may explain this associa-tion observed in the current study between co-expressed gene network and the circulating concentrations of Pro, Glu and Gln. Almost no significant relationships were however identified between interconnected gene modules related to FCR and the profile in other AA concentrations in the whole blood of the pigs. On the opposite, networks of co-expressed genes related to %loin (such as **brown** and **green**) displayed significant correlations with the metabolomic profile. This suggests direct relationships between AA metabolism and the molecular regulation of muscle growth in the pigs. Indeed, protein (amino acids) metabolism is essential for optimizing efficiency of nutrient absorption and metabolism and to enhance growth performance, notably
because various AA are the building blocks of proteins. Moreover, the darkgreen module which was significantly associated with ADFI, was highly correlated the circulating concentrations in some amino acids (Gln, Glu and Pro) and hydroxy-butyrate.

## Declarations

### Ethics approval and consent to participate

This study was based on previous published studies. The original publications have included a statement on ethics approval to use animals into genetics and feeding experiments. In these publications, the care and use of pigs were performed in compliance with the European Union legislation (directive 2010/63/EU). All animals were reared in compliance with the French national regulations on ethics in animal experimentation (Décret n° 2013-118 du 1er février 2013 relatif à la protection des animaux utilisés à des fins scientifiques), and killed according to procedures approved by the French veterinary services. The protocols were approved by an Ethics Committee in Animal Experiment (Comité Rennais d’Ethique en matière d’Expérimentation Animale, agreement N°07-2012; Ministère de l’Enseignement, de la Recherche et de l’Innovation) from which INRAE is affiliated. All methods have been reported in accordance with ARRIVE guidelines (Animal research: Reporting of in vivo Experiments). The present study follows the internal charts on ethics and animal experimentation at INRAE since we only reused data and did not include additional living animals.

### Consent for publication

Not applicable.

### Availability of data and materials

Transcriptomic dataset was obtained from NCBI's Gene Expression Omnibus (GEO) Subserie accession number GSE70838 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70838). The metabolomic dataset was retrieved in Jegou et al. [14]. Phenotypic traits and fatty acid composition in plasma are are deposited in the publicly available repository at https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/TM2ANC.

The adjacency matrix is available at https://data-access.cesgo.org/index.php/s/YPz0J2ItxIEuN5M.

We provide R notebooks detailing the analysis https://github.com/cjuigne/multiomics_and feed_ efficiency.

### Competing interests

The authors declare that they have no competing interests.

### Funding

Camille Juigné was supported by a doctoral grant from INRAE (Metaprogramme DigitBio) and Région Bretagne (France).

### Authors’ contributions

FG and EB: conceived the study; FG: provided datasets; CJ: implemented the analysis and wrote the initial draft; CJ, FG, EB: discussed the analyses and drafted the manuscript; All authors read and approved the
Acknowledgements

The authors are grateful to Cécile Perrier (PEGASE, INRAE) for gas chromatography analysis to comple-
ment the transcriptomics and metabolomics datasets.

Abbreviations

ADG: average daily gain; ADFI: average daily feed intake; DE: differentially-expressed; FCR: feed conversion ratio; GS: gene significance; ME: module eigengene; PCA: principal component analysis

References

[1] Koch, R.M., Swiger, L.A., Chambers, D., Gregory, K.E.: Efficiency of Feed Use in Beef Cattle. Journal of Animal Science **22**(2), 486–494 (1963). doi:10.2527/jas1963.222486x. Accessed 2023-02-13

[2] Gilbert, H., Billon, Y., Brossard, L., Faure, J., Gatellier, P., Gondret, F., Labussière, E., Lebret, B., Lefaucheur, L., Floch, N.L., Louveau, I., Merlot, E., Meunier-Salaün, M.-C., Montagne, L., Mormede, P., Renaudeau, D., Riquet, J., Rogel-Gaillard, C., Milgen, J.v., Vincent, A., Noblet, J.: Review: divergent selection for residual feed intake in the growing pig. animal **11**(9), 1427–1439 (2017). doi:10.1017/S175173111600286X. Accessed 2022-12-20

[3] Wang, Z., He, Y., Tan, Z.: Transcription Analysis of Liver and Muscle Tissues from Landrace Finishing Pigs with Different Feed Conversion Ratios. Genes **13**(11), 2067 (2022). doi:10.3390/genes13112067. Accessed 2022-12-20

[4] Cantalapiedra-Hijar, G., Abo-Ismail, M., Carstens, G.E., Guan, L.L., Hegarty, R., Kenny, D.A., McGee, M., Plastow, G., Relling, A., Ortigues-Marty, I.: Review: Biological determinants of between-animal variation in feed efficiency of growing beef cattle. animal **12**(s2), 321–335 (2018). doi:10.1017/S1751731118001489. Accessed 2023-02-13

[5] Taiwo, G.A., Idowu, M., Denvir, J., Cervantes, A.P., Ogunade, I.M.: Identification of Key Pathways Associated With Residual Feed Intake of Beef Cattle Based on Whole Blood Transcriptome Data Analyzed Using Gene Set Enrichment Analysis. Frontiers in Veterinary Science **9** (2022). Accessed 2022-12-21

[6] Schmidt, M., Hopp, L., Arakelyan, A., Kirsten, H., Engel, C., Wirkner, K., Krohn, K., Burkhardt, R., Thiery, J., Loeffler, M., Loeffler-Wirth, H., Binder, H.: The Human Blood Transcriptome in a Large Population Cohort and Its Relation to Aging and Health. Frontiers in Big Data **3**, 548873 (2020). doi:10.3389/fdata.2020.548873

[7] Jégou, M., Gondret, F., Vincent, A., Tréfeu, C., Gilbert, H., Louveau, I.: Whole Blood Transcriptomics Is Relevant to Identify Molecular Changes in Response to Genetic Selection for Feed Efficiency and Nutritional Status in the Pig. PLOS ONE **11**(1), 0146550 (2016). doi:10.1371/journal.pone.0146550. Accessed 2022-12-20

[8] Liu, H., Nguyen, Y.T., Nettleton, D., Dekkers, J.C.M., Tuggle, C.K.: Post-weaning blood transcriptomic differences between Yorkshire pigs divergently selected for residual feed intake. BMC Genomics **17**(1), 73 (2016). doi:10.1186/s12864-016-2395-x. Accessed 2023-03-10

[9] Messad, F., Louveau, I., Renaudeau, D., Gilbert, H., Gondret, F.: Analysis of merged whole blood transcriptomic datasets to identify circulating molecular biomarkers of feed efficiency in growing pigs. BMC Genomics **22**(1), 501 (2021). doi:10.1186/s12864-021-07843-4. Accessed 2022-12-21
[10] Carmelo, V.A.O., Banerjee, P., da Silva Diniz, W.J., Kadarmideen, H.N.: Metabolomic networks and pathways associated with feed efficiency and related-traits in Duroc and Landrace pigs. Scientific Reports 10(1), 255 (2020). doi:10.1038/s41598-019-57182-4. Accessed 2022-12-21

[11] Langfelder, P., Horvath, S.: WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9(1), 559 (2008). doi:10.1186/1471-2105-9-559. Accessed 2022-03-23

[12] Salleh, S.M., Mazzoni, G., Løvendahl, P., Kadarmideen, H.N.: Gene co-expression networks from RNA sequencing of dairy cattle identifies genes and pathways affecting feed efficiency. BMC Bioinformatics 19(1), 513 (2018). doi:10.1186/s12859-018-2553-z. Accessed 2022-12-27

[13] Banerjee, P., Carmelo, V.A.O., Kadarmideen, H.N.: Integrative Analysis of Metabolomic and Transcriptomic Profiles Uncovers Biological Pathways of Feed Efficiency in Pigs. Metabolites 10(7), 275 (2020). doi:10.3390/metabo10070275. Accessed 2023-03-24

[14] Jégou, M., Gondret, F., Lalande-Martin, J., Tea, I., Baéza, E., Louveau, I.: NMR-based metabolomics highlights differences in plasma metabolites in pigs exhibiting diet-induced differences in adiposity. European Journal of Nutrition 55(3), 1189-1199 (2016). doi:10.1007/s00394-015-0932-z

[15] Gondret, F., Louveau, I., Mourot, J., Duclos, M.J., Lagarrigue, S., Gilbert, H., van Milgen, J.: Dietary energy sources affect the partition of body lipids and the hierarchy of energy metabolic pathways in growing pigs differing in feed efficiency1,2. Journal of Animal Science 92(11), 4865-4877 (2014). doi:10.2527/jas.2014-7995. Accessed 2023-03-24

[16] Ghazalpour, A., Doss, S., Zhang, B., Wang, S., Plaisier, C., Castellanos, R., Brozell, A., Schadt, E.E., Drake, T.A., Lusis, A.J., Horvath, S.: Integrating Genetic and Network Analysis to Characterize Genes Related to Mouse Weight. PLoS Genetics 2(8), 130 (2006). doi:10.1371/journal.pgen.0020130. Accessed 2023-03-13

[17] Zhang, B., Horvath, S.: A general framework for weighted gene co-expression network analysis. Statistical Applications in Genetics and Molecular Biology 4, 17 (2005). doi:10.2202/1544-6115.1128

[18] Yip, A.M., Horvath, S.: Gene network interconnectedness and the generalized topological overlap measure. BMC Bioinformatics 8(1), 22 (2007). doi:10.1186/1471-2105-8-22. Accessed 2021-11-24

[19] Lê, S., Josse, J., Husson, F.: Factominer: an r package for multivariate analysis. Journal of statistical software 25, 1-18 (2008)

[20] Kassambara, A., Mundt, F.: Package 'factoextra'. Extract and visualize the results of multivariate data analyses 76(2) (2017)

[21] Christensen, O.F., Börner, V., Varona, L., Legarra, A.: Genetic evaluation including intermediate omics features. Genetics 219(2) (2021). doi:10.1093/genetics/iyab130. iyab130. https://academic.oup.com/genetics/article-pdf/219/2/iyab130/41280101/iyab130.pdf

[22] Taiwo, G., Idowu, M., Collins, S., Sidney, T., Wilson, M., Pech-Cervantes, A., Ogunade, I.M.: Chemical Group-Based Metabolome Analysis Identifies Candidate Plasma Biomarkers Associated With Residual Feed Intake in Beef Steers. Frontiers in Animal Science 2 (2022). Accessed 2023-03-10

[23] Goldansaz, S.A., Markus, S., Berjanskii, M., Rout, M., Guo, A.C., Wang, Z., Plastow, G., Wishart, D.S.: Candidate serum metabolite biomarkers of residual feed intake and carcass merit in sheep. Journal of Animal Science 98(10), 298 (2020). doi:10.1093/jas/skaa298. Accessed 2023-03-10

[24] Hudson, N.J., Dalrymple, B.P., Reverter, A.: Beyond differential expression: the quest for causal mutations and effector molecules. BMC Genomics 13(1), 356 (2012). doi:10.1186/1471-2164-13-356. Accessed 2023-03-10
[25] Alexandre, P.A., Kogelman, L.J.A., Santana, M.H.A., Passarelli, D., Pulz, L.H., Fantinato-Neto, P., Silva, P.L., Leme, P.R., Strefezzi, R.F., Coutinho, L.L., Ferraz, J.B.S., Eler, J.P., Kadarmideen, H.N., Fukumasu, H.: Liver transcriptomic networks reveal main biological processes associated with feed efficiency in beef cattle. BMC Genomics **16**(1), 1073 (2015). doi:10.1186/s12864-015-2292-8. Accessed 2023-03-10

[26] Cho, D.-Y., Kim, Y.-A., Przytycka, T.M.: Chapter 5: Network Biology Approach to Complex Diseases. PLOS Computational Biology **8**(12), 1002820 (2012). doi:10.1371/journal.pcbi.1002820. Accessed 2023-03-10

[27] Horodyska, J., Oster, M., Reyer, H., Mullen, A.M., Lawlor, P.G., Wimmers, K., Hamill, R.M.: Analysis of meat quality traits and gene expression profiling of pigs divergent in residual feed intake. Meat Science **137**, 265-274 (2018). doi:10.1016/j.meatsci.2017.11.021

[28] Gondret, F., Vincent, A., Houée-Bigot, M., Siegel, A., Lagarrigue, S., Causeur, D., Gilbert, H., Louveau, I.: A transcriptome multi-tissue analysis identifies biological pathways and genes associated with variations in feed efficiency of growing pigs. BMC Genomics **18**(1), 244 (2017). doi:10.1186/s12864-017-3639-0. Accessed 2023-04-05

[29] Cao, S., Tang, W., Diao, H., Li, S., Yan, H., Liu, J.: Reduced Meal Frequency Decreases Fat Deposition and Improves Feed Efficiency of Growing-Finishing Pigs. Animals: an open access journal from MDPI **12**(19), 2557 (2022). doi:10.3390/ani12192557

[30] Labussière, E., Dubois, S., Gilbert, H., Thibault, J.N., Le Floc'h, N., Noblet, J., van Milgen, J.: Effect of inflammation stimulation on energy and nutrient utilization in piglets selected for low and high residual feed intake. Animal: An International Journal of Animal Bioscience **9**(10), 1653-1661 (2015). doi:10.1017/S1751731115000932

[31] Rodrigues, L.A., Ferreira, F.N.A., Costa, M.O., Wellington, M.O., Columbus, D.A.: Factors affecting performance response of pigs exposed to different challenge models: a multivariate approach. Journal of Animal Science **99**(6), 035 (2021). doi:10.1093/jas/sab035. Accessed 2023-04-05

[32] Widmann, P., Reverter, A., Weikard, R., Suhre, K., Hammon, H.M., Albrecht, E., Kuehn, C.: Systems biology analysis merging phenotype, metabolomic and genomic data identifies Non-SMC Condensin I Complex, Subunit G (NCAPG) and cellular maintenance processes as major contributors to genetic variability in bovine feed efficiency. PloS One **10**(4), 0124574 (2015). doi:10.1371/journal.pone.0124574

[33] Gutiérrez, S., Svahn, S.L., Johansson, M.E.: Effects of Omega-3 Fatty Acids on Immune Cells. International Journal of Molecular Sciences **20**(20), 5028 (2019). doi:10.3390/ijms20205028. Accessed 2023-03-10

[34] Christ, A., Lauterbach, M., Latz, E.: Western Diet and the Immune System: An Inflammatory Connection. Immunity **51**(5), 794-811 (2019). doi:10.1016/j.immuni.2019.09.020

[35] Morell, P., Fiszman, S.: Revisiting the role of protein-induced satiation and satiety. Food Hydrocolloids **68**, 199-210 (2017). doi:10.1016/j.foodhyd.2016.08.003. Accessed 2023-04-05
<figure><img src="image_3.png" /><figcaption>Figure 1: Heatmap of correlations between module eigengenes and animal phenotypic traits.</figcaption></figure>

Module eigengene (ME) was the representative of gene expression profile in the module of co-expressed molecular probes elicited from the weighted gene correlated network analysis (WGCNA) from microarray data in the whole blood of 47 growing pigs. Animal phenotypic traits were recorded during a test period of 58 days. The heatmap indicates the Pearson correlation coefficient between ME and the phenotypic trait together with the statistical significance (*P*value).
<figure><img src="image_4.png" /><figcaption>Figure 2: Hierarchical clustering of module eigengenes.</figcaption></figure>

The modules that were found highly correlated with feed conversion ratio (FCR) are enlightened (*** $P \le 0.001$, * $P \le 0.05$, and † $0.05 < P \le 0.10$).
<figure><img src="image_5.png" /><figcaption>Figure 3: Graphical representation of the royalblue module.</figcaption></figure>

The network was constructed from the adjacency matrix of the *royalblue* module using Cytoscape. The nodes are the molecular probes, labeled with their annotation and connected by purple edges that represent the correlation coefficients greater than the 95th percentile and by green edges, for probes that were not sufficiently correlated to other probes, to the annotated probe to which they were the most correlated. The size of the nodes is a function of their degree.
<figure><img src="image_6.png" /><figcaption>Figure 4: Correlation circle of the principal component analysis summarizing the profiles of circulating metabolites.</figcaption></figure>

1H-NMR spectra were obtained in the plasma prepared from the whole blood of 47 growing pigs. The matrix of correlations was calculated from 94 individual variables corresponding to the different annotated spectra.
<figure><img src="image_7.png" /><figcaption>Figure 5: Correlation circle of the principal component analysis summarizing fatty acid composition in blood.</figcaption></figure>

The fatty acid composition (in percentage) was obtained in the plasma prepared from the whole blood of 47 growing pigs by using gas chromatography. Some of the individual FA were grouped in biologically relevant families (saturated FA with 14 carbons or less, omega-6 sum of *n* − 6 and omega-3 sum of *n* − 3), whereas the other fatty acids were kept as these.
<figure><img src="image_8.png" /><figcaption>Figure 6: Heatmap of correlations between molecular modules and profiles of circulating molecules.</figcaption></figure>

Modules of co-expressed probes were obtained from a weighted gene correlation network analysis (WGCNA) from microarray data in the whole blood of 47 growing pigs. The eigengene of each module (ME) was considered as a mathematical representative of the expression levels of the molecular probes within the module. Circulating biochemical molecules were analyzed and the data were summarized by weighted linear correlation using principal component analysis (PCA). The first two dimensions of the PCA were called dim1 met and dim2 met for the metabolites obtained by 1H-NMR high throughput method and dim1 FA and dim2 FA for the fatty acids analyzed by target gas chromatography, respectively.