The animal material used in this study was previously reported in Ballester et al. (7). In brief, we used a total of 432 weaned piglets (217 males and 215 females) from a commercial Duroc pig line. The pigs were distributed in six batches and belonged to 134 litters obtained from 132 sows and 22 boars. All animals were raised in the same farm and fed ad libitum with a commercial cereal-based diet. At the time of sample collection (60 ± 8 days of age), all animals were apparently healthy, without any sign of infection.

Details of sampling and laboratory determination of immunological, haematological and stress parameters are described in Ballester et al. (7). In brief, total concentrations of immunoglobulins IgA, IgG and IgM in plasma, and IgA in saliva (IgAsal), were measured by ELISA with commercial kits (Bethyl laboratories Inc., Bionova, Spain), following manufacturer’s instructions. Concentration of acute phase proteins: C-reactive protein (CRP) and haptoglobin (HP) were measured in serum by ELISA kit (Abcam Plc., Spain) and colorimetric assay (Tridelta Development Limited, Ireland), respectively, following manufacturer’s instructions. Total concentration of nitric oxide (NO) in serum was measured by colorimetric assay (Termo Fisher Scientifc, Spain). To determine the percentage of gamma-delta T cells (γδ T cells), peripheral blood mononuclear cells (PBMCs) were separated from heparinized peripheral blood by density-gradient centrifugation with Histopaque-1077 (Sigma, Spain) and 10⁶ PBMCs were stained using the monoclonal antibody APC Rat Anti-Pig γδ T Lymphocytes (MAC320 clone, BD Pharmigen, Spain) and the APC Rat IgG2a κ isotype control (R35-95 clone, BD Pharmigen, Spain). The cells were analysed by flow cytometry using the MACSQuant Analyzer 10 Flow cytometer (Miltenyi Biotec GmbH, Bergisch Gladbach, Germany) and the MACSQuantify sofware v2.6 (Miltenyi Biotec GmbH, Bergisch Gladbach, Germany). Phagocytosis assay was carried out in heparinized whole blood samples incubated with fluorescein (FITC)-labelled opsonized E. coli bacteria by using the Phagotest kit (BD Pharmigen, Spain) according to manufacturer’s instructions. The following percentages of phagocytosis traits in blood: percentage of total phagocytic cells (PHAGO_%); percentage of phagocytic cells among granulocytes (GRANU_PHAGO_%), monocytes (MON_PHAGO_%) and lymphocytes (LYM_PHAGO_%); mean fluorescence in fluorescein isothiocyanate (FITC) among the total phagocytic cells (PHAGO_FITC); and mean fluorescence in FITC among the granulocytes (GRANU_ PHAGO_FITC), monocytes (MON_PHAGO_FITC) and lymphocytes (LYM_PHAGO_FITC) that phagocyte were analysed by flow cytometry (MACSQuant Analyzer 10 Flow cytometer; Miltenyi Biotec GmbH, Bergisch Gladbach, Germany) and the MACSQuantify sofware v2.6 (Miltenyi Biotec GmbH, Bergisch Gladbach, Germany). Haematological traits: haematocrit (HCT), haemoglobin (HB), mean corpuscular volume (MCV), mean corpuscular haemoglobin (MCH), mean corpuscular haemoglobin concentration (MCHC), total number of leukocytes (LEU), eosinophils (EO), lymphocytes (LYM), monocytes (MON), neutrophils (NEU), erythrocytes (ERY) and platelets (PLA) were determined from blood samples by haemogram. Finally, the stress indicators analysed were the neutrophil to lymphocyte ratio (NLR) and the cortisol (CORT) levels in hair, which were measured by ELISA kit (Cusabio Technology LLC., Bionova, Spain). Table S1 summarizes the descriptive statistics of the analysed traits.

The 432 animals were genotyped with the GGP Porcine HD Array (Illumina, San Diego, CA) using the Infinium HD Assay Ultra protocol (Illumina). The quality control of the 68,516 SNPs was performed using Plink software (24), removing SNPs with a minor allele frequency lower than 5%, SNPs with more than 10% missing genotype data, and SNPs that did not map to the porcine reference genome (Sscrofa11.1 assembly). Thereafter, a subset of 42,641 SNPs was retained for subsequent analysis.

GWAS for the 30 health-related traits was performed using the following mixed linear model with the Genome-wide Complex Trait Analysis (GCTA) software tool (25):

