Blood samples from swine were obtained from a local abattoir. Prior to blood sampling, animals were anesthetized electrically and sacrificed by exsanguination in accordance with Austrian Animal Welfare Slaughter Regulation. PBMCs were isolated from fresh heparinized blood of six animals of approximately six months of age by density gradient centrifugation (Pancoll human, density: 1.077 g/ml, PAN-Biotech, Aidenbach, Germany; 30 min at 920 x g).

CD8⁺ T cells were enriched by positive selection of CD8β-labeled PBMCs using magnetic-activated cell sorting (MACS, Miltenyi Biotec, Bergisch Gladbach, Germany). For enrichment of CD8β⁺ T cells, freshly isolated PBMCs (1 x 10⁹) were stained with an in-house produced primary monoclonal anti-CD8β antibody (clone PPT23, IgG1) for 20 min on ice. Subsequently, cells were washed once with MACS buffer [PBS w/o Ca/Mg + 2% (v/v) FCS (both Gibco™, Thermo Fisher Scientific) + 2mM EDTA (Carl Roth)], resuspended in 1,5 mL MACS buffer and incubated with magnetically labeled secondary antibody (rat-anti mouse IgG1, Miltenyi Biotec) for 30 min on ice. After a further washing step, cells were resuspended in 3 mL MACS buffer and loaded on pre-wetted LS columns (Miltenyi Biotec). The columns were applied to a magnetic field and unlabeled cells were removed by extensive washing. For final elution of the positive fraction, columns were removed from the magnetic field and CD8β⁺ T cells were eluted in 5 mL MACS buffer. Finally, sorted cells were resuspended in cold culture medium (RPMI 1640 + 100 IU/mL penicillin + 0.1 mg/mL streptomycin (all PAN Biotech) + 10% (v/v) FCS), centrifuged and counted with a Cell Counter (XP-300 Hematology Analyzer, Sysmex Europe GmbH, Norderstedt). Purity of the positively sorted cells was over 90% (FACSCanto™II, BD Biosciences, San Jose, CA, USA).

In order to further separate MACS-enriched CD8β⁺ cells into subpopulations, CD8β⁺ cells were FACS sorted based on surface expression of CD27 and CD11a ( Supplementary Figure S1 ).

Upon magnetic-activated cell sorting, CD8β⁺ cells were washed once with FACS buffer (RPMI 1640 + 100 IU/mL penicillin + 0.1 mg/mL streptomycin + 5% FCS + 5% porcine plasma (in-house preparation) + 2 mM EDTA) and then labeled with a goat anti-mouse IgG1-PE secondary antibody to stain residual CD8β⁺ cells (Southern Biotech, Birmingham, AL, USA).

Free binding sites of the PE-labeled antibody were blocked with whole mouse IgG molecules (2 μg per sample, ChromPure, Jackson ImmunoResearch, West Grove, PA, USA). Afterwards, cells were incubated with directly labeled primary antibodies: CD27-Alexa647 (b30c7, mouse IgG1, in-house preparation and labeling with Alexa Fluor-647 Protein Labeling Kit, Thermo Fisher Scientific) and CD11a-FITC (BL1H8, mouse IgG2b, BioRad, Hercules, CA, USA).

Cell sorting was performed on a FACSAria (BD Biosciences) and CD8⁺ T-cell subsets were defined as follows: naïve (CD8β⁺CD27⁺CD11a^low), intermediate differentiated (CD8β⁺CD27^dimCD11a⁺), and terminally differentiated cells (CD8β⁺CD27^-CD11a^high). Subsets were sorted with an average purity greater than 96%.

To identify transcriptomic differences between the CD8⁺ T-cell subsets as well as between ex vivo and stimulated cells within the same CD8⁺ T-cell subset, cells from each sorted subpopulations with at least 5 x 10⁵ sorted cells were cultivated at 37°C and 5% CO₂ under following conditions: (i) 16 hours, unstimulated in culture medium, (ii) cultivation in culture medium for 14 hours followed by stimulation for two hours with phorbol 12-myristate 13-acetate (PMA, 50 ng/mL, Sigma-Aldrich, Schnelldorf, Germany) and ionomycin (500 ng/mL, Sigma-Aldrich), (iii) stimulated with concanavalin A (ConA) (5 μg/mL, Amersham Biosciences, Uppsala, Sweden) for 16 hours. Both stimulation protocols are established in our laboratory and used as high controls for proliferation experiments and cytokine induction in ELISpot assays (ConA) and as positive control for intracellular cytokine staining in flow cytometry (PMA/ionomycin). Furthermore, each CD8⁺ T-cell subset with 5 x 10⁵ was used immediately after sorting for RNA isolation without any further cell culture (ex vivo). Altogether four different conditions for each CTL subset were applied: cultivation in medium, stimulation with PMA/ionomycin or ConA and ex vivo isolation. Therefore, 72 samples (3 subsets x 4 conditions x 6 animals) were generated.

