A complete description of the experimental animals, phenotypes, sample extraction, and RNA-sequencing of muscle and liver tissues are described in Almeida et al. (2021), Fanalli et al. (2022), and Freitas et al. (2024). Briefly, a total of 72 immunocastrated Large White male pigs (28.44 ± 2.95 kg) were used in a 98-day experimental period. All animals had ad libitum access to feed and water throughout the experimental period (98 days). Four days prior to their slaughter, blood samples were taken from all pigs for determination of glucose (GLU; mg/dL), aspartate aminotransferase (AST; U/L), total proteins (TP; g/dL), albumin (ALB; g/dL), globulin (GLOB; g/dL), triglycerides (TG; mg/dL), cholesterol (CHOL; mg/dL), high-density lipoprotein (HDL; mg/dL), low-density lipoprotein (LDL; mg/dL), and very low-density lipoprotein (VLDL; mg/dL).

After 98 days, all pigs were slaughtered (average final body weight of 133.9 ± 9.4 kg), and skeletal muscle (Longissimus lumborum) between the 10th and 11th ribs and right lobe of the liver samples were collected within a maximum of 30 min after bleeding. The tissue samples were quickly collected, snap-frozen in liquid nitrogen, and then stored at −80^oC until further analyses. At slaughter, meat and carcass quality traits were also measured, including slaughter weight (SW; in kg), cold carcass yield as a percentage of the slaughter weight (CCY; %), loin eye area measured by ultrasound (LEA; cm²), backfit thickness measured by ultrasound (BFT; cm²), intramuscular fat content (IMF, %), and liver fat content (LFC; in %). Furthermore, tissue samples were used to measure the abundance of cytokine levels in skeletal muscle, liver, and blood for Interleukin-10 (IL-10; MFI), interferon-gamma (IFN-γ; MFI), interleukin-1 beta (IL-1β; MFI), interleukin-6 (IL-6; MFI), interleukin-18 (IL-18; MFI), and tumor necrosis factor-alpha (TNF-α; MFI).

Skeletal muscle (Longissimus lumborum) and right lobe liver tissue samples were collected after slaughter for RNA-seq. Total RNA was extracted from the frozen tissue samples, and the RNA integrity was verified based on RNA integrity number (RIN). Only samples with RIN higher than seven were used. Sequencing adaptors and low-complexity reads were removed using the Trim Galore 0.6.5 software (Krueger et al., 2019), and reads longer than or equal to 70 bases and a Phred score threshold greater than 33 were kept for further analyses. After this quality filtering step, alignment and mapping were performed using the current reference pig genome (Sscrofa 11.1) (Warr et al., 2020), generating Genomic Variant Call Format (GVCF) (McKenna et al., 2010; Van der Auwera et al., 2013; Franke and Crowgey, 2020) files for each sample from the liver and skeletal muscle tissues.

The Genome Analysis Toolkit (GATK, v. 4.1.9.0) was used in the GVCF format (Van der Auwera et al., 2013), adopting the HaplotypeCaller algorithm (Van der Auwera et al., 2013) for individually calling the variants for each sample. The output data files with the annotated variants were merged by tissues using the CombineGVCF tool (Van der Auwera et al., 2013; Poplin et al., 2018), and the joint genotyping analysis was performed using the GenotypeGVCF. Subsequently, a VCF file with all samples genotyped for each tissue was obtained and the variants’ annotation and functional consequences were predicted using the Ensembl Variant Effect Predictor tool v. 101 (VEP) (McLaren et al., 2016). The SNP data from muscle, liver, and GGP50K (GeneSeek Genomic Profiler -GGP Porcine 50K, a medium-density SNP chip array with 50,915 SNPs) were merged into a single dataset for subsequent analyses. When the alleles of SNPs between the datasets were different, is the SNP was considered as missing, and both were removed. The variants were filtered based on variant quality scores equal to or greater than 30 (QUAL) and total coverage depth (DP) greater than 10, using BCFtools v. 1.9 (Danecek et al., 2011; 2021; Li, 2011). Moreover, SNP with call rate lower than 95%, minor allele frequency (MAF) lower than 5%, SNPs located in non-autosomal chromosomes, SNPs with extreme departure from the Hardy-Weinberg equilibrium (P < 10⁻⁶), and non-biallelic markers were removed from the genomic dataset. Finally, linkage disequilibrium (LD) pruning was applied based on a r ² threshold of 0.8 within a 100 kb window (Freitas et al., 2024). LD pruning was incorporated into the quality control process to minimize false positives and remove redundant markers. Furthermore, as demonstrated by Freitas et al. (2024), pruning SNPs in linkage disequilibrium can substantially improve the detection of relevant eQTLs in complex traits by reducing confounding effects.