where y_ijk corresponds to the phenotypic trait (either log-transformed or raw data) of the i^th individual of sex j in the k^th batch; sex_j corresponds to the j^th sex effect (2 levels); b_k corresponds to the k^th batch effect (6 levels) for most traits but for phagocytosis related traits, for which the data of laboratory analysis (12 levels, two by batch) was considered instead; u_i is the infinitesimal genetic effect of individual i, with u∼N(0,G σ u 2 ), where G is the genomic relationship matrix calculated using the filtered autosomal SNPs based on the methodology of (25) and σ u 2 is the additive genetic variance; s_li is the genotype (coded as 0,1,2) for the l^th SNP, and a_l is the allele substitution effect of the SNP on the trait under study; and e_ijk is the residual term.

The GCTA software tool (25) was also used to estimate the heritability value for each one of the 30 phenotypes by performing a restricted maximum likelihood (REML) analysis with sex and batch as fixed effects.

The results from GWAS analysis were subjected to an AWM approach (11, 12). For the initial step, only those SNPs located in the coding region or within 5 kb of an annotated gene based on the Sscrofa11.1 reference genome assembly were kept. Then, γδ T cells was chosen as key phenotype and those γδ T cells-associated SNPs (P<0.05) were selected. In addition, we also retained those SNPs that were associated (P<0.05) with three or more of the 30 health-related phenotypes and referred to as pleiotropic SNPs. For assessing the gene-gene interactions, the standardized SNP effects across phenotypes were used to infer the gene co-association network using the Partial Correlation and Information Theory (PCIT) approaches (26). In the network, every node represents a gene (or a SNP), whereas every edge connecting two nodes represents a significant gene–gene interaction (based on SNP–SNP co-association). In order to quantify the number of connections as well as to prioritize potential regulators of the gene network we applied an information lossless approach (27) that explored the connectivity of all regulators (transcription factors, miRNAs and lnRNAs) and its target genes. The list of the regulators that was cross-matched with the list of 3,636 genes corresponds to the census of transcription factors reported by Vaquerizas et al. (28). The Cytoscape software (29) was used to visualize the gene network and calculate node centrality values.

When only the Ensembl ID was available for some of the associated novel genes, they were searched against the database for human orthologues (GRCh38.p13) to extract their putative gene name using the BioMart web-based tool (30).

ClueGO plugin (31) was used to identify the over-represented gene ontology terms, KEGG pathways, and immune functions. The cut-off for considering a significant over-representation (corrected p-value ≤ 0.05) was established by Benjamini and Hochberg multiple-test correction (32). Further categorization of candidate genes was performed using information from the Mouse Genome Database (www.informatics.jax.org) (33) and GeneCards (www.genecards.org) (34).

A comparative analysis between our results and previous published data was performed by retrieving all pig and human quantitative trait loci (QTLs) and association data on Sscrofa11.1 and on GRCh38.p13 (dbSNP Build 153), for health traits from the pigQTL database (35) and from the NHGRI-EBI GWAS Catalog (v1.0.2) (36), respectively.

From the 42,641 SNPs used in the GWAS analysis, 2,230 SNPs were directly associated (P<0.05) with the key phenotype (γδ T cells), while 8,421 SNPs were associated (P<0.05) with three or more phenotypes. Thus, a total of 9,830 SNPs were obtained after merging both lists of associated SNPs. Two SNPs located in Sus scrofa chromosome (SSC) 13 (rs323669145 and rs81448146) were associated with the maximum number of phenotypes, 13 out of the 30 traits (CRP, HP, NO, CORT, MCHC, LYM, NLR, PHAGO_%, PHAGO_FITC, GRANU_PHAGO_%, GRANU_ PHAGO_FITC, MON_PHAGO_FITC, LYM_PHAGO_%). Those two SNPs were located in intronic regions of the fibroblast growth factor 12 (FGF12) gene. Among other SNPs associated with a high number of phenotypes (n≥10) we also identified rs81284215 on SSC9. Remarkably, this SNP was located in the haemoglobin subunit beta (HBB) gene and was associated with haematological-related phenotypes such as HCT, HB and ERY, as well as with γδ T cells and six out of the eight phagocytosis phenotypes, excluding those related to lymphocytes.

