Single-cell analysis reveals lasting immunological consequences of influenza infection and respiratory immunisation in the pig lung

The pig is a natural host for influenza viruses and integrally involved in virus evolution through interspecies transmissions between humans and swine. Swine have many physiological, anatomical, and immunological similarities to humans, and are an excellent model for human influenza. Here, we employed single RNA-sequencing (scRNA-seq) and flow cytometry to characterize the major leucocyte subsets in bronchoalveolar lavage (BAL), twenty-one days after H1N1pdm09 infection or respiratory immunization with an adenoviral vector vaccine expressing haemagglutinin and nucleoprotein with or without IL-1β. Mapping scRNA-seq clusters from BAL onto those previously described in peripheral blood facilitated annotation and highlighted differences between tissue resident and circulating immune cells. ScRNA-seq data and functional assays revealed lasting impacts of immune challenge on BAL populations. First, mucosal administration of IL-1β reduced the number of functionally active Treg. Second, influenza infection upregulated IFI6 in BAL cells, decreasing their susceptibility to virus replication in vitro. Our data provides a reference map of porcine BAL cells and reveals lasting immunological consequences of influenza infection and respiratory immunisation in a highly relevant large animal model for respiratory virus infection. Author Summary Pigs and humans have a similar anatomy and physiology. In humans, cells from lung-washes are used to study immune responses and it was shown that these cells are crucial in protection against respiratory diseases such as influenza and COVID-19. To better understand lung immunity, we compared genes expressed in cells of pig lung-wash to white blood cells, providing an atlas for future studies of immunity in the lung. We also tested a vaccine given to the lung containing IL-1β, a strong immune activator that protects mice against influenza virus infection. However, although IL-1β increased pig immune responses it did not protect pigs against infection. We also showed that the number of immune cells that dampen immune responses (regulatory T cells) is reduced. In addition, we demonstrated increased expression of a protein, IFI6, 21 days after infection showing that while immune cells in the lung have common properties, the invading organisms influence them significantly. Our study elucidates why some vaccines fail despite inducing powerful immune responses, emphasizes the need for caution when applying results from small animals like mice to humans, and indicates the importance of the pig as a model to study disease in humans and livestock.


