Genetic architecture of innate and adaptive immune cells in pigs

Pig industry is facing new challenges that make necessary to reorient breeding programs to produce more robust and resilient pig populations. The aim of the present work was to study the genetic determinism of lymphocyte subpopulations in the peripheral blood of pigs and identify genomic regions and biomarkers associated to them. For this purpose, we stained peripheral blood mononuclear cells to measure ten immune-cell-related traits including the relative abundance of different populations of lymphocytes, the proportions of CD4⁺ T cells and CD8⁺ T cells, and the ratio of CD4⁺/CD8⁺ T cells from 391 healthy Duroc piglets aged 8 weeks. Medium to high heritabilities were observed for the ten immune-cell-related traits and significant genetic correlations were obtained between the proportion of some lymphocytes populations. A genome-wide association study pointed out 32 SNPs located at four chromosomal regions on pig chromosomes SSC3, SSC5, SSC8, and SSCX as significantly associated to T-helper cells, memory T-helper cells and γδ T cells. Several genes previously identified in human association studies for the same or related traits were located in the associated regions, and were proposed as candidate genes to explain the variation of T cell populations such as CD4, CD8A, CD8B, KLRC2, RMND5A and VPS24. The transcriptome analysis of whole blood samples from animals with extreme proportions of γδ T, T-helper and memory T-helper cells identified differentially expressed genes (CAPG, TCF7L1, KLRD1 and CD4) located into the associated regions. In addition, differentially expressed genes specific of different T cells subpopulations were identified such as SOX13 and WC1 genes for γδ T cells. Our results enhance the knowledge about the genetic control of lymphocyte traits that could be considered to optimize the induction of immune responses to vaccines 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.

1 Introduction

Pork is one of the most consumed meats worldwide. The industrialization and specialization of the swine production coupled with the ever-increasing movement of pigs and pork products at national and international scales facilitate the appearance and spread 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 bacteria, and the threats of climate change, stresses the need for guiding breeding programs to produce more robust and disease resistant pig populations (1). Along with conventional methods such as vaccination and strict biosecurity measures, improving the overall immunocompetence of healthy animals is a good approach to prevent and control infectious pathogens (2, 3).

In this context, immunity traits (ITs) are considered biologically relevant parameters to measure immunocompetence (4). ITs may be divided in the two major components of the immune system, the 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 pathogens through a restricted set of pattern recognition receptors. Cell types involved in innate immunity include Natural Killer (NK) cells, dendritic cells, monocytes or γδ T cells (5). In contrast, specialized cells of the acquired immune system can provide immunological memory after an initial response mediated by a very large number of specific receptors for antigens present in pathogens (6). While T lymphocytes are the effectors of adaptive cellular immune responses, antibody-producing cells, B lymphocytes, mediate adaptive humoral immunity (6).

The immune system has a high degree of interindividual variation that is controlled by genetic and environmental factors (7). However, leukocyte subsets, mostly the adaptive ones, are highly affected by genetics (8–10). Genetic studies have estimated moderate to high heritabilities for both innate (e.g. NK cells and γδ T cells) and adaptive (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). Furthermore, the development of high-density genotyping arrays and in-depth immune phenotyping techniques have allowed the identification of candidate genes and genetic variants associated with the phenotypic variation of immune cells in humans (9, 14) and pigs (15, 16). However, genome-wide association studies (GWAS) for innate and adaptive immune cells in pigs are still scarce. The few studies on this topic have been focused on identifying genetic variants that affect T-cell subpopulations in animals vaccinated at 21 days of age with live classical swine fever (CSF) vaccine (15, 16).

The objective of this work was to study the genetic determinism of both innate and adaptive immune cells traits in healthy pigs of a Duroc commercial line by estimating their genetic parameters and identifying genomic regions, candidate genes and genetic markers associated to their phenotypic variation. Finally, knowing the relevance of the pig as a human biomodel, a comparison between both species was performed with the significantly associated regions found in pigs.

2 Materials and methods

2.1 Ethics approval

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).

2.2 Animal material

A total of 391 healthy piglets aged 60 ± 8 days (males and females) belonging to a commercial Duroc line were used in this study. All animals were vaccinated at weaning (aged ~26 days old) against PCV2 with Porcilis^® PCV. The pigs came from six batches (72 ± 1 animals per batch) and were raised in the same farm and ad libitum fed with a commercial cereal-based diet. All animals were apparently healthy, without any observable signs of infection.

Blood was collected via the external jugular vein into vacutainer tubes with anticoagulants (Sangüesa S.A., Spain), which required the restraint of the animals but not their sedation. Genomic DNA was extracted from blood samples using the NucleoSpin^® Blood (Macherey-Nagel, Germany). DNA concentration and purity were measured in a Nanodrop ND-1000 spectrophotometer.

2.3 Immunophenotyping

Peripheral blood mononuclear cells (PBMCs) isolated from heparinized peripheral blood by density-gradient centrifugation with Histopaque-1077 (Sigma, Spain) and previously stored at -80°C from a study for γδ T-lymphocyte quantification (12), were herein employed for immune phenotyping. Frozen PBMCs were thawed and resuspended in RPMI 1640 medium supplemented with 5% Foetal Bovine Serum (FBS) (Sigma, Spain), 1% Penicillin-Streptomycin (10,000 U/ml - 10 mg/mL) and 1% L-Glutamine (200 mM) (Cultek, Spain) after centrifugation at 450 g at 4°C for 10 min.

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, CA, USA) diluted in PBS for 25 min at 4°C following the manufacturer’s indications. After two washes, PBMCs were labelled with the primary-conjugated antibody mixture in PBS-1%FBS for 25 min at 4°C (Supplementary Table 1).

After two washes with 1xPBS-1% FBS at 4°C, cells were resuspended in 1xPBS-1% FBS and analysed by flow cytometry using the MACSQuant Analyzer 10 Flow cytometer (Miltenyi Biotec GmbH, Bergisch Gladbach, Germany) and the Flowlogic™ software v7.3 (Inivai Technologies, Melbourne, Australia). Beside test samples, unstained cells, stained samples with sole viability marker, isotypes for each antibody subclass, and Fluorescence Minus One (FMO) stained samples were included as controls to adjust the analysis and discard false positive results.

The gating strategy per each animal used for this calculation is shown in Figure S1. Absolute number of events were used to calculate percentages of lymphocyte populations depending on PBMCs gate. After discarding the dead cells gated by the use of the viability staining, the following cell subsets were quantified: B lymphocytes (CD21⁺), T lymphocytes (CD3⁺), natural killer (NK) cells (CD3^-CD21^-CD8⁺), cytotoxic (CTL) T cells (CD3⁺CD4^-CD8⁺), T helper cells (CD3⁺CD4⁺CD8^-), memory T-helper cells (CD3⁺CD4⁺CD8⁺) and naïve T cells (CD3⁺CD4^-CD8^-). In addition, we quantified the total proportion of CD3⁺CD4⁺, CD3⁺CD8⁺ and the ratio of CD4⁺ to CD8⁺ (CD4⁺/CD8⁺).

