TYPE Original Research PUBLISHED 06 February 2023 DOI 10.3389/fimmu.2023.1058346 Genetic architecture of innate OPEN ACCESS and adaptive immune cells in pigs EDITED BY John C. Schwartz, The Pirbright Institute, United Kingdom Maria Ballester 1*, Teodor Jové-Juncà 1, Afra Pascual 1, REVIEWED BY Sergi López-Serrano 2,3, Daniel Crespo-Piazuelo 1, Jianzhong Zhu, Yangzhou University, China Carles Hernández-Banqué 1, Olga González-Rodrı́guez 1, Haifei Wang, Yuliaxis Ramayo-Caldas 1 and Raquel Quintanilla 1* Yangzhou University, China Rayner Gonzalez-Prendes, 1 Animal Breeding and Genetics Program, Institute of Agrifood Research and Technology (IRTA), Torre Universitat de Lleida, Spain Marimon, Caldes de Montbui, Spain, 2 Unitat mixta d’Investigació IRTA-UAB en Sanitat Animal, Centre de Recerca en Sanitat Animal (CReSA), Universitat Autònoma de Barcelona (UAB), Bellaterra, *CORRESPONDENCE Catalonia, Spain, 3 Institute of Agrifood Research and Technology (IRTA), Programa de Sanitat Animal, Maria Ballester Centre de Recerca en Sanitat Animal (CReSA), Universitat Autònoma de Barcelona (UAB), Bellaterra, maria.ballester@irta.cat Catalonia, Spain Raquel Quintanilla raquel.quintanilla@irta.cat SPECIALTY SECTION This article was submitted to Pig industry is facing new challenges that make necessary to reorient breeding Comparative Immunology, a section of the journal programs to produce more robust and resilient pig populations. The aim of the Frontiers in Immunology present work was to study the genetic determinism of lymphocyte subpopulations RECEIVED 30 September 2022 in the peripheral blood of pigs and identify genomic regions and biomarkers ACCEPTED 16 January 2023 associated to them. For this purpose, we stained peripheral blood mononuclear PUBLISHED 06 February 2023 cells to measure ten immune-cell-related traits including the relative abundance CITATION of different populations of lymphocytes, the proportions of CD4+ T cells and CD8+ Ballester M, Jové-Juncà T, Pascual A, López-Serrano S, Crespo-Piazuelo D, T cells, and the ratio of CD4+/CD8+ T cells from 391 healthy Duroc piglets aged 8 Hernández-Banqué C, weeks. Medium to high heritabilities were observed for the ten immune-cell- González-Rodrı´guez O, Ramayo-Caldas Y and Quintanilla R (2023) Genetic related traits and significant genetic correlations were obtained between the architecture of innate and adaptive proportion of some lymphocytes populations. A genome-wide association study immune cells in pigs. pointed out 32 SNPs located at four chromosomal regions on pig chromosomes Front. Immunol. 14:1058346. doi: 10.3389/fimmu.2023.1058346 SSC3, SSC5, SSC8, and SSCX as significantly associated to T-helper cells, memory COPYRIGHT T-helper cells and gd T cells. Several genes previously identified in human © 2023 Ballester, Jové-Juncà, Pascual, association studies for the same or related traits were located in the associated López-Serrano, Crespo-Piazuelo, regions, and were proposed as candidate genes to explain the variation of T cell Hernández-Banqué, González-Rodrı´guez, Ramayo-Caldas and Quintanilla. This is an populations such as CD4, CD8A, CD8B, KLRC2, RMND5A and VPS24. The open-access article distributed under the transcriptome analysis of whole blood samples from animals with extreme terms of the Creative Commons Attribution License (CC BY). The use, distribution or proportions of gd T, T-helper and memory T-helper cells identified differentially reproduction in other forums is permitted, expressed genes (CAPG, TCF7L1, KLRD1 and CD4) located into the associated provided the original author(s) and the regions. In addition, differentially expressed genes specific of different T cells copyright owner(s) are credited and that the original publication in this journal is subpopulations were identified such as SOX13 and WC1 genes for gd T cells. Our cited, in accordance with accepted results enhance the knowledge about the genetic control of lymphocyte traits that academic practice. No use, distribution or could be considered to optimize the induction of immune responses to vaccines reproduction is permitted which does not comply with these terms. against pathogens. Furthermore, they open the possibility of applying effective selection programs for improving immunocompetence in pigs and support the use of the pig as a very reliable human biomedical model. KEYWORDS immunocompentence, immune cells, RNA-Seq, biomodel, gd T cells, pig Frontiers in Immunology 01 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 1 Introduction Duroc commercial line by estimating their genetic parameters and identifying genomic regions, candidate genes and genetic markers Pork is one of the most consumed meats worldwide. The associated to their phenotypic variation. Finally, knowing the industrialization and specialization of the swine production coupled relevance of the pig as a human biomodel, a comparison between with the ever-increasing movement of pigs and pork products at both species was performed with the significantly associated regions national and international scales facilitate the appearance and spread found in pigs. of endemic and emerging pathogens which causes significant economic losses to the swine industry. Furthermore, the restriction in the use of antimicrobials, due to the appearance of drug-resistant 2 Materials and methods bacteria, and the threats of climate change, stresses the need for guiding breeding programs to produce more robust and disease 2.1 Ethics approval resistant pig populations (1). Along with conventional methods such as vaccination and strict biosecurity measures, improving the All experimental procedures with pigs were performed according overall immunocompetence of healthy animals is a good approach to to the Spanish Policy for Animal Protection RD 53/2013, which meets prevent and control infectious pathogens (2, 3). the European Union Directive 2010/63/EU about the protection of In this context, immunity traits (ITs) are considered biologically animals used in experimentation. The experimental protocol was relevant parameters to measure immunocompetence (4). ITs may be approved by the Ethical Committee of the Institut de Recerca i divided in the two major components of the immune system, the Tecnologia Agroalimentàries (IRTA). innate or natural immunity and the acquired or adaptive immunity. The innate immunity is composed by non-specific components that constitute the first line of host defence and recognize a wide range of 2.2 Animal material pathogens through a restricted set of pattern recognition receptors. Cell types involved in innate immunity include Natural Killer (NK) A total of 391 healthy piglets aged 60 ± 8 days (males and females) cells, dendritic cells, monocytes or gd T cells (5). In contrast, belonging to a commercial Duroc line were used in this study. All specialized cells of the acquired immune system can provide animals were vaccinated at weaning (aged ~26 days old) against immunological memory after an initial response mediated by a very PCV2 with Porcilis® PCV. The pigs came from six batches (72 ± 1 large number of specific receptors for antigens present in pathogens animals per batch) and were raised in the same farm and ad libitum (6). While T lymphocytes are the effectors of adaptive cellular fed with a commercial cereal-based diet. All animals were apparently immune responses, antibody-producing cells, B lymphocytes, healthy, without any observable signs of infection. mediate adaptive humoral immunity (6). Blood was collected via the external jugular vein into vacutainer The immune system has a high degree of interindividual variation tubes with anticoagulants (Sangüesa S.A., Spain), which required the that is controlled by genetic and environmental factors (7). However, restraint of the animals but not their sedation. Genomic DNA was leukocyte subsets, mostly the adaptive ones, are highly affected by extracted from blood samples using the NucleoSpin ® Blood genetics (8–10). Genetic studies have estimated moderate to high (Macherey-Nagel, Germany). DNA concentration and purity were heritabilities for both innate (e.g. NK cells and gd T cells) and adaptive measured in a Nanodrop ND-1000 spectrophotometer. (e.g. CD4+ T lymphocytes, CD8+ T lymphocytes and B lymphocytes) immune cells under different conditions in pigs (10–12), suggesting that selection for those traits is a plausible strategy for the swine industry (13). 2.3 Immunophenotyping Furthermore, the development of high-density genotyping arrays and in- depth immune phenotyping techniques have allowed the identification Peripheral blood mononuclear cells (PBMCs) isolated from of candidate genes and genetic variants associated with the phenotypic heparinized peripheral blood by density-gradient centrifugation with variation of immune cells in humans (9, 14) and pigs (15, 16). However, Histopaque-1077 (Sigma, Spain) and previously stored at -80°C from a genome-wide association studies (GWAS) for innate and adaptive study for gd T-lymphocyte quantification (12), were herein employed immune cells in pigs are still scarce. The few studies on this topic have for immune phenotyping. Frozen PBMCs were thawed and been focused on identifying genetic variants that affect T-cell resuspended in RPMI 1640 medium supplemented with 5% Foetal subpopulations in animals vaccinated at 21 days of age with live Bovine Serum (FBS) (Sigma, Spain), 1% Penicillin-Streptomycin classical swine fever (CSF) vaccine (15, 16). (10,000 U/ml - 10 mg/mL) and 1% L-Glutamine (200 mM) (Cultek, The objective of this work was to study the genetic determinism of Spain) after centrifugation at 450 g at 4°C for 10 min. both innate and adaptive immune cells traits in healthy pigs of a Cell concentration was adjusted using RPMI to obtain around 800,000 cells per each sample. Once adjusted, all samples in duplicate were viability stained (Live-or-dyeTM Aqua 405, Biotium, Fremont, Abbreviations: CTL, cytotoxic T cells; DE, differentially expressed; FBS, foetal CA, USA) diluted in PBS for 25 min at 4°C following the bovine serum; FC, fold change; FDR, false discovery rate; GCTA, genome-wide manufacturer’s indications. After two washes, PBMCs were labelled complex trait analysis; GO, gene ontologies; GWAS, genome-wide association with the primary-conjugated antibody mixture in PBS-1%FBS for studies; h2, heritability; H, high; IT, immunity traits; L, low; MAF, minor allele 25 min at 4°C (Supplementary Table 1). frequency; NK, natural killer; PBMCs, peripheral blood mononuclear cells; Q-Q, After two washes with 1xPBS-1% FBS at 4°C, cells were quantile-quantile; Rg, genetic correlations; Rp, phenotypic correlations. resuspended in 1xPBS-1% FBS and analysed by flow cytometry Frontiers in Immunology 02 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 using the MACSQuant Analyzer 10 Flow cytometer (Miltenyi Biotec 2.5 Estimation of genetic parameters GmbH, Bergisch Gladbach, Germany) and the Flowlogic™ software v7.3 (Inivai Technologies, Melbourne, Australia). Beside test samples, The heritability (h2) , i.e. the proportion of phenotypic variance unstained cells, stained samples with sole viability marker, isotypes attributable to additive genetic effects, was estimated for the for each antibody subclass, and Fluorescence Minus One (FMO) percentage of the different lymphocyte subsets of cells showed in stained samples were included as controls to adjust the analysis and Table 1. Variance components and the corresponding h2 were discard false positive results. estimated from an univariate animal model as follows: The gating strategy per each animal used for this calculation is Y=Xb+Zu+e shown in Figure S1. Absolute number of events were used to calculate where Y is the vector of phenotypes, the percentage of lymphocyte percentages of lymphocyte populations depending on PBMCs gate. subpopulations (either transformed or raw data) of all individuals; b is After discarding the dead cells gated by the use of the viability the vector of systematic (fixed) effects on the trait, including when staining, the following cell subsets were quantified: B lymphocytes significant the effects of sex (2 levels) and/or batch (6 levels); u is the (CD21+), T lymphocytes (CD3+), natural killer (NK) cells (CD3- vector of animal’s genetic additive (random) effects; X and Z are the CD21-CD8+), cytotoxic (CTL) T cells (CD3+CD4-CD8+), T helper corresponding incidence matrices for b and u ; and e is the vector of cells (CD3+CD4+CD8-), memory T-helper cells (CD3+CD4+CD8+) random residual terms. The assumed distribution of additive genetic and naïve T cells (CD3+CD4-CD8-). In addition, we quantified the effects was u∼(0, Asu2 ), where A is the numerator relationship matrix total proportion of CD3+CD4+, CD3+CD8+ and the ratio of CD4+ to computed on the basis of pedigree and su2 is the additive genetic CD8+ (CD4+/CD8+). variance; random errors were distributed as e∼N(0,Ise2 ). Estimation of the model variance components and the corresponding heritability (h2 = su2 =(su2 + se2 )) for each lymphocytes subset was performed by REML using the aireml program included in the BGF90 package (17). 2.4 Statistics of PBMCs subpopulations The standard errors (SE) of the heritability estimates were also computed, thus obtaining the corresponding confidence intervals at Basic descriptive statistics of absolute counts and proportion of 95% (CI95). the different analysed lymphocyte subsets of cells among PBMCs Subsequently, pairwise genetic correlations (for each two traits are shown in Table 1. Normality of each PBMC subpopulation data combination) were estimated in a two-traits animal model described was checked through Shapiro-Wilk test in R. For some of the as follows: variables, logarithm or square root transformations were applied " # " #" # " #" # " # to reach normal distribution of residuals (p-value > 0.05). Y t1 X t1 0 b t1 Z t1 0 ut1 et1 A linear regression model was then applied on transformed and = + + Y t2 0 X t2 b t2 0 Z t2 ut2 et2 raw data using the lm() function in R with sex and batch as fixed effects to test their significance on each cell subset through likelihood where Yt1 and Yt2 are the vectors of phenotypic observations for ratio tests. When significant, sex and/or batch effects were considered trait 1 and trait 2, respectively; bt1 and bt2 are the vectors of systematic for subsequent analyses. (fixed) effects on each trait previously described, and Xt1 and Xt2 the Pairwise phenotypic correlations (rp) among all analysed correspondent incidence matrices; ut1 and ut2 are the vectors of phenotypes were computed after adjusting for significant animal genetic additive effects on trait 1 or trait 2 (random effects), systematic factors. and Zt1 and Zt2 the corresponding incidence matrices; finally et1 and TABLE 1 Descriptive statistics of absolute numbers and proportions of peripheral blood mononuclear cells (PBMCs). Cell counts % of PBMCs Trait Mean SD CV Mean SD CV T lymphocytes 10327.01 2466.48 0.24 43.81 9.19 0.21 CTL cells 4563.98 2039.04 0.45 19.41 8.52 0.44 Memory T cells 748.85 671.54 0.90 3.18 2.79 0.88 T helper cells 428.55 216.16 0.50 1.84 0.95 0.52 Naïve T cells 4585.62 2035.98 0.44 19.37 8.23 0.43 Natural killer (NK) cells 3346.56 1479.30 0.44 14.16 6.12 0.43 B lymphocytes 5048.28 1938.44 0.38 21.22 7.84 0.37 + CD4 T cells 1177.41 776.91 0.66 5.02 3.25 0.65 + CD8 T cells 5312.84 2294.21 0.43 22.59 9.49 0.42 CD4+/CD8+ 0.12 0.09 0.78 – – – standard deviation (SD), coefficient of variation (CV). Frontiers in Immunology 03 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 et2 are the vectors of residual errors for each trait. The (co)variance using the Ensembl Genes 106 Database of Scrofa11.1 reference matrix of random genetic effects was defined as: assembly. Gene features and gene ontologies (GO) were retrieved. " # " # For those pig genes without a name, orthologous human gene names ut1 Asu1 2 Asu1,u2 were retrieved. Furthermore, information from Genecards (24) and Var = ut2 Asu1,u2 Asu1 2 Mouse Genome Database (25) was used to identify gene functions affecting the analysed phenotypes. where su12 and su2 2 are the additive genetic variance of traits 1 and 2, respectively, su1,u2 is the genetic covariance between the traits, and A is the numerator relationship matrix as defined above. Estimation of the 2.8 RNA-Seq of blood (co)variance components of each pairwise analysis was also performed by REML using the aireml program included in the BGF90 package (17).  RNA-Seq of whole blood was performed for 30 pigs classified as Genetic correlations (rg) between traits were obtained as rg= (su1,u2 su1 su2 low (L) or high (H) percentage of T gd (L=12, H=12), T helper and ), and the SE of the genetic correlation estimates were also computed. memory T-helper cells (L=15, H=15), which were selected from the 391 pigs of the Duroc population by a principal component analysis 2.6 SNP genotyping (PCA; Figure S2). PC2 was considered to classify 30 animals extreme for T helper (H group: 3.49 ± 0.25; L group: 0.76 ± 0.06) and memory The GGP Porcine HD Array (Illumina, San Diego, CA) was T helper cells (H group: 9.92 ± 1.16; L group: 1 ± 0.18). Since these employed to genotype 68,516 single nucleotide polymorphisms two populations showed the same direction in the PCA, we (SNPs) in the 391 animals using the Infinium HD Assay Ultra considered the same animals per groups. Subsequently, these 30 protocol (Illumina). Plink software (18) was used to remove those animals were also classified for extreme proportions of T gd cells SNPs with a minor allele frequency (MAF) lesser than 5%, SNPs with considering the phenotype by (12), where PBMCs were stained with a more than 10% missing genotypes, and SNPs that did not map to the specific antibody for T gd cells (APC Rat Anti-Pig gd T Lymphocytes, porcine reference genome (Sscrofa11.1 assembly). Finally, a subset of clone MAC320, BD Biosciences). For this analysis, we selected 12 42,641 SNPs remained for further analysis. animals per group to maximize the differences in the proportions of T gd cells (H group: 16.36 ± 1.86; L group: 3.64 ± 0.51). Total RNA from blood was isolated and purified with a Tempus™ 2.7 Genome wide association study (GWAS) Spin RNA Isolation Kit (Thermo Fisher Scientific, Spain), following the manufacturer’s instructions. Total RNA concentration was GWAS was carried out between the 42,641 filtered SNPs and the measured in a Nanodrop One spectrophotometer (Thermo Fisher 10 traits related to the percentages of PBMCs subpopulations, Scientific, Spain) and RNA purity and integrity was checked by using described in Table 1. To this end, the genome-wide complex trait a Fragment Analyzer equipment (Agilent Technologies, INC., Santa analysis (GCTA) software (19) was used following an additive model Clara, CA). Libraries were prepared using a Stranded total RNA Prep for each trait across all SNPs: with Ribo-Zero Plus kit (Illumina Inc., CA) and were paired-end sequenced (2 × 100 bp) on a NovaSeq 6000 platform (Illumina Inc., Y = Xb + Zg + Sl al + e CA) at the National Centre for Genomic Analysis (CNAG, Barcelona, where Y corresponds to the vector of phenotypic observations Spain). More than 50 million of paired-end reads were obtained for all (log or sqrt transformed or raw data) for each analysed trait; b is the samples. The quality of the raw sequenced reads in the FASTQ files vector of systematic effects on the trait described above; g is the vector was analysed with the FASTQC software (Babraham Bioinformatics, of infinitesimal genetic effects of each individual, with distribution g http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads ∼N(0,Gsg2 ), where G is the genomic relationship matrix (GRM) were mapped to the reference pig genome Sscrofa11.1 and the calculated using the filtered autosomal SNPs based on the annotation database Ensembl Genes 106 by using STAR v. 2.75.3a. methodology of (19) and sg2 is the additive genetic variance; Sl is Transcript quantification was performed with RSEM v. 1.3.0. The the vector of genotypes (coded as 0, 1, 2) of each individual for the lth EdgeR R package (26) was used to identify differentially expressed SNP, and al is the allele substitution effect of the lth SNP on the trait (DE) genes between the two divergent groups (L vs H). After under study. The FDR method of multiple testing described by correcting for multiple testing, genes with a Fold Change (FC) Benjamini and Hochberg (20) was used to measure the statistical between groups greater than 1.2 (i.e. |log2FC| > 0.26) and an FDR< significance for association studies at genome-wide level with the 0.1 were classified as DE. Enriched GO terms and pathway analysis p.adjust function of R. The significant association threshold was set at (padj<0.1) of the DE genes was performed with the ClueGO v2.5.8 FDR ≤ 0.1. Manhattan plots based on the significance of the plug-in of Cytoscape v3.9.1 (27). associations across the whole genome were generated using the qqman R package (21). Quantile-quantile (Q-Q) plots of the p- value distribution were generated with the CMplot R package 3 Results available at: https://github.com/YinLiLin/CMplot. VEP software (22) was used to locate all the significant SNPs and 3.1 Characterization of immune phenotypes determine their functional effects. Gene annotation for 1Mb downstream/upstream of genomic intervals around the most In the present study, we stained PBMCs to measure ten immune- significant SNPs was performed with BIOMART software (23) cell-related traits including the relative abundance of different Frontiers in Immunology 04 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 populations of lymphocytes, the proportions of CD4+ T cells and immunity traits previously analysed in the same individuals (12), CD8+ T cells, and the ratio of CD4+/CD8+ T cells from 391 including plasma immunoglobulins and acute phase proteins individuals of a commercial Duroc pig line. All proportions and concentrations, gd T lymphocytes, haemogram figures and absolute numbers of cell subsets with descriptive statistics are phagocytosis capacity, were also computed (Supplementary reported in Table 1. The T cells were the most abundant Table 2). Among these traits, only the gd T lymphocytes showed a population of lymphocytes, representing 43.8% of PBMCs. Among strong and highly significant phenotypic association with naïve T cells T lymphocytes subpopulations, the CTL and naïve T cells were the (rp = 0.81, p-value = 2.2×10-16), indicating that naïve T cells (CD3+ most abundant (~19% of PBMCs each one), whereas the memory plus CD4- CD8-) corresponded mainly to gd T cells. The gd T cells helper T cells (i.e. CD4+ T cells) represented 5% of PBMCs. The mean proportion also showed a moderate correlation with T cells and proportion of B lymphocytes (21.2% of PBMCs) was less than half certain antagonism with NK and CTL cells proportions. Finally, the that of T lymphocytes, whereas the NK cells represented 14.16% of percentage of phagocytic cells and the phagocytosis capacity of PBMCs. A high variability between individuals was observed in the lymphocytes showed low but significant positive associations with relative abundance of some PBMC types, the CV ranging from 0.21 to the percentage of B cells, jointly with low negative correlations with T 0.89. The percentage of T lymphocytes showed the lowest dispersion cells (CTL and Memory T cells) and NK relative abundances. (CV=0.21), although within T cells, the percentage of memory T cells showed the highest variability (CV=0.88) among the analysed PBMCs subtypes. 3.2 Genetic parameters of T and Phenotypic correlations were also computed to analyse the B-cell populations association between the relative abundance of the different lymphocyte subsets (Supplementary Table 2). Positive correlations The heritability of the different cell immunity subsets was were observed between proportions of T cells and its different estimated to ascertain the genetic determinism of cellular subpopulations, ranging from 0.34 to 0.51, but for T helper cells, immunocompetence traits. Medium to high heritabilities were that exhibited a low correlation with T cells percentage. Among T observed for the ten immune-cell-related traits comprising the cells subtypes, the relative abundance of memory T-helper cells relative abundance among PBMCs of different lymphocyte correlated positively with both T-helper and CTL cells (rP= 0.39 populations (Table 2). The heritability estimates ranged between and 0.32), whereas an antagonism between the proportions of CTL 0.361 to 0.841, being the proportion of memory T-helper cells (h2 = and naïve T cells was observed (rP= -0.50). The relative abundance of 0.841) the most heritable trait, followed by the total proportion of T B lymphocytes showed a strong negative correlation with total T cells (h2 = 0.771). In contrast, the percentage of B cells exhibited the lymphocytes and CTL proportions of PBMCs (rP= -0.56 and -0.50). lowest heritability value (h2 = 0.361); despite moderate, the 95% NK cells proportion was positively correlated to CTL abundance but confidence interval for its heritability did not encompasses zero showed negative correlations with the relative abundances of naïve T (Table 2). The NK cells abundance and the proportion of different cells and B lymphocytes (rP= -0.51 and -0.45, respectively). As T cells subtypes except for memory T-helper cells showed medium expected, the proportion of CD4+ T cells was highly correlated with heritabilities between 0.470 and 0.558. T helper and memory T-helper cells. The proportion of CD8+ T cells When we inferred the genetic correlations between each pairwise was highly correlated with CTL T cells and moderately correlated combination of traits, positive and negative genetic associations were with T cells abundance, whereas it showed antagonism with Naïve T identified between the proportion (among PBMCs) of some cells and B lymphocytes proportions. lymphocytes populations (Figure 1 and Supplementary Table 3). In Additionally, phenotypic correlations of the different general the genetic correlation pattern matched with that observed lymphocytes subtypes abundance with a plethora of innate with phenotypic correlations, but stronger associations were observed TABLE 2 Heritability estimates (h2) for the relative abundance (among PBMCs) of different cell immunity subsets, plus standard errors (SE) and confidence intervals at 95% (CI95) of the h2 estimates. Trait h2 SE CI95 T lymphocytes (% of PBMCs) 0.771 0.152 0.473 - 1.070 CTL cells (% of PBMCs) 0.563 0.155 0.259 - 0.868 Memory T cells (% of PBMCs) 0.841 0.161 0.526 - 1.157 T helper cells (% of PBMCs) 0.470 0.148 0.181 - 0.759 Naïve T cells (% of PBMCs) 0.481 0.144 0.198 - 0.764 B lymphocytes (% of PBMCs) 0.361 0.165 0.036 - 0.685 Natural killer cells (% of PBMCs) 0.558 0.168 0.228 - 0.887 + CD4 T cells (% of PBMCs) 0.725 0.159 0.413 - 1.036 CD8+ T cells (% of PBMCs) 0.487 0.158 0.178 - 0.796 CD4+/CD8+ 0.462 0.141 0.186 - 0.738 Frontiers in Immunology 05 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 at genetic level despite its significance being limited in a number of were mainly gd T cells, jointly with a genetic antagonism between the cases due to the sample size of the population. Positive and significant proportion of gd T and T helper lymphocytes subpopulation (rg = genetic correlations were obtained between the proportion of T cells -0.568). The CTL abundance showed positive genetic correlation with and the different T cells subpopulations (rg from 0.50 to 0.68) with the the number of lymphocytes and other white blood cell counts in the exception of T helper cells. Among T cells subpopulations, a positive haemogram (rg from 0.622 to 0.702) as well as with haematocrit (rg = genetic association of memory T-helper cells with T helper (rg = 0.58) 0.741). Besides, negative genetic correlations of the proportion of B and CTL (rg = 0.37) cells was estimated, whereas naïve T cells were lymphocytes with both the phagocytic capacity and the percentage of negatively correlated with CTL and T helper cells (rg = -0.37 and phagocytic cells (rg = -0.803 and -0.733, respectively) were obtained, -0.36). The proportion of B lymphocytes exhibited negative genetic in contrast with the positive correlation with the percentage of correlations with the rest of lymphocytes subsets, particularly with the phagocytic lymphocytes (rg = 0.544). total percentage of T cells (rg = -0.78), CD4+ T cells (rg = -0.56) and CTL T cells (rg = -0.41). NK cells proportion were negatively correlated with memory T-helper cells (rg = -0.44) and to a lesser 3.3 Genomic regions and candidate genes extent with B lymphocytes (rg = -0.32) proportion. associated with immune-cell traits As far as genetic correlation of the analysed PBMCs subsets with innate immunity and other health-related traits (Supplementary To identify genomic regions associated with lymphocytic traits, a Table 3), we identified a strong genetic correlation between the GWAS was performed using the ten immune-cell-related traits and proportion of gd T lymphocytes and naïve T cells (rg = 0.903) the genotypes of 42,641 SNPs of the Porcine GGPSNP70 BeadChip subpopulation, which corroborates that CD3+ CD4- CD8- T cells (Illumina) in 391 Duroc pigs. Significant associations at whole- genome level (FDR ≤ 0.1) were detected for T helper cells, memory T-helper cells and naïve T cells (Table 3). A total of 32 significantly associated SNPs located at four chromosomal regions on pig chromosomes SSC3, SSC5, SSC8 and SSCX were identified (Supplementary Table 4). All the associated SNPs, with their predicted consequences, are shown in Supplementary Table 5. In addition, the identified genomic region on SSC5 was also associated with the total proportion of CD4+ T cells at FDR threshold< 0.2 (Table 3 and Figure S3). Graphical representation displayed in Manhattan plots of the GWAS results for T helper (A), memory T-helper (B) and naïve T (C) cells are shown in Figure 2. Furthermore, Q-Q plots of the data represented in the Manhattan plots are also shown in Figure 2. In SSC3, a region comprising 55.60 to 58.34 Mb was declared to be associated with the proportion of T helper cells (Figure 2A). Genes related to immune response were found annotated in this region: the family of immunoglobulin kappa variable region genes, T-cell surface glycoproteins (CD8A and CD8B), Required for Meiotic Nuclear Division 5 Homolog A (RMND5A), Charged Multivesicular Body FIGURE 1 Heatmap of genetic correlations estimated by pairwise combination Protein 3 (CHMP3), SET and MYND Domain Containing 1 among cell immunity traits. (SMYD1), Zeta Chain of T cell Receptor Associated Protein Kinase 70 (ZAP70) and Granulysin (GNLY). The three top significant TABLE 3 Description of the four chromosomal regions associated with cell immunity traits and the annotated candidate genes. Region Chr Start End N Top SNPs MAF p- FDR Trait Candidate genes (Mbp) (Mbp) SNPs value rs81340900, T helper cells CD8A, CD8B, CHMP3, GNLY, IGKV, 1 3 55.60 58.34 14 rs81293514, 0.49 2.18E-06 3.10E-02 (% of PBMCs) RMND5A, SMYD1, ZAP70 rs81336780 Memory T cells AICDA, CD4, CD69, CD163, MFAP5, 2 5 61.62 62.44 10 rs326461238 0.26 1.18E-05 8.70E-02 (% of PBMCs) PHC1, STYK1, TAS2Rs, KLRs, C-type CD4+ T cells lectins rs81406759, T helper cells 3 8 20.51 20.57 2 0.40 3.45E-05 9.80E-02 RBPJ, STIM2 rs81406807 (% of PBMCs) Naïve T cells 4 X 33.36 33.63 6 rs342772739 0.13 1.18E-06 1.20E-02 CYBB, ssc-mir-9786-1 (% of PBMCs) Frontiers in Immunology 06 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 A B C FIGURE 2 Manhattan plots and quantile-quantile plots representing the p-values profiles corresponding to the association analysis between immunity traits and SNPs, (A) for Helper T cells (B) for Memory T cells and (C) for Naïve T cells. Red line indicates those SNPs that are below the genome-wide significance threshold (FDR ≤ 0.1). associated SNPs (rs81340900, rs81293514 and rs81336780) were performed. Remarkably, we found several candidates genes previously located in an intron of the RMND5A gene. In addition, two reported as candidates in human GWAS studies for the same or significant SNPs at 20.51-20.57 Mb of SSC8 were also associated similar phenotypes such as CD8A, CD8B, RMND5A and VPS24 (also with the variation of T helper cells (Figure 2A and Table 3). In this known as CHMP3) for CD4+ proliferating, CD4+CD8dim T cell region, Recombination Signal Binding Protein for Immunoglobulin Absolute Counts, CD4 + CD8 dim T cell % lymphocytes and Kappa J region (RBPJ) and Stromal Interaction Molecule 2 (STIM2) CD4+CD8dim T cell % leukocyte traits (9, 14, 28) and CD4 and genes were annotated. KLRC2 for CD4+CD8+ and effector memory CD4+ T cell % T cell In SSC5, 10 significant SNPs at 61.97 to 62.29 Mb were associated traits (9, 28, 29). In pigs, two overlapping or very close genomic with the proportion of memory T-helper cells (Figure 2B), being regions for T helper cells and memory T-helper cells were identified. rs326461238 the most significant one (Table 2). In this region, a high Regarding T helper cells trait, overlapping regions for leukocyte number of immune-related genes (AICDA, CD4, CD69, CD163, percentages of CD8+, CD8- and CD3-CD8- traits of an F2 Duroc × MFAP5, PHC1, STYK1, a family of TAS2Rs, a family of KLR such Erhualian piglets vaccinated with a live CSF vaccine were identified as KLRB1, KLRC1, KLRK1 and KLRD1, and a family of C-type lectins (15). A region at 55 Mbp in SSC5 for CD4+CD8+ percentage of such as CLECL1, CLEC4D and CLEC4E) were identified. leukocytes and located less than 7 Mb away from the region identified Furthermore, the most associated SNP was found inside a long- in our analysis for CD4+CD8+ memory T-helper cells was identified noncoding RNA (lncRNA, LOC100524679). in a population of animals consisting of Landrace, Yorkshire, and Finally, six SNPs at 33.36 to 33.63 Mb in SSCX were associated Songlio Black pigs (16). with naïve T cells (Figure 2C). The most significant SNP (rs342772739) was found in the locus for LANCL3 gene. Despite that, no reports have been found to indicate that the gene is directly 3.5 Genes differentially expressed between associated with immunity phenotypes. In addition, CYBB gene, also groups diverging in percentage of T gd, T related to immunity functions, was found in this region. helper and T memory cells The transcriptome of whole blood samples from animals with 3.4 Comparison with other GWAS studies extreme proportions of T gd, T helper and memory T-helper cells was analysed by RNAseq to identify genes DE between groups. To identify overlaps between our QTL regions and those A total of 18 genes were identified as DE (FC>1.2; FDR<0.1) in previously identified in pigs and humans for immune cellular traits blood between the two groups of animals with high (H) and low (L) T a comparison between genomic regions and candidate genes was helper and T memory cell percentages. From them, 5 genes were Frontiers in Immunology 07 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 downregulated and 13 were upregulated in the H group when 3.6 Coincident association between GWAS compared to the L group (Figure 3A and Supplementary Table 6). and DE analyses The most up- and downregulated genes in the H group were DNA nucleotidylexotransferase (DNTT; FC = 17.65, P-value = 8.19 x 10-05) In order to identify potential candidate genes whose differences in and Cyclin dependent kinase 20 (CDK20; FC−1 = 5.67, P-value = 1.24 gene expression may be modulating our cellular phenotypic traits, we x 10-04), respectively. overlapped the lists of DE genes and candidate genes located in The functional annotation of DE genes identified CD4 as being associated genomic regions for T helper, memory T-helper and gd T involved in “helper T cell enhancement of adaptive immune cells traits. We identified four DE genes between groups of animals response” GO immune system process, among other immune- with extreme proportions of T helper and memory T-helper cells: related functions (Supplementary Table 7). Other genes related to CAPG and TCF7L1 on SSC3, and KLRD1 and CD4 on SSC5, at or very immune functions were KLRD1, NOTCH3 and TCF7L1 close to the genomic regions associated with T helper cells and (Supplementary Table 7). Furthermore, the following genes CAPG, memory T-helper cell, respectively. No genes colocalizing between DNTT, KCNIP3, and P2RY14 have also been reported in the literature GWAS and DE studies were identified for dg T cells. to be related to immunity (30–32). Regarding blood transcriptomic changes between H and L groups for percentages of gd T cells, 246 DE genes were identified (FC>1.2; 4 Discussion FDR<0.1; Figure 3B and Supplementary Table 8). Of these, 102 genes were upregulated and 144 were downregulated in the H group when Nowadays, the pig industry is facing new challenges that make compared to the L group. The most up- and downregulated genes in the necessary to reorient breeding programs to produce more robust and H group were the novel genes ENSSSCG00000046089 (FC = 35.25, P- resilient pig populations while maintaining production efficiency. It is value = 2.48 x 10-04) and ENSSSCG00000038561 (FC−1 = 43.43, P-value known that the immune system and particularly lymphocytes play an = 3.41 x 10-04), respectively. Out of the DE genes, 215 were identified as important role in controlling many porcine pathogens. In this study, protein coding genes, 24 as lncRNAs, and 7 as pseudogenes. we focused our analysis on peripheral blood mononuclear cell traits, After performing a functional analysis with the list of DE genes, mainly related to adaptive immunity, but also including cell we mainly identified GO immune system processes associated with a- populations such as dg T cells with functions associated with both b T cell differentiation and activation from the list of downregulated arms of immunity. genes in H group, gathering animals with high proportion of gd T After PBMC immunophenotyping, the relative abundance among cells, compared to L group. Other processes identified were negative PBMCs of total T and B lymphocytes in the analysed Duroc population regulation of lymphocyte activation and proliferation. Some genes was similar to figures reported in literature for other pig populations identified in these processes were CD274, HLX, IL27, LGALS9B, MYB, (33–35). However, differences with other studies were observed in the RORC, RSAD2 and SOCS1 (Supplementary Table 9). Remarkably, we proportion of some T-cells subsets, with CD4+ cells showing lower identified in the list of upregulated genes in H group a transcriptional relative abundance than that reported before (36). Although we could factor (SOX13) associated with positive regulation of gd T cell not discard an effect on immunophenotyping associated to handling activation and differentiation. Immune processes also identified and freezing processes, other factors such as the age of animals, from the list of up-regulated genes were regulation of B cell environmental factors or the lack of standardization of laboratory receptor signalling pathway (BLK and GCSAML), regulation of tests could also have generated discrepancies between studies. It is erythrocyte (ACVR2A, RHEX), and dendritic, monocyte and also plausible that differences in the genetic background of the porcine macrophage (ZBTB46) differentiation, neutrophil-mediated populations analysed in these studies could be behind the differences in functions (DPP4, CTSG and ADAM8) and T cell differentiation in the abundance of certain lymphocytes subsets. For instance, genetic thymus (ADAM8). variation of CD4 has been associated to a different reactivity to the CD4 A B FIGURE 3 Volcano plot displaying DE genes in blood between H and L groups for (A) T helper and memory T-helper cell and (B) dg T cells. The vertical axis (y-axis) corresponds to the -log10 (P-value), and the horizontal axis (x-axis) displays the log2 fold change (logFC) value. The vertical lines mark the thresholds located to FC = 1.5 and FC = -1.5. Blue dots represent downregulated genes in HvsL groups with FC< -1.5 and FDR< 0.1. Cyan dots represent downregulated genes in HvsL groups with -1.5< FC< -1.2 and FDR< 0.1. Red dots represent upregulated genes in HvsL groups with FC > 1.5 and FDR< 0.1. Frontiers in Immunology 08 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 antibody (36). In our study, CD4 was one of the DE genes between between the proportion of B lymphocytes and the phagocytic animals with extreme phenotypes for T helper and memory T-helper capacity, and between some T cells subpopulations (particularly cells, and was located in a genomic region associated with proportion of CTL and memory T-cells) and serum immunoglobulin CD4+ and memory T cells, evidencing a genetic contribution to the concentrations. As a whole, the genetic interaction map among variation of these traits in our Duroc population. Therefore, changes in innate and adaptive immunity traits obtained in our Duroc CD4 mRNA levels among animals may be the most plausible population partially support the possibility of applying an effective explanation for the differences in CD4 T cells percentages observed multi-trait selection for both innate and adaptive immunity. In this between studies. scenario of complex genetic relationships, the possibility of A particularly high heritability was estimated for both memory identifying genetic markers, biomarkers and candidate genes and CD4+ T-cells relative abundances, whereas medium to high associated to the phenotypic variation that could allow a more heritabilities were obtained for the rest of immune-cell-related traits directed selection of some of these traits, without impairing other in our Duroc population. These heritability values are in line with elements of the immune systems, takes particular relevance. those obtained in previous human and pig studies, reaffirming A genetic correlation close to unity was observed between naïve genetics as an important determinant of adaptive cell immune traits and dg T cells quantified with a specific monoclonal antibody (8, 10, 28, 37). The lowest heritability among PBMCs subpopulations (MAC320) in a previous study in the same Duroc population (12), relative abundances was observed for B cells. This is concordant with showing that most of the CD4-CD8- T cells were in fact gd T results in humans (29), showing that B cell counts were more lymphocytes. Consequently, GWAS study identified the same environmentally influenced while T cells were more strongly driven genomic region on SSCX using both phenotypic datasets (12). In by genetic factors. the present study, we annotated a new promising candidate gene, The map of genetic correlations among PBMC subtypes relative CYBB encoding for one of the two chains of cytochrome b (CYB). abundances reported positive and negative genetic associations CYB is a component of the microbicidal oxidase system expressed by between these immune-cell-related traits, but also reflected the neutrophils and T cells, a critical process of innate immunity for interdependence between some of these compositional variables and defence against bacterial and fungal infections. Mutations in human should be interpreted with caution. The antagonism between T and B CYBB have been previously associated to immune system dysfunction cells proportions in blood observed in our Duroc population might and can lead to chronic granulomatous disease (39). Since dg T cells have some physiological implications deserving further research. constitute a high proportion of lymphocytes in porcine peripheral However, we cannot discard this antagonism to be the consequence blood and could have the potential to combine conventional adaptive of dealing with compositional data. Conversely, we can affirm that our and innate-like responses (35, 40), the activation of WC1-bearing dg results confirm a positive genetic association between memory and T- T cells has emerged as a potential alternative to optimize new vaccine helper cells in blood composition. In this regard, two significant strategies (41). In this sense, an increased recruitment of dg T cells genomic regions associated with the proportion of T helper and when animals are vaccinated earlier in life appears to correlate with memory T-helper cells in blood were identified. This may allow improved efficacy of current vaccines (reviewed in (42)). Indeed, selection of animals to investigate abundance of these cells in blood in TRDC-knockout pigs, defective in dg T cells, showed lower relation to immune responses toward different stimuli. neutralizing antibody titres compared to syngenic controls when Selecting animals for producing higher immune responses should they were vaccinated against classical swine fever virus using a probably consider both innate and adaptive immunity traits in highly attenuated live virus (43). Thus, our results indicate the conjunction with qualitative responses. The map of genetic possibility to select animals for higher dg T cells counts in blood to correlations of PBMCs subpopulation proportions with innate check for potential enhancement of vaccine efficacy. immunity traits might help to evaluate the consequences/ Another promising candidate genes were located in the SSC3, SSC5 convenience of selecting for both innate and adaptive immunity and SSC8 genomic regions associated with the variation of T helper and traits. Most PBMCs subset abundances (but naïve T cells) showed a memory T-helper cells and proportion of CD4+ cells. CD8A, CD8B, GNLY negligible genetic association with the concentration of acute phase and RBPJ were some of the genes associated with the proportion of T helper proteins (CRP and Hp) in blood. Similar results were shown by (10), cells. CD8a and CD8b are integral membrane glycoproteins found on the whereas a positive genetic correlation between the concentration of surface of many immune-cell types and act as co-receptors during antigen CRP and B cell counts was described by (38) in animals tested at 90kg recognition by TCR (44). GNLY is an antimicrobial secretion protein of live weight. The genetic independence of acute phase proteins present in CTLs and NK cells (45). It also functions as an alarmin, concentration versus most lymphocytes subsets abundance observed activating dendritic cells and promoting antigen-specific immune in our Duroc population would support the possibility of selecting responses by using TLR4 (46). Finally, RBPJ encodes a transcription some immune-cell-related traits independently of some innate factor acting through the NOTCH signalling pathway to protect T- immune pathways governing CRP synthesis and then avoiding helper cells from apoptosis (47). Also associated with the percentage of adverse effects due to inflammation. memory T-helper and CD4+ cells we identified CD4 and CD69. While Genetic correlations of the different PBMC subpopulations CD4 is the main membrane marker for helper (naïve and memory) T-cells, proportions with other immunity traits also pinpoint to other it serves as a co-receptor of the TCR, enhancing recognition of the MHC interesting results. The relative abundance of CTL T cells was class II complex (48). Of interest CD69 expressed in resident memory T positively associated at the genetic level with haematological traits, cells and gd T cells, is an activation marker determining migration- particularly with leucocyte counts and with the haematocrit. retention of T cells from tissues (49). Furthermore, we also identified Conversely, our results also suggest certain genetic antagonism CLECL1 and the killer cell lectin-like receptors genes KLRC1 and KLRD1 Frontiers in Immunology 09 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 which are highly expressed in memory T cells and are involved in their molecular mechanisms underlying lymphocyte traits with special activation (CLECL1) and killing function (KLRs expressed by human interest to WC1+ gd T cells to explore further optimization of CD8+ T effector memory cells) (50). Indeed, a particularly remarkable current and new vaccines. result arising in this study was the overlap of genes and genomic regions for the same or similar lymphocyte-related phenotypes with previous studies reported in pig and humans (9, 14–16, 28, 29). To the best of our Data availability statement knowledge, this is the first study to report overlapping between humans and pig genes (CD4, CD8A, CD8B, KLRC2, RMND5A and VPS24) The original contributions presented in the study are included in associated with the variation of T-cell populations, further supporting the article/Supplementary Material. Further inquiries can be directed the use of pig as a very reliable biomodel for human infectious diseases, to the corresponding authors. vaccine development and T-cell research (51–53). Also, another remarkable result obtained in this study was the overlap between GWAS and DE analyses, pointing out genes (CAPG, TCF7L1, KLRD1 and CD4) Ethics statement that may be modulating T lymphocytes traits through variation of its gene expression. Since transcriptional and post-transcriptional regulation of All experimental procedures with pigs were performed according immune gene expression is critical for modulating immune responses (54, to the Spanish Policy for Animal Protection RD 53/2013, which meets 55), eGWAS studies are required to identify potential cis and/or trans the European Union Directive 2010/63/EU about the protection of regulatory factors regulating the expression of these genes. However, we animals used in experimentation. The experimental protocol was cannot rule out that the differences in the expression of some of these genes approved by the Ethical Committee of the Institut de Recerca i may be related to the different proportions of these populations in blood. Tecnologia Agroalimentàries (IRTA). For example, CAPG is a gene mainly expressed in monocytes, B cells and dendritic cells (Human Protein Atlas proteinatlas.org (56),), and was downregulated in the group with higher proportions of T helper and T Author contributions memory cells. Furthermore, no overlapped genes were identified between GWAS and DE analyses for dg T cells. However, we identified an MB designed the study. MB supervised the generation of the interesting pattern of gene expression with genes related to immune material animal used in this work. MB, OG-R, RQ and YR-C processes of ab T cell differentiation and activation downregulated in performed the sampling. AP, MB, SL-S and TJ-J carried out the the group with higher proportion of gd T cells. This result is in agreement laboratory analyses. AP, CH-B, DC-P, RQ, SL-S, TJ-J, YR-C and MB with the upregulation of SOX13 in animals with higher proportion of gd T analysed the data and interpreted the results. MB and RQ wrote the cells as this gene encodes a transcription factor essential for normal gd T cell manuscript. All authors contributed to the article and approved the development (57). It has been reported a higher expression of SOX13 in submitted version. bovine WC1+ gd T cells relative to WC1- gd T cells (58). Remarkably, we also identified three upregulated genes ENSSSCG00000034914 (CD163L1), ENSSSCG00000039217 (cow WC1 orthologue) and Funding ENSSSCG00000031085 (antigen WC1.1) in animals with higher proportions of gd T cells, consistent with the transcription of WC1 co- The study was funded by grants AGL2016-75432-R and receptors by circulating gd T cells (42). WC1+ cells can be divided in PID2020-112677RB-C21 awarded by the Spanish Ministry of different subpopulations with different ability to respond to particular Economy and Competitiveness (MINECO). D. Crespo-Piazuelo was pathogens and cytokine responses (42). A model in cows has been supported by the GENE-SWitCH project (https://www.gene-switch. proposed in which a differential occupancy on these WC1 gene loci by eu), which is funded by the European Union’s Horizon 2020 research SOX13 influence the development of these subpopulations (58). However, and innovation programme under the grant agreement n° 817998. T. still there is a lack of information about these subpopulations in swine (59). Jové -Juncà was supported by an IRTA fellowship (CPI1221) and C. Thus, further analysis is required to determine if pigs used in the present Herná ndez-Banqué was funded with a FPI grant (PRE2021-097825) report have different proportions of WC1+ subpopulations. Besides, it from the PID2020-112677RB-C21 project. YR-C was financially should be noted that the proportion of some lymphocyte subsets may vary supported by a Ramon y Cajal contract (RYC2019-027244-I) from with the age of animals (33, 34, 36), so it would be worthy to study the the Spanish Ministry of Science and Innovation. The authors evolution of these lymphocytes traits throughout the life of the animals. belonged to a Consolidated Research Group AGAUR, ref. In conclusion this study confirmed the genetic determinism of 2017SGR-1719. lymphocytes traits in pigs and provide a list of promising genes and biomarkers associated to the variation of these traits. The genetic correlations described for innate and adaptive traits and the Acknowledgments molecular markers identified for T helper, memory T-helper and gd T cells support the possibility of applying multi-trait selection for We gratefully acknowledge the technical staff from IRTA and immunocompetence in pigs. Finally, the overlap between genes for Selección Batallé S.A., and specially to Josep Reixach, for their the same T cell traits in humans and pigs reinforces the use of pig as a collaboration in the experimental protocols at farm and biomodel for vaccine development and the study of human infectious slaughterhouse, as well as to Albert Bensaid for valuable discussions diseases. Overall, our results provide new knowledge about the and comments on the study. Frontiers in Immunology 10 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 Conflict of interest organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or The authors declare that the research was conducted in the claim that may be made by its manufacturer, is not guaranteed or absence of any commercial or financial relationships that could be endorsed by the publisher. construed as a potential conflict of interest. Supplementary material Publisher’s note The Supplementary Material for this article can be found online at: All claims expressed in this article are solely those of the authors https://www.frontiersin.org/articles/10.3389/fimmu.2023.1058346/ and do not necessarily represent those of their affiliated full#supplementary-material References 1. Millet S, Maertens L. The European ban on antibiotic growth promoters in animal 20. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and feed: from challenges to opportunities. Vet J (2011) 187:143–4. doi: 10.1016/ powerful approach to multiple testing. J R Stat Soc Ser B (1995) 57:289–300. doi: 10.2307/ J.TVJL.2010.05.001 2346101 2. Bai X, Plastow GS. Breeding for disease resilience: opportunities to manage 21. Turner SD. Qqman: an r package for visualizing GWAS results using q-q and polymicrobial challenge and improve commercial performance in the pig industry. manhattan plots. J Open Source Softw (2018) 3:731. doi: 10.21105/JOSS.00731 CABI Agric Biosci (2022) 3:6. doi: 10.1186/S43170-022-00073-Y 22. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The 3. Mellencamp MA, Galina-Pantoja L, Gladney CD, Torremorell M. Improving pig ensembl variant effect predictor. Genome Biol (2016) 17:122. doi: 10.1186/s13059-016- health through genomics: A view from the industry. Dev Biol (Basel) (2008) 132:35–41. 0974-4 doi: 10.1159/000317142 23. Smedley D, Haider S, Durinck S, Pandini L, Provero P, Allen J, et al. The BioMart 4. Visscher AH, Janss LL, Niewold TA, de Greef KH. Disease incidence and community portal: an innovative alternative to large, centralized data repositories. Nucleic immunological traits for the selection of healthy pigs. A review Vet Q (2002) 24:29–34. Acids Res (2015) 43:W589–98. doi: 10.1093/nar/gkv350 doi: 10.1080/01652176.2002.9695121 24. Safran M, Solomon I, Shmueli O, Lapidot M, Shen-Orr S, Adato A, et al. 5. Mair KH, Sedlak C, Käser T, Pasternak A, Levast B, Gerner W, et al. The porcine GeneCardsTM 2002: towards a complete, object-oriented, human gene compendium. innate immune system: An update. Dev Comp Immunol (2014) 45:321–43. doi: 10.1016/ Bioinformatics (2002) 18:1542–3. doi: 10.1093/bioinformatics/18.11.1542 j.dci.2014.03.022 25. Eppig JT, Smith CL, Blake JA, Ringwald M, Kadin JA, Richardson JE, et al. Mouse 6. Bonilla FA, Oettgen HC. Adaptive immunity. J Allergy Clin Immunol (2010) 125: genome informatics (MGI): Resources for mining mouse genetic, genomic, and biological S33–40. doi: 10.1016/J.JACI.2009.09.017 data in support of primary and translational research. Methods Mol Biol (Clifton N.J.) 7. Liston A, Carr EJ, Linterman MA. Shaping variation in the human immune system. (2017), 47–73. doi: 10.1007/978-1-4939-6427-7_3 Trends Immunol (2016) 37:637–46. doi: 10.1016/J.IT.2016.08.002 26. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for 8. Mangino M, Roederer M, Beddall MH, Nestle FO, Spector TD. Innate and adaptive differential expression analysis of digital gene expression data. Bioinformatics (2010) immune traits are differentially affected by genetic and environmental factors. Nat 26:139–40. doi: 10.1093/bioinformatics/btp616 Commun (2017) 8:13850. doi: 10.1038/NCOMMS13850 27. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. 9. Orrù V, Steri M, Sole G, Sidore C, Virdis F, Dei M, et al. Genetic variants regulating ClueGO: a cytoscape plug-in to decipher functionally grouped gene ontology and pathway immune cell levels in health and disease. Cell (2013) 155:242. doi: 10.1016/ annotation networks. Bioinformatics (2009) 25:1091–3. doi: 10.1093/bioinformatics/ J.CELL.2013.08.041 btp101 10. Flori L, Gao Y, Laloë D, Lemonnier G, Leplat J-J, Teillaud A, et al. Immunity traits 28. Orrù V, Steri M, Sidore C, Marongiu M, Serra V, Olla S, et al. Complex genetic in pigs: substantial genetic variation and limited covariation. PloS One (2011) 6:e22717. signatures in immune cells underlie autoimmunity and inform therapy. Nat Genet (2020) doi: 10.1371/journal.pone.0022717 52:1036–45. doi: 10.1038/s41588-020-0684-4 11. Clapperton M, Glass EJ, Bishop SC. Pig peripheral blood mononuclear leucocyte 29. Aguirre-Gamboa R, Joosten I, Urbano PCM, van der Molen RG, van Rijssen E, van subsets are heritable and genetically correlated with performance. Animal (2008) 2:1575– Cranenbroek B, et al. Differential effects of environmental and genetic factors on T and b 84. doi: 10.1017/S1751731108002929 cell immune traits. Cell Rep (2016) 17:2474–87. doi: 10.1016/j.celrep.2016.10.053 12. Ballester M, Ramayo-Caldas Y, Gonzá lez-Rodrı́guez O, Pascual M, Reixach J, Dı́az 30. Aono A, Enomoto H, Yoshida N, Yoshizaki K, Kishimoto T, Komori T. Forced M, et al. Genetic parameters and associated genomic regions for global expression of terminal deoxynucleotidyl transferase in fetal thymus resulted in a decrease immunocompetence and other health-related traits in pigs. Sci Rep (2020) 10:18462. in gammadelta T cells and random dissemination of Vgamma3Vdelta1 T cells in skin of doi: 10.1038/s41598-020-75417-7 newborn but not adult mice. Immunology (2000) 99:489–97. doi: 10.1046/J.1365- 2567.2000.00987.X 13. Knap PW, Bishop SC. Relationships between genetic change and infectious disease in domestic livestock. Occ Publ Br Soc Anim Sci (2000) 27:65–80. doi: 10.1017/ 31. Li Q, Xu L, Li Y, Yang R, Qiao Q, Wang Y, et al. P2RY14 is a potential biomarker S1463981500040553 of tumor microenvironment immunomodulation and favorable prognosis in patients with head and neck cancer. Front Genet (2021) 12:670746. doi: 10.3389/ 14. Lagou V, Garcia-Perez JE, Smets I, Van Horebeek L, Vandebergh M, Chen L, et al. FGENE.2021.670746 Genetic architecture of adaptive immune system identifies key immune regulators. Cell Rep (2018) 25:798–810.e6. doi: 10.1016/J.CELREP.2018.09.048 32. Parikh SS, Litherland SA, Clare-Salzler MJ, Li W, Gulig PA, Southwick FS. CapG (-/-) mice have specific host defense defects that render them more susceptible than CapG 15. Zhang J, Chen JH, Liu XD, Wang HY, Liu XL, Li XY, et al. Genomewide association (+/+) mice to listeria monocytogenes infection but not to salmonella enterica serovar studies for hematological traits and T lymphocyte subpopulations in a duroc × erhualian F2 typhimurium infection. Infect Immun (2003) 71:6582–90. doi: 10.1128/IAI.71.11.6582- resource population. J Anim Sci (2016) 94:5028–41. doi: 10.2527/jas.2016-0924 6590.2003 16. Lu X, Fu WX, Luo YR, Ding XD, Zhou JP, Liu Y, et al. Genome-wide association 33. Yang H, Parkhouse RME. Phenotypic classification of porcine lymphocyte study for T lymphocyte subpopulations in swine. BMC Genomics (2012) 13:488. subpopulations in blood and lymphoid tissues. Immunology (1996) 89:76–83. doi: 10.1186/1471-2164-13-488 doi: 10.1046/j.1365-2567.1996.d01-705.x 17. Misztal I, Tsuruta S, Strabel T, Auvray B, Druet T, Lee D. 7th world congress on 34. Pietrasina O, Miller J, Rzas̨ a A. Differences in the relative counts of peripheral genetics applied to livestock production, August 19-23, 2002, Montpellier, France. 7th blood lymphocyte subsets in various age groups of pigs(2020) (Accessed September 14, World Congress Genet Appl to Livestock Production (2002). 2022). 18. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. 35. Takamatsu HH, Denyer MS, Stirling C, Cox S, Aggarwal N, Dash P, et al. Porcine PLINK: a tool set for whole-genome association and population-based linkage analyses. gammadelta T cells: possible roles on the innate and adaptive immune responses Am J Hum Genet (2007) 81:559–75. doi: 10.1086/519795 following virus infection. Vet Immunol Immunopathol (2006) 112:49–61. doi: 10.1016/ 19. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: A tool for genome-wide J.VETIMM.2006.03.011 complex trait analysis. Am J Hum Genet (2011) 88(1):76–82. doi: 10.1016/ 36. Blanc F, Pré vost-Blondel A, Piton G, Bouguyon E, Leplat JJ, André oletti F, et al. j.ajhg.2010.11.011 The composition of circulating leukocytes varies with age and melanoma onset in the Frontiers in Immunology 11 frontiersin.org Ballester et al. 10.3389/fimmu.2023.1058346 MeLiM pig biomedical model. Front Immunol (2020) 11:291. doi: 10.3389/ cellular longevity. Proc Natl Acad Sci U.S.A. (2012) 109:9041–6. doi: 10.1073/PNAS. FIMMU.2020.00291 1206044109 37. Clapperton M, Glass EJ, Bishop SC. Pig peripheral blood mononuclear leucocyte 48. Gascoigne NRJ, Zal T, Yachi PP, Hoerter JAH. Co-Receptors and recognition of subsets are heritable and genetically correlated with performance. Animal (2008) 2:1575– self at the immunological synapse. Curr Top Microbiol Immunol (2010) 340:171–89. 84. doi: 10.1017/S1751731108002929 doi: 10.1007/978-3-642-03858-7_9 38. Clapperton M, Diack AB, Matika O, Glass EJ, Gladney CD, Mellencamp MA, et al. 49. Cibriá n D, Sá nchez-Madrid F. CD69: from activation marker to metabolic Traits associated with innate and adaptive immunity in pigs: heritability and associations gatekeeper. Eur J Immunol (2017) 47:946–53. doi: 10.1002/eji.201646837 with performance under different health status conditions. Genet Sel Evol (2009) 41:54. 50. Weng NP, Araki Y, Subedi K. The molecular basis of the memory T cell response: doi: 10.1186/1297-9686-41-54 differential gene expression and its epigenetic regulation. Nat Rev Immunol (2012) 39. Rae J, Newburger PE, Dinauer MC, Noack D, Hopkins PJ, Kuruto R, et al. X- 12:306–15. doi: 10.1038/NRI3173 Linked chronic granulomatous disease: Mutations in the CYBB gene encoding the gp91- 51. Käser T. Swine as biomedical animal model for T-cell research-success and phox component of respiratory-burst oxidase. Am J Hum Genet (1998) 62:1320–31. potential for transmittable and non-transmittable human diseases. Mol Immunol doi: 10.1086/301874 (2021) 135:95–115. doi: 10.1016/J.MOLIMM.2021.04.004 40. Vantourout P, Hayday A. Six-of-the-best: unique contributions of gd T cells to immunology. Nat Rev Immunol (2013) 13:88–100. doi: 10.1038/NRI3384 52. Meurens F, Summerfield A, Nauwynck H, Saif L, Gerdts V. The pig: a model for human infectious diseases. Trends Microbiol (2011) 20:50–7. doi: 10.1016/ 41. Teixeira AF, Gillespie A, Yirsaw A, Britton E, Telfer JC, Nascimento ALTO, et al. j.tim.2011.11.002 Identification of leptospiral protein antigens recognized by WC1 + gd T cell subsets as target for development of recombinant vaccines. Infect Immun (2022) 90(1):e0049221. 53. Pabst R. The pig as a model for immunology research. Cell Tissue Res (2020) doi: 10.1128/IAI.00492-21 380:287–304. doi: 10.1007/S00441-020-03206-9 42. Le Page L, Baldwin CL, Telfer JC. gd T cells in artiodactyls: Focus on swine. Dev 54. Smale ST. Transcriptional regulation in the immune system: A status report. Comp Immunol (2022) 128:104334. doi: 10.1016/J.DCI.2021.104334 Trends Immunol (2014) 35:190. doi: 10.1016/J.IT.2014.03.003 43. Petersen B, Kammerer R, Frenzel A, Hassel P, Dau TH, Becker R, et al. Generation 55. Schwerk J, Savan R. Translating the untranslated region. J Immunol (2015) and first characterization of TRDC-knockout pigs lacking gd T cells. Sci Rep (2021) 195:2963. doi: 10.4049/JIMMUNOL.1500756 11:14965. doi: 10.1038/S41598-021-94017-7 56. Karlsson M, Zhang C, Mé ar L, Zhong W, Digre A, Katona B, et al. A single-cell 44. Cole DK, Laugel B, Clement M, Price DA, Wooldridge L, Sewell AK. The type transcriptomics map of human tissues. Sci Adv (2021) 7: eabh2169. doi: 10.1126/ molecular determinants of CD8 co-receptor function. Immunology (2012) 137:139–48. SCIADV.ABH2169 doi: 10.1111/J.1365-2567.2012.03625.X 57. Melichar HJ, Narayan K, Der SO, Hiraoka Y, Gardiol N, Jeannet G, et al. 45. Crespo  C, Mulik S, Dotiwala F, Ansara JA, Sen Santara S, Ingersoll K, et al. Regulation of gammadelta versus alphabeta T lymphocyte differentiation by the Decidual NK cells transfer granulysin to selectively kill bacteria in trophoblasts. Cell transcription factor SOX13. Science (2007) 315:230–3. doi: 10.1126/SCIENCE.1135344 (2020) 182:1125–1139.e18. doi: 10.1016/J.CELL.2020.07.019 58. Damani-Yokota P, Zhang F, Gillespie A, Park H, Burnside A, Telfer JC, et al. 46. Tewary P, Yang D, de la Rosa G, Li Y, Finn MW, Krensky AM, et al. Granulysin Transcriptional programming and gene regulation in WC1 + gd T cell subpopulations. activates antigen-presenting cells through TLR4 and acts as an immune alarmin. Blood Mol Immunol (2022) 142:50–62. doi: 10.1016/J.MOLIMM.2021.12.016 (2010) 116:3465–74. doi: 10.1182/BLOOD-2010-03-273953 59. Le Page L, Gillespie A, Schwartz JC, Prawits LM, Schlerka A, Farrell CP, et al. 47. Helbig C, Gentek R, Backer RA, De Souza Y, Derks IAM, Eldering E, Subpopulations of swine gd T cells defined by TCRg and WC1 gene expression. Dev Comp et al. Notch controls the magnitude of T helper cell responses by promoting Immunol (2021) 125:104214. doi: 10.1016/J.DCI.2021.104214 Frontiers in Immunology 12 frontiersin.org