The cis- and trans-eQTL were evaluated using an additive linear model implemented in the Matrix eQTL package (Shabalin, 2012). Principal components (PCs) were fitted as covariates in the models to correct for potential population stratification, along with sire information, dummy categories for treatment effect, and initial body weight as a linear covariate. The expression levels in muscle and liver tissues were tested for association using the combined SNP dataset described above. Gene expression levels were normalized using the average method while preserving rank, a method endorsed by the GTEx consortium (Aguet et al., 2020). Cis-eQTL, defined as local effects, were considered if they were within 1 Mb upstream or downstream of the genes (first and final base pair positions in the gene map). Trans-eQTL were defined as those with a distance greater than 1 Mb from the genes. The model fitted can be defined as: g = α + γ   x + β s + ε , where g is the gene expression level in transcripts per million (TPM); α is the mean term; γ is the slope coefficient of x , which represents the kth covariates, including the top 10 principal components (which explained approximately 28% of the structural population variance), initial body weight in kg as a covariate, sire dummy variables, and categorical effect of treatment dummy variables; β is the slope coefficient of s , which represents the genotype coded as 0 (homozygous for the reference allele), 1 (heterozygous), and 2 (homozygous for the alternative allele), and ε is the residual of the model. We adopted the False Discovery Rate (FDR) correction method for p-values (Benjamini and Hochberg, 1995; Huang et al., 2018) with a significance threshold of 0.05.

After the eQTL analyses, significant eQTL were used for association with the biochemical and cytokine profiles in pigs. The traits analyzed included GLU, AST, TP, ALB, GLOB, TG, CHOL, HDL, LDL, VLDL, SW, CCY, LEA, BFT, IMF, LFC, IL-10, IFNg, IL-1β, IL-6, IL-18, and TNF-α in skeletal muscle or liver tissues and blood serum of pigs. The association analysis was conducted using the GCTA software (v.1.94.1) (Yang et al., 2011), employing the mixed linear model (MLMA-LOCO) with the genomic relationship matrix of the animals (GRM). Phenotypic traits were previously adjusted to treatment and block effects. The model used was: y = a + bx + g − + ε ,

where y represents the trait, a is the observed mean, b is the additive effect (fixed effect) of the SNP or eQTL candidate being tested for association, x is the indicator variable of the SNP genotype coded as 0 (homozygous for the reference allele), 1 (heterozygous), and 2 (homozygous for the alternative allele), g − is the polygenic effect (random effect), representing the cumulative effect of all SNPs, except those located on the chromosome of the candidate SNP (the variance is re-estimated each time when a chromosome is excluded), and ε is the residual effect. The p-values resulting from the association analysis were corrected for multiple testing using the FDR method (Benjamini and Hochberg, 1995; Huang et al., 2018) and a significance threshold of 0.05. To illustrate the results, we used the CMplot (Yin et al., 2021) package of the R, to generate the density plot, Q-Q plot, and Manhattan plot in the R environment.

The functional genomic analyses were performed following Freitas et al. (2024). Briefly, the GALLO R package (Fonseca et al., 2020) was used to perform the QTL annotation and enrichment of the eQTL associated with biochemical parameters and cytokine profiles. The eQTL annotation and enrichment were performed using known QTL data obtained from the PigQTLdb database (Release 53 – Sscrofa11.1, 28 Apr 2024), considering a genomic window of up to 500 kb downstream and upstream of the genomic coordinates of the eQTL. The QTL enrichment analyses were performed using a hypergeometric test to reduce the bias of overrepresented traits.