2.4 Statistics of PBMCs subpopulations

Basic descriptive statistics of absolute counts and proportion of the different analysed lymphocyte subsets of cells among PBMCs are shown in Table 1. Normality of each PBMC subpopulation data was checked through Shapiro-Wilk test in R. For some of the variables, logarithm or square root transformations were applied to reach normal distribution of residuals (p-value > 0.05).

TABLE 1

Table 1 Descriptive statistics of absolute numbers and proportions of peripheral blood mononuclear cells (PBMCs).

A linear regression model was then applied on transformed and 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 ratio tests. When significant, sex and/or batch effects were considered for subsequent analyses.

Pairwise phenotypic correlations (rp) among all analysed phenotypes were computed after adjusting for significant systematic factors.

2.5 Estimation of genetic parameters

The heritability (h²) , i.e. the proportion of phenotypic variance attributable to additive genetic effects, was estimated for the percentage of the different lymphocyte subsets of cells showed in Table 1. Variance components and the corresponding h² were estimated from an univariate animal model as follows:

Y=Xβ+Zu+e

where Y is the vector of phenotypes, the percentage of lymphocyte subpopulations (either transformed or raw data) of all individuals; β is the vector of systematic (fixed) effects on the trait, including when significant the effects of sex (2 levels) and/or batch (6 levels); u is the vector of animal’s genetic additive (random) effects; X and Z are the corresponding incidence matrices for β and u ; and e is the vector of random residual terms. The assumed distribution of additive genetic effects was u∼(0, A${\sigma }_u^2$), where A is the numerator relationship matrix computed on the basis of pedigree and ${\sigma }_u^2$ is the additive genetic variance; random errors were distributed as e∼N(0,I${\sigma }_e^2$). Estimation of the model variance components and the corresponding heritability $\left(h^2={\sigma }_u^2/\left({\sigma }_u^2+{\sigma }_e^2\right)\right)$ for each lymphocytes subset was performed by REML using the aireml program included in the BGF90 package (17). The standard errors (SE) of the heritability estimates were also computed, thus obtaining the corresponding confidence intervals at 95% (CI95).

Subsequently, pairwise genetic correlations (for each two traits combination) were estimated in a two-traits animal model described as follows:

$$\left[\begin{array}{c}{\mathit{Y}}_{t1}\\ {\mathit{Y}}_{t2}\end{array}\right]=\left[\begin{array}{cc}{\mathit{X}}_{t1}& 0\\ 0& {\mathit{X}}_{t2}\end{array}\right]\left[\begin{array}{c}{\mathit{\beta }}_{t1}\\ {\mathit{\beta }}_{t2}\end{array}\right]+\left[\begin{array}{cc}{\mathit{Z}}_{t1}& 0\\ 0& {\mathit{Z}}_{t2}\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_{t1}\\ {\mathit{u}}_{t2}\end{array}\right]+\left[\begin{array}{c}{\mathit{e}}_{t1}\\ {\mathit{e}}_{t2}\end{array}\right]$$

where Y_t1  and Y_t2 are the vectors of phenotypic observations for trait 1 and trait 2, respectively; β_t1 and β_t2 are the vectors of systematic (fixed) effects on each trait previously described, and X_t1 and X_t2 the correspondent incidence matrices; u_t1 and u_t2 are the vectors of animal genetic additive effects on trait 1 or trait 2 (random effects), and Z_t1 and Z_t2 the corresponding incidence matrices; finally e_t1 and e_t2 are the vectors of residual errors for each trait. The (co)variance matrix of random genetic effects was defined as:

$$Var\left[\begin{array}{c}{\mathit{u}}_{t1}\\ {\mathit{u}}_{t2}\end{array}\right]=\left[\begin{array}{cc}\mathit{A}{\sigma }_{u1}^2& \mathit{A}{\sigma }_{u1,u2}\\ \mathit{A}{\sigma }_{u1,u2}& \mathit{A}{\sigma }_{u1}^2\end{array}\right]$$

where ${\sigma }_{u1}^2$ and ${\sigma }_{u2}^2$ are the additive genetic variance of traits 1 and 2, respectively, σ_u1,u2 is the genetic covariance between the traits, and A is the numerator relationship matrix as defined above. Estimation of the (co)variance components of each pairwise analysis was also performed by REML using the aireml program included in the BGF90 package (17). Genetic correlations (rg) between traits were obtained as rg= $\left(\raisebox{1ex}{${\sigma }_{u1,u2}$}\!\left/ \!\raisebox{-1ex}{${\sigma }_{u1}{\sigma }_{u2}$}\right.\right)$, and the SE of the genetic correlation estimates were also computed.

2.6 SNP genotyping

The GGP Porcine HD Array (Illumina, San Diego, CA) was employed to genotype 68,516 single nucleotide polymorphisms (SNPs) in the 391 animals using the Infinium HD Assay Ultra protocol (Illumina). Plink software (18) was used to remove those SNPs with a minor allele frequency (MAF) lesser than 5%, SNPs with more than 10% missing genotypes, and SNPs that did not map to the porcine reference genome (Sscrofa11.1 assembly). Finally, a subset of 42,641 SNPs remained for further analysis.

2.7 Genome wide association study (GWAS)

GWAS was carried out between the 42,641 filtered SNPs and the 10 traits related to the percentages of PBMCs subpopulations, described in Table 1. To this end, the genome-wide complex trait analysis (GCTA) software (19) was used following an additive model for each trait across all SNPs:

$$\mathit{Y}=\mathit{X}\mathit{\beta }+\mathit{Z}\mathit{g}+{\mathit{S}}_{\mathit{l}} a_l+\mathit{e}$$

where Y corresponds to the vector of phenotypic observations (log or sqrt transformed or raw data) for each analysed trait; β is the vector of systematic effects on the trait described above; g is the vector of infinitesimal genetic effects of each individual, with distribution g ∼N(0,G${\sigma }_g^2$), where G is the genomic relationship matrix (GRM) calculated using the filtered autosomal SNPs based on the methodology of (19) and ${\sigma }_g^2$ is the additive genetic variance; S_l is the vector of genotypes (coded as 0, 1, 2) of each individual for the l^th SNP, and a_l is the allele substitution effect of the l^th SNP on the trait under study. The FDR method of multiple testing described by Benjamini and Hochberg (20) was used to measure the statistical significance for association studies at genome-wide level with the p.adjust function of R. The significant association threshold was set at FDR ≤ 0.1. Manhattan plots based on the significance of the 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 available at: https://github.com/YinLiLin/CMplot.