Total RNA was isolated from the samples mentioned above using RNeasy Mini Kit with on-column DNase treatment using the RNAse-Free DNase Set (both Qiagen, Hilden, Germany), following manufacturer’s protocol. Quantification and quality control of isolated RNA were assessed with both Qubit 3.0 fluorometer (RNA HS assay kit, ThermoFisher, Massachusetts, MA, USA) and Agilent 2100 Bioanalyzer (Agilent RNA 6000 Pico Kit, Agilent Technologies, Palo Alto, CA, USA). Samples with both a final yield comprised between 0.03 – 1.25 ng/µl and a RIN of 9 were prepared for sequencing with the SMARTer Stranded Total RNA-Seq v2 – Pico Input Mammalian Kit (Takara Bio Inc., Shiga, Japan). Fully automated library preparation was performed on a Microlab Star Hamilton robotic station (Hamilton Company, Reno, NV, USA). Briefly, 8 µl per sample were used for the cDNA synthesis via the SMART^® technology (SMART technology, Clontech, USA). Thereafter, each sample was amplified to generate Illumina-compatible libraries according to the manufacturer’s guidance. Libraries were validated using the Agilent 2100 Bioanalyzer (Agilent High Sensitivity DNA Kit, Agilent Technologies, Palo Alto, CA) and the Qubit 3.0 fluorometer (DNA HS assay kit, ThermoFisher, Massachusetts, MA, USA). Libraries were paired-end sequenced on two SP flow cell on NovaSeq 6000 system (Illumina Inc., San Diego, CA, USA).

Standard raw sequencing data in BCL format was converted to FASTQ files using the software bcl2fastq v2.19.1.403. After importing the FASTQ files into CLC Genomics Workbench 21.0.3 (Qiagen, Aarhus, Denmark), the reads were adapter- and quality trimmed. Prior to mapping, sequence reads were trimmed using quality score (Phred score ≤ 25) and with maximum number of 2 ambiguous nucleotides allowed. Next, the adapter sequences were trimmed off according to the Illumina Adapter List. Reads shorter than 35 and longer than 75 nucleotides were discarded.

The filtered reads were mapped to the Sus scrofa 11.1 reference genome from NCBI database (GCA_000003025.6) using default parameters of CLC Genomics RNA-Seq Analysis tool (mismatch cost = 2, insertion cost = 3, deletion cost = 3, length fraction = 0.8 and similarity fraction = 0.8). For principal component analysis (PCA), mapped reads were TMM normalized, log CPM values calculated and Z-normalization performed. For the ex vivo condition, differential gene expression test for differences between all pairs of CD8⁺ T-cell subsets using Wald test was performed. Therefore, three pairwise comparisons were made: (i) naïve vs. terminally differentiated, (ii) intermediate vs. terminally differentiated, and (iii) naïve vs. intermediate differentiated. To assess the effect of stimulation on gene expression profiles of CD8⁺ T-cell subsets, Wald test with medium condition as control group was used. Correspondingly, that yielded two pairwise comparisons for each CD8⁺ T-cell subset: (i) ConA stimulation vs. medium and (ii) PMA/ionomycin stimulation vs. medium. As criteria to define differentially expressed genes (DEGs), fold-change > |2|, maximum of the average reads per kilobase per million mapped reads (RPKM’s) > 2 and a false discovery rate corrected p-value < 0.01 (FDR) were used. Venn diagram and heat map visualization of DEGs were constructed using ggvenn and pheatmap packages in R software version 4.0.2 (R Core Team, GNU General Public License). Bar charts were visualized with Tableau Desktop 2020.3 (Tableau Software Inc.).

