A Co-Association Network Analysis Reveals Putative Regulators for Health-Related Traits in Pigs

In recent years, the increase in awareness of antimicrobial resistance together with the societal demand of healthier meat products have driven attention to health-related traits in livestock production. Previous studies have reported medium to high heritabilities for these traits and described genomic regions associated with them. Despite its genetic component, health- and immunity-related traits are complex and its study by association analysis with genomic markers may be missing some information. To analyse multiple phenotypes and gene-by-gene interactions, systems biology approaches, such as the association weight matrix (AWM), allows combining genome wide association study results with network inference algorithms. The present study aimed to identify gene networks, key regulators and candidate genes associated to immunocompetence in pigs by integrating multiple health-related traits, enriched for innate immune phenotypes, using the AWM approach. The co-association network analysis unveiled a network comprised of 3,636 nodes (genes) and 451,407 edges (interactions), including a total of 246 regulators. From these, five genes (ARNT2, BRMS1L, MED12L, SUPT3H and TRIM25) were selected as key regulators as they were associated with the maximum number of genes with the minimum overlapping (1,827 genes in total). The five regulators were involved in pathways related to immunity such as lymphocyte differentiation and activation, platelet activation and degranulation, megakaryocyte differentiation, FcγR-mediated phagocytosis and response to nitric oxide, among others, but also in immunometabolism. Furthermore, we identified genes co-associated with the key regulators previously reported as candidate genes (e.g., ANGPT1, CD4, CD36, DOCK1, PDE4B, PRKCE, PTPRC and SH2B3) for immunity traits in humans and pigs, but also new candidate ones (e.g., ACSL3, CXADR, HBB, MMP12, PTPN6, WLS) that were not previously described. The co-association analysis revealed new regulators associated with health-related traits in pigs. This approach also identified gene-by-gene interactions and candidate genes involved in pathways related to cell fate and metabolic and immune functions. Our results shed new light in the regulatory mechanisms involved in pig immunity and reinforce the use of the pig as biomedical model.


INTRODUCTION
Health-related traits are becoming more and more relevant due to the great challenges that the pig sector has faced in recent years. The current emergence of antibiotic resistance, with pig production as one of the largest consumers of antimicrobials, and the societal demand for healthier livestock products and more sustainable production systems, make necessary to include new traits such as those related to health in the selection programs (1). Thus, health-related traits that have traditionally played a minor role in breeding programs are currently considered relevant traits to take into account.
Several studies have determined an important genetic contribution to the phenotypic variability of health-related traits, describing medium to high heritabilities for some of these traits (2)(3)(4)(5)(6)(7) and identifying genomic regions with candidate genes associated to immunocompetence (7,8). However, although genome wide association studies (GWAS) have identified single nucleotide polymorphisms (SNPs), or genes, associated with health-related traits, these SNPs generally had a little effect on the phenotypic variation of the traits under study due to their complex nature. Besides, some interesting SNPs could be overlooked due to stringent significance thresholds after multiple test correction (9). Furthermore, these studies do not consider gene-by-gene interactions or the contribution of genetic correlations between related traits.
During the last years, the use of systems biology approaches has allowed the integration of multiple sources of biological information, providing a more complete picture of the functional mechanisms involved in the determination of complex traits, and allowing the identification of candidate genes and genetic variants implicated in their phenotypic diversity (10). One of these approaches, the association weight matrix (AWM), simultaneously considers multiple phenotypes and gene-by-gene interactions by combining GWAS results with network inference algorithms (11,12). To date, multiple studies in livestock species have implemented this approach to identify gene networks and key regulatory elements implicated in the genetic architecture of meat quality (13,14), growth (15), gut microbial communities (16), reproductive (11,17), feed efficiency (18,19), and milk-related (20)(21)(22)(23) traits, among others. However, to the best of our knowledge, no studies have applied this methodology for health-related (mainly immunological and haematological) traits in livestock or other animal species.
The present study aimed to identify gene networks, key regulators and candidate genes associated to immuno-competence in pigs by integrating 30 health-related traits, enriched for innate immune phenotypes, using the AWM co-association analysis.

Animal Material and Phenotypic Parameters
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.

SNP Genotyping and Genome Wide Association Analysis
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 logtransformed 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,Gs 2 u ), where G is the genomic relationship matrix calculated using the filtered autosomal SNPs based on the methodology of (25) and s 2 u 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.