VEP software (22) was used to locate all the significant SNPs and determine their functional effects. Gene annotation for 1Mb downstream/upstream of genomic intervals around the most significant SNPs was performed with BIOMART software (23) using the Ensembl Genes 106 Database of Scrofa11.1 reference assembly. Gene features and gene ontologies (GO) were retrieved. For those pig genes without a name, orthologous human gene names were retrieved. Furthermore, information from Genecards (24) and Mouse Genome Database (25) was used to identify gene functions affecting the analysed phenotypes.

2.8 RNA-Seq of blood

RNA-Seq of whole blood was performed for 30 pigs classified as low (L) or high (H) percentage of T γδ (L=12, H=12), T helper and memory T-helper cells (L=15, H=15), which were selected from the 391 pigs of the Duroc population by a principal component analysis (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 T helper cells (H group: 9.92 ± 1.16; L group: 1 ± 0.18). Since these two populations showed the same direction in the PCA, we considered the same animals per groups. Subsequently, these 30 animals were also classified for extreme proportions of T γδ cells considering the phenotype by (12), where PBMCs were stained with a specific antibody for T γδ cells (APC Rat Anti-Pig γδ T Lymphocytes, clone MAC320, BD Biosciences). For this analysis, we selected 12 animals per group to maximize the differences in the proportions of T γδ cells (H group: 16.36 ± 1.86; L group: 3.64 ± 0.51).

Total RNA from blood was isolated and purified with a Tempus™ Spin RNA Isolation Kit (Thermo Fisher Scientific, Spain), following the manufacturer’s instructions. Total RNA concentration was measured in a Nanodrop One spectrophotometer (Thermo Fisher Scientific, Spain) and RNA purity and integrity was checked by using a Fragment Analyzer equipment (Agilent Technologies, INC., Santa Clara, CA). Libraries were prepared using a Stranded total RNA Prep with Ribo-Zero Plus kit (Illumina Inc., CA) and were paired-end sequenced (2 × 100 bp) on a NovaSeq 6000 platform (Illumina Inc., CA) at the National Centre for Genomic Analysis (CNAG, Barcelona, Spain). More than 50 million of paired-end reads were obtained for all samples. The quality of the raw sequenced reads in the FASTQ files was analysed with the FASTQC software (Babraham Bioinformatics, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were mapped to the reference pig genome Sscrofa11.1 and the annotation database Ensembl Genes 106 by using STAR v. 2.75.3a. Transcript quantification was performed with RSEM v. 1.3.0. The EdgeR R package (26) was used to identify differentially expressed (DE) genes between the two divergent groups (L vs H). After correcting for multiple testing, genes with a Fold Change (FC) between groups greater than 1.2 (i.e. |log2FC| > 0.26) and an FDR< 0.1 were classified as DE. Enriched GO terms and pathway analysis (padj<0.1) of the DE genes was performed with the ClueGO v2.5.8 plug-in of Cytoscape v3.9.1 (27).

3 Results

3.1 Characterization of immune phenotypes

In the present study, we stained PBMCs to measure ten immune-cell-related traits including the relative abundance of different populations of lymphocytes, the proportions of CD4⁺ T cells and CD8⁺ T cells, and the ratio of CD4⁺/CD8⁺ T cells from 391 individuals of a commercial Duroc pig line. All proportions and absolute numbers of cell subsets with descriptive statistics are reported in Table 1. The T cells were the most abundant population of lymphocytes, representing 43.8% of PBMCs. Among T lymphocytes subpopulations, the CTL and naïve T cells were the most abundant (~19% of PBMCs each one), whereas the memory plus helper T cells (i.e. CD4⁺ T cells) represented 5% of PBMCs. The mean proportion of B lymphocytes (21.2% of PBMCs) was less than half that of T lymphocytes, whereas the NK cells represented 14.16% of PBMCs. A high variability between individuals was observed in the relative abundance of some PBMC types, the CV ranging from 0.21 to 0.89. The percentage of T lymphocytes showed the lowest dispersion (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.

Phenotypic correlations were also computed to analyse the association between the relative abundance of the different lymphocyte subsets (Supplementary Table 2). Positive correlations were observed between proportions of T cells and its different subpopulations, ranging from 0.34 to 0.51, but for T helper cells, that exhibited a low correlation with T cells percentage. Among T cells subtypes, the relative abundance of memory T-helper cells correlated positively with both T-helper and CTL cells (r_P= 0.39 and 0.32), whereas an antagonism between the proportions of CTL and naïve T cells was observed (r_P= -0.50). The relative abundance of B lymphocytes showed a strong negative correlation with total T lymphocytes and CTL proportions of PBMCs (r_P= -0.56 and -0.50). NK cells proportion was positively correlated to CTL abundance but showed negative correlations with the relative abundances of naïve T cells and B lymphocytes (r_P= -0.51 and -0.45, respectively). As expected, the proportion of CD4⁺ T cells was highly correlated with T helper and memory T-helper cells. The proportion of CD8⁺ T cells was highly correlated with CTL T cells and moderately correlated with T cells abundance, whereas it showed antagonism with Naïve T cells and B lymphocytes proportions.

Additionally, phenotypic correlations of the different lymphocytes subtypes abundance with a plethora of innate immunity traits previously analysed in the same individuals (12), including plasma immunoglobulins and acute phase proteins concentrations, γδ T lymphocytes, haemogram figures and phagocytosis capacity, were also computed (Supplementary Table 2). Among these traits, only the γδ T lymphocytes showed a strong and highly significant phenotypic association with naïve T cells (r_p = 0.81, p-value = 2.2×10^-16), indicating that naïve T cells (CD3⁺ CD4^- CD8^-) corresponded mainly to γδ T cells. The γδ T cells proportion also showed a moderate correlation with T cells and certain antagonism with NK and CTL cells proportions. Finally, the percentage of phagocytic cells and the phagocytosis capacity of lymphocytes showed low but significant positive associations with the percentage of B cells, jointly with low negative correlations with T cells (CTL and Memory T cells) and NK relative abundances.

3.2 Genetic parameters of T and B-cell populations

The heritability of the different cell immunity subsets was estimated to ascertain the genetic determinism of cellular immunocompetence traits. Medium to high heritabilities were observed for the ten immune-cell-related traits comprising the relative abundance among PBMCs of different lymphocyte populations (Table 2). The heritability estimates ranged between 0.361 to 0.841, being the proportion of memory T-helper cells (h^{2 =} 0.841) the most heritable trait, followed by the total proportion of T cells (h^{2 =} 0.771). In contrast, the percentage of B cells exhibited the lowest heritability value (h^{2 =} 0.361); despite moderate, the 95% confidence interval for its heritability did not encompasses zero (Table 2). The NK cells abundance and the proportion of different T cells subtypes except for memory T-helper cells showed medium heritabilities between 0.470 and 0.558.

TABLE 2

Table 2 Heritability estimates (h²) for the relative abundance (among PBMCs) of different cell immunity subsets, plus standard errors (SE) and confidence intervals at 95% (CI95) of the h² estimates.

When we inferred the genetic correlations between each pairwise combination of traits, positive and negative genetic associations were identified between the proportion (among PBMCs) of some lymphocytes populations (Figure 1 and Supplementary Table 3). In general the genetic correlation pattern matched with that observed with phenotypic correlations, but stronger associations were observed at genetic level despite its significance being limited in a number of cases due to the sample size of the population. Positive and significant genetic correlations were obtained between the proportion of T cells and the different T cells subpopulations (r_g from 0.50 to 0.68) with the exception of T helper cells. Among T cells subpopulations, a positive genetic association of memory T-helper cells with T helper (r_g = 0.58) and CTL (r_g = 0.37) cells was estimated, whereas naïve T cells were negatively correlated with CTL and T helper cells (r_g = -0.37 and -0.36). The proportion of B lymphocytes exhibited negative genetic correlations with the rest of lymphocytes subsets, particularly with the total percentage of T cells (r_g = -0.78), CD4⁺ T cells (r_g = -0.56) and CTL T cells (r_g = -0.41). NK cells proportion were negatively correlated with memory T-helper cells (r_g = -0.44) and to a lesser extent with B lymphocytes (r_g = -0.32) proportion.

FIGURE 1

Figure 1 Heatmap of genetic correlations estimated by pairwise combination among cell immunity traits.

As far as genetic correlation of the analysed PBMCs subsets with innate immunity and other health-related traits (Supplementary Table 3), we identified a strong genetic correlation between the proportion of γδ T lymphocytes and naïve T cells (r_g = 0.903) subpopulation, which corroborates that CD3⁺ CD4^- CD8^- T cells were mainly γδ T cells, jointly with a genetic antagonism between the proportion of γδ T and T helper lymphocytes subpopulation (r_g = -0.568). The CTL abundance showed positive genetic correlation with the number of lymphocytes and other white blood cell counts in the haemogram (r_g from 0.622 to 0.702) as well as with haematocrit (r_g = 0.741). Besides, negative genetic correlations of the proportion of B lymphocytes with both the phagocytic capacity and the percentage of phagocytic cells (r_g = -0.803 and -0.733, respectively) were obtained, in contrast with the positive correlation with the percentage of phagocytic lymphocytes (r_g = 0.544).

3.3 Genomic regions and candidate genes associated with immune-cell traits

To identify genomic regions associated with lymphocytic traits, a GWAS was performed using the ten immune-cell-related traits and the genotypes of 42,641 SNPs of the Porcine GGPSNP70 BeadChip (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).

TABLE 3

Table 3 Description of the four chromosomal regions associated with cell immunity traits and the annotated candidate genes.

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.

FIGURE 2

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).

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 Protein 3 (CHMP3), SET and MYND Domain Containing 1 (SMYD1), Zeta Chain of T cell Receptor Associated Protein Kinase 70 (ZAP70) and Granulysin (GNLY). The three top significant associated SNPs (rs81340900, rs81293514 and rs81336780) were located in an intron of the RMND5A gene. In addition, two significant SNPs at 20.51-20.57 Mb of SSC8 were also associated with the variation of T helper cells (Figure 2A and Table 3). In this region, Recombination Signal Binding Protein for Immunoglobulin Kappa J region (RBPJ) and Stromal Interaction Molecule 2 (STIM2) genes were annotated.

In SSC5, 10 significant SNPs at 61.97 to 62.29 Mb were associated with the proportion of memory T-helper cells (Figure 2B), being rs326461238 the most significant one (Table 2). In this region, a high number of immune-related genes (AICDA, CD4, CD69, CD163, MFAP5, PHC1, STYK1, a family of TAS2Rs, a family of KLR such as KLRB1, KLRC1, KLRK1 and KLRD1, and a family of C-type lectins such as CLECL1, CLEC4D and CLEC4E) were identified. Furthermore, the most associated SNP was found inside a long-noncoding RNA (lncRNA, LOC100524679).

Finally, six SNPs at 33.36 to 33.63 Mb in SSCX were associated 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 associated with immunity phenotypes. In addition, CYBB gene, also related to immunity functions, was found in this region.

3.4 Comparison with other GWAS studies

To identify overlaps between our QTL regions and those previously identified in pigs and humans for immune cellular traits a comparison between genomic regions and candidate genes was performed. Remarkably, we found several candidates genes previously reported as candidates in human GWAS studies for the same or similar phenotypes such as CD8A, CD8B, RMND5A and VPS24 (also known as CHMP3) for CD4⁺ proliferating, CD4⁺CD8^dim T cell Absolute Counts, CD4⁺CD8^dim T cell % lymphocytes and CD4⁺CD8^dim T cell % leukocyte traits (9, 14, 28) and CD4 and KLRC2 for CD4⁺CD8⁺ and effector memory CD4⁺ T cell % T cell traits (9, 28, 29). In pigs, two overlapping or very close genomic regions for T helper cells and memory T-helper cells were identified. Regarding T helper cells trait, overlapping regions for leukocyte percentages of CD8⁺, CD8^- and CD3^-CD8^- traits of an F2 Duroc × Erhualian piglets vaccinated with a live CSF vaccine were identified (15). A region at 55 Mbp in SSC5 for CD4⁺CD8⁺ percentage of leukocytes and located less than 7 Mb away from the region identified in our analysis for CD4⁺CD8⁺ memory T-helper cells was identified in a population of animals consisting of Landrace, Yorkshire, and Songlio Black pigs (16).

3.5 Genes differentially expressed between groups diverging in percentage of T γδ, T helper and T memory cells

The transcriptome of whole blood samples from animals with extreme proportions of T γδ, T helper and memory T-helper cells was analysed by RNAseq to identify genes DE between groups.

A total of 18 genes were identified as DE (FC>1.2; FDR<0.1) in blood between the two groups of animals with high (H) and low (L) T helper and T memory cell percentages. From them, 5 genes were downregulated and 13 were upregulated in the H group when compared to the L group (Figure 3A and Supplementary Table 6). The most up- and downregulated genes in the H group were DNA nucleotidylexotransferase (DNTT; FC = 17.65, P-value = 8.19 x 10^-05) and Cyclin dependent kinase 20 (CDK20; FC⁻¹ = 5.67, P-value = 1.24 x 10^-04), respectively.

FIGURE 3

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) δγ T cells. The vertical axis (y-axis) corresponds to the -log10 (P-value), and the horizontal axis (x-axis) displays the log₂ 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.

The functional annotation of DE genes identified CD4 as being involved in “helper T cell enhancement of adaptive immune response” GO immune system process, among other immune-related functions (Supplementary Table 7). Other genes related to immune functions were KLRD1, NOTCH3 and TCF7L1 (Supplementary Table 7). Furthermore, the following genes CAPG, DNTT, KCNIP3, and P2RY14 have also been reported in the literature to be related to immunity (30–32).

Regarding blood transcriptomic changes between H and L groups for percentages of γδ T cells, 246 DE genes were identified (FC>1.2; FDR<0.1; Figure 3B and Supplementary Table 8). Of these, 102 genes were upregulated and 144 were downregulated in the H group when compared to the L group. The most up- and downregulated genes in the H group were the novel genes ENSSSCG00000046089 (FC = 35.25, P-value = 2.48 x 10^-04) and ENSSSCG00000038561 (FC⁻¹ = 43.43, P-value = 3.41 x 10^-04), respectively. Out of the DE genes, 215 were identified as protein coding genes, 24 as lncRNAs, and 7 as pseudogenes.

After performing a functional analysis with the list of DE genes, we mainly identified GO immune system processes associated with α-β T cell differentiation and activation from the list of downregulated genes in H group, gathering animals with high proportion of γδ T cells, compared to L group. Other processes identified were negative regulation of lymphocyte activation and proliferation. Some genes identified in these processes were CD274, HLX, IL27, LGALS9B, MYB, RORC, RSAD2 and SOCS1 (Supplementary Table 9). Remarkably, we identified in the list of upregulated genes in H group a transcriptional factor (SOX13) associated with positive regulation of γδ T cell activation and differentiation. Immune processes also identified from the list of up-regulated genes were regulation of B cell receptor signalling pathway (BLK and GCSAML), regulation of erythrocyte (ACVR2A, RHEX), and dendritic, monocyte and macrophage (ZBTB46) differentiation, neutrophil-mediated functions (DPP4, CTSG and ADAM8) and T cell differentiation in thymus (ADAM8).

3.6 Coincident association between GWAS and DE analyses

In order to identify potential candidate genes whose differences in gene expression may be modulating our cellular phenotypic traits, we overlapped the lists of DE genes and candidate genes located in associated genomic regions for T helper, memory T-helper and γδ T cells traits. We identified four DE genes between groups of animals with extreme proportions of T helper and memory T-helper cells: CAPG and TCF7L1 on SSC3, and KLRD1 and CD4 on SSC5, at or very close to the genomic regions associated with T helper cells and memory T-helper cell, respectively. No genes colocalizing between GWAS and DE studies were identified for δγ T cells.

4 Discussion

Nowadays, the pig industry is facing new challenges that make necessary to reorient breeding programs to produce more robust and resilient pig populations while maintaining production efficiency. It is known that the immune system and particularly lymphocytes play an important role in controlling many porcine pathogens. In this study, we focused our analysis on peripheral blood mononuclear cell traits, mainly related to adaptive immunity, but also including cell populations such as δγ T cells with functions associated with both arms of immunity.

After PBMC immunophenotyping, the relative abundance among PBMCs of total T and B lymphocytes in the analysed Duroc population was similar to figures reported in literature for other pig populations (33–35). However, differences with other studies were observed in the proportion of some T-cells subsets, with CD4⁺ cells showing lower relative abundance than that reported before (36). Although we could not discard an effect on immunophenotyping associated to handling and freezing processes, other factors such as the age of animals, environmental factors or the lack of standardization of laboratory tests could also have generated discrepancies between studies. It is also plausible that differences in the genetic background of the porcine populations analysed in these studies could be behind the differences in the abundance of certain lymphocytes subsets. For instance, genetic variation of CD4 has been associated to a different reactivity to the CD4 antibody (36). In our study, CD4 was one of the DE genes between animals with extreme phenotypes for T helper and memory T-helper cells, and was located in a genomic region associated with proportion of CD4⁺ and memory T cells, evidencing a genetic contribution to the variation of these traits in our Duroc population. Therefore, changes in CD4 mRNA levels among animals may be the most plausible explanation for the differences in CD4 T cells percentages observed between studies.

A particularly high heritability was estimated for both memory and CD4⁺ T-cells relative abundances, whereas medium to high heritabilities were obtained for the rest of immune-cell-related traits in our Duroc population. These heritability values are in line with those obtained in previous human and pig studies, reaffirming genetics as an important determinant of adaptive cell immune traits (8, 10, 28, 37). The lowest heritability among PBMCs subpopulations relative abundances was observed for B cells. This is concordant with results in humans (29), showing that B cell counts were more environmentally influenced while T cells were more strongly driven by genetic factors.

The map of genetic correlations among PBMC subtypes relative abundances reported positive and negative genetic associations between these immune-cell-related traits, but also reflected the interdependence between some of these compositional variables and should be interpreted with caution. The antagonism between T and B cells proportions in blood observed in our Duroc population might have some physiological implications deserving further research. However, we cannot discard this antagonism to be the consequence of dealing with compositional data. Conversely, we can affirm that our results confirm a positive genetic association between memory and T-helper cells in blood composition. In this regard, two significant genomic regions associated with the proportion of T helper and memory T-helper cells in blood were identified. This may allow selection of animals to investigate abundance of these cells in blood in relation to immune responses toward different stimuli.

Selecting animals for producing higher immune responses should probably consider both innate and adaptive immunity traits in conjunction with qualitative responses. The map of genetic correlations of PBMCs subpopulation proportions with innate immunity traits might help to evaluate the consequences/convenience of selecting for both innate and adaptive immunity traits. Most PBMCs subset abundances (but naïve T cells) showed a negligible genetic association with the concentration of acute phase proteins (CRP and Hp) in blood. Similar results were shown by (10), whereas a positive genetic correlation between the concentration of CRP and B cell counts was described by (38) in animals tested at 90kg of live weight. The genetic independence of acute phase proteins concentration versus most lymphocytes subsets abundance observed in our Duroc population would support the possibility of selecting some immune-cell-related traits independently of some innate immune pathways governing CRP synthesis and then avoiding adverse effects due to inflammation.

Genetic correlations of the different PBMC subpopulations proportions with other immunity traits also pinpoint to other interesting results. The relative abundance of CTL T cells was positively associated at the genetic level with haematological traits, particularly with leucocyte counts and with the haematocrit. Conversely, our results also suggest certain genetic antagonism between the proportion of B lymphocytes and the phagocytic capacity, and between some T cells subpopulations (particularly CTL and memory T-cells) and serum immunoglobulin concentrations. As a whole, the genetic interaction map among innate and adaptive immunity traits obtained in our Duroc population partially support the possibility of applying an effective multi-trait selection for both innate and adaptive immunity. In this scenario of complex genetic relationships, the possibility of identifying genetic markers, biomarkers and candidate genes associated to the phenotypic variation that could allow a more directed selection of some of these traits, without impairing other elements of the immune systems, takes particular relevance.

A genetic correlation close to unity was observed between naïve and δγ T cells quantified with a specific monoclonal antibody (MAC320) in a previous study in the same Duroc population (12), showing that most of the CD4^-CD8^- T cells were in fact γδ T lymphocytes. Consequently, GWAS study identified the same genomic region on SSCX using both phenotypic datasets (12). In the present study, we annotated a new promising candidate gene, CYBB encoding for one of the two chains of cytochrome b (CYB). CYB is a component of the microbicidal oxidase system expressed by neutrophils and T cells, a critical process of innate immunity for defence against bacterial and fungal infections. Mutations in human CYBB have been previously associated to immune system dysfunction and can lead to chronic granulomatous disease (39). Since δγ T cells constitute a high proportion of lymphocytes in porcine peripheral blood and could have the potential to combine conventional adaptive and innate-like responses (35, 40), the activation of WC1-bearing δγ T cells has emerged as a potential alternative to optimize new vaccine strategies (41). In this sense, an increased recruitment of δγ T cells when animals are vaccinated earlier in life appears to correlate with improved efficacy of current vaccines (reviewed in (42)). Indeed, TRDC-knockout pigs, defective in δγ T cells, showed lower neutralizing antibody titres compared to syngenic controls when they were vaccinated against classical swine fever virus using a highly attenuated live virus (43). Thus, our results indicate the possibility to select animals for higher δγ T cells counts in blood to check for potential enhancement of vaccine efficacy.

Another promising candidate genes were located in the SSC3, SSC5 and SSC8 genomic regions associated with the variation of T helper and memory T-helper cells and proportion of CD4⁺ cells. CD8A, CD8B, GNLY and RBPJ were some of the genes associated with the proportion of T helper cells. CD8a and CD8b are integral membrane glycoproteins found on the surface of many immune-cell types and act as co-receptors during antigen recognition by TCR (44). GNLY is an antimicrobial secretion protein present in CTLs and NK cells (45). It also functions as an alarmin, activating dendritic cells and promoting antigen-specific immune responses by using TLR4 (46). Finally, RBPJ encodes a transcription factor acting through the NOTCH signalling pathway to protect T-helper cells from apoptosis (47). Also associated with the percentage of memory T-helper and CD4⁺ cells we identified CD4 and CD69. While CD4 is the main membrane marker for helper (naïve and memory) T-cells, it serves as a co-receptor of the TCR, enhancing recognition of the MHC class II complex (48). Of interest CD69 expressed in resident memory T cells and γδ T cells, is an activation marker determining migration-retention of T cells from tissues (49). Furthermore, we also identified CLECL1 and the killer cell lectin-like receptors genes KLRC1 and KLRD1 which are highly expressed in memory T cells and are involved in their activation (CLECL1) and killing function (KLRs expressed by human CD8⁺ T effector memory cells) (50). Indeed, a particularly remarkable 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 knowledge, this is the first study to report overlapping between humans and pig genes (CD4, CD8A, CD8B, KLRC2, RMND5A and VPS24) associated with the variation of T-cell populations, further supporting the use of pig as a very reliable biomodel for human infectious diseases, 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) that may be modulating T lymphocytes traits through variation of its gene expression. Since transcriptional and post-transcriptional regulation of immune gene expression is critical for modulating immune responses (54, 55), eGWAS studies are required to identify potential cis and/or trans regulatory factors regulating the expression of these genes. However, we cannot rule out that the differences in the expression of some of these genes may be related to the different proportions of these populations in blood. 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 memory cells. Furthermore, no overlapped genes were identified between GWAS and DE analyses for δγ T cells. However, we identified an interesting pattern of gene expression with genes related to immune processes of αβ T cell differentiation and activation downregulated in the group with higher proportion of γδ T cells. This result is in agreement with the upregulation of SOX13 in animals with higher proportion of γδ T cells as this gene encodes a transcription factor essential for normal γδ T cell development (57). It has been reported a higher expression of SOX13 in bovine WC1⁺ γδ T cells relative to WC1^- γδ T cells (58). Remarkably, we also identified three upregulated genes ENSSSCG00000034914 (CD163L1), ENSSSCG00000039217 (cow WC1 orthologue) and ENSSSCG00000031085 (antigen WC1.1) in animals with higher proportions of γδ T cells, consistent with the transcription of WC1 co-receptors by circulating γδ T cells (42). WC1⁺ cells can be divided in different subpopulations with different ability to respond to particular pathogens and cytokine responses (42). A model in cows has been proposed in which a differential occupancy on these WC1 gene loci by SOX13 influence the development of these subpopulations (58). However, still there is a lack of information about these subpopulations in swine (59). Thus, further analysis is required to determine if pigs used in the present report have different proportions of WC1⁺ subpopulations. Besides, it should be noted that the proportion of some lymphocyte subsets may vary with the age of animals (33, 34, 36), so it would be worthy to study the evolution of these lymphocytes traits throughout the life of the animals.