Introduction
Respiratory viruses pose significant global health threats, with influenza and coronaviruses responsible for major human epidemics and pandemics.The immune responses to these respiratory pathogens are initiated within the respiratory tract.Alveolar macrophages and dendritic cells continuously survey and sample the respiratory lumen, contributing to the induction of immune responses by secretion of cytokines and chemokines.These inflammatory mediators recruit innate cells including neutrophils and natural killer cells, which possess the ability to kill virus-infected cells.Adaptive immune responses initiated in the local draining lymph nodes are mediated by CD4 and CD8 T cells that enter the lung.These cells recognize relatively conserved internal viral proteins and can mediate cross protection against infections with serologically distinct virus strains (1,2).A proportion of these effector T cells remain in the lung to become tissue resident memory cells (TRM) cells (3,4).TRM can mount a rapid protective immune response on antigen encounter and have been shown to be critical in heterotypic cross protection against influenza in mice (5,6).In contrast, the best correlate of type specific adaptive immunity is antibody, and multiple lung-resident B cell (BRM) subsets exhibiting spatial and phenotypic diversity have been described (7)(8)(9).
Understanding how an immune response is induced and maintained in the lung is crucial for the design of more effective vaccines that provide long-lasting immunity.
Respiratory immune responses are highly localized, and direct investigation of lung cells and fluids, rather than assessment of correlates in blood, is essential (10).Mucosal immune responses following antigen exposure have been extensively characterized in mouse models of infection, helping to define phenotypic, functional and transcriptional features of lung immune cells (11)(12)(13)(14).However, investigating respiratory immune responses in humans is more challenging, although techniques allowing direct sampling of the respiratory tract are constantly improving.One common approach is the collection of bronchoalveolar lavage (BAL), containing cells and fluid from the airways and alveoli.BAL is enriched with TRM, and it has been shown that the presence of respiratory syncytial virus-specific CD8 TRM in the BAL correlates with reduced symptoms and viral load in human challenge studies (15).High dimensional phenotypic, transcriptomic, and functional profiling of immune responses in paired BAL and blood samples from COVID-19 patients revealed key roles for airway T cells, monocytes and macrophages associated with disease outcome (16).Single-cell transcriptome profiling (scRNA-seq) has enabled high resolution mapping of cellular heterogeneity, developmental trajectories, and activation status of BAL cells, revealing how these responses differ across species and between different infections (17)(18)(19)(20)(21)(22)(23)(24).These studies offer insight for the development of targeted therapies and vaccines.
In recent years it has become clear that targeting TRM is a promising vaccine strategy in inducing long term protection against respiratory viruses.Intranasal administration of a liveattenuated influenza vaccine generated long-term virus-specific TRM in the lungs of mice, which mediated cross protective heterotypic immunity in the absence of neutralizing antibodies, suggesting that for optimal induction of heterotypic immunity, virus infection of the lungs or immunisation of the respiratory tract is required (14,(25)(26)(27).Recent immunization strategies for achieving heterotypic immunity via TRM include pulmonary delivery of influenza vaccine candidates, targeting conserved internal virus proteins such as the nucleoprotein (NP).These strategies prevent infection, limit viral replication or reduce disease severity in the absence of cross-neutralising antibodies in pre-clinical models (5,(28)(29)(30).
Pigs are an important livestock species and serve as a significant biomedical model for studying human disease (31).Pigs have a longer life span, and are genetically, immunologically, physiologically, and anatomically more like humans than small laboratory animals (32,33).Pigs are also an important, natural, large animal host for influenza A viruses and are infected by the same subtypes of H1N1 and H3N2 viruses (34).Pigs can also be a source of new viruses potentially capable of initiating pandemics (35,36).Pigs exhibit similar clinical manifestations and pathogenesis to humans when infected with influenza viruses making them an excellent model to study immunity to influenza (37).This similarity extends to the lobar and bronchial anatomy, as well as the histological structure of pig lungs which closely resembles that in humans.Additionally, pigs exhibit a comparable distribution of sialic acid receptors in the respiratory tract (38).We have established a robust influenza challenge pig model to study immunity to influenza and evaluate the efficacy of novel vaccines and therapeutics (39)(40)(41)(42).We have further shown that pig immune responses following infection with H1N1pdm09 influenza virus are similar to those induced in humans (43).
Here we used scRNA-seq to analyse alveolar macrophages, CD4 T cells, CD8 T cells and B cells from the BAL of inbred Babraham pigs to define cellular states of homeostasis and activation following H1N1pdm09 (pH1N1) infection.We also characterized these cells following respiratory immunization with an adenoviral vector expressing hemagglutinin and nucleoprotein in the presence or absence of IL-1β, which is a potent mucosal adjuvant that significantly increases antibody and T cell responses.This high-resolution analysis of BAL cells provides insights into porcine resident leucocytes in the lung environment and the changes they undergo following infection and immunization.

Experimental design and cell sorting of BAL.
In the present study, we utilized samples from a previous animal experiment described in (44).The experiment involved four groups of inbred Babraham pigs: one group (n=5) was infected intranasally with pH1N1 influenza virus, another group (n=5) received intranasal immunisation with a recombinant adenoviral vector expressing hemagglutinin (HA) and nucleoprotein (NP) Ad-HA/NP, and the third group (n=5) was immunized intranasally with Ad-HA/NP and a recombinant adenoviral vector encoding porcine IL-1β (Ad-HA/NP+Ad-IL-1β).Controls (n=3) received PBS intranasally (Figure 1A).
The animals were culled 21 days post immunization or infection and tissues harvested for immunological and virological analysis.In the previous study, immune responses of these groups were extensively analyzed and showed that IL-1β acts as a potent mucosal adjuvant, increasing antibody responses, proportions of plasma and activated B cells, cytokine production of T cells and NP tetramer-specific CD8 T cells in the Ad-HA/NP+Ad-IL-1β group (44).
Given the well-defined nature of these samples we chose to perform scRNA-seq analysis on leukocyte populations sorted from the BAL.We were particularly interested in the BAL as T cells in this niche are inaccessible to intravenous CD3 antibody and have a phenotype characteristic of TRM in other species (45,46).As the cellular composition of BAL samples is strongly dominated by alveolar macrophages, we applied cell sorting to isolate macrophages, CD4 T, CD8 T and B cells from 3 animals per group and subjected equal proportions of these cells (12,000 cells each) to scRNA-seq (Figure 1B, Figure S1A).
Pig BAL scRNA-seq identifies typical immune cell populations.Graph-based clustering of the processed transcriptomic dataset identified 21 distinct cell clusters.Of these, three clusters were excluded for being of low quality or under suspicion of containing doublets.Three clusters were additionally split into subclusters to separate CD4 from CD8 T cells or split T cells from non-T cells, creating a final dataset with 22 cell clusters for annotation.To annotate each cluster with conventional cell-type names, marker genes defining each cluster were extracted and each cluster was screened for expression of a list of typical defining marker genes (Table S1, Figure S1B).From this information, each cluster was manually annotated with an established cell-type name (Figure 1C).
Most clusters were found to be comprised of well-defined cell-types and were easily annotated, such as macrophages, conventional dendritic cells (cDCs), monocytes, plasma cells, B cells and T cells, though particular subclusters presented more niche and complex phenotypes.These included 18a and 18b, which were annotated as 'mitotic T cells' based on their strong expression of E2F1 and E2F2, and subclusters 9a, 9b and 9c, which were found to have varying degrees of expression of CD8 T cell and NK cell related genes, including CD3E, CD8A, CD8B, KLRK1, KLRB1 and KIT.We also identified very small clusters of CD4 and CD8 T cells (20a and 20b, with 18 and 8 cells, respectively) that differed from the major CD4 and CD8 T cell clusters (7 and 6, respectively) by expressing higher levels of genes associated with naïve and/or central memory T cell states (e.g.BACH2, SELL, LEF1) and lower levels of genes encoding cytoskeletal and adhesion molecules (e.g.LGALS1, S100A4, VIM) that were previously found to characterize lung-localized human T cells (47).
Cluster 14 expressed T cell-related genes in an atypical pattern: strong expression of CD28 and CD5, with low expression of CD4, CD8A and CD8B, bifurcated expression of CD3D and CD3E, and very high expression of EOMES.Pending further investigation, cluster 14 has been annotated as 'T cell like'.
Cluster 21 is suspected to be T cell-related, due to adjacent clustering to identified T cell clusters in both UMAP (Figure 1C) and TSNE (Figure S1D) space.However, cluster 21 does not express any of the conventional αβ T cell markers: CD28, CD5, CD3D, CD3E, CD4, CD8A, or CD8B.Expressed T cell-related genes include GATA3, IL1RL1, IL7R and KIT, suggesting that these cells may be innate lymphoid cells, but without corroborating evidence, we have currently classified these cells as 'unknown'.
Cluster 19 showed strong expression of ACTG1, GPR37, NAPSA and SFTPA1, typically expressed by mucus producing epithelial cells.Cluster 9 and its subclusters were enriched for genes typically associated with NK cells such as KLRK1, KLRB1, KIT and NCR1.As our cell sorting strategy was not intended to capture NK cells, NK cell-like T cells or epithelial cells, clusters 9 and 19 risked presenting a biased representation of these larger cell populations.We therefore excluded clusters 9 and 19, leaving 18 clusters for further analysis.
Confirmation of scRNA-seq derived cell types by flow cytometry.We aimed to confirm cell types identified by scRNA-seq by flow cytometry (Figure 2).By using an antibody panel designed for phenotyping of myeloid cells (Table S2) and gating on FSC-A high SSC-A high cells (Figure S2A), we identified a prominent population of CD172a high CD163 + CD14 -/dim CD16 + macrophages (Figure 2A), which represent scRNA-seq clusters 1, 4 and 8. Within CD163 - cells, we identified a CD14 + CD16 + subset, representing putative monocytes identified in scRNA-seq cluster 15.CD14 -cells were further analysed for CD172a/CADM1 co-expression which allowed the identification of CD172 dim CADM1 high subset and a CD172a high CADM1 + subset, representing cDC1 and cDC2 cells, respectively (48), and identified in cluster 13 by scRNA-seq.CD4 + CD172a + CADM1 -pDCs were not identified in BAL by flow cytometry (data not shown) and no pDC phenotype was found by scRNA-seq.
The presence of B cells and plasma cells was also confirmed by flow cytometry (Figure 2B) using a specifically designed B cell panel (Table S2).CD79α + cells contained a prominent subset of IRF4 high Blimp-1 + plasma cells, which were also Pax-5 -and thus resembled clusters 2 and 3 of the sc-RNAseq data.CD79α + IRF4 dim Blimp-1 -cells had variable levels of Pax-5 but overall likely represent B cell clusters 12 and 16 from the scRNA-seq data.
CD4 and CD8 T cells were investigated by a dedicated T cell panel (Table S2).
After exclusion of γδ T cells, the majority of remaining CD8 T cells were CD161 -NKp46 - perforin -, resembling CD8 T cell cluster 6 identified by scRNA-seq.Within this phenotype, some cells expressed Ki-67, which may represent "mitotic" CD8 T cells found in cluster 18a.
Our sorting strategy for scRNA-seq excluded most NK and γδ T cells present in BAL.
However, a flow cytometry panel designed for the characterisation of nonconventional T cells (Table S2) indicated the presence of a substantial proportion of γδ T cells, which were largely CD8α + CD16 -NKp46 -(Figure S2C).Within the CD3 -population, CD161/CD8α defined putative NK cell subsets, which exhibited substantial phenotypic heterogeneity in terms of CD16, NKp46 and perforin expression.This phenotypic analysis confirmed cell types identified by scRNA-seq and indicated further heterogeneity within NK cells and nonconventional T cells.
Comparison with PBMCs delineates similar and distinct cell-types.To validate our celltype annotations and to further examine features of BAL tissue resident populations, we compared our pig BAL scRNA-seq dataset with a publicly available pig PBMC scRNA-seq dataset, published by Herrera-Uribe et al (23).Both datasets were restricted to shared genes, batch corrected and merged using a mutual nearest neighbors approach (Figure 3).When plotted as a UMAP, several cell types such as monocytes, cDCs, plasma cells (antibody secreting cells, ASCs in the PBMC dataset), B cells and some T cells formed integrated and/or adjacent clusters between BAL and PBMC (Figures 3A-B), implying similar patterns of overall gene expression.Cell-types without near neighbors between the datasets, i.e.BAL macrophages, the BAL 'T cell like' cluster 14 and PBMC γδ T cells formed distant and distinct clusters (Figures 3A-B).Looking more closely at T cell populations (Figure 3C) revealed two clusters that partially merge: BAL conventional T cells with PBMC αβ T cells (both CD4 + and CD8 + ) and BAL mitotic T cells with a separate cluster of PBMC αβ T cells.
To further delineate how individual clusters compared between BAL and PBMC, the top 20 marker genes for each cluster in the PBMC dataset were used to generate an index for projecting the annotated BAL cells onto the PBMC clusters.Figure 3D provides a visual representation of this mapping, identifying the proportion of cells in each BAL cluster that showed gene expression similarity to each PBMC cluster.To examine which genes were driving or inhibiting mapping between clusters, all clusters and all genes used for mapping were plotted as a heatmap (Figure S3A).Subsets of B cells and T cells were separately plotted (Figures 3E-F) using gene clusters identified in Figure S3A alongside defining genes from Table S1; additional subset heatmaps are also provided in Figure S3B-G.Broadly, cell-types in the BAL mapped to the same or closest equivalent cell types in the PBMC dataset.BAL macrophages and monocytes mapped to PBMC monocytes (there are no macrophages in PBMC), BAL cDCs to PBMC cDCs, BAL plasma cells to PBMC ASCs and BAL T cells to PBMC T cells.In some cases, mapping at the subset level was ambivalent.For example, the BAL macrophage clusters mapped to PBMC monocyte clusters 20 and 27 (Figure 3D), suggesting minimal differences in the patterns of top marker gene expression between these PBMC monocyte clusters.At the same time, this example also demonstrated preferential mapping such that BAL macrophages mapped to PBMC monocyte clusters 20 and 27 but not to monocyte clusters 13, 19 or 25.Such preferential mapping was apparent across most cell types analysed, with only a subset of PBMC clusters identifying as close matches for BAL cells.This implies that the heterogeneity in gene expression that drove the clustering in the PBMC dataset also embodied some of the divergence between BAL and PBMC cell types.Of the 33 PBMC clusters, BAL cells only substantially mapped to 12, suggesting that many of the PBMC cell types/states are distinct from those found in the BAL or relate to cell-types not included in our sorting strategy.
There are also exceptions to expected mappings.In the BAL, both clusters 12 and 16 were annotated as B cells, but only cluster 12 mapped to PBMC B cells, while cluster 16 failed to map (Figure 3D). Figure 3E illustrates the gene expression of these B cell clusters, with BAL cluster 16 uniquely showing reduced expression of numerous ribosomal genes (RPL*, RPS*).Likewise, BAL T cells mapped to PBMC T cells, but there was heterogeneity in mapping of T cell subsets.BAL clusters 18 and 20 were previously subclustered into separate mitotic CD4 + and CD8 + populations.However, BAL mitotic T cells (cluster 18) mapped primarily to PBMC CD8 + αβ T cell cluster 14, and BAL cluster 20 mapped primarily to PBMC CD4 + αβ T cell cluster 0 (Figure 3D).For the largest conventional BAL T cell clusters (6 and 7) the primary mapping was to PBMC CD2 -γδ T cell cluster 21 and, to a lesser extent, PBMC CD4 + αβ T cell cluster 4 (Figure 3D,F).Furthermore, the 'unknown' BAL cluster 21 also mapped strongly to PBMC CD2 -γδ T cell cluster 21.Conceptually it is unsurprising that a TRM-dominant tissue like the BAL would exhibit distinct gene expression from the circulating T cells that inhabit the peripheral blood.These examples also demonstrate how gene expression representing cellular functions such as cell cycle regulators or cytokine regulation (Figure 3F,     S3) can outweigh the relatively stable cell surface markers commonly used to distinguish different cell types when performing mapping between datasets using the whole transcriptome.
Frequencies of Treg cells.We next tested the scRNA-seq data for differential abundance of cells within the annotated cell types between treatment conditions (Table S3).This analysis identified that pigs immunized with Ad-HA/NP+Ad-IL1b have greatly reduced proportions of cells in cluster 10, annotated as Tregs, within the BAL.Among CD4 + T cells, Tregs make up an average of 0.02% in Ad-HA/NP+Ad-IL1b pigs versus an average of 0.09% in the other conditions (adjusted p=7.67e -8 , Figure 4A).This apparent loss of Tregs was not mirrored by any of the other CD4 T cell clusters and was consistent between samples within each condition (Figure S4A).
Based on this significant reduction of Tregs in the BAL of IL-1β treated animals seen by scRNA-seq, we further analysed Treg frequency by flow cytometry in BAL and tracheobronchial lymph nodes from the same pigs, plus two additional pigs from immunized and infected groups.Interestingly, there was a significant reduction of Tregs, defined as CD3 + CD4 + CD25 high FOXP3 + in the tracheobronchial lymph nodes of the Ad-HA/NP+Ad-IL-1β immunized group compared to the other three groups (Figure 4B).However, the reduction in Tregs that was observed in BAL scRNA-seq data was not replicated by the BAL flow analysis.This is likely because the Treg populations isolated by scRNA-seq and by flow are, whilst overlapping, not the same.In both scRNA-seq and flow, Tregs are defined as CD3 + CD4 + FOXP3 + cells, however in the scRNA-seq analysis, cells are grouped into clusters based on their top 5000 highly variable genes (HVGs) prior to being identified according to their expression of defining genes.Hence scRNA-seq Tregs (cluster 10, 799 cells), represent a population of uniformly CD3 + CD4 + FOXP3 + T cells that are also phenotypically distinct from the other CD4 T cell clusters (based on 5000 HVGs).FOXP3 + cells also exist within CD4 T cell cluster 7 and the flow analysis would include these FOXP3 + cells alongside cluster 10 as a single Treg population.Hence scRNA-seq Tregs and flow Tregs represent two overlapping but not identical populations.These data therefore imply that the BAL Treg population that is diminished during IL-1β immunization is a phenotypically distinct subset of BAL Treg and that others are not diminished.
We previously showed that IL-1β is a potent mucosal adjuvant and that intranasal immunization with Ad-HA/NP+Ad-IL-1β increased T cell cytokine production in Babraham pigs (44).Based on our observations of altered Treg abundance in this condition, we next investigated whether Tregs modulate cytokine production in the BAL.To this end, sorted CD4 + CD25 high Tregs and CD4 + CD25 neg control cells from the BAL of Ad-HA/NP immunized Babraham pigs were co-cultured with BAL cells from Ad-HA/NP+Ad-IL-1β immunized animals in the presence of pH1N1.The expression of IFNγ and TNFα was then assayed by intracellular staining.The addition of Tregs (CD4 + CD25 high ) to the BAL nearly abolished cytokine production, while a minor effect was detected when CD4 + CD25 neg control cells were introduced (Figure 4C).Together with scRNA-seq and flow cytometry quantification of Treg populations, these data suggest that IL-1β reduces the proportion of functionally active Tregs, contributing to increased cytokine production, in agreement with previous studies demonstrating the suppressive role of porcine Tregs (49).

Differential expression analysis reveals condition-and cell-specific phenotypes. To investigate differential gene expression between our experimental conditions we used a
NEgative Binomial mixed model Using a Large-sample Approximation (NEBULA) (50) run discretely on all cells, B cells, CD4 T cells, CD8 T cells and macrophages (Figure S5).To distinguish if differentially expressed (DE) genes were distinct or shared between cell types for each experimental condition, the identified significant genes were visualized as Venn diagrams (Figures 5A-C) and tested for Gene Ontology (GO) enrichment to approximate biological function (Figures 5D-E, Table S4).For all experimental conditions, most DE genes were distinct to specific cell types, with only a minority being shared (Figures 5A-C).For Ad-HA/NP versus PBS (Figure 5A), only macrophages and CD8 T cells presented any significant DE genes (adjusted p value less than 0.01), all of which were downregulated (Figure 5A).GO enrichment for these comparisons returned fewer than 5 DE genes per GO term and is therefore not shown.Ad-HA/NP+Ad-IL1β versus PBS identified considerably more DE genes and significant GO terms, with 11 DE genes shared between cell types, 7 of which were shared between CD4 and CD8 T cells (Figure 5B).DE genes in both CD4 and CD8 T cell populations were enriched for lymphocyte activation and apoptosis-related GO terms, enrichment testing in macrophages returned only broad GO terms, and analyses in B cells returned fewer than 5 DE genes per GO term.The enriched GO terms in T cells implied that Ad-HA/NP+Ad-IL1β elicited a strong immunogenic response, potentially overstimulation given the presence of apoptotic terms.Comparing pH1H1 infection with PBS, in addition to numerous distinctly expressed DE genes, NEBULA identified 10 DE genes shared by cell-types Figure 5D).IFI6 was particularly notable as being universally upregulated across leukocyte subsets.GO-term enrichment of DE genes in macrophages, CD4 T cells and CD8 T cells presented a uniform trend towards anti-viral response mechanisms, whereas B cells again returned fewer than 5 DE genes per GO term (Figure 5E).
Overall, these comparisons suggest that both Ad-HA/NP+Ad-IL1β and pH1N1 invoke a greater immune response than Ad-HA/NP, resulting in greater differential gene expression.
IFI6 marks lasting impacts of pH1N1 infection on BAL leukocytes.The differential expression analyses identified IFI6 as strongly upregulated across multiple cell types in pH1N1 infected pigs.Further investigation revealed that IFI6 was upregulated in all scRNA-seq cell clusters in pH1N1 infected pigs (Figure 6A).To assess if IFI6 expression was associated with the expression of other genes in different cell types, we used the mixed model of NEBULA to perform a separate differential expression analysis with IFI6 expression itself as the dependent variable, asking which genes were associated with IFI6 expression.Genes associated with IFI6 expression in each cell type were compared in Figure 6B and Figure S6 A-D, identifying numerous IFI6 related genes, dominantly in macrophages, CD4 T cells and B cells.GO enrichment (Figure 6C) of these genes in CD4 and CD8 T cells identified terms related to antiviral and effector responses, overlapping, as would be expected, with the GO terms enriched in pH1N1 infected pigs.For macrophages and B cells, IFI6-associated genes were enriched for metabolism-related GO terms.
We also compared the antiviral state of porcine BAL cells between pH1N1-infected and control pigs through a functional interferon assay.To this end, we infected BAL cells with recombinant virus-like particles of vesicular stomatitis virus in which glycoprotein had been replaced with GFP (VSV∆G-GFP) (51,52).The kinetics of GFP expression from VSV∆G-GFP was used as a measure of virus replication.The mean GFP expression observed in BAL cells from pH1N1 infected pigs was significantly lower than that in BAL cells from the control group, confirming the anti-viral state of pH1N1 BAL cells at 21 days post infection (Figure 6D, Figure S6E-F).This analysis indicates that although it would be expected that cells in BAL always represent a population that had migrated to this site in response to an antigenic stimulus in the respiratory tract, local factors modify gene expression of BAL cells, as demonstrated by the ubiquitous expression of IFI6 following influenza infection.

Discussion
Studies in animal models have highlighted the unique nature of local lung immunity and the critical role this compartmentalization plays in facilitating protective immune responses to respiratory pathogens.Here we present for the first time a single-cell transcriptomic map of steady-state airways and airways following influenza infection or pulmonary immunisation in the pig, a highly relevant large natural host model for influenza.We performed our study in inbred Babraham pigs, which have identical MHC (swine leucocyte antigen) and allow for fine-grained immune response analysis (53).The similar size and anatomy of pig and human organs make this model particularly valuable for translational research.
Within the BAL, the overwhelming majority of cells exhibit the characteristics of tissue residency.Intravenous CD3 administration, which is the most commonly used method to define tissue residency, does not stain porcine BAL T cells (45,46).BAL is thus a particularly useful source of information about tissue-resident immune responses as it reflects the site of respiratory infections and is accessible during bronchoscopy in humans.Single-cell RNA-seq studies of BAL have become more common in recent years, particularly since the COVID-19 pandemic revealed the utility of monitoring the immune response at this accessible site of infection.Studies in humans (18,54), ferrets (19), horses (20), dogs (21) and mice (55) have characterized the various immune and epithelial cells from the BAL, revealing extensive heterogeneity among these populations, substantial differences in their abundances as compared with PBMC (54) and respiratory infection/inflammation-associated changes (18,19,55).
Our study examined BAL extracted from the lung 21 days after vaccine or influenza challenge to understand the lasting impacts of induced immune responses.While porcine BAL cells are dominated by macrophages, our sorting strategy allowed us to make key observations about relatively rare lymphocyte populations (e.g.Treg) while still capturing information about the phenotype and abundance of more prominent macrophage populations and finding shared gene expression changes across immune cell subsets.Leveraging available porcine PBMC single-cell RNAseq data (23) enabled validation of our cell type annotations and highlighted likely phenotypic differences between tissue and circulating immune cell populations, particularly among T and B lymphocytes.
Finally, our successful use of frozen BAL samples from a previous immunisation study (44) corroborates conclusions drawn in a recent equine study (20) that frozen BAL is a useful source of single-cell transcriptomic information for large-animal model immune responses, opening the door to further re-use of samples from key respiratory immunology studies.
By utilizing samples from influenza infected and immunized animals, we had the opportunity to investigate parameters that change under these circumstances.Our data indicate that IL-1β administration resulted in a significant reduction of Treg in the BAL as determined by scRNA-seq, although there were no changes in Treg frequency in BAL when identified as CD3 + CD4 + , CD25 high Foxp3 + on the protein level by flow cytometry in the same samples.This highlights the power of transcriptomic profiling, allowing higher dimensional clustering of immune cell subsets compared to the limited marker sets used in flow cytometry.
Mucosal production of IL-1β has been shown to attract innate and adaptive immune cells by the induction of cytokines and adhesion molecules and to have a direct proliferative effect on CD4 and CD8 T cells (56)(57)(58).Intranasal delivery of IL-1β, together with HA-and NP-encoding adenoviral vectors (Ad-HA/NP+Ad-IL-1β), significantly increased the immunogenicity and heterotypic protection against the influenza virus in mice (30).However, in pigs although IL-1β increased antibody and T cell immune responses, this did not reduce lung pathology or viral load following H3N2 influenza challenge (44).Tregs are responsible for maintaining immune homeostasis by supressing excessive immune responses (49,59).Of note, effector T cells in the BAL were still susceptible to Treg inhibition when cocultured with Tregs from Ad-HA/NP immunized animals.Hence, reduced numbers of functionally active Tregs can lead to imbalance of cytokine production, and this proinflammatory environment can trigger inflammation and tissue damage.Indeed, outbred pigs immunized with AdHA/NP+Ad-IL-1β showed increased lung pathology (44), confirming that maintaining a balance of Tregs is crucial for the proper functioning of the immune system.
Another interesting observation from this study is that pH1N1 infection drives upregulation of IFN alpha-inducible protein 6 (IFI6) in most BAL cell subsets.The expression of IFI6 is tightly regulated in response to IFN alpha stimulation and plays a crucial role in restricting viral infections including West Nile Virus, Dengue Virus, Yellow Fever, Hepatitis B and Hepatitis C viruses through interference with viral RNA synthesis and translation (60)(61)(62).In agreement with these studies, we observed reduced GFP expression in VSV∆G-GFP infected BAL cells from pH1N1 infected pigs compared to control pigs, confirming an inhibitory role of IFI6 following influenza infection.This is in contrast with a recent report showing that IFI6 can positively affect virus infection and that IFI6 overexpression results in increased influenza replication in cell cultures and in the murine model (63).Irrespective of the mechanism of action and whether this reflects species differences between mice and pigs, or between VSV and influenza infection, the sustained upregulation of IFI6 suggests that the immune system is still actively engaged in combating the virus or resolving inflammation and tissue damage at day 21 post infection.Our results also indicate that although all tissue resident lymphocytes must have a common gene expression programme enabling them to migrate to and survive in tissue sites, local pathogen-related factors can profoundly modify these programs.These modifications are detectable for a surprisingly long time, 21 days, after infection or immunization.
In summary our findings provide a comprehensive high dimensional cellular reference for porcine airways in different immunological contexts, paving the way for a deeper understanding of respiratory physiology and immune responses in pigs.

Materials and methods
Infection and immunization animal studies.Samples obtained from a previous experiment conducted on inbred Babraham pigs were utilized (44).Briefly three groups of 5 pigs each, were treated as follows: 1) infected intranasally with 8x10 6 PFU of A/swine/England/1353/2009 (pH1N1) using a mucosal atomization device (MAD); 2) immunized intranasally with 1x10 9 particles of recombinant adenoviral vectors expressing NP from H1N1/PR8/34 virus and HA from pH1N1/Texas/05/2009 (Ad-HA/NP) using a MAD and 3) immunized intranasally with a mix of 1x10 9 particles of Ad-HA/NP and 1x10 9 particles recombinant adenoviral vector encoding porcine IL-1β (Ad-HA/NP+Ad-IL-1β) using a MAD.
Control pigs (n=3) received PBS intranasally.The animals were culled three weeks later, and bronchoalveolar lavage (BAL), tracheobronchial lymph nodes and blood collected.Virus shedding and immune responses in respiratory tissues and peripheral blood were described in (44).
Representative phenotyping data is shown in Figure S1A.Cells from each fraction were pooled and subjected to partitioning and barcoding on a 10x Genomics Chromium iX Controller.scRNA-seq library preparation and sequencing.For 10X, single cells were processed with the Chromium iX controller, using the Next GEM chip G, and single cell 3' kit v3.1 library preparation kit (10X Genomics, Pleasanton, CA).An estimated 10,000 cells were targeted for recovery during partitioning.Samples were randomly grouped and run with two libraries per sequencing run on three NextSeq High Output 150 cycle reagent kits (Illumina, San Diego, CA), with 1% PhiX, targeting approximately 200 million reads per sample.Bcl files were demultiplexed using 'cellranger mkfastq' and aligned/counted using 'cellranger count' (cellranger-7.0.0) using default parameters and Sus scrofa genome (genome assembly 11.1, Ensembl release 107).
scRNA-seq QC and pre-processing.Data analysis was run using R (version 4.1.1)and RStudio (build 382).Empty droplets, barcode-swapped droplets and ambient RNA were removed from the unfiltered cellranger output using DropletUtils (1.18.1) (64,65).Genes with no counts and cells with outlying mitochondrial/nuclear gene log-ratios (> 3 median absolute deviations) were removed.One of twelve samples was excluded as a strong outlier due to low sequencing saturation and excessive detected cell number.All remaining samples were normalized by pooled factors (66) and combined to form the analysis dataset.Batch effects were corrected between sequencing runs using fastmnn (batchelor 1.14.0)implementation of mutual nearest neighbors batch correction (67).