Association Weight Matrix
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, gd T cells was chosen as key phenotype and those gd 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.

Functional and Comparative Analyses
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.

Gene Co-Association Network Description
From the 42,641 SNPs used in the GWAS analysis, 2,230 SNPs were directly associated (P<0.05) with the key phenotype (gd 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 gd 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 2 >0.7) were observed for the haematological traits (MCH, MCV, HB, HCT), while CRP, MON, and IgAsal showed low heritabilities (h 2 <0.3). However, the other immunoglobulinrelated 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 2 >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 (gd 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 subclustered 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 gd 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 gd T cell number in the latter (37).

Key Transcription Factors Regulating Health-Related Traits
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.

Functional Classification and Identification of Candidate Genes Associated to Immunocompetence
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 coassociated genes were involved (Table S5). Figure 2 shows a simplified co-association network representing those immunityassociated 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 gd T cell differentiation and activation that was co-associated with ARNT2, BRMS1L and MED12L.

Comparison With Literature Reported Candidate Genes and QTLs for Health-Related Traits
When we studied in detail the top five regulators and their coassociated 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)(39)(40)(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)(39)(40). MYH9 was another candidate gene associated with haematological parameters in humans (38)(39)(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)(39)(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 gd 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, Tolllike receptor (TLR)-2 and TLR-9 levels, and basophil number, among others (Table S6), in other pig populations (35).

DISCUSSION
In the present study, a systems genetic approach using the AWM methodology was used to identify potential regulators, gene-gene interactions and pathways for immunocompetence in pigs considering 30 health-related phenotypes enriched for innate FIGURE 2 | Simplification of the full gene co-association network formed by 171 nodes and 832 interactions representing the immunity-associated genes (orange) that were correlated with the five top regulator genes (green; ARNT2, BRMS1L, MED12L, SUPT3H, and TRIM25). immune phenotypes. We selected gd T cells as the key phenotype, since these cells play important roles at the interface between innate and adaptive immunity and had the potential to impact other phenotypes analysed in our population. In fact, in a previous study in the same population (7), positive and negative genetic correlations were identified between this phenotype and acute phase proteins and platelets (positive correlations), and IgAsal, PHAGO-FITC and cortisol levels (negative correlations). Similar functions such as response to pathogens and tumours, wound healing, tissue homeostasis, stressed cell clearance, and antigen presentation, among others, have been described for gd T cells in several species [reviewed in Holderness et al. (52)]. However, compared to humans and mice, this population of cells is strongly enriched in pig blood (53), which also suggest specific immune responses and functions.
After GWAS analysis and the SNP annotation step, only 3,544 SNPs were retained for further analysis. When compared with the heritability values obtained with the whole dataset of SNPs, this subset of 3,544 SNPs helped to explain a greater proportion of the additive variance of the analysed traits. As expected, the heritability of our key phenotype (gd T cells) was greatly improved from moderate to high after using the reduced dataset, despite other phenotypes such as HB and HCT had a greater improvement and a higher heritability value. This may also be due to the fact that they were the two phenotypes with the highest number of associated SNPs after the key phenotype. Among the associated SNPs with HB and HCT, we identified a pleiotropic SNP (rs81284215) that was associated with more than 10 phenotypes and was located in the HBB gene. It is well established the relationship between HB levels and growth on the early-stages of life (54,55). Similarly, a positive association between HB and HCT and average daily weight gain has also been described in pigs (56). Therefore, this SNP could be a useful marker to determine the haematological status of our animals, although further analysis will be required to determine the relationship of this polymorphism with weight gain in piglets.
Interestingly, the AWM additive effects correlations of the 3,544 SNPs primarily clustered health-related traits according to their nature and function and were informative enough to capture the genetic associations among the health-related traits in a similar way than the published genetic correlations calculated for the same population (7). Of note, there was a lack of genetic correlation among our key phenotype and the rest of the health-related phenotypes with 225 SNPs exclusively associated with gd T cells. However, we also identified 534 SNPs associated with gd T cells and one or more health-related phenotypes, highlighting the genetic interplay between several immune pathways and functions. A possible explanation for these results is the multitasking role of gd T cells in the immune system with a conventional adaptive role, but with increasing relevance in their influences on innate immunity (57,58). In fact, it has been proposed a direct but also and indirect role of gd T cells in promoting monocyte and macrophage differentiation and macrophage recruitment during inflammatory response, as well as in macrophage elimination during the resolution of inflammation (57). Additionally, a specialized role of epithelial gd T cells in wound repair by promoting platelet and neutrophil accumulation in the limbal vessels has been proposed (59). Overall, these functions are the reflection of the results found during this work that will be discussed below.
In this study, AWM co-association analysis allowed the identification of five key regulators (ARNT2, BRMS1L, MED12L, SUPT3H and TRIM25) for immunocompetence in pigs. It is relevant to note that for most of these regulators (BRMS1L, MED12L, SUPT3H) there is still little information in the literature about their relationship with immunity-related functions, despite all regulators were located within QTL regions associated to immunity phenotypes (pigQTLdb). Among the identified top regulators, BRMS1L and TRIM25 gathered the highest number of associations with other genes. While BRMS1L has been described as a metastatic suppressor gene (60), TRIM25 has been involved in numerous cellular processes including regulation of innate immune response against viral infection (61,62). Specifically, TRIM25 is involved in the RIG-I/interferon pathway in response to viral RNAs (63). Among the list of genes regulated by TRIM25, we identified two genes (MMP12 and PTPN6) related to regulation of type I interferon-mediated signaling pathway. Another key regulator that also plays a role in the innate immune system was ARNT2. This transcription factor forms heterodimers with the aryl hydrocarbon receptor (AHR), that is expressed in all gd T cell lymphocyte subsets, being essential for some aspects of their functions [reviewed in Stange and Veldhoen (64)]. This complex regulates the transcription of multiple immune-related genes and the activity and differentiation of phagocytic cells and lymphocytes (65). Accordingly, this transcription factor was associated with leukocyte, lymphocyte and monocyte counts in our population and was co-associated with the PTPRC gene that was related to gd T cell differentiation. Furthermore, a polymorphism in ARNT2 has been associated with an impaired fungicidal activity by depleting the phagocytic activity of macrophages (66). Notably, in our gene-network we identified ARNT2-regulated genes involved in FcgR-mediated phagocytosis (ACTR3C, AKT2, DOCK1, GAB2, PRKCE, PTPRC) and response to nitric oxide (AIFM1, CCR7, MTR, TXN2) processes. Therefore, these genes are interesting candidate genes to contribute to the phenotypic differences observed for these health-related traits. Consistently with our results, some of these genes (AKT2, PRKCE, PTPRC) were previously associated with health-related traits in humans and/ or pigs, and seven of the ten genes have been at least associated with a phagocytosis-related trait in our population. MED12L is a paralog of the MED12 gene, which is required to activate innate immunity genes in Drosophila (67). A recent GWAS study identified a variant within the MED12L gene associated with viral infection response (Epstein-Barr virus (EBV) nuclear antigen (EBNA)) in humans (68). In the network, the MED12L gene was also interacting with genes associated to gd T cell differentiation (PTPRC) and activation (CXADR and PTPRC). Finally, a genetic study found an association between a polymorphism in the SUPT3H gene with the ADAMTS13 activity, a protein with antithrombotic properties (69). Except for some genes directly related to the key phenotype (gd T cell), the five regulators and their 1,828 co-associated genes were involved in a wide variety of biological functions, some of them directly related to immunity but also to immunometabolism. Among them, there were functions related to lymphocytes (differentiation and activation), platelets (activation, degranulation and megakaryocyte differentiation), phagocytosis (FcgR-mediated phagocytosis), NO (response to nitric oxide), and to a lesser extent other functions associated with other immune phenotypes (such as those related to haematological related traits) were identified. Remarkably, for some of these phenotypes (platelets count, NO in serum, and some phagocytosis traits) we did not identify any significantly associated genomic region after performing GWAS (7). Here, it is relevant to highlight that for some of these traits, we identified more than one of the top transcription factors coassociated with genes related to the same functions. Therefore, a plausible explanation for our results would be that due to the complex nature of immunity traits with the involvement of multiple genes and regulators, GWAS analyses are not able to capture all the significant regions associated with their phenotypic variation, since many of the biologically relevant polymorphisms are lost after multiple test corrections. For instance, it is worthy to note the overrepresentation of platelet associated functions, which were commonly found for all the key identified regulators. Apart from being important for coagulation and fibrinolysis processes, platelets are also activators of the immune system and mediators of inflammatory response [reviewed in Trzeciak-Ryczek et al. (70)], therefore participating in multiple functions. In fact, in the hierarchical analysis, platelets phenotype grouped within a heterologous group composed by stress indicators, acute phase proteins and immunoglobulins among others, which may be indicative of its multifunctional nature. Another overrepresented pathway regulated by all the key transcription factors was phospholipase D (PLD) signaling pathway. PLD is an enzyme that catalyses the hydrolysis of phosphatidylcholine, to produce the signal molecule phosphatidic acid (PA). PLD participates in multiple functions of the immune system mediating phagocytosis by macrophages, activation of NADPH oxidation in neutrophils, mast cell degranulation and T cell activation [reviewed in Zhu et al. (71)]. Finally, it is worth mentioning some specific immune functions for the poorly characterized BRMS1L and SUPT3H regulators; for instance, the regulation of the NOTCH signaling pathway identified for BRMS1L. NOTCH signaling plays a pivotal role in T cell fate decision, influencing the differentiation between ab T cells versus gd T cells (72,73). Furthermore, SUPT3H seems to have an active role in immunometabolism by controlling signaling pathways such as AMPK, mTOR and PI3K-Akt. In recent years, immunometabolism has emerged as an important mechanism affecting innate and adaptive immune system. In fact, AMPK and PI3K/Akt/mTOR regulate metabolic pathways important for cell fate and function in the immune system (74)(75)(76). Overall, these results underline the relevance of the five key transcription factors identified in our network regulating important functions of the immune system in pigs. Furthermore, and taking into account the plethora of different functions that these pathways/molecules perform, we could hypothesize, that depending on the function, either one or the interaction of several key transcription factors will be required.
Among the genes co-associated with the top five regulators, we also identified candidate genes that have been previously described as candidate genes for health-related traits in humans and pigs, some of them associated with similar immunological and haematological traits (ANGPT1, PTPRC, CD4, KLRK1, and FTO). In addition, there were also candidate genes with functions directly related to the associated traits. Indeed, PRKCE, SH2B3 and CD36 were associated in our study with phagocytosis-related traits. PRKCE plays a role in FcgR-mediated phagocytosis (77), while CD36 is a phagocytic receptor that mediates the uptake of different substrates (78,79). SH2B3 is involved in macrophage function (80), whereas MYH9, which encodes a conventional non-muscle myosin and is required for normal T lymphocyte migration (81), was associated with the quantity of lymphocytes. ANGPT1, a candidate gene in humans, belongs to the angiopoietin family regulating vascular integrity and angiogenesis through endothelium-derived NO release (82), and accordingly, it was also found in association with NO and some haematological traits in our pig population. Also, it is relevant to highlight that the candidate gene PTPRC was associated with the quantity of lymphocytes in blood. This member of the protein tyrosine phosphatase (PTP) family is expressed in all nucleated haematopoietic cells and plays an essential role in T and B cell antigen receptor signaling, and it has been associated to many autoimmune disorders [reviewed in Ulloa-Aguirre et al. (83)]. Finally, PDE4B, which was associated with haematological and phagocytosis-related traits in our study, is essential for the LPS-activated immune response of phagocytic cells (82,84).
In conclusion, the AWM gene co-association analysis applied in the present study represented an important step to better characterize the genetic architecture of complex traits for immunity in pigs. The network-based approach allowed the identification of five key transcription factors (ARNT2, BRMS1L, MED12L, SUPT3H and TRIM25) regulating pathways directly involved in cell fate and in metabolic and immune functions. Furthermore, we have identified gene-gene interactions and candidate genes, some of which have been previously associated to immunity traits in humans and pigs. This study has been the first one to use AWM for analysing health-related traits in pigs and supports the effectiveness of using network-based approaches to identify potential regulators and candidate genes, representing a major step in understanding the regulatory mechanisms involved in porcine immunity, and reinforcing the use of the pig as biomedical model.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
All experimental procedures with pigs were performed according to the Spanish Policy for Animal Protection RD 53/2013, which meets the European Union Directive 2010/63/EU about the protection of animals used in experimentation. The experimental protocol was approved by the Ethical Committee of the Institut de Recerca i Tecnologia Agroalimentàries (IRTA).

AUTHOR CONTRIBUTIONS
RQ and MB designed the study. MB supervised the generation of the animal material used in this work. YR-C, OG-R, MP, RQ and MB performed the sampling. OG-R and MB carried out the laboratory analyses. DC-P, YR-C, RQ and MB analysed the data and interpreted the results. DC-P and MB wrote the manuscript. All authors contributed to the article and approved the submitted version.