In conclusion this study confirmed the genetic determinism of 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 molecular markers identified for T helper, memory T-helper and γδ T cells support the possibility of applying multi-trait selection for immunocompetence in pigs. Finally, the overlap between genes for the same T cell traits in humans and pigs reinforces the use of pig as a biomodel for vaccine development and the study of human infectious diseases. Overall, our results provide new knowledge about the molecular mechanisms underlying lymphocyte traits with special interest to WC1⁺ γδ T cells to explore further optimization of current and new vaccines.

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

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

Funding

The study was funded by grants AGL2016-75432-R and PID2020-112677RB-C21 awarded by the Spanish Ministry of Economy and Competitiveness (MINECO). D. Crespo-Piazuelo was supported by the GENE-SWitCH project (https://www.gene-switch.eu), which is funded by the European Union’s Horizon 2020 research and innovation programme under the grant agreement n° 817998. T. Jové-Juncà was supported by an IRTA fellowship (CPI1221) and C. Hernández-Banqué was funded with a FPI grant (PRE2021-097825) from the PID2020-112677RB-C21 project. YR-C was financially supported by a Ramon y Cajal contract (RYC2019-027244-I) from the Spanish Ministry of Science and Innovation. The authors belonged to a Consolidated Research Group AGAUR, ref. 2017SGR-1719.

Acknowledgments

We gratefully acknowledge the technical staff from IRTA and Selección Batallé S.A., and specially to Josep Reixach, for their collaboration in the experimental protocols at farm and slaughterhouse, as well as to Albert Bensaid for valuable discussions and comments on the study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1058346/full#supplementary-material

Abbreviations

CTL, cytotoxic T cells; DE, differentially expressed; FBS, foetal bovine serum; FC, fold change; FDR, false discovery rate; GCTA, genome-wide complex trait analysis; GO, gene ontologies; GWAS, genome-wide association studies; h², heritability; H, high; IT, immunity traits; L, low; MAF, minor allele frequency; NK, natural killer; PBMCs, peripheral blood mononuclear cells; Q-Q, quantile-quantile; Rg, genetic correlations; Rp, phenotypic correlations.

References

1. Millet S, Maertens L. The European ban on antibiotic growth promoters in animal feed: from challenges to opportunities. Vet J (2011) 187:143–4. doi: 10.1016/J.TVJL.2010.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bai X, Plastow GS. Breeding for disease resilience: opportunities to manage polymicrobial challenge and improve commercial performance in the pig industry. CABI Agric Biosci (2022) 3:6. doi: 10.1186/S43170-022-00073-Y

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Mellencamp MA, Galina-Pantoja L, Gladney CD, Torremorell M. Improving pig health through genomics: A view from the industry. Dev Biol (Basel) (2008) 132:35–41. doi: 10.1159/000317142

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Visscher AH, Janss LL, Niewold TA, de Greef KH. Disease incidence and immunological traits for the selection of healthy pigs. A review Vet Q (2002) 24:29–34. doi: 10.1080/01652176.2002.9695121

CrossRef Full Text | Google Scholar

5. Mair KH, Sedlak C, Käser T, Pasternak A, Levast B, Gerner W, et al. The porcine innate immune system: An update. Dev Comp Immunol (2014) 45:321–43. doi: 10.1016/j.dci.2014.03.022

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bonilla FA, Oettgen HC. Adaptive immunity. J Allergy Clin Immunol (2010) 125:S33–40. doi: 10.1016/J.JACI.2009.09.017

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liston A, Carr EJ, Linterman MA. Shaping variation in the human immune system. Trends Immunol (2016) 37:637–46. doi: 10.1016/J.IT.2016.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Mangino M, Roederer M, Beddall MH, Nestle FO, Spector TD. Innate and adaptive immune traits are differentially affected by genetic and environmental factors. Nat Commun (2017) 8:13850. doi: 10.1038/NCOMMS13850

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Orrù V, Steri M, Sole G, Sidore C, Virdis F, Dei M, et al. Genetic variants regulating immune cell levels in health and disease. Cell (2013) 155:242. doi: 10.1016/J.CELL.2013.08.041

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Flori L, Gao Y, Laloë D, Lemonnier G, Leplat J-J, Teillaud A, et al. Immunity traits in pigs: substantial genetic variation and limited covariation. PloS One (2011) 6:e22717. doi: 10.1371/journal.pone.0022717

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Clapperton M, Glass EJ, Bishop SC. Pig peripheral blood mononuclear leucocyte subsets are heritable and genetically correlated with performance. Animal (2008) 2:1575–84. doi: 10.1017/S1751731108002929

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Ballester M, Ramayo-Caldas Y, González-Rodríguez O, Pascual M, Reixach J, Díaz M, et al. Genetic parameters and associated genomic regions for global immunocompetence and other health-related traits in pigs. Sci Rep (2020) 10:18462. doi: 10.1038/s41598-020-75417-7

PubMed Abstract | CrossRef Full Text | Google Scholar

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/S1463981500040553

CrossRef Full Text | Google Scholar

14. Lagou V, Garcia-Perez JE, Smets I, Van Horebeek L, Vandebergh M, Chen L, et al. 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

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhang J, Chen JH, Liu XD, Wang HY, Liu XL, Li XY, et al. Genomewide association studies for hematological traits and T lymphocyte subpopulations in a duroc × erhualian F2 resource population. J Anim Sci (2016) 94:5028–41. doi: 10.2527/jas.2016-0924

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Lu X, Fu WX, Luo YR, Ding XD, Zhou JP, Liu Y, et al. Genome-wide association study for T lymphocyte subpopulations in swine. BMC Genomics (2012) 13:488. doi: 10.1186/1471-2164-13-488

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Misztal I, Tsuruta S, Strabel T, Auvray B, Druet T, Lee D. 7th world congress on genetics applied to livestock production, August 19-23, 2002, Montpellier, France. 7th World Congress Genet Appl to Livestock Production (2002).

Google Scholar

18. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet (2007) 81:559–75. doi: 10.1086/519795

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: A tool for genome-wide complex trait analysis. Am J Hum Genet (2011) 88(1):76–82. doi: 10.1016/j.ajhg.2010.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc Ser B (1995) 57:289–300. doi: 10.2307/2346101

CrossRef Full Text | Google Scholar

21. Turner SD. Qqman: an r package for visualizing GWAS results using q-q and manhattan plots. J Open Source Softw (2018) 3:731. doi: 10.21105/JOSS.00731

CrossRef Full Text | Google Scholar

22. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The ensembl variant effect predictor. Genome Biol (2016) 17:122. doi: 10.1186/s13059-016-0974-4

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Smedley D, Haider S, Durinck S, Pandini L, Provero P, Allen J, et al. The BioMart community portal: an innovative alternative to large, centralized data repositories. Nucleic Acids Res (2015) 43:W589–98. doi: 10.1093/nar/gkv350

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Safran M, Solomon I, Shmueli O, Lapidot M, Shen-Orr S, Adato A, et al. GeneCardsTM 2002: towards a complete, object-oriented, human gene compendium. Bioinformatics (2002) 18:1542–3. doi: 10.1093/bioinformatics/18.11.1542

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Eppig JT, Smith CL, Blake JA, Ringwald M, Kadin JA, Richardson JE, et al. Mouse genome informatics (MGI): Resources for mining mouse genetic, genomic, and biological data in support of primary and translational research. Methods Mol Biol (Clifton N.J.) (2017), 47–73. doi: 10.1007/978-1-4939-6427-7_3

CrossRef Full Text | Google Scholar

26. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, Kirilovsky A, et al. ClueGO: a cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics (2009) 25:1091–3. doi: 10.1093/bioinformatics/btp101

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Orrù V, Steri M, Sidore C, Marongiu M, Serra V, Olla S, et al. Complex genetic signatures in immune cells underlie autoimmunity and inform therapy. Nat Genet (2020) 52:1036–45. doi: 10.1038/s41588-020-0684-4

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Aguirre-Gamboa R, Joosten I, Urbano PCM, van der Molen RG, van Rijssen E, van Cranenbroek B, et al. Differential effects of environmental and genetic factors on T and b cell immune traits. Cell Rep (2016) 17:2474–87. doi: 10.1016/j.celrep.2016.10.053

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Aono A, Enomoto H, Yoshida N, Yoshizaki K, Kishimoto T, Komori T. Forced expression of terminal deoxynucleotidyl transferase in fetal thymus resulted in a decrease in gammadelta T cells and random dissemination of Vgamma3Vdelta1 T cells in skin of newborn but not adult mice. Immunology (2000) 99:489–97. doi: 10.1046/J.1365-2567.2000.00987.X

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Li Q, Xu L, Li Y, Yang R, Qiao Q, Wang Y, et al. P2RY14 is a potential biomarker of tumor microenvironment immunomodulation and favorable prognosis in patients with head and neck cancer. Front Genet (2021) 12:670746. doi: 10.3389/FGENE.2021.670746

PubMed Abstract | CrossRef Full Text | Google Scholar

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(+/+) mice to listeria monocytogenes infection but not to salmonella enterica serovar typhimurium infection. Infect Immun (2003) 71:6582–90. doi: 10.1128/IAI.71.11.6582-6590.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yang H, Parkhouse RME. Phenotypic classification of porcine lymphocyte subpopulations in blood and lymphoid tissues. Immunology (1996) 89:76–83. doi: 10.1046/j.1365-2567.1996.d01-705.x

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Pietrasina O, Miller J, Rząsa A. Differences in the relative counts of peripheral blood lymphocyte subsets in various age groups of pigs(2020) (Accessed September 14, 2022).

Google Scholar

35. Takamatsu HH, Denyer MS, Stirling C, Cox S, Aggarwal N, Dash P, et al. Porcine gammadelta T cells: possible roles on the innate and adaptive immune responses following virus infection. Vet Immunol Immunopathol (2006) 112:49–61. doi: 10.1016/J.VETIMM.2006.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Blanc F, Prévost-Blondel A, Piton G, Bouguyon E, Leplat JJ, Andréoletti F, et al. The composition of circulating leukocytes varies with age and melanoma onset in the MeLiM pig biomedical model. Front Immunol (2020) 11:291. doi: 10.3389/FIMMU.2020.00291

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Clapperton M, Glass EJ, Bishop SC. Pig peripheral blood mononuclear leucocyte subsets are heritable and genetically correlated with performance. Animal (2008) 2:1575–84. doi: 10.1017/S1751731108002929

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Clapperton M, Diack AB, Matika O, Glass EJ, Gladney CD, Mellencamp MA, et al. Traits associated with innate and adaptive immunity in pigs: heritability and associations with performance under different health status conditions. Genet Sel Evol (2009) 41:54. doi: 10.1186/1297-9686-41-54

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Rae J, Newburger PE, Dinauer MC, Noack D, Hopkins PJ, Kuruto R, et al. X-Linked chronic granulomatous disease: Mutations in the CYBB gene encoding the gp91-phox component of respiratory-burst oxidase. Am J Hum Genet (1998) 62:1320–31. doi: 10.1086/301874

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Vantourout P, Hayday A. Six-of-the-best: unique contributions of γδ T cells to immunology. Nat Rev Immunol (2013) 13:88–100. doi: 10.1038/NRI3384

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Teixeira AF, Gillespie A, Yirsaw A, Britton E, Telfer JC, Nascimento ALTO, et al. Identification of leptospiral protein antigens recognized by WC1 + γδ T cell subsets as target for development of recombinant vaccines. Infect Immun (2022) 90(1):e0049221. doi: 10.1128/IAI.00492-21

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Le Page L, Baldwin CL, Telfer JC. γδ T cells in artiodactyls: Focus on swine. Dev Comp Immunol (2022) 128:104334. doi: 10.1016/J.DCI.2021.104334

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Petersen B, Kammerer R, Frenzel A, Hassel P, Dau TH, Becker R, et al. Generation and first characterization of TRDC-knockout pigs lacking γδ T cells. Sci Rep (2021) 11:14965. doi: 10.1038/S41598-021-94017-7

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Cole DK, Laugel B, Clement M, Price DA, Wooldridge L, Sewell AK. The molecular determinants of CD8 co-receptor function. Immunology (2012) 137:139–48. doi: 10.1111/J.1365-2567.2012.03625.X

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Crespo ÂC, Mulik S, Dotiwala F, Ansara JA, Sen Santara S, Ingersoll K, et al. Decidual NK cells transfer granulysin to selectively kill bacteria in trophoblasts. Cell (2020) 182:1125–1139.e18. doi: 10.1016/J.CELL.2020.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Tewary P, Yang D, de la Rosa G, Li Y, Finn MW, Krensky AM, et al. Granulysin activates antigen-presenting cells through TLR4 and acts as an immune alarmin. Blood (2010) 116:3465–74. doi: 10.1182/BLOOD-2010-03-273953

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Helbig C, Gentek R, Backer RA, De Souza Y, Derks IAM, Eldering E, et al. Notch controls the magnitude of T helper cell responses by promoting cellular longevity. Proc Natl Acad Sci U.S.A. (2012) 109:9041–6. doi: 10.1073/PNAS.1206044109

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Gascoigne NRJ, Zal T, Yachi PP, Hoerter JAH. Co-Receptors and recognition of self at the immunological synapse. Curr Top Microbiol Immunol (2010) 340:171–89. doi: 10.1007/978-3-642-03858-7_9

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Cibrián D, Sánchez-Madrid F. CD69: from activation marker to metabolic gatekeeper. Eur J Immunol (2017) 47:946–53. doi: 10.1002/eji.201646837

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Weng NP, Araki Y, Subedi K. The molecular basis of the memory T cell response: differential gene expression and its epigenetic regulation. Nat Rev Immunol (2012) 12:306–15. doi: 10.1038/NRI3173

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Käser T. Swine as biomedical animal model for T-cell research-success and potential for transmittable and non-transmittable human diseases. Mol Immunol (2021) 135:95–115. doi: 10.1016/J.MOLIMM.2021.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

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/j.tim.2011.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Pabst R. The pig as a model for immunology research. Cell Tissue Res (2020) 380:287–304. doi: 10.1007/S00441-020-03206-9

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Smale ST. Transcriptional regulation in the immune system: A status report. Trends Immunol (2014) 35:190. doi: 10.1016/J.IT.2014.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Schwerk J, Savan R. Translating the untranslated region. J Immunol (2015) 195:2963. doi: 10.4049/JIMMUNOL.1500756

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Karlsson M, Zhang C, Méar L, Zhong W, Digre A, Katona B, et al. A single-cell type transcriptomics map of human tissues. Sci Adv (2021) 7: eabh2169. doi: 10.1126/SCIADV.ABH2169

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Melichar HJ, Narayan K, Der SO, Hiraoka Y, Gardiol N, Jeannet G, et al. Regulation of gammadelta versus alphabeta T lymphocyte differentiation by the transcription factor SOX13. Science (2007) 315:230–3. doi: 10.1126/SCIENCE.1135344

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Damani-Yokota P, Zhang F, Gillespie A, Park H, Burnside A, Telfer JC, et al. Transcriptional programming and gene regulation in WC1 + γδ T cell subpopulations. Mol Immunol (2022) 142:50–62. doi: 10.1016/J.MOLIMM.2021.12.016

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Le Page L, Gillespie A, Schwartz JC, Prawits LM, Schlerka A, Farrell CP, et al. Subpopulations of swine γδ T cells defined by TCRγ and WC1 gene expression. Dev Comp Immunol (2021) 125:104214. doi: 10.1016/J.DCI.2021.104214

PubMed Abstract | CrossRef Full Text | Google Scholar