The gene enrichment analysis was performed using the Over-Representative Analysis (ORA) method on the WEB-based Gene Set Analysis Toolkit (Elizarraras et al., 2024; Liao et al., 2019; Elizarraras et al., 2024). The Gene Ontology terms include Biological Processes, Cellular Components (non-redundant), Molecular Functions (non-redundant), and Biological Pathways. The gene list was based on annotations, considering a window of 500 kb up and downstream of significantly associated eQTL genomic coordinates. The gene data annotation of the Sus scrofa (Assembly Sscrofa11.1; genome-build-accession GCA_000003025.6; available at: https://ftp.ensembl.org/pub/release-112/gtf/sus_scrofa/) were extracted from the Ensembl platform (Ensembl release 112 - August 2024) (Aken et al., 2016) in the General Transfer Format (“. gtf” format). Finally, multiple protein-protein interaction (PPI) analyses were performed using the STRING 12.0 package (version: 26 July 2023, https://string-db.org/). We explored protein-protein interactions using the gene list with a focus on the Sus scrofa species. For that, the same genes annotated around eQTL were used as input. Furthermore, functional genomic annotations were obtained by consulting databases such as UniProt (www.uniprot.org), Ensembl (www.ensembl.org/Sus_scrofa/), and National Center for Biotechnology Information (NCBI, www.ncbi.nlm.nih.gov/).

Initially, 84,809 SNPs were identified in liver samples with a genotyping rate of 0.947. Upon merging with the GGP-50k dataset, the number of SNPs increased to 122,325, with a combined genotyping rate of 0.952. Further merging with 75,447 muscle tissue SNPs resulted in 146,344 unique SNPs, and the final genotype rate was 0.926. Quality filtering was then applied to the combined SNP dataset. After removing SNPs with a missing genotype rate greater than 5% (67,861 variants removed), an extreme departure from the Hardy-Weinberg equilibrium test (p-value <10⁻⁶, 539 variants removed), and MAF below 5% (7,251 variants removed), 70,693 SNPs remained for further analyses. Additional 28,934 SNPs were removed during the LD pruning step and 41,759 SNPs from seventy-two animals remained for further analyses, with a final genotyping rate of 0.989. The detailed process and results of SNP data combination, quality filtering, and LD pruning are provided in Supplementary Tables S1–S3.

The SNPs (n = 41,759) were tested for association with the expression level of 15,093 and 15,516 genes from muscle and liver, respectively. The number of significant eQTL and regulated genes found in muscle and liver tissue are shown in Figure 1, and a summary of eQTL analysis is provided in Supplementary Material 1.

Figure 2 shows the number of eQTLs for each 1 Mb window, dispersed across the chromosomes used to test associations with the traits (1,199 eQTLs). The descriptive statistics of the performance traits, blood biochemical parameters, and cytokine profiles of pigs are presented in Supplementary Table S4.

The Q-Q plot (Figure 3) shows the expected distribution of -log10(p) values versus the observed distribution for the phenotypes in blood serum AST, HDL, IL-6, and liver tissue IL-6 and IL-18. A summary of the genomic inflation factors of the p-values resulting from the GWAS analysis with eQTL are presented in Supplementary Table S5.

Figures 4, 5 illustrate the distribution of eQTLs across the genome for the phenotypes in liver tissue IL-6 and IL-18 (Figure 4), as well as blood serum, AST, HDL, and IL-6 (Figure 5). The highlighted points in the plots represent the significant eQTL. Specific SNPs such as rs345667860 (3′UTR), rs695637860 (Downstream), and rs337362164 (Missense) are indicated, with effects predicted using the Variant Effect Predictor (VEP) from Ensembl.

The Manhattan and QQ plots of all other traits are in Supplementary Material 2 (.zip). Furthermore, the summary statistics table of GWAS, which includes all genomic inflation values is presented in Supplementary Table S5.

The QTL types enriched with the significant eQTLs (FDR <0.05) located at SSC12:39,493,883, 8:39,107,307, and 5:88,678,346 (chromosome:base-pair) was “Health” followed by “Production” and Meat and Carcass.” Top significant traits in Production, Health, Meat, and Carcass enrichment analyses around eQTL associated with AST, HDL, IL-6 in pig serum, and IL-18 and IL-6 in pig liver tissue are shown in Figure 6.

The gene enrichment type analysis for the genes around significant eQTL included 34 mapped genes from unique Entrez gene IDs. Parameters included a minimum of three gene IDs per category, with an enrichment significance level adjusted to FDR <0.05. The bar charts in Figure 7 show the distribution of these Gene Ontology categories (Biological Processes–BP, Molecular Functions–MF, and Cellular Components–CC).

The significant GO (Gene Ontology) and MP (Molecular Pathway) terms are shown in Figure 8, where the x-axis corresponds to the log₂ enrichment ratio relative to the -log₁₀ FDR on the y-axis for genes annotated around the eQTL associated with biochemical parameters from pig blood and liver tissue.

For better understanding the enrichment analysis, a directed acyclic graph (DAG, output report for ORA from WebGestalt, the complete results can be accessed on https://2024.webgestalt.org/results/1732400034/) was generated to describe the hierarchical relationships between the enriched biological processes, and the full version can be viewed in Supplementary Figure 1 (.png). We highlighted the nodes corresponding to “Cytokine Response–GO:0034097” in Figure 9.

A PPI network analysis for genes annotated around eQTLs associated with blood serum biochemical parameters and cytokines in Sus scrofa, which was performed to investigate their relationship with the annotated genes around eQTLs associated with biochemical parameters from pig blood serum and liver tissue are shown in Figure 10. The PPI network displayed an average node degree of 1.79 and a local clustering coefficient of 0.455, suggesting moderate connectivity among the nodes (medium confidence). The PPI enrichment p-value <10⁻¹⁶ indicates that the observed interactions are significantly more frequent than expected by chance, suggesting potential functional relationships among these genes. Clusters within the network were identified using Markov Cluster Algorithm (MCL) clustering with an inflation parameter of 3, highlighting groups of genes involved in processes such as chemotaxis and cytokine signaling pathways. The nodes without links were hidden, but the complete analysis can be accessed on https://version-12-0.string-db.org/cgi/network?networkId=bE85LK5REtmx and Supplementary Material 3.