scRNA-seq clustering and annotation. Shared-nearest-neighbors graph-based clustering was
performed on batch-corrected data using clusterCells (bluster 1.8.0) with the top 5000 highly variable genes (HVGs), specifically using Jaccard index to weight edges and the Louvain method for community detection, with a k value of 20.This generated 21 initial clusters.Cluster 5 exhibited particularly low counts, limited numbers of genes detected and lack of any positive defining features and was therefore defined as poor-quality and excluded from further analysis.
Co-expression based doublet scoring by cxds (scds 1.14.0)indicated that cluster 17 was likely comprized of a large proportion of doublets.Cluster 11 also had a substantial proportion of doublets, and manual inspection revealed co-expressed macrophage and B cell marker genes, suggesting that this cluster was also likely doublets.These clusters were therefore excluded from further analysis.Clusters 18 and 20 were found to express a dichotomy of CD4 and CD8 expression and so were split into CD4 + and CD8 + subclusters.This was achieved by identifying HVGs within the cluster of interest, finding those that correlated best with CD4 and CD8B expression using a zero-inflated Kendall's Tau correlation metric implemented in scHOT (v1.10.0), and then re-running batch correction and subclustering (Jaccard index to weight edges and the fast greedy modularity optimization algorithm for community detection) on cells in the cluster of interest.Varied expression of genes including CD3E, CD8A, CD8B, RORC and KLRB1 within cluster 9 also prompted sub-clustering.This was done in an unsupervized manner by identifying the top 2000 HVGs in cluster 9 cells and then re-running batch correction and clustering (Jaccard index to weight edges and the fast greedy modularity optimization algorithm for community detection).This split the population into three subclusters.
Clusters were manually annotated with conventional names (cell-types) by identifying the positive and negative expression of typical defining markers, as listed in Table S1 and illustrated in Figure 1C and Figure S1B, and looking at top marker genes within each cluster as defined by the scoreMarkers function (scran 1.26.0,Table S5, provided in online repository).
Cluster 19 showed strong expression of ACTG1, GPR37, NAPSA and SFTPA1, typically expressed by mucus producing epithelial cells.Cluster 9 and its subclusters were enriched for NK cell typical genes such as KLRK1, KLRB1, KIT and NCR1.As our cell sorting strategy was not intended to capture NK cells, NK-cell-like T cells or epithelial cells, clusters 9 and 19 risked presenting a biased representation of these larger cell populations and were therefore removed from further analysis.
scRNA-seq analysis, differential abundance analyses.To compare the abundance of cells within each cluster across conditions, we performed differential abundance analyses using a negative binomial generalized linear model with empirical Bayes quasi-likelihood F-tests (68,69), as implemented in edgeR (v 3.40.0).We included sex as a covariate and fitted the model with all conditions to test for differential abundance in any condition.The number of cells per cluster was normalized by sample cell number within the comparison, and trend estimation was turned off when estimating dispersion and quasi-likelihood dispersion.To account for the sort and pool strategy used to balance cell types among the sequenced cells, clusters were divided according to their membership of the 4 sorted cell types (macrophages, CD4 T cells, CD8 T cells, B cells) before running differential abundance analyses within each.Those clusters that did not obviously fall into one of these cell types were excluded in order to keep the universes for differential abundance analysis equivalent to or smaller than the sorted cell types, thereby avoiding sorting biases in intra-cell-type differential abundance analyses.5 and 6) were performed using a NEgative Binomial mixed model Using a Large-sample Approximation (NEBULA) (50), using PBS treated pigs as a baseline, blocking by sample and accounting for cell type when multiple cell types were present.Final p values were adjusted using the Benjamini-Hochberg procedure (70) as implemented in p.adjust (stats v4.1.1).Gene Ontology Enrichment.Lists of differential expressed genes were passed to topGO (2.44.0) using the Genome wide annotation for pig (DOI: 10.18129/B9.bioc.org.Ss.eg.db, release 3.13.0),biological process subontology and a minimum node size of 5. Three sets of enrichment tests were performed, classicFisher, classicKS and elimKS with all results available in Table S4.This PBMC dataset was then trimmed to include only genes shared between the BAL and PBMC datasets (13050 shared genes); clustering and cell type annotation in the PBMC dataset was left unchanged.The top 20 differentially expressed genes for each cluster within each dataset (55 total clusters) were identified using findMarkers (scran 1.20.1)and combined (402 final genes after removing duplicates) to provide a data subset for mutual nearest neighbors correction (67).mnnCorrect (batchelor 1.8.1) was used to correct for batch effects in each dataset and then to merge the datasets prior to plotting of joint UMAPs.