After the annotation step, only 3,544 out of the 9,830 SNPs were retained as they were located in the coding region or 5 kb upstream or downstream of a gene. However, as some SNPs were annotated to more than one gene, a total of 3,650 unique genes were identified. To explore the genetic determinism of the 3,544 SNPs selected from the AWM approach over the 30 phenotypes, their heritability was calculated using these SNPs and compared to the heritabilities obtained using the whole dataset comprised by 42,641 SNPs ( Table 1 ). Heritability calculated with the SNPs subset ranged between 0.238 and 0.791. The highest heritabilities (h²>0.7) were observed for the haematological traits (MCH, MCV, HB, HCT), while CRP, MON, and IgAsal showed low heritabilities (h²<0.3). However, the other immunoglobulin-related traits (IgG, IgA and IgM) had heritabilities above 0.5. In addition, most of the traits (24 out of 30) showed moderate to high heritabilities (h²>0.4). When compared with the heritabilities obtained with the whole dataset, the heritabilities of the SNPs subset showed higher values in all comparisons except one (LYM_PHAGO_FITC), whose values were almost similar (0.374 vs. 0.379). The greatest improvement in the heritability value was observed for HB and HCT, while others, such as the values for CRP, IgG, and the aforementioned LYM_PHAGO_FITC, did not change much. Of note, the heritability of the key phenotype (γδ T cells) was also greatly improved from moderate in the whole dataset (0.439) to high in the subset (0.688). Furthermore, to validate the heritabilities obtained by the reduced dataset of 3,544 SNPs, the heritability values of the 30 phenotypes were calculated by randomly selecting the same number of SNPs from the whole dataset in 100 different iterations. Thus, for each phenotype, the heritability values of the reduced dataset were always higher than the heritability values randomly obtained with each iteration with the exception, again, of the LYM_PHAGO_FITC phenotype.

Genetic relationships among phenotypes based on the normalized additive values of the 3,544 SNPs are shown in Figure 1 . The phagocytosis capacity phenotypes (measured through FITC) and the phenotypes for the proportion of phagocytic cells were highly correlated ( Figure 1A and Table S2 ), and clustered together, jointly with NEU and NLR ( Figure 1B ). Haematological erythrocyte-related measurements had also a high positive correlation, similar to the ones reported for the white blood cell types ( Table S2 ). These phenotypes were separated in a second cluster, with HB, HCT and ERY being sub-clustered together and in close proximity to MCH and MCV, and all leukocyte-related counts but NEU (i.e., LEU, EO, MON and LYM) were grouped together in a differentiated branch ( Figure 1B ). Associated to them, the last sub-cluster gathered most of the innate immunity traits such as positive correlated acute-phase proteins (CRP and HP) or positive correlated immunoglobulins concentrations (IgA, IgG and IgM) in plasma, but also other phenotypes as CORT, PLA or MCHC. The latter showed a high correlation with IgA measured in saliva ( Figure 1A and Table S2 ). The percentage of γδ T cells, the key phenotype used in the AWM approach, was also enclosed in this last sub-cluster, but far from the rest of phenotypes, since it was not correlated with any phenotype ( Figure 1A ).

The gene co-association network was comprised of 3,636 nodes and 451,407 interactions ( Figure S1 ). The topology of the network reflected an average degree of 248.299, a clustering coefficient of 0.347 and an average path length of 2.219, which manifest a great level of connection between the nodes of the network. Of note, the three genes with the highest degree (i.e., the highest number of connections) were REPS2, ACSL3 and WLS, with 893, 858 and 837, respectively. ACSL3 and WLS knockout mice showed abnormal hematopoietic and immune responses, including decreased γδ T cell number in the latter (37).

After PCIT/network analysis, we identified a total of 246 potential regulators, including three microRNAs, within the 3,636 co-associated genes for the 30 health-related traits. The information lossless approach (12) allowed identifying trios of regulators with maximum connectivity (i.e., interactions) within the network; the top 10 trios are shown in Table 2 . Among them, those composed by BRMS1L, MED12L and SUPT3H, and by ARNT2, SUPT3H and TRIM25, were among the top trios that spanned most of the network topology with highest connectivity (1,482 and 1,472 unique interactions, respectively) and minimum redundancy ( Table 2 ). BRMS1L and TRIM25 were the regulators with the highest number of co-associated genes (n=603), followed by SUPT3H (n=580), ARNT2 (n=552), and MED12L (n=466). All these protein-coding genes are transcription factors with repressor (BRMS1L), co-activator (MED12L), or activator functions. For each significant association between a SNP of the five key transcription factors and a trait, the allele substitution effect of the SNP on the respective phenotype is summarised in Table S3 .

