Abstract Single-cell RNA sequencing (scRNAseq) has accelerated characterizing cellular phenotypes in pigs under healthy and diseased conditions. To pair scRNAseq with immune receptor profiling, we developed porcine-specific T cell receptor (TCR) and B cell receptor (BCR) enrichment primers that are compatible with the 10 × Genomics VDJ sequencing protocol. Using these assays, we profiled the immune repertoire of cryopreserved lung cells from CD1D-expressing and CD1D-deficient pigs after one or two infections with influenza A virus (IAV) to examine whether natural killer T (NKT) cells influence pulmonary TCR and BCR receptor repertoires. We also profiled T cells longitudinally sampled from the lung fluid of IAV-vaccinated and -infected pigs to track clonal expansion. While all pigs presented highly diverse repertoires, pigs re-exposed to IAV had more expanded T cell clonotypes with activated phenotypes, suggesting potential IAV-reactive clones. Our results demonstrate the utility of high throughput single cell TCR and BCR sequencing in pigs. Subject terms: Adaptive immunity, High-throughput screening __________________________________________________________________ We developed porcine-specific primers for TCR and BCR, enabling paired single-cell RNA and receptor sequencing of pig T and B cells. Using this assay, we tracked treatment-induced changes in TCR and BCR repertoires in pig lung tissue. Introduction High throughput single-cell RNA sequencing (scRNAseq) technology has greatly increased our understanding of the phenotypic diversity and plasticity of immune cell types as well as their cellular interactions in complex tissues, in both healthy and disease states^[43]1,[44]2. In addition to revealing cellular heterogeneity through RNA expression, scRNAseq can be coupled with additional assays to enhance cellular phenotyping, such as enrichment of T cell receptor (TCR) and B cell receptor (BCR) repertoires using primers that target the C regions in mRNA transcripts of TCR and BCR chains and isotypes. This allows construction of expressed TCRs and BCRs of individual cells including genomic rearrangements between the variable (V), diversity (D), and joining (J) regions of the TCR and BCR intervals responsible for creating diversity in receptor binding surfaces. Pairing TCR/BCR sequencing with scRNAseq provides a powerful approach for studying the relationship between immune repertoire and many types of immune responses. Additionally, V(D)J recombination at the TCR and BCR loci can be used as endogenous barcodes to trace T and B cell clonotypes as they expand or transition through different states, including within the same individual over time^[45]3,[46]4. Domestic pigs (Sus scrofa) are an important agricultural species that are intensively farmed making them vulnerable to many infectious pathogens. Therefore, a thorough understanding of the porcine immune system is needed to optimize vaccine and drug design and to identify immune targets for increased disease resistance through selective breeding and genetic engineering. Because of their many similarities to humans, swine are increasingly used in place of non-human primate models^[47]5,[48]6. However, a limitation that prevents fully exploiting these pig models is our incomplete understanding of the porcine immune system, which is due in part to a scarcity of immune profiling reagents for pigs. scRNAseq which does not require marker-based sorting of cell subsets is helping to address this gap. In the current work, we developed porcine-specific TCR and BCR primers compatible with the droplet-based protocols of the 10 × Genomics Next GEM Single Cell 5’ sequencing protocol. These assays were used to compare cryopreserved lung cells from pigs genetically engineered to lack CD1D, which encodes the antigen presenting molecule CD1d that is required for the development of natural killer T (NKT) cells, a subset of innate-like T cell that accumulates in barrier organs such as the lungs^[49]7–[50]12. Pigs were analyzed after one or two infections with influenza A virus (IAV) to interrogate the evolution of the TCR and BCR repertoires after primary or secondary infection and to determine if NKT cells exert T helper cell functions that influence receptor diversity. Additionally, we applied our protocol to a longitudinal assessment of T cells recovered from the lung lavage fluid of infant pigs exposed to IAV vaccination and infection. Here, we were able to track clonal expansion and monitor changes in the frequency of clonotypes within the same pig. Together, these results demonstrate the potential of TCR/BCR profiling to better understand a wide variety T and B cell-related immune responses in pigs. Results Single-cell RNA sequencing analysis of influenza virus-infected genetically edited pigs Cryopreserved cells liberated from the enzyme digested lungs of mixed-breed pigs carrying an inactive form of the CD1D (CD1D−/−) gene and littermates carrying one inactivated copy (CD1D−/+) were subjected to paired scRNAseq and scTCR/BCRseq using our custom primer sets. One cohort of pigs was necropsied five days after a single infection (1X) with H1N1 A/Missouri/CS20N08/2020 (MO20) (Fig.[51]1A and Table [52]S1). To compare heterologous adaptive immune responses, a second cohort of pigs was infected with MO20 virus two weeks after an initial infection with H1N1 A/California/04/2009 (pdmH1N1) and necropsied 5 days later (2X). While 1X pigs presented an increase in body temperature and shed more virus compared to 2X pigs, there was no difference in how the different genotypes responded to infection (Fig. [53]S1). Single-cell sequencing was performed on 12 pigs (3 per group), totaling 45,850 cells. A dimensionality reduction analysis identified 28 clusters by Uniform Manifold Approximation and Projection (UMAP) that we annotated according to established lineage markers (Figs. [54]1B, [55]S2, [56]S3A). Clusters 1 and 2 are characterized by the naïve T cells signatures (SELL, CCR7, LEF1, and TCF7), as well as the memory-associated and tissue residency marker CD69, suggesting their potential for tissue residency, as reported in human non-lymphoid tissues^[57]13–[58]16. The frequencies of CD8^+ tissue resident memory T cells (TRM – cluster 4) and cycling T cells (cluster 16) were significantly higher in 2X CD1D−/+ pigs than 1X CD1D−/+ pigs (Figs. [59]1C, [60]S3B). The cytotoxic CD8^+ T cell frequency (cluster 5) was also higher in 2X CD1D−/+ pigs, but the difference was not significant. The same cell types were increased in 2X CD1D−/− compared to 1X CD1D−/− pigs, but the differences were smaller than in CD1D−/+ pigs and not significant. Fig. 1. Single-cell transcriptomic analysis of IAV-infected CD1D−/− and CD1D−/+ pig lungs. [61]Fig. 1 [62]Open in a new tab A Overview of experiment setup. 3 CD1D−/− and 3 CD1D−/+ pigs were infected with pdmH1N1 IAV. Fourteen days later, the same 6 pigs were infected with H1N1 A/Missouri/CS20N08/2020 (MO20) IAV (designated 2X pigs) along with an additional 3 CD1D−/− and 3 CD1D−/+ pigs (designated 1X pigs). Necropsies were performed 5 days after the MO20 infection to collect lung tissue for single-cell immune profiling. Created with BioRender. B Uniform manifold approximation and projection (UMAP) visualization of lung leukocyte populations colored by cell clusters. Clusters were identified using the graph-based Louvain algorithm at a resolution of 0.5. C The average frequency of each cell type is presented for each treatment. D Bar graphs displaying the number of upregulated and downregulated differentially expressed genes (DEGs) in 2X compared to 1X CD1D−/+ and CD1D−/− pigs. E Bar graphs displaying the number of DEGs in CD1D−/+ versus CD1D−/− pigs after 1X or 2X infections. The DEGs were selected using a threshold of adjusted p-value < 0.05 (Data file [63]1). F Ingenuity Pathway Analysis (IPA) using 2X versus 1X NK cell (cluster 11) DEGs in CD1D−/+ and CD1D−/− pigs. G IPA analysis using CD1D−/+ versus CD1D−/− monocyte/macrophage (cluster 17) DEGs in 1X and 2X infected pigs. A z score of −2.0 < Z  >  2.0 and p < 0.05 were considered significant. The y-axis displays the top 20 pathways. The dot size displays significance [−log10 (P value)] of gene sets. Dot color saturation represents the z score. Data file [64]1 contains a complete list of significantly enriched cluster 11 and 17 IPA pathways. Next, we compared differentially expressed genes (DEGs) within individual immune cell types between 2X and 1X pigs by genotype (Fig.[65]1D, Data file [66]1) or between CD1D−/+ and CD1D−/− pigs by number of infections (Fig. [67]1E, Data file [68]1). Overall, we detected substantially more upregulated genes in 2X versus 1X pigs in the CD1D−/+ group (3241 genes) compared to their CD1D−/− counterparts (1162 genes), especially within cytotoxic T cells (cluster 5), CD2^- γδ T cells (cluster 9), NK cells (clusters 11 and 12), and B cells (cluster 14). This suggests that stimuli originating from NKT cells during a primary influenza infection prime multiple lung immune cell types to react more robustly upon re-infection. For example, an Ingenuity Pathway Analysis (IPA) of DEGs from NK cells (cluster 11) revealed that, compared to 2X CD1D−/− NK cells, 2X CD1D−/+ NK cells upregulated pathways related to C-type lectin receptors, PTEN regulation, interleukin-1 family signaling, and nuclear factor-κB (NF-kB) pathway signaling (Fig.[69]1F, Data file [70]1). Within genotype (CD1D−/+ versus CD1D−/−), monocytes and macrophages in cluster 17 had a markedly higher number of DEGs in CD1D−/+ compared to CD1D−/− pigs under both 1X (558 DEGs) and 2X (282 DEGs) infection conditions. An IPA of cluster 17 DEGs revealed that, compared to 1X CD1D−/− pigs, 1X CD1D−/+ monocytes and macrophages upregulated a diverse array of biological pathways, including pathways related to interferon alpha/beta signaling, interleukin-1 family signaling, NF-κB pathway signaling, and class I major histocompatibility complex (MHC) mediated antigen processing and presentation (Fig.[71]1G, Data file [72]1). This agrees with prior reports that activated NKT cells induce the downstream maturation of antigen presenting cells^[73]17. A subset of these pathways was upregulated by 2X CD1D−/+ versus 2X CD1D−/− monocytes and macrophages along with other pathways that were unique to the 2X group, including several cell cycling-associated pathways. TCR repertoire of lung T lymphocytes We analyzed the pig lung tissue cells using our pig-specific V(D)J primers designed for compatibility with the 10 × Genomics Chromium Next GEM Single-cell 5’ kit (Table [74]S2). Over the 12 samples, we obtained 927,503,456 sequence reads with an average of 77,291,955 reads per library (Data file [75]2). A de novo assembly of raw sequencing reads produced unannotated porcine V(D)J contigs which were identified using the sequencing primers and thereafter aligned to cells in the gene expression clusters in Fig. [76]1. Across all samples, approximately 60% and 70% of identified TRA or TRB contigs were respectively aligned to single cells (Fig. [77]2A, Data file [78]2). Assembled VDJ sequences were blasted against the international ImMunoGeneTics (IMGT) germline TRBV, TRBD, TRBJ, and TRAJ databases^[79]18. Approximately 70% of TRB contigs in a pairing with TRA contigs mapped to an annotated TRVB or TRJB gene (Data file 2). Because Vα genes are not annotated in IMGT, Vα sequences were assigned according to pig TRAV sequences from our previous publication^[80]19 (Fig. [81]S4, Data file [82]2). Only cells with IMGT-annotated TRB genes that were paired with productive TRA were used for further analysis (Fig. [83]2B). These cells accounted for between 23% and 38% of cells within αβ T cell clusters (clusters 1-6) across samples. Fig. 2. Characterization of T cell clonotypes from lung tissue. [84]Fig. 2 [85]Open in a new tab A UMAP plot of cells expressing one (tra, trb) or both (tratrb) TCR α and TCR β chains. B UMAP plot of cells with paired TCRα and TCRβ chains that were used for downstream analysis. C Relationship between TRBV and TRBJ usage in T cell receptor rearrangements by treatment. Cell barcode counts for TRBV and TRBJ gene segments were normalized by cell numbers across treatments and scaled by TRBV segments. D Principal component analysis of TRBV and TRBJ gene usages by individual pig. E, F Heatmaps showing row-scaled mean expression of the 10 highest differentially expressed SLA class I (E) and SLA class II (F) genes per pig. G Heatmap of overlapping clonotypes defined by identical CDR3β region sequences between pigs. H Proportion of clonotypes by treatment with 1, 2, ≥ 3 cells per clonotype. I Abundance of clonotypes by cluster with 1, 2, ≥ 3 cells per clonotype in the combined dataset. J Expression of naive and activation T cell markers in expanded and non-expanded clonotypes in the combined dataset. K UMAP plot displaying the five most expanded clonotypes in the combined dataset. L, M Number of cells in each of the five most expanded clonotypes by treatment (L) and cluster (M). The percentage values displayed in the (M) figure legend represent the proportion of each clonotype relative to the total number of TCR+ cells used in the analysis, as shown in (B). See Fig. [86]1B for cluster annotations for (I) and (M). The expression of Vβ/Jβ and Vα/Jα combinations was analyzed to determine if favored rearrangements were correlated with CD1D genotype or the number of IAV infections (Figs. [87]2C and Fig. [88]S4). Overall, TRBJ1-1*01, TRBJ1-2*01, TRBJ3-3*01, TRBJ2-3*01, and TRBJ2-6*01 were used in many rearrangements, while pairings with TRBJ1-5*01, TRBJ2-5*01, TRBJ3-7*01, TRBJ1-6*01, and TRBJ2-1*01 were relatively rare (Fig. [89]2C). This is similar to our previous analysis of pig peripheral blood T cells that used a bulk RNA sequencing approach^[90]19. Although very rare, only CD1D-expressing CD1D−/+ pigs had cells bearing the invariant NKT cell TCR α-chain composed of pTRAV10 and TRAJ18*01^[91]19 (Fig. [92]S5, Data file [93]2). We observed that 2X CD1D−/+ pigs had Vβ/Jβ recombinations that were more similar to each other than to other pigs (Fig. [94]2D). This was also the case for pigs s303, s305, s307, and s308, which were from the same litter. Since structurally rearranged TCRs are selected by major histocompatibility complex (MHC) molecules, we analyzed swine leukocyte antigen (SLA) class I and II expression in individual pigs by aligning scRNAseq transcripts across clusters (Fig. [95]S6) to the IPD-MHC database that provides a repository of MHC sequences for a number of species, including swine^[96]18. To visualize SLA usage, a row-scaled mean expression of the 10 highest class I (Fig. [97]2E) and class II (Fig. [98]2F) genes in each pig was performed. Taking the expression of both SLA classes together, we could group pigs into 6 distinct MHC haplotype inheritance patterns (designated A-F) as follows: A – s403, s406, s408; B – s305, s307, s308, s407, s412; C – s303; D – s402; E – s411; F – s404 (Table [99]S1). Vβ/Jβ pairings in the five “B” inheritance pattern pigs clustered together (Fig. [100]2D), suggesting a link between MHC inheritance and VDJ selection. Next, since TRB gene segments are fully annotated in IMGT whereas TRAV segments are not, we used identical CDR3β region sequences to analyze our samples for expanded clonotypes (Data file [101]2). Some clones were found in more than one pig, with some of the highest numbers of shared clones among the 2X infected pigs (Fig. [102]2G). Additionally, 2X pigs, but especially the CD1D−/+ group, had more expanded clones than the 1X groups (Fig. [103]2H). The highest concentrations of clonally expanded T cells were among CD4^+ and CD8^+ TRMs, cytotoxic CD8^+ T cells, and cycling T cells while comparatively few expanded clones were present among naïve-like T cells in clusters 1 and 2 (Fig. [104]2I). Expanded clones were enriched for immune activation/effector related genes compared to unexpanded clones (Fig. [105]2J). 2X CD1D−/+ pigs contributed four of the five most expanded clonotypes, all of which included more than twelve cells (Figs. [106]2K, [107]L). One of these clonotypes originated from CD4^+ TRMs while the remaining four were from CD8^+ TRMs (Fig. [108]2M). Interestingly, 2X pigs, but especially the CD1D−/+ group, had shorter average CDR3β lengths than the 1X groups (Fig. [109]S7). This phenomenon has been associated with converging TCR motif signatures in a number of diseases^[110]20,[111]21. BCR repertoire of lung B lymphocytes Using a similar approach to our TCR profiling, we analyzed the lung tissue cells using our pig-specific primers for the IGK, IGL, IGHM/IGHD, IGHA, and IGHG genes (Fig. [112]3A, Table [113]S2), which respectively encode the immunoglobulin κ and λ light chains and IgM/IgD, IgA, and IgG heavy chains. Over the 12 samples, we obtained 1,395,728,419 sequence reads with an average of 11,631,0701 reads per library (Data file [114]2). Assembled V(D)J sequences were blasted against the IMGT germline IGL, IGK, and IGH databases^[115]18. Only cells with IMGT-annotated IGH genes that were paired with productive IGL or IGK chains were used for further analysis. Pig IGHM and IGHD CH1 sequences are almost identical making it difficult to design chain-specific primers^[116]22. Instead, we designed a primer for a CH1 region common to both chains and used chain-specific single nucleotide polymorphisms to distinguish IGHM and IGHD contigs (Fig. [117]S8). The percentage of B cells expressing different light and heavy chain contigs are as follows: 21% IGK, 40% IGL, 45% IGHM, 7% IGHD, 5% IGHG, and 3% IGHA. These cells accounted for between 53% and 70% of B cells among the 12 samples. Fig. 3. B cell receptor repertoire profiling of lung tissue B cells. [118]Fig. 3 [119]Open in a new tab A UMAP plots of cells expressing IGK and IGL chains and IGHM, IGHD, IGHG, and IGHA chains. B Proportion of cells expressing IGK and IGL in each B cell cluster. C Percentage of cells expressing IGHM, IGHD, IGHG, and IGHA in each B cell cluster. D Percentage of IGHG+ cells expressing different porcine IGHG subclasses in each B cell cluster. E Number of cells expressing light and heavy chain V(D)J gene segments. F Percentage of clonotypes by cluster with 1, 2, ≥ 3 cells per clonotype. G Percentage of clonotypes by IGH chain with 1, 2, ≥ 3 cells per clonotype. H Proportion of cells expressing IGHM, IGHD, IGHG, and IGHA by treatment. I Percentage of clonotypes by treatment with 1, 2, ≥ 3 cells per clonotype. J Percentage of IGHG+ cells expressing different porcine IGHG subclasses by treatment. K Length of light and heavy chain CDR3 sequences by treatment. The overall IGL to IGK ratio was 1.8:1 (Fig. [120]3B). Prior studies have reported that IGL:IGK usage in circulating immature pig B cells is ~1:1^[121]23–[122]25. Plasma B cells (cluster 15) had the highest proportions of IGHG+ and IGHA+ cells consistent with B cells that have undergone diversification for mucosal antibody secretion (Fig. [123]3C). Since our IGHG primers targeted a conserved constant region of the eight genomic IgG constant region gene sequences cataloged in IMGT (IGHG1, IGHG2, IGHG3, IGHG4, IGHG5-1, IGHG5-2, IGHG6-1, and IGHG6-2)^[124]26–[125]28, we mapped IGHG sequences to the IMGT database. IGHG1 dominated all three clusters (Fig. [126]3D). IGHG5-2 was also detected in every B cell subset. However, only antigen-experienced memory B cells and plasma cells in clusters 14 and 15 had detectable IGHG2, IGHG3, and IGHG5-1 subclasses. IGHG4 and IGHG6 could not be distinguished from the other subclasses. Next, we analyzed V(D)J gene segment usage (Figs. [127]3E, [128]S9, [129]S10). Unlike humans or laboratory rodents, pigs use a limited number of V, D, and J genes in VDJ heavy chains to form most of their BCR repertoire, with >90% of the VDJ repertoire made up of seven V segments, two D segments, and a single J segment^[130]29–[131]31. Consistent with these studies, we found six V segments (V1-15, V1S2, V1-4, V1-8, V1S5, and V1-6) accounted for ~70% of V usage (Figs. [132]3E, [133]S10) and that D and J usage was largely restricted to D1 or D2 and J5. We analyzed our samples for expanded clonotypes defined by the recombination of identical IG light and heavy chain CDR3 sequences (Data file [134]2). This revealed far fewer expanded B cells compared to the T cell compartment, with no single clonotype containing more than four cells. Almost all expanded clonotypes were among the small number of IgG or IgA secreting plasma cells (cluster 15) (Figs. [135]3F, [136]G). No striking differences in IGH usage (Fig. [137]3H) and clonotype expansion (Fig. [138]3I) were observed among the four treatment groups. However, the ratio of IGHG5-2/IGHG1 was greater in 2X than 1X pigs (Fig. [139]3J) while heavy chain CDR3 length tended to be greater in 1X than 2X pigs (Fig. [140]3K), which mirrors the TCR CDR3β length differences we observed between these two groups. Longitudinal assessment of T cells from lung lavage fluid Clonotype tracking is a powerful approach to monitor changes in the frequency of clonotypes of interest in cancer and vaccine immunology within a single individual^[141]32–[142]36. Using our swine TCR α and β chain primers, we analyzed T cells for V(D)J clonotypes and TCR α and β chain gene usage from the lungs of three specific pathogen free raised piglets (Fig. [143]4A). Two pigs (FLU1 and FLU2) were vaccinated with a combination of inactivated pdmH1N1 virus and adjuvant at 28 days of age and boosted 15 days later. Both pigs were infected with live pdmH1N virus 2 weeks after the booster. Lung fluid was collected 3 days before infection (T1) and 7 days after infection (T2). The third pig (NAIVE) was lavaged at 28 (T1) and 33 (T2) days of age. Between 942 and 4205 (average = 2469) CD3D+ cells, sorted from each lung fluid sample using an anti-CD3 antibody (Fig. [144]4B), were subjected to paired single-cell RNA and TCR sequencing. The combined dataset separated into 15 clusters (Fig. [145]4C, [146]D). Clusters 10 and 11 were γδ T cells that are relatively abundant in pigs. Clusters 1, 2, 3, 4, 5, 6, 7, 11, 12, 13, and 14 upregulated tissue residency markers, while clusters 8, 9, and 10 expressed circulating T cell markers (Figs. [147]4E, S[148]11). Based on their upregulation of CCR7, SELL, and CXCR3, cluster 8 cells are likely CD8^+ central memory (TCM) T cells (Fig. [149]S11), whereas cluster 9 cells are probably circulating naïve T cells. Tissue resident memory T cells were composed of both CD4^+ (cluster 1) and CD8^+ cells (clusters 2, 4, 5, 6, and 12) as well as proliferating cells that contained a mixture of CD4^+ and CD8^+ TRMs (clusters 13 and 14). In response to IAV infection, both FLU1 and FLU2 pigs presented an increase in the frequency of proliferating CD4^+ and CD8^+ TRMs, CD2^- γδ T cells, CD4^+ TRMs, naïve-like CD4^+ T cells, CD8αα T cells (clusters 6 and 7), and CD8^+ tissue effector memory (TEM) clusters, and a decrease in the frequency of CD8^+ TRMs (Fig. [150]4F). We identified 10 modules of co-regulated genes and regulatory networks across cell types (Fig. [151]4G, [152]H, Data file [153]2) among which Module 5 in clusters 7, 8, 9, and 10 harbored circulating T cell genes, Module 6 in clusters 2, 3, 4, 5, 6, and 11 harbored tissue residency and cytotoxic genes, and Module 10 in clusters 1, 2, and 11 harbored interferon stimulated genes (ISG). Fig. 4. Single-cell transcriptomic analysis of T cells isolated from lung lavage fluid. [154]Fig. 4 [155]Open in a new tab A Overview of experiment setup. T cells were FACS-sorted from the lung lavage fluid of 3 infant pigs: 2 pdmH1N1-vaccinated and -infected pigs (FLU1 and FLU2) and 1 naïve pig (NAÏVE). FLU pigs were sampled 3 days before IAV infection (T1) and 7 days after infection (T2). The NAÏVE pig was sampled at 28 (T1) and 33 (T2) days of age. Created with BioRender. B FACS plot showing acquisition of lung lavage T cells. C UMAP plot of the combined T cell datasets. D UMAP plots displaying individual samples. E Examples of genes used to identify resident (BHLHE40 and CXCR3) and circulating (SELL and S1PR1) T cells. F Proportions of T cell subsets in individual pigs at T1 and T2 timepoints. G Heatmap of 10 gene modules whose genes had a similar expression pattern across cell clusters. H UMAPs showing select genes from each module. Approximately 60% of αβ T cells recovered had paired TCRα and TCRβ chains (Fig. [156]5A). Sequence identity was used to map the phylogenetic relatedness of TCRα and TCRβ chains in individual samples (Fig. [157]5B). Many clones, identified by CDR3β sequences, were present at both T1 and T2 in the same animal, especially the FLU2 pig (Fig. [158]5C). A high proportion of expanded clones were CD8^+ TRMs (Fig. [159]5D). The two vaccinated and infected pigs presented the five most expanded clones, especially the FLU2 pig (Fig. [160]5E, [161]F). The low number of expanded clonotypes present in the NAÏVE pig samples is due partly to the lower number of cells we were able to collect from this pig. Next, we examined a larger collection of expanded clonotypes in FLU1 and FLU2 pigs, before and after infection, to identify potential influenza-reactive T cells (Fig. [162]5G). While the frequency of several clones decreased from T1 to T2, some increased, including CSAGERSNYEQIF, CASSVRSYPLNDLHF, CASSFGGVHTGQLYF, CAWSTTGTVTGQLYF, and CSAGEGGFGDTCFF. We also analyzed specificity groups within the CDR3β repertoire using immunarch^[163]37, a program that enables clustering of TCRs with an increased probability of sharing antigen specificity due to conserved motifs of CDR3 sequences. Two CDR3β 4-mer motif patterns (ASSL and SSLV) were found to be enriched in FLU1 and FLU2 pig samples (Fig. [164]5H). Interestingly, a few of the CDR3β sequences in expanded clones were identical or differed by a single amino acid from curated human TCR CDR3β sequences which recognize IAV epitopes in M1, NP, PB1, and PB2^[165]38. For example, the seventh most expanded CDR3β sequence in all pigs, ASSPGQGYEQ, matches a human CDR3β sequence that recognizes the immunodominant IAV Matrix protein 1 epitope GILGFVFTL when presented by human HLA-A*0201 (Fig. [166]5I)^[167]39,[168]40. Fig. 5. T cell clonotype tracking using VDJ recombination at the TRB locus. [169]Fig. 5 [170]Open in a new tab A UMAP plot of cells expressing one (tra, trb) or both (tratrb) TCR α and TCR β chains in the combined dataset. B Phylogenetic trees of TCRα and TCRβ sequences in the FLU1_T1 sample. C Heatmap of overlapping clonotypes between samples. D Abundance of clonotypes by cluster with 1, 2, ≥ 3 cells per clonotype in the combined dataset. See Fig. [171]4C for cluster annotations. E UMAP plot displaying the five most expanded clonotypes in the combined dataset. F Number of cells in each of the five most expanded clonotypes by sample. G Proportion of the most abundant clonotypes in FLU1 and FLU2 pigs by sample time. H Heatmap displaying the most abundant CDR3β 4-mer motifs in each pig normalized by cell numbers in each sample. I T cell clones expressing the CDR3β sequence ASSPGQGYEQ that matches a human CDR3β sequence recognizing the human HLA-A*0201 M1 epitope GILGFVFTL. Discussion The current study describes assays for single-cell TCR and BCR sequencing in pigs which presents a useful tool for enhancing effective vaccine and therapeutic design to protect swine health and to increase the potential of pigs as biomedical models for studying human immune-related physiological processes and diseases. We examined the utility of the assay in the setting of the pulmonary immune response against IAV infection since influenza is a respiratory pathogen of major importance for both humans and swine^[172]41. One set of samples was from cryopreserved lung sections of CD1D-expressing and -deficient pigs after one IAV infection or two infections with heterologous IAVs. This is of interest since NKT cell effector responses are important for anti-IAV immunity in mice, as evidence by reports that NKT cell-deficient mouse strains are significantly more susceptible to IAV infections than standard mice^[173]42–[174]45. Moreover, NKT cells have been shown capable of supplying a wide array of helper function reminiscent of CD4 T helper cells, which have the capacity to elicit wide-ranging cellular and humoral responses that can substantially boost the quality and durability of immune responses, including against heterologous and heterosubtypic IAV infections^[175]46. Among the downstream events regulated by activated NKT cells is the maturation of professional antigen presenting cells that subsequently induce various adaptive immune responses^[176]17,[177]47. This may account for why monocytes and macrophages in NKT sufficient CD1D−/+ pigs had a large number of upregulated genes and pathways compared to NKT cell deficient CD1D−/− pigs after both 1X and 2X IAV infections. Among the 12 samples in this dataset, we detected almost all VDJ gene segments annotated for TCR β chains, including some that were annotated as pseudogenes. In addition, most of the Vα genes detected overlapped with TCR α chains that we identified in a previous analysis of TCR chain usage using bulk RNAseq^[178]19. While we found that cells which expressed both α and β TCR chains were mostly within αβ T cell clusters (clusters 1, 2, 3, 4, and 5), we also detected cells expressing unpaired α and β chains in some non-αβ T cell clusters, such as CD2^- γδ T cells that expressed TCR β chains and NK cells and B cells that expressed TCR α chains. This is consistent with prior reports that a significant proportion of peripheral γδ T cells express TCR β chains^[179]48 and that NK cells express germline TCR transcripts^[180]49. Indeed, we found that virtually all NK cell TRA transcripts matched to TRAJ segments but that very few matched to TRAV segments. Examination of Vβ/Jβ combinations confirmed previous studies showing that certain rearrangements are favored in pigs^[181]19,[182]50. In several pigs, VDJ rearrangements clustered by MHC inheritance pattern. We also found that VDJ rearrangements preferred by 2X CD1D−/+ pig clonotypes clustered together which could mean that IAV infection influenced Vβ/Jβ selection in NKT cell-intact pigs. The most expanded clones were among CD4^+ and CD8^+ TRMs, cytotoxic CD8^+ T cells, and proliferating T cells, consistent with reports that these populations harbor antigen experienced T cells that are poised for rapid responses during an IAV infection^[183]51. As expected, 2X pigs had more expanded clones than 1X pigs. This was less apparent in CD1D−/− compared to CD1D−/+ pigs, which might suggest that induction of IAV-specific lung T cells is reduced in the absence of NKT cell helper functions. Lung tissue B cells consisted of naïve B cells, memory B cells, and a small population of plasma cells. Plasma cells displayed the highest diversity in immunoglobulin heavy chain usage and clonally expanded B cells. Consistent with prior reports^[184]23–[185]26,[186]28,[187]29,[188]31,[189]52,[190]53, we found (i) a mixture of IGL+ and IGK + B cells, (ii) IgG heavy chain usage was dominated by the IGHG1 subclass, and (iii) a limited number of IGHV, IGHD, and IGHJ genes formed most of the BCR repertoire. Pigs are interesting in so far as their antibody repertoire has evolved somewhat differently from mice and humans, including that they possess a highly streamlined IGH gene complex, which contains IGHV genes that all belong to a single ancestral IGHV3 family, and that only one IGHJ segment is functional^[191]26,[192]28,[193]29,[194]31,[195]52. Because of the small number of IGHV, IGHD, and IGHJ segments used, the combinatorial diversity in pigs is comprised of a mere ~14 possibilities compared to ~9,000 in humans^[196]31. T cell receptor transcript capture was more efficient in longitudinally collected lavage fluid T cell samples than in the CD1D lung tissue samples, probably owing to the former consisting entirely of FACS-purified T cells. As in lung tissue cells, most αβ T cell clusters exhibited paired α and β TCR chains, whereas a significant fraction of CD2^− γδ T cells expressed TCR β chains and CD2^+ γδ T cells and CD8αα T cells were enriched for TCR α chains. Clusters with the highest numbers of expanded clones were again tissue resident memory T cell populations. Our ability to detect a substantial number of the same T cell clones at two different timepoints in the same individual shows the utility of TCR profiling for monitoring changes in clonotypes of interest in pigs. The higher number of expanded clonotypes within FLU1 and FLU2 pig samples, some of which overlapped, compared to the NAÏVE pig, may be the result of IAV exposure by vaccination and infection leading to an increase in antigen experienced T cells in lung tissue and/or because fewer cells were collected in the NAÏVE pig. While there is a lack of reagents to identify influenza-specific T cell clones in commercial pig breeds, this information can to some extent be inferred by studying clonotypes in vaccinated pigs that increase after virus exposure. Using this approach, we identified CDR3 sequences that are potentially reactive to IAV antigens, including a clone with a CDR3β sequence identical to a human motif that recognizes an immunodominant epitope from Matrix protein 1^[197]54. This may arise because peptide binding motifs of some common swine SLA molecules partly overlap with the binding motifs of human HLA molecules^[198]54. In summary, the assays presented in this study can easily be applied to 5’ 10 × Genomics protocols for use in swine. Our protocols can be employed to profile TCRs and BCRs in the same sample which enhances the utility of the method as most adaptive immune responses involve both cellular and humoral responses. Accordingly, the combined protocol could shed light on acquired immunity that develops in response to vaccination and infection in production pigs, as well as in the growing number of immune-related pig models being developed for biomedical use. Limitations of the study Although this study demonstrates that single cell immune profiling can be performed in swine by customizing the 10 × Genomics VDJ sequencing protocol, conclusions on how an NKT cell deficiency as well as pre-existing immunity impacts TCR and BCR selection are based on a limited number of samples. In future, it will be important to expand our dataset to capture more individuals and timepoints. Additionally, tissue processing steps, including enzymatic digestion and tissue dispersion may introduce biases in cell recovery and viability, particularly affecting sensitive or rare immune subsets. Thus, validation of our findings is required using additional datasets as they become available. Methods Pigs The National Swine Resource and Research Center (NSRRC) at the University of Missouri bred a boar homozygous for a CD1D gene deletion with two sows heterozygous for the CD1D deletion that were full sisters to produce piglets that were homozygous for the CD1D deletion (CD1D−/−) as well as heterozygous segregants (CD1D−/+). Our CD1D breeding herd^[199]55 is on a commercial Large White crossbred background and maintained under specific pathogen free conditions. The CD1D genotypes of pigs were determined by PCR and flow cytometry^[200]55,[201]56. Piglets used for longitudinal assessment of T cells from lung lavage fluid were commercial Large White crossbred background pigs provided by the NSRRC. Virus infection and sample collection Eight CD1D−/− and eight CD1D−/+ pigs were transferred to biocontainment rooms at 4 weeks of age after being confirmed seronegative for IAV nucleoprotein antibodies by an ELISA developed by at the Veterinary Diagnostic Laboratory at Iowa State University. At day 0, 4 CD1D−/− and 3 CD1D−/+ pigs were intratracheally infected with 1 × 10^6 tissue culture infectious dose (TCID[50]) of H1N1 A/California/04/2009 (pdmH1N1) in 2 mL of DMEM (Gibco, Brooklyn, NY) after sedation with an intramuscularly delivered combination of midazolam, butorphanol, and xylazine. Fourteen days later, all 16 pigs were sedated and infected with 1 × 10^6 TCID[50] of H1N1 A/Missouri/CS20N08/2020 (MO20) IAV. Pigs were measured for clinical disease signs and daily nasal swab virus titers using previously established protocols^[202]57. Five days later, all pigs were sedated, euthanized by pentobarbital sodium intracardiac injections (70 mg/kg of body weight), and exsanguinated. At necropsy, the lungs were removed from the thoracic cavity for tissue collection. Cells were isolated for scRNAseq and receptor profiling from 3 randomly selected animals of each genotype that were infected once (1X) or twice (2X). Litter, sex, CD1D genotype, and MHC inheritance pattern of the piglets used for sequencing are described in Table [203]S1. For the study to collect T cells through successive lung lavages, two pigs (FLU1 and FLU2) were intramuscularly vaccinated with a combination of 1 × 10^6 TCID[50] of ultraviolet-inactivated pdmH1N1 virus and an oil-in water adjuvant (Emulsigen, 1:5 vaccine volume) at 28 days of age and boosted 15 days later. Both pigs were intratracheally infected with 1 × 10^6 TCID[50] live pdmH1N1 virus 2 weeks after the booster. Lung fluid was collected 3 days before and 7 days after infection. A third unvaccinated, uninfected pig (NAÏVE), housed in a separate room from infected pigs, was lavaged at 28 and 33 days of age. Pigs were intramuscularly sedated with midazolam, butorphanol, and xylazine to perform lung lavages and infections. Lavages involved inserting a size 10 French catheter attached to a syringe into the lung after which the lung was flushed twice with 5 ml of sterile saline solution. Recovered lung fluid was ejected into 10 ml phosphate buffered solution containing 10% fetal bovine serum (FBS). We have complied with all relevant ethical regulations for animal use. The studies were in accordance with the University of Missouri’s Institutional Animal Care and Use Committee (protocol number 34343) and Institutional Biosafety Committee (protocol number 17320). Tissue sampling and cell isolation Approximately 1 g of tissue collected from the left cranial, middle, and caudal lung lobes was combined and then digested with 2.5 mg/mL of Liberase TL (Roche, Indianapolis, IN) in Dulbecco’s Modified Eagle Medium (Thermo Fisher, Waltham, MA) at 37 °C for 45 min. The digested tissue was dispersed into single cells using a previously established protocol^[204]58, and then passed through a 70 μm cell strainer (Thermo Fisher, Waltham, MA). Cells were immediately cryopreserved in freezing media [90% FBS, 10% dimethylsulfoxide (DMSO)] in temperature-controlled freezing containers at 3 × 10^7 cells per/mL and stored at −80 °C until use. Samples were thawed in thawing media (RPMI –HyClone, Logan, UT, - 20% FBS), resuspended in RPMI with 10% FBS, filtered through a 40 μm filter (Bel-Art SP Scienceware, Wayne, NJ), and counted for viable cells using a Countess 3 automated cell counter (Thermo Fisher, Waltham, MA). To obtain T cells from lung washes, lavage fluid was filtered first through a 70 μm cell strainer (Thermo Fisher, Waltham, MA) and then a 40 μm pipette tip Flowmi® cell strainer (SP Scienceware, Warminster, PA), washed in RPMI with 10% FBS, counted, stained with PE-Cy7-conjugated anti-porcine CD3 antibody (clone BB23-8E6-8C8; BD Biosciences, San Jose, CA) and propidium iodide viability dye, and sorted for live CD3 positive cells using a BD FACSMelody Cell Sorter (BD Biosciences, San Jose, CA). Recovered live cells from cryopreserved lung tissue samples and lung lavage fluid were loaded onto the 10 × Chromium controller (10 × Genomics, Pleasanton, CA) for single-cell encapsulation. Primer design Pig-specific V(D)J primers were designed according to guidelines from the 10 × Genomics Chromium Next GEM Single-cell 5’ V2 user guide, which is well described in a recent publication on scTCRseq in dogs^[205]59. This protocol employs a nested PCR design which involves two rounds of V(D)J amplification using the same forward primer that primes off the 5’ Illumina adapter sequence, which is annealed to the 10 × barcoded cDNA during the conversion from mRNA. The first round of amplification uses the 5’ forward primer in combination with a 3’ outer reverse primer that matches the C region of the targeted chain. The second round uses the same forward primer in combination with a second reverse primer that primes the C region at an inner 5’ region from the outer reverse primer. To design the 3’ reverse primer sets, C-gene sequences were identified for TRA, TRB, IGH, and IGL and IGK transcripts. Inner and outer primers were respectively designed to target regions between 50–200 and 200–300 base pairs away from the 5’ region of the C region. Primer candidates were selected based on the percent of GC content, similar melting temperature range, and minimal predicted interaction with other regions in the pig genome. The final primer sequences as well as accession numbers for the DNA and rearranged mRNA transcripts used to identify the C regions are listed in Table [206]S2. Single cell processing Libraries were constructed by following the manufacturer’s protocol with reagents supplied in the Chromium Next GEM Single Cell 5′ Kit v2 (10 × Genomics). Briefly, cell suspension concentration and viability were measured with a Cellometer K2 (Nexcelom Biosciences, Lawrence, MA) stained with an acridine orange/propidium iodine dye mix (Invitrogen, Waltham, MA). Cell suspension combined with reverse transcription master mix was loaded on a Chromium Next GEM chip K along with gel beads and partitioning oil to generate gel emulsions (GEMs). GEMs were transferred to a PCR strip tube and reverse transcription performed on a Veriti thermal cycler (Applied Biosystems, Waltham, MA) at 53 °C for 45 min. All samples underwent 11 cycles of cDNA amplification upon which cDNA concentration and quality were assessed using a Fragment Analyzer 5200 (Agilent, Santa Clara, CA). For the gene expression library, up to 50 ng of the cDNA was fragmentated, end-repaired, A-tail added, and ligation of sequencing adaptors was performed according to manufacturer specifications. V(D)J libraries were constructed using pig-specific primers for two successive enrichments of TCR and BCR transcripts. The TCR and BCR assay used primer pools containing 1.43 μM per gene specific primer and 1.43 μM of the 10 × forward primer. Nested PCR amplification of the TCR and BCR sequences was performed using adapted mouse and human PCR protocols. This involved amplifying 2 μL of cDNA in 100 μL total reaction volume using 50 μL Amp Mix (10 × Genomics). Steps involved in the first reaction were 98 °C for 45 s for initial denaturation followed by 9 cycles of 98 °C for 20 s, 65 °C for 30 s, and 72 °C for 60 s. PCR product was purified using AxyPrep MagPCR Clean-up beads (Axygen, Union City, CA) and subjected to a second round of amplification using the same cycling conditions to the first except a total of 8 cycles was used. The amplicons were purified using AxyPrep MagPCR Clean-up beads (Axygen, Union City, CA). Libraries were constructed from PCR amplicons according to manufacturer specifications. The concentrations for all libraries were measured with the Qubit HS DNA kit (Invitrogen, Waltham, MA) and fragment sizes were determined on a Fragment Analyzer 5200 (Agilent, Santa Clara, CA). Libraries were pooled and sequenced on a NovaSeq 6000 (Illumina, San Diego, CA) to obtain paired end reads. The minimum sequencing depths targeted were 40,000 reads per cell for 5’ GEX libraries and 10,000 reads per cell for V(D)J libraries. Single-cell RNA sequencing data analysis The Sscrofa 11.1 genome assembly was used to align sequencing reads to generate gene matrix data by Cell Ranger (v8.0.0). Clustering analyses were performed using Seurat (v.4.4.0)^[207]60. To filter out low-quality genes and cells, only genes expressed in more than 3 cells and cells with more than 200 genes and less than 10% mitochondrial reads were included in the analysis. Afterwards, we followed a standard integration workflow to integrate samples. Briefly, transcript counts were log normalized, and the top 2,000 most variable genes in each dataset were identified using the FindVariableFeatures function. Then, the SelectIntegrationFeatures function was applied to genes that were consistently variable across datasets. Next, the FindIntegrationAnchors function identified a set of anchors between datasets using the top 30 dimensions from the canonical correlation analysis to specify the neighbor search space. Next, an integrated dataset was created by running the IntegrateData function. Then, the clustering analysis workflow was performed using RunPCA, FindNeighbours, FindClusters, and RunUMAP. Cell types were assigned based on the expression of known cell type markers (Fig. [208]S2). The differentially expressed genes (DEGs) between treatments in each cluster were computed using FindMarkers function with the Wilcoxon test (logfc.threshold = 0.1). The pathway enrichment analysis was performed using Ingenuity Pathway Analysis (IPA) (QIAGEN Inc.). The Core Analysis function was used to identify canonical pathways significantly associated with the DEGs using an adjusted p-value threshold of 0.05. Single-cell V(D)J data analysis Single-cell TCR sequencing reads were assembled into contigs using cellranger vdj (10 × Genomics) pipeline in denovo mode rather than reference-based mode due to the incomplete annotation of the pig germline Vα chain sequences. To identify the V(D)J chains, we searched assembled contigs against inner-enrichment primers using the usearch_global command. Primer matched TCR αβ and BCR IGH and IGL, and IGK chains were selected and integrated with the above cellular gene expression profiles. The TCR β chains with matched TCR α chains were selected for downstream analysis. TRBV, TRBD, TRBJ, and CDR3 sequences were mapped to the pig TRB reference in IMGT using the IMGT/V-QUEST sequence alignment software. Cells with duplicated or multiple contigs were removed. Immunarch (v1.0.0)^[209]37 was used to track clones across samples and identify k-mers from CDR3 sequencing. Scirpy (v.0.12.0) was used to analyze TCR β repertoires for clonal expansion (each unique CDR3) using the scirpy.pl.clonal_expansion command^[210]61. In addition, CDR3 sequences were searched within the Immune Epitope Database (IEDB)^[211]38 to find matches predicted to recognize epitope specificity using the TCRMatch T cell epitope prediction tool. The unique contigs of TCR α and β chains in each cell were aligned to each other using [212]CLUSTALW, after which phylogenetic trees were generated to infer the similarities among contigs using the package ape (v5.7-1)^[213]62. Using MMseqs2^[214]63, TRAJ sequences were mapped to the pig TRAJ reference in IMGT while TRAV segments were annotated according to TRAV sequences deposited in GeneBank^[215]64 that we previously named according to similar human TRAV genes^[216]19. Similarly, BCR repertoires were annotated using the available pig IGH and IGK/IGL reference in IMGT using IMGT/V-QUEST^[217]65. The reference mapped IGH chains and IGK/L chains were selected for downstream analysis. The IGH constant segments were mapped to IMGT/GENE-DB database. Scirpy was used to perform clonal (the recombination of IGH and IGL CDR3s) expansion, gene usage, and CDR3 sequencing length analyses. SLA alleles annotation from scRNA-seq data The Cellranger mkgtf was used to build the SLA alleles reference using the SLA FASTA files downloaded from the IPD-MHC database^[218]18. Afterwards, SLA reads were quantified using Cellranger count command, SLA counts were log normalized, and then the average expression of each SLA allele in each sample was computed with the AverageExpression function in Seurat. The data were visualized in heatmap using ComplexHeatmap (v2.25.1)^[219]66 and scaled for PCA analysis. Statistics and reproducibility Statistical analyses were performed using GraphPad Prism version 9.0.1 (GraphPad Software, La Jolla, CA). One-way analysis of variance (ANOVA) with Tukey’s post hoc test for multiple comparisons was used to assess differences in body temperature changes, virus titers, cell frequencies, and CDR3 lengths. A p-value of < 0.05 was considered statistically significant. The number of replicates and sample sizes are noted and explained in each related figure legend. Reporting summary Further information on research design is available in the [220]Nature Portfolio Reporting Summary linked to this article. Supplementary information [221]Supplementary Information^ (3.8MB, pdf) [222]42003_2025_8507_MOESM2_ESM.pdf^ (30.5KB, pdf) Description of Additional Supplementary Files [223]Supplementary Data 1^ (32.1MB, xls) [224]Supplementary Data 2^ (1.2MB, xlsx) [225]Reporting Summary^ (3.8MB, pdf) Acknowledgements