For DEGs, gene ontology (GO) and enrichment analysis for immune system processes were executed using the ClueGO v.2.5.8 plug-in in the bioinformatic software Cytoscape 3.8.2. version (https://cytoscape.org). The analysis was performed for upregulated genes between CD8⁺ T-cell subsets and based on GO data for Sus scrofa. Following cut-off thresholds were set: at least 3 genes per GO term, two-sided hypergeometric statistical testing corrected with the Bonferroni step-down method (p < 0.05) and a Kappa score of 0.4. Moreover, organism-specific pathway analysis of DEGs were constructed by using KEGG mapper based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database with KEGG Orthology (KO) assignment.

Based on our hypothesis that within the CD8β⁺ T-cell subpopulation three subsets with distinct differentiation stages can be defined, we analyzed the presumable naïve (T_n; CD8β⁺CD27⁺CD11a^low), intermediate differentiated (T_inter; CD8β⁺CD27^dimCD11a⁺), and terminally differentiated cells (T_term; CD8β⁺CD27^-CD11a^high).

In total 3.59 billion paired-end reads were generated by sequencing 72 libraries. Overall, the percentage of mapping reads to the reference genome was between 90.44% and 94.87% (mean = 93.1%) with approximately 50 million paired-end reads per sample. PCA of gene expression data from all ex vivo CD8⁺ T-cell subsets revealed distinguishable differences between CTL subsets as PCA plot clustered data into three distinct groups ( Figure 1A ).

For hierarchical cluster analysis, we selected 1439 genes, which were significantly expressed in at least one pairwise comparison between CD8⁺ T-cell subsets (as defined in the Methods section). Afterwards, a heat map based on their gene expression values was generated. Notably, the hierarchical clustering of selected genes identified three well-defined groups of samples. The first contained all T_n samples, the second all T_term samples and the third all T_inter samples ( Figure 1B ). Genes highly expressed in T_n were downregulated in T_term and vice versa. This clear separation regarding gene expression could indicate transcriptional switch that CD8⁺ T cells undergo while differentiating from naïve to terminally differentiated CD8⁺ T cells. In comparison to T_n and T_term, T_inter showed upregulation of genes expressed in both groups. However, Venn diagram analysis showed that T_inter and T_term share more DEGs (n=386) than T_inter and T_n (n=130) ( Figure 1C ). In contrast, only one upregulated DEGs was shared between the T_n and T_term when compared to T_inter. Next, using Wald test for pairwise comparison, 575 and 709 DEGs were identified as upregulated in T_n and T_term, respectively ( Supplementary Table S1 ). The number of upregulated DEGs was smaller in T_inter vs. T_term comparison than T_n vs. T_inter comparison. A higher number of upregulated DEGs (n = 492) was observed in T_inter compared to T_n (n = 215) CD8⁺ T cells. Also, higher numbers of upregulated DEGs were discovered in T_inter (n = 208) than in T_term (n = 132). To obtain further information about each stage of CD8 T-cell differentiation, gene expression profiles were compared between T_n, T_inter and T_term. We found that genes related to early stages of CD8 T-cell differentiation were highly expressed in the T_n ( Table 1 ). Expression of several genes encoding transcription factors associated with naïve lymphocytes (27), including LEF1, BACH2, TCF7 (TCF1), SATB1, ZEB1 and BCL2 were markedly increased in the T_n. In contrast, genes encoding transcription factors associated with terminally differentiated effector cells, such as TBX21 (T-bet), PRDM1 (Blimp-1), ZEB2, ZNF683 (Hobit), BATF, EZH2 and ID2 were highly upregulated in the T_term.

Furthermore, T_term showed high expression of several genes involved in cell adhesion and migration including CX3CR1, CCR5, CCL4 and CCL5. Moreover, higher expression of adhesion genes ITGAM (CD11b) and ITGAL (CD11a) (28) was observed among T_term compared with T_n and T_inter. Expression of ITGA4 (CD49d), which together with CD44 is expressed in effector T cells and effector memory T cells (13), was upregulated in T_inter and T_term. In addition, expression of CD44 was increased in both T_inter and T_term but not in the T_n ( Supplementary Table S1 ). Conversely, genes encoding lymph node homing receptor molecules such as CCR7, SELL (CD62L) and CCR9 were highly upregulated in the T_n. Sphingosine-1-Phosphate Receptor 1 (S1PR1), important for lymphocyte trafficking and upregulated in human naïve T cells (29), was also increased in the porcine T_n. Also, T_n showed high expression of genes encoding CD27 and CD28 molecules, the former in accordance with cell surface expression used for the sorting strategy.

We observed that several genes involved in T-cell effector functions and cytolytic killing, including GNLY (Granulysin), PRF1 (Perforin), GZMB (Granzyme B), FAS, FASL, IFNG and TNF, were highly increased in T_term in comparison to T_n or T_inter. Notably, T_term expressed the GNLY 1457-fold higher in comparison to T_n. Moreover, T_term showed high expression of KLRG1, KLRD1 and KLRK1, whereas T_n displayed high mRNA levels of IL-7R (CD127). In mouse a selective expression of IL-7R (CD127) is used for the discrimination between MPEC and SLEC, with the high expression specific for MPEC (8). In addition to the high expression of IL-7R, human MPEC show low expression of KLRG1, while SLEC show upregulation of KLRG1 and low expression of IL7R (13).

We found higher expression of co-inhibitory molecule PDCD1 (PD-1) in T_inter and T_term when compared to the T_n. Previous research suggests that high expression of PDCD1 (PD-1) is specific for SLEC formation, whereas low PDCD1 expression contributes to the T effector memory generation (30). Several genes encoding cytokine receptors associated with effector T cells were increased in the T_term, including IL2RB (CD122), IL2RG (CD132), IL12RB1 and IL12RB2. In comparison to the T_n, T_inter and T_term showed high expression of IRF8, which supports the transition from naïve to effector CD8⁺ T cells in independent matter to T-bet and Eomes (31). Furthermore, upregulation in transcript levels of ITGB2 (CD18) and ANXA2, known to be increased in CD8⁺ effector T cells (32), as well as LGALS1, which is expressed only on activated CD8⁺ effector T cells but not resting CD8⁺ T cells (33), were observed in T_inter and T_term. Additionally, genes strongly linked to cytotoxic T cells such as S1PR5 and ADGRG1 were substantially upregulated in the T_term. By contrast, T_n showed high expression of genes, which enforce quiescence state of naïve T cells (MYB, FOXP1, KLF9 and SOCS3). In comparison to T_n, we found other members of SOCS family, namely SOCS1 and SOCS7, highly expressed in T_term. Furthermore, expression of MKI67, encoding proliferation marker Ki-67, was upregulated in T_inter and T_term. Both TNFRSF1A (TNFR1) and TNFRSF1B (TNFR2) were upregulated in T_inter and T_term. While transcripts of TNFSF12 (TWEAK) and TNFSF10 (TRAIL) were upregulated in T_inter and T_term, expression of costimulatory TNFRSF25 (DR3) was highly induced in the T_n. Only T_term expressed high levels of GZMM and STAT4, on the other hand T_inter showed upregulation of GZMA2. Notably, the expression of RUNX3, which is important for the acquisition and maintenance of cytolytic functions of CD8⁺ effector T cells (34), was upregulated in T_inter and T_term. Furthermore, T_inter and T_term showed increased levels of TNFAIP2 and NFKBIE. Regarding genes involved in metabolism, we observed high expression of ARNTL and SLC1A5 in more differentiated CTL subsets, whereas expression of HIF1A was upregulated in the T_n. Interestingly, when compared to T_n and T_term, T_inter shared a more comparable gene expression profile associated with later stages of T-cell differentiation. In particular, most of genes highly expressed in T_term were also upregulated in T_inter. However, the difference in expression of genes related to early stages of T-cell differentiation was substantially smaller between T_n and T_inter than T_n and T_term. Also, those genes were higher expressed in T_inter than T_term.

For better understanding of the relationship between porcine, human and mouse CD8⁺ T cells we assessed the orthology of their genes expressed in CD8⁺ T cells in corresponding subsets publicly available on GEO Data sets (NCBI) under GDS3834 and GDS592. Here, we have focused on the analysis of DEGs in T_n and T_term, which cover the vast majority of DEGs generated from porcine CD8⁺ T-cell subsets. When compared to DEGs in porcine T_n, we found 495 (86.1%) orthologs in human and 226 (39.3%) in mouse data set. In case of T_term, out of 709 DEGs, 612 (86.3%) were recorded in human, and 277 (39.1%) in mouse data set ( Figure 1D ).

In order to further highlight the heterogeneity in gene expression between CD8⁺ T-cell subsets, cells were analyzed upon stimulation with ConA and PMA/ionomycin and compared to cells cultured in medium control. Overall, a substantially higher number of upregulated DEGs in all three CD8⁺ T-cell subsets was observed in response to PMA/ionomycin compared to ConA stimulation. Additionally, gene expressions of all PMA/ionomycin-stimulated CD8⁺ T-cell subsets are clearly distinct from all other CD8⁺ T-cell subsets as showed in PCA plot ( Figure 2A ). Further investigation of CD8⁺ T-cell subsets stimulated with PMA/ionomycin revealed the highest number of upregulated DEGs in the T_term (1717), followed by the T_inter (1667) and the T_n (1383). Conversely, in CD8⁺ T-cell subsets stimulated with ConA, the highest number of DEGs was found in the T_n (100), followed by the T_inter (35) and the T_term (36) ( Supplementary Table S1 ). In order to obtain a more detailed view upon PMA/ionomycin stimulation, PCA was performed additionally only on PMA/ionomycin-stimulated CD8⁺ T-cell subsets ( Figure 2B ). Interestingly, CD8⁺ T-cell subsets clustering is unaltered to PMA/ionomycin stimulation, resulting again in the three distinct groups of T_n, T_inter and T_term. Despite this separate clustering, Venn diagram analysis revealed high number of DEGs shared between PMA/ionomycin-stimulated CD8⁺ T-cell subsets (903), indicating that all three subsets acquire more similar cell properties following PMA/ionomycin stimulation. In case of ConA-stimulated CD8⁺ T-cell subsets, we found much smaller number of shared DEGs ( Figure 2C, D ).

Several genes encoding cytokines involved in T-cell response were differentially expressed between CD8⁺ T-cell subsets ( Figure 3A ). Looking at the expression of IFNG (IFN-γ) and TNF, we observed overexpression in all three CD8⁺ T-cell subsets following PMA/ionomycin stimulation, but only moderate expression in T_inter with ConA stimulation. Porcine T_n and T_inter showed high expression of IL2 and its receptor chains IL2RA (CD25) and IL2RG (CD132) as well as IRF7 when stimulated with PMA/ionomycin. It was reported that IL2 and its receptor chains IL2RA (CD25) and IL2RG (CD132) are involved in terminal effector differentiation but also in memory development of CD8⁺ T cells (37). Moreover, expression of IL4, IL17A, IL18RAP and IL22 was induced only in the T_inter stimulated with PMA/ionomycin. In contrary, IL12RB1, IL27RA and ILF3 were only expressed by the T_term. Nevertheless, both CD8⁺ T-cell subsets showed high expression of IL10 and IRF2BP2 after PMA/ionomycin stimulation. Although the expression of IRF4 was upregulated in all three CD8⁺ T-cell subsets following PMA/ionomycin and ConA stimulations, the highest expression was induced by PMA/ionomycin-stimulated T_n followed by T_inter and T_term. Consistently, studies in mice showed that IRF4 contributes to expansion and maintenance of effector functions of CTL as well as to memory formation of CTL (36). In case of IRF8 transcript, we found similar expression between CD8⁺ T-cell subsets stimulated with PMA/ionomycin. In comparison, ConA stimulation induced a much smaller extent expression of IRF8 in T_inter and T_term. Expression of both genes, IL6ST and ILF2 were similarly increased in all three CD8⁺ T-cell subsets upon PMA/ionomycin stimulation. Remarkably, the highest expression of IL4R, IL15RA and IRF1 was recorded in the T_n followed by T_inter and T_term. Of interest, we found three genes of TNF-induced proteins, namely TNFAIP2, TNFAIP3 and TNFAIP8 highly expressed in CD8⁺ T-cell subsets following PMA/ionomycin stimulation. Moreover, expression of TNFAIP2 and TNFAIP3, both inhibiting canonical NF-kB signaling pathway and thus negatively effecting cytokine production, was highest in the T_n (38, 39). Apart from its aforementioned functions, TNFAIP3, also highly expressed on naïve T cells, restricts MAP kinases and CD8⁺ T cell proliferation (40).

Chemokines and chemokine receptors play a pivotal role in attracting and guiding the naïve and effector T cells to lymph nodes and sites of inflammation, respectively (41). Overall, in all three CD8⁺ T-cell subsets the PMA/ionomycin stimulation induced stronger expression of genes associated with chemokines than after ConA stimulation ( Figure 3B ). Expression of CCL4 (MIP-1ß) and XCL1 (ATAC/lymphotactin), the inflammatory chemokines secreted by activated CD8⁺ T cell (42), was induced in all three CD8⁺ T-cell subsets upon both stimulations, although with significantly higher increase in PMA/ionomycin-stimulated CD8⁺ T-cell subsets. All three PMA/ionomycin-stimulated CTL subsets showed similar expression of CCL3L1, while following ConA stimulation it was increased in T_inter and T_term. Interestingly, only T_inter stimulated with PMA/ionomycin showed significant increase in CCL20, CXCL8 and CXCL10 expression, later known as one of interferon-inducible ligands of CXCR3 (43). The PMA/ionomycin stimulation induced also transcriptional upregulation of CCL5 (RANTES) and CXCL16 in all three CD8⁺ T-cell subsets. Furthermore, transcription of CCL1 was increased in T_inter and T_term upon PMA/ionomycin stimulation.

Transition from naïve T cell to activated effector T cell is accompanied by metabolic adjustment necessary for specific cellular functions (44). Overall, the PMA/ionomycin stimulation induced stronger expression of genes linked to T-cell metabolism in comparison to the ConA stimulation. T_n and T_inter upregulated BCAT1 and GCLC upon PMA/ionomycin stimulation, while T_term were enriched in transcripts for LDHA and TPI1 gene. Both T_inter and T_term induced high expression of PDPK1 and SLC2A1, whereas FASN, GLS and TPP2 were similarly expressed by all three CD8⁺ T-cell subsets following PMA/ionomycin stimulation. In comparison to T_n, we recorded higher expression of HIF1A and SLC7A5 in T_inter and T_term. Upon PMA/ionomycin stimulation, the highest levels of SLC1A5, HK2 and MYC were expressed in T_n, T_inter and T_term, respectively. The ConA stimulation induced upregulation of ID2 expression in T_n and SLC1A5 in T_inter only. Contrary, the PMA/ionomycin stimulation induced upregulation of ID2 in all three CD8⁺ T-cell subsets ( Figure 3B ).

Next, we examined the impact of stimulation with PMA/ionomycin and ConA on expression of transcription factor genes in CD8⁺ T-cell subsets ( Figure 4A ). Several genes encoding transcription factors associated with terminally differentiated effector cells, including BATF, BATF3, EZH2, MYC and TBX21 were upregulated in all three CD8⁺ T-cell subsets upon PMA/ionomycin stimulation but not after ConA stimulation. While the highest expression of BATF3, EZH2 and MYC was observed in the T_term, highest upregulation of TBX21 which encodes the T-bet, the master regulator of cytotoxic T-cell development (45), was observed in the T_n compared to T_inter and T_term. Interestingly, upregulation of EOMES and ID3 was limited only to the ConA-stimulated T_n. T_inter and T_term displayed upregulation of FOXO1, FOXP1, PRDM1 (Blimp-1), SATB1 and SREBF2 following PMA/ionomycin stimulation. Recent studies in mice and human have shown that Blimp-1, encoded by PRDM1, enhances formation of SLECs, production of IL-10 and cytotoxic functions of CD8⁺ T cells (46). Along the PMA/ionomycin-stimulated differentiation subsets, we observed a gradual increase of expression of EGR family of zinc-finger transcription factors, including EGR1, EGR2 and EGR3, which are upregulated upon TCR activation. Similar expression was recorded in case of NAB2, a coactivator and corepressor of T-cell function (47). Also, transcriptions of NR4A2 and NR4A3, two members of the Nuclear receptor 4A (NR4A) family known for their important role during acute and chronic CD8⁺ T cell response (48), were highly expressed along the differentiation subsets. Moreover, stimulation with PMA/ionomycin induced the highest expression of both genes in T_term, followed by T_inter and T_n. In contrary, the highest expression of NR4A2 and NR4A3 in ConA-stimulated subsets was recorded in T_n. Recently it has been reported that NR4A3 increases early expression of transcription factors involved in the SLEC differentiation and its absence favors differentiation of MPEC and central memory CD8⁺ T cells (49). Notably, T_n and T_inter but not T_term expressed high levels of BCL2 upon PMA/ionomycin and ConA stimulation. These results are in accordance with the recent findings which show that naïve T cells highly express BCL2 and are more dependent on it for survival than effector and memory T cells (50). The transcription factor MYB promotes formation of stem-like memory cell and restrains terminal effector differentiation by inducing expression of BCL-2 and TCF7 as well as inhibition of ZEB2 (51). While T_n strongly expressed MYB following PMA/ionomycin and ConA stimulation, no upregulation was induced by T_inter or T_term. Also, expression of BACH2, described as transcriptional repressor of terminal differentiation that restrains formation of short-lived effector cells (52, 53), was upregulated in T_n and T_inter but not T_term after PMA/ionomycin stimulation. Furthermore, the PMA/ionomycin stimulation induced the expression of ZEB1 and TCF3 in T_inter and T_term, respectively. Studies in mice showed that STAT1 and STAT4 are important transcription factors for the clonal expansion and promotion of antigen-activated CD8⁺ T cells. Whereas STAT1 effects type I IFN-dependent clonal expansion of CD8⁺ T cells, STAT4 contributes to proliferation and effector maturation of CD8⁺ T cells triggered by IL-12-mediated signaling (54, 55). In fact, the PMA/ionomycin stimulation induced high expression of STAT1 in all three CD8⁺ T-cell subsets, while the ConA stimulation induced the upregulation only in the T_inter. In case of STAT4 expression, we observed upregulation in the T_inter stimulated with ConA and T_term stimulated with PMA/ionomycin. Notably, while BCL6 was highly expressed in T_n and T_term following PMA/ionomycin stimulation, expression of ID2, a transcriptional regulator upregulated by activated CD8⁺ T cells late in effector phase which can also influence their differentiation into memory cells (56), was upregulated in all three CD8⁺ T-cell subsets stimulated with PMA/ionomycin. Moreover, expression of ID2 was also increased in the T_n following ConA stimulation. Transcription of KLF9 in the T_n was increased with ConA and PMA/ionomycin stimulation as well as in PMA/ionomycin-stimulated T_inter.

Looking at co-stimulatory genes, we found that expression of CD27 (TNFRSF7), expressed mostly on naïve T cells and also required for T-cell memory in mice (35), was upregulated only in T_inter stimulated with PMA/ionomycin ( Figure 4B ). Also, another co-stimulatory gene CD28, which is absent from human effector CTLs (14), was upregulated in T_n and T_inter stimulated with PMA/ionomycin as well as in T_n stimulated with ConA. Furthermore, upon PMA/ionomycin stimulation all three CD8⁺ T-cell subsets expressed ITGAL (CD11a), a β2 integrin reported to be important for homing of T cells and generation of antigen-specific T cells (9). In all three CD8⁺ T-cell subsets PMA/ionomycin stimulation induced high expression of CD40LG, a member of the tumor necrosis factor superfamily transiently expressed on activated CD8⁺ T cells that promotes expansion and differentiation in a cell-extrinsic manner (57), and CD83, a member of the immunoglobulin superfamily. Expression of CD69, an early activation marker, was highly expressed in all three CD8⁺ T-cell subsets stimulated with PMA/ionomycin, and to a much smaller extent in T_n and T_inter upon ConA stimulation. These differences in transcripts of CD69 concerning different stimulations can be explained with the fact that the CD69 expression is upregulated already after 30 to 60 minutes after activation and declines promptly after 4-6 hours (58). Furthermore, transcription of the inducible T cell co-stimulator (ICOS), a member of the immunoglobulin family structurally close to CD28 and rapidly expressed on activated T cells (59), was highly increased in all three CD8⁺ T-cell subsets following PMA/ionomycin and to a lesser extent in T_n and T_inter after ConA stimulation. In addition, all three CTL subsets stimulated with PMA/ionomycin showed high increase of the lymphotoxin alpha (LTA), described to positively affect antigen-specific T-cell response during an acute LCMV infection through increase of IFN-γ production (60). Two members of the signaling lymphocytic activation molecule family (SLAMF), namely, SLAMF1 and SLAMF6, were induced in T_n and T_inter after PMA/ionomycin stimulation. Highest expression of SLAMF1 and SLAMF6 has been reported on central memory and effector memory subsets of CD8⁺ T cells (61).

In case of co-inhibitory genes, known to inhibit T-cell activation, cytolytic function and cytokine production (62), we observed that expression of PDCD1 (PD-1) was induced in all three CD8⁺ T-cell subsets stimulated with PMA/ionomycin, whereas its ligand CD274 (PD-L1) was only expressed on T_inter. Moreover, both T_n and T_inter showed upregulation of cytotoxic T lymphocyte antigen-4 (CTLA4) upon stimulations, with PMA/ionomycin stimulation inducing stronger expression. It has been shown that CTLA4 is closely related to CD28, binds to the same ligands (CD80 and CD86) and inhibits T cell response (63). Next, we found that expression of lymphocyte activation gene-3 (LAG3) was induced in all three CD8⁺ T-cell subsets after both stimulations. This is in accordance with previous research in mice suggesting that naïve CD8⁺ T cell show low expression of LAG3, but increase its expression in response to stimulation (64). Notably, expression of HAVCR2, which encodes T-cell immunoglobulin and mucin domain-3 (Tim-3) inhibitory molecule, was upregulated only in T_inter and T_term. The tumor necrosis factor superfamily (TNFSF) and its corresponding receptor superfamily (TNFRSF) were differently expressed among porcine CD8⁺ T-cell subsets.

In case of TNFSFs, PMA/ionomycin stimulation induced expression of TNFSF8 (CD30L) and TNFSF11 (RANKL) in T_n and T_inter, whereas transcript of TNFSF14 (LIGHT) was upregulated in T_inter and T_term. The highest expression of TNFSF9 (4-1BBL) could be observed in T_term, followed by T_inter and T_n. In addition, the ConA stimulation induced its expression only in T_term. Also, only T_term showed increased upregulation of TNFSF10 (TRAIL) upon PMA/ionomycin stimulation. Regarding TNFRSFs, transcript of TNFRSF1A (TNFR1) was enriched in T_term, while TNFRSF1B (TNFR2) was expressed in all three CD8⁺ T-cell subsets following PMA/ionomycin stimulation. Both T_n and T_inter expressed TNFRSF6B (DCR3) and TNFRSF25 (DR3) after PMA/ionomycin stimulations. Notably, expression of TNFRSF9 (4-1BB) was strongly induced in all CD8⁺ T-cell subsets upon both stimulations. Similarly, the PMA/ionomycin stimulation induced high expression of TNFRSF18 (GITR) and TNFRSF4 (OX40), an intermediate activation marker, in all CD8⁺ T-cell subsets, while stimulation with ConA induced the upregulation of these two genes in T_inter but only TNFRSF18 in T_term.

Genes associated with effector functions of CTLs were primarily highly expressed by T_term, followed by T_inter and in just few cases by T_n following PMA/ionomycin stimulation ( Figure 4C ). Moreover, the ConA stimulation had almost no effect on upregulation of those genes in CD8⁺ T-cell subsets. Several genes linked to cytolytic activity, including GZMA1, PRF1 (Perforin), FASLG, JUN, MCL1 and HSP90B1, were upregulated only in T_inter and T_term following PMA/ionomycin stimulation. Also, the highest expression of genes belonging to Jun (JUN, JUNB) and Fos (FOS, FOSB) families was detected in the T_term, followed by T_inter and T_n. Similarly, expression of several other genes involved in effector function and apoptosis, including BIRC3, CFLAR, CRTAM, HSPA9, NFATC1, NFKB1 and NFKBIE was induced among CD8⁺ T-cell subsets. Notably, only T_n and T_inter showed upregulation of BBC3, HSPD1 and HSPE1. While killer cell lectin like receptor k1 (KLRK1) gene was upregulated in the T_term, transcript of FADD was induced only in T_n. Compared to T_inter and T_term, we found higher expression of GZMA2 and PMAIP in T_n stimulated with PMA/ionomycin. Interestingly, both ConA and PMA/ionomycin stimulation induced upregulation of BCL2A1 in all CD8⁺ T-cell subsets, although with markedly stronger expression after PMA/ionomycin stimulation and in T_term. What is more, the highest expression of TRAF1 and RELB was observed in the T_term, followed by T_inter and T_n. In case of T_term, expression of TRAF1 and RELB was also induced by the ConA stimulation, whereas T_inter showed only the upregulation of RELB transcript. Expression of another IkB family gene linked to apoptosis (65), namely the NFKBIA, was highest in the T_inter stimulated with PMA/ionomycin.

To extend the understanding of the immunological roles and functions of genes across different ex vivo CD8⁺ T-cell subsets, we performed GO term enrichment analysis. For GO terms related to immune system, most of the upregulated DEGs in the T_term were assigned to lymphocyte activation involved in immune response (42.39%) ( Figure 5A ). In contrast, upregulated DEGs in the T_n were mostly associated with T-cell differentiation (42.31%), T-cell receptor signaling pathway (23.08%) and V(D)J recombination (11.54%). The majority of upregulated DEGs in the T_inter compared to naïve subsets were related to T-cell differentiation involved in immune response (27.69%), T-cell cytokine production (27.69%) and alpha-beta T-cell differentiation (23.08%). On the other side, upregulated DEGs in the T_inter compared to T_term were mainly linked to the regulation of T-cell differentiation (90.0%). When compared to T_inter, most of the upregulated DEGs in T_n were enriched for V(D)J recombination (55.17%), regulation of T-cell receptor signaling pathway (17.24%) and T-cell differentiation (17.24%), whereas within the T_term they were associated with regulation of lymphocyte differentiation (70.0%) ( Supplementary Table S2 , and Figure S2 ).

KEGG pathway analysis revealed that upregulated DEGs in the T_term compared to T_n were assigned to 272 pathways, including 21 pathways related to the immune system. Within immune-related pathways, the highest number of upregulated DEGs were enriched in chemokine signaling and T-cell receptor signaling pathways. Additionally, upregulated DEGs were linked to metabolic, MAPK signaling and cytokine-cytokine receptor interaction pathways. Compared to T_term, DEGs within T_n were associated with 278 pathways, with 20 immune-related pathways, as well as metabolic and MAPK signaling pathways. When compared to T_inter, DEGs of T_n and T_term were enriched in 192 and 167 pathways, respectively. Furthermore, DEGs of both T_n and T_term were enriched in 15 immune-related pathways. In comparison to T_n and T_term, DEGs within T_inter were enriched in 268 and 235 pathways, respectively. All pathways including corresponding genes retrieved through KEGG pathway analysis were listed in Supplementary Table S3 .

To further explore DEGs of in vitro stimulated CD8⁺ T-cell subsets, we conducted GO term enrichment analysis for immune system processes using previously mentioned bioinformatic software and plug-in package. The results demonstrated that the DEGs of PMA/ionomycin-stimulated T_n compared to medium control were mostly associated with the regulation of T-cell activation (60.23%) ( Figure 5B ). In contrast, the DEGs in PMA/ionomycin-stimulated T_term were differently linked to the leukocyte differentiation (55.77%). For the DEGs of PMA/ionomycin-stimulated T_inter, we found enrichment in GO terms associated with the leukocyte differentiation (31.37%) and the regulation of lymphocyte activation (22.14%). Next, we investigated GO terms for immune processes in different CD8⁺ T-cell subsets stimulated with ConA, showing that the lymphocyte differentiation and the regulation of T-cell activation were more related with the T_n upon stimulation, while response to interferon-gamma term was typically associated with T_inter and T_term ( Supplementary Table S2 and Supplementary Figure S2 ).

We next performed KEGG pathway analysis of each CD8⁺ T-cell subset upon PMA/ionomycin and ConA stimulation as described above. Upregulated DEGs within PMA/ionomycin-stimulated T_n were assigned to 302 pathways in KEGG pathway database. Similar observations were made with T_inter (318 pathways) and T_term (305 pathways) CD8⁺ T-cell subsets. Although similar number of pathways related to immune system were observed among PMA/ionomycin-stimulated CD8⁺ T-cell subsets, a much higher number of DEGs was found from T_inter and T_term than T_n. Furthermore, T_inter showed the highest number of DEGs enriched in T-cell receptor signaling pathway (36), followed by T_term (33) and T_n (25). Looking at the chemokine signaling pathway, T_inter and T_term showed same number of DEGs enriched in the pathway (26), whereas T_n had only 18 DEGs involved. Besides immune-related pathways, DEGs from all PMA/ionomycin-stimulated CD8⁺ T-cell subsets were highly enriched in metabolic, MAPK-signaling and cytokine-cytokine receptor interaction pathways. Interestingly, number of DEGs enriched in MAPK-signaling pathway was gradually increased along the CD8⁺ T-cell subsets. Based on upregulated DEGs upon ConA stimulation, we recorded 175 pathways in case of T_n, 120 pathways for T_inter and only 43 pathways for T_term. Also, upregulated DEGs within the T_n showed the highest number of immune-related pathways. In contrast to PMA/ionomycin-stimulated CD8⁺ T-cell subsets, ConA stimulation induced a limited number of genes associated with metabolic and MAPK-signaling pathways. Moreover, in all CD8⁺ T-cell subsets only few DEGs were enriched in T-cell receptor signaling, chemokine signaling and cytokine-cytokine receptor interaction pathways after ConA stimulation ( Supplementary Table S3 ).