scRNA-seq analysis, differential expression analyses. Differential expression analyses for all cells, macrophages, CD4 T cells, CD8 T cells and B cells (Figures
To directly compare gene expression between cell types and clusters, the top 20 differentially expressed genes for each cluster in the PBMC dataset (36 clusters) were identified using findMarkers (scran 1.20.1),duplicates removed (268 final genes) and used as selected features to index the clusters in the batch corrected PBMC dataset via scmap (scmap 1.14.1).
The clusters in the batch corrected BAL dataset were then projected onto the indexed PBMC dataset via scmap, using the default threshold of 0.7 (i.e.70%, equal to 188 genes), and the results plotted as a Sankey plot using networkD3 (networkD3 0.4).
To identify the sets of genes driving the mapping between BAL and PBMC clusters, all genes and all clusters inputted into scmap were plotted as a heatmap (Figure S3A), with clusters of genes manually annotated (guided by Gene Ontology terms and the literature).
Selected subsets of cell clusters and gene clusters were then plotted as separate heatmaps for Figures 3E-F and Figures S3B-G.
Flow cytometry (FCM) of BAL cells.BAL cells were thawed, washed in PBS, plated into 96well round bottom plates for flow cytometry staining and labelled with antibodies grouped into four staining panels (Table S2).Some antibodies were labelled with isotype or species-specific secondary antibodies conjugated to fluorochromes.Biotinylated antibodies were labelled with streptavidin conjugates.Free binding sites of secondary antibodies were blocked with mouse IgG (1μg per sample, Jackson) prior to incubation with further, directly conjugated antibodies.
Live/dead discrimination was achieved by labelling of samples with VDeFluor780 (ThermoFisher).Following labelling of surface markers, cells were fixed and permeabilized with eBioscience™ Foxp3 / Transcription Factor Staining Buffer Set (ThermoFisher) according to manufacturer's instructions, followed by labelling with antibodies against intracellular targets.Incubation steps for extracellular markers lasted for 20 minutes, incubation with antibodies for intracellular targets was performed for 30 minutes.Cells were kept at 4°C at all times.Cells from BAL and TBLN were analysed on a Cytek Aurora spectral flow cytometer, equipped with 5 lasers (UV 355nm, Violet 405nm, Blue 488nm, Yellow-Green 561nm, Red 640nm) and 64 fluorescence detection channels UV: 16, Violet: 16, Blue: 14, Yellow-Green: 10, Red: 8. Per sample, at least 2 x 10 5 cells were acquired.Data was acquired using Cytek's Spectro Flo software and later analysed with FlowJo version 10 (BD Biosciences).
CD4 + CD25 high and CD4 + CD25 neg cells were sorted on a FACSAria UIII (Figure S4B).To determine the effect of Tregs, 4 x 10 4 cells of the sorted cell populations were added to 2 x 10 6 BAL cells from animals immunized with Ad-HA/NP+Ad-IL-1β.These cocultures were stimulated with H1N1pdm09 (MOI of 1) or medium control overnight.Intracellular staining (ICS) was performed to analyse IFNγ and TNFα production by CD4 and CD8b T cells.
Additionally, Near-Infrared Fixable LIVE/DEAD stain (Invitrogen) was used for identification of live cells.
Virus infectivity assays to determine antiviral state of pig BAL.Pig BAL cells from pH1N1 infected or control groups were challenged with infectious recombinant VLPs of vesicular stomatitis virus (VSV-ΔG-GFP) in which the glycoprotein was replaced with GFP which was used as a marker to monitor the replication of VSV-ΔG-GFP.The procedure for infecting cells with VSV-ΔG-GFP was previously described (51,52,71) with slight modifications in which we measured the kinetics of GFP expression from VSV-ΔG-GFP in a live cell imager (IncuCyte-S3 ® ).

Figure 1 :
Figure 1: Experimental design and UMAP clustering of BAL immune cells (A) Design of pig immunisation/ infection experiment.Total number of pigs per treatment is given.For scRNA-seq 3 pigs were randomly selected for treatments 1, 2 and 3. (B) Graphical overview of applied sorting strategy.(C) Uniform manifold approximation and projection (UMAP) of all viable cells from pig BAL scRNA-seq analysis, coloured by cell-type and labelled with cluster numbers.

Figure 2 :
Figure 2: Confirmation of scRNA-seq cluster annotation by flow cytometry Annotation of clusters in scRNA-seq was confirmed by flow cytometry using antibody panels addressing key marker combinations.Colours used in UMAP clustering (left) match with coloured gates and histograms in flow cytometry to indicate the same cell type (right).(A) Different subsets of myeloid cells.(B) B cell and plasma cell subsets.(C) CD4 T cell subsets.(D) CD8 T cell subsets.

Figure 3 :
Figure 3: Comparison of pig BAL and PBMC scRNA-seq datasets (A-C) Uniform manifold approximation and projection (UMAP) of BAL and PBMC cells.(A) All cells from BAL and PBMC scRNA-seq datasets coloured by cell type.(B) All cells from BAL and PBMC scRNA-seq datasets, coloured by tissue.(C) T cells from BAL and PBMC scRNA-seq datasets, coloured by cell type.(D) Mapping of BAL cells (left) to PBMC clusters (right) via scMAP, coloured by cell type.Each node represents one cluster and is labelled with the cell type and cluster number.Unmapped PBMC nodes and links constituting less than 15% of the cells in a single BAL cluster were omitted for legibility.(E-F) Heatmaps of BAL and PBMC cluster gene expression.Each column represents one cluster, labelled with the tissue of origin, cell type and cluster number.Colour scale is generated from mean logcounts.(E) Heatmap of B cell clusters and (F) Heatmap of T cell clusters from BAL and PBMC scRNAseq datasets.

Figure 4 :
Figure 4: Abundance and suppressive function of Tregs (A) Tregs identified by scRNA-seq (cluster 10) as a proportion of all CD4 T cells for each treatment.(B) Tregs identified by flow cytometry as percentage within CD3 + CD4 + T cells for each treatment in cell preparations from BAL and tracheobronchial lymph nodes (TBLN).(C) Suppression of IFNγ and TNFα production in pH1N1 re-stimulated CD4 (top panel) and CD8

Figure 5 :
Figure 5: Differential gene expression with immune challenge (A-C) Venn diagrams presenting the overlapping and distinct differentially expressed genes for macrophages, CD4 + T cells, CD8 + T cells and B cells under each experimental condition.Red gene names indicate gene upregulation, blue indicates gene downregulation.Genes lacking conventional symbols (i.e.only ensembl IDs) have been omitted for legibility.All genes shown have an adjusted p value (Benjamini-Hochberg) of less than 0.01 when comparing gene expression in each condition versus PBS.(D-E) Bubbleplots presenting GO term enrichment for macrophages, CD4 T cells and CD8 T cells under each experimental condition.Based on all genes with an adjusted P value of less than 0.01.GO terms annotated with fewer than five differentially expressed genes have been omitted.(A) Ad-HA/NP treated pigs versus PBS, (B, D) Ad-HA/NP+IL1β treated pigs versus PBS and (C, E) pH1N1 infected pigs versus PBS.

Figure 6 :
Figure 6: IFI6 abundance, associated genes and functional consequences of IFI6 expression.(A) IFI6 expression across clusters in each experimental condition.(B) Venn diagram of genes differentially expressed as IFI6 is upregulated for macrophages, CD4 T cells, CD8 T cells and B cells.Red gene names indicate gene upregulation, blue indicates gene downregulation.All genes shown have an adjusted p value (Benjamini-Hochberg) of less than 0.05 and a logFC that is either more than 0.6 or less than -0.6 (C) Bubbleplot presenting GO term enrichment for genes differentially expressed as IFI6 is upregulated for macrophages, CD4 T cells, CD8 T cells and B cells.Based on all genes with an adjusted P value of less than 0.05.GO terms annotated with fewer than five differentially expressed genes have been omitted.(D) GFP counts of BAL cells following infection with VSV∆G-GFP 7hrs post-infection.Comparison of BAL cells from pH1N1 and control pigs.Asterisk indicates significant difference (p<0.05).

Figure legends supplementary figuresFigure
Figure legends supplementary figures

Figure S2 :
Figure S2: Gating strategy for BAL samples and phenotyping of unconventional T cells, related to Figure 2. (A) For phenotyping of myeloid cells, FSC-A high SSC-A high cells were gated, followed by exclusion of doublets and dead cells (low fluorescence intensity for live/dead dye).(B) As in (A) but for phenotyping of lymphocyte subsets, initially FSC-A low SSC-A low cells were gated, followed by exclusion of doublets and dead cells.(C) Unconventional T cells, which had been largely excluded from scRNA-seq by the applied sorting strategy, were investigated by flow cytometry.After exclusion of dead cells and doublets as shown in (B), γδ T cells (CD3 + TCR-γδ + ), non-γδ T cells (CD3 + TCR-γδ -) and non-T cells (CD3 -TCR-γδ -) were gated and analysed for expression of CD8α, CD161, CD16, NKp46 (CD335) and perforin.

Figure S3 :
Figure S3: Comparison of pig BAL and PBMC scRNA-seq datasets, related to Figure 3 (A) Heatmap of all cell clusters and genes used in the BAL to PBMC scMAP.(B-G) Additional heatmaps of cluster and gene subsets used in the BAL to PBMC scMAP.

Figure S4 :
Figure S4: Treg abundance and sorting of BAL Tregs, related to Figure 4. (A) Proportions of UMAP clusters (7, 10, 18b, 20a) representing CD4 T cells across all samples subjected to scRNA-seq.Different treatments are indicated at the bottom of the graph.(B).Gating hierarchy applied for sorting of CD4 + CD25 -(control) and CD4 + CD25 high (Treg) cells from BAL of Ad-HA/NP treated pig.

Figure S6 :
Figure S6: Differentially expressed genes associated with IFI6 upregulation and VSV infection of BAL cells, related to Figure 6.(A-D) Volcano plots of genes differentially expressed as IFI6 is upregulated in macrophages, CD4 T cells, CD8 T cells and B cells.Horizontal red line is at -log10(0.05).Vertical red lines are at -0.3 and 0.3.Adjusted p values were calculated using the Benjamini-Hochberg procedure.Adjusted p values below 1e-100 were capped at 1e-100.(A) Macrophages.(B) CD4 T cells.(C) CD8 T cells.(D) B cells.(E-F) GFP expression from infectious recombinant VLPs of vesicular stomatitis virus expressing GFP (VSV∆G-GFP) in porcine BAL cells at 7 hrs postinfection captured in live cell imager (IncuCyte-S3).(E) GFP expression from VSV∆G-GFP at 7hrs post-infection in H1N1 infected porcine BAL cells, (F) GFP expression from VSV∆G-GFP at 7hrs post-infection in naïve (control) porcine BAL cells.