A functional classification considering all the genes (n=1,828) co-associated with the five top regulators was performed to identify biological processes, metabolic pathways and candidate genes directly related to immunocompetence and study their role in immune system processes. The list of immune biological processes and pathways identified is shown in Table S4 . Overall, genes for immune cells activation, differentiation and proliferation, regulation of innate and adaptive immune responses, phagocytosis, immunoglobulin and cytokine production and regulation of inflammatory response processes were identified. It is worth to highlight the JAK/STAT pathway, the ERK1/2 and MAPK cascades, and the NOTCH and Phospholipase D signaling pathways among those processes. Furthermore, we also identified other metabolic pathways directly related to immunometabolism such as the AMPK, mTOR, PI3K-Akt and PPAR signaling pathways, which control oxidative phosphorylation, carbohydrate, lipid and peptide metabolic processes, among others. This functional analysis revealed that a total of 589 candidate genes belonged to immune-related functions ( Table S4 ).

In addition, we studied in detail for each top transcription factor those immune processes in which a high number of co-associated genes were involved ( Table S5 ). Figure 2 shows a simplified co-association network representing those immunity-associated genes interacting with the top 5 regulators. Table 3 summarizes the main immune and metabolic processes in which the top transcription factors were associated. It is worthy to note that all the regulators were co-associated with genes related to platelet functions and phospholipase D signaling. Furthermore, among the top five regulators SUPT3H was the transcription factor co-associated with the highest number of pathways and processes related to immunometabolism and haemopoiesis and immune development, and BRMS1L was the only transcription factor co-associated with genes related to the NOTCH signaling pathway. Finally, we identified a specific gene, protein tyrosine phosphatase receptor type C (PTPRC), involved in γδ T cell differentiation and activation that was co-associated with ARNT2, BRMS1L and MED12L.

When we studied in detail the top five regulators and their co-associated genes related to immune functions, we identified interesting candidate genes enclosed in genomic regions previously associated with immune phenotypes in studies with humans and pigs ( Table S6 ). For instance, PRKCE, a candidate gene for haematological parameters in humans (38–41) and different leukocyte subpopulations in pigs (42), was associated in our study with phagocytosis-related traits. Other genes also associated with phagocytosis-related traits in our population were CD36 and SH2B3. In humans, CD36 has been associated with haematological traits (38, 39), while SH2B3 has been associated with a plethora of immune phenotypes, mainly with leukocyte levels and their types (38–40). MYH9 was another candidate gene associated with haematological parameters in humans (38–40) together with ANGPT1 (38, 39, 43), where the latter was also associated with the levels of exhaled nitric oxide and IgG glycosylation in serum (44, 45). In our analysis, MYH9 was associated with LYM and phagocytosis-related traits while ANGPT1 was also associated with NO and some haematological traits. The previously mentioned PTPRC (also known as CD45) was associated in humans with leukocyte, lymphocyte and monocyte counts (38, 39, 43) and memory T cell differentiation (46). In our population this gene was associated with LYM, among other traits. Other candidate genes commonly associated with immune cell types in humans were: ATXN1 and CLEC16A (38–40, 43), CD4 (46), FTO and KLRK1 (47); while in pigs were AKT2 (48), CLEC16A (42), JMJD1C and RBFOX1 (8), PRKG2 and PAQR3 (49), and SLC30A10 (50). CD4 and KLRK1 are mainly expressed in T cells and Natural Killer cells and they were associated with γδ T cells in our population. Another interesting gene described in association with CRP in humans (51) and with platelet levels in pigs (50) was PDE4B. This gene was associated with haematological and phagocytosis-related traits in our study.

Furthermore, while none of the top five regulators were located in QTL regions associated with health-related traits in our Duroc pig population (7), they were located in QTLs for C3c concentration, HCT, HB, Interleukin (IL)-2 and IL-10 levels, Toll-like receptor (TLR)-2 and TLR-9 levels, and basophil number, among others ( Table S6 ), in other pig populations (35).