Establishment of the infection model. The work flow of a typical experiment is shown in Fig. 1A. In order to optimize RNA quality for RNAseq, we compared two methods for RNA isolation. In the first method, the lung tissue pieces were incubated with RNAlater for 24 h at 4°C and then stored at -80°C in RNAlater prior to extraction. In the second method, they were snap frozen in liquid nitrogen and then stored at -80°C. The liquid nitrogen method yielded significantly higher RNA Integrity Numbers (Fig. 1B) and was thus used throughout. We measured LDH release in uninfected and IAV-infected HLTEs in order to assess spontaneous degradation of tissue during culturing and potential virus-induced cytopathic effects. LDH was released throughout the 72 h time course, reaching maximum values of about 7% and 12% in uninfected and infected HLTEs, respectively, indicating a mild-moderate cytopathic effect (Fig. 1C). To assess whether the HLTE model supports release of active IAV particles, we employed the immuno-foci assay to quantify replication-competent virus particles in the supernatant and measured viral RNA transcription through HA mRNA quantification. An apparent decline of replication-competent virus particles in the supernatant was noted over a 72 h period (Figure S1A). This was likely due to the fact that the initial viral inoculum was not replaced by virus-free medium, and it was thus concluded that viral titers in supernatant could not be used to measure viral replication. On the other hand, viral HA mRNA expression increased steadily and peaked by 48 h post infection (p.i.) (Fig. 1D). However, significant inter- and intra-donor variation was observed, likely due to differences in tissue integrity and cell type proportions between samples. To verify whether the tissues could mount the expected differential responses to IAV and the two bacterial pathogens (P. aeruginosa and BCG) in the first 24 h of infection, we measured several classic cytokines in supernatants and the corresponding mRNAs in tissue lysates (Fig. 1E). As expected, CXCL10 and IFNγ levels were highest in influenza infection, whereas IL6, IL10, IL1β, and PROK2 (which is highly inducible in macrophages by LPS and other TLR2/4 agonists (21)) reached higher levels in the bacterial infections. Taken together, these results suggested that HLTEs from emphysema patients constituted a valid model to study early transcriptomic tissue responses to viral and bacterial infections ex vivo.
Modest reprogramming of small noncoding RNA (sncRNA) populations. Small (< 50 nucleotides long) RNA sequencing of total RNA extracted 24 h p.i. revealed the expected sequencing depth: a mean 8.3*107 reads were obtained per sample; a mean 4.4*107 of these passed the length filter, of which 2.8*107 could be mapped to the human genome (hg38) (Fig. 2A). By sncRNA class, most reads mapped to miRNAs and substantially fewer to PIWI-interacting RNA (piRNA) and small nucleolar RNA (snoRNA) (Fig. 2B). Very few reads mapped to small nuclear RNA (snRNA) and ribosomal RNA (rRNA), which was consistent with the high quality of the input RNA. After filtering out lowly expressed sncRNA (< 20 reads in the four experimental groups combined), piRNA contributed the highest number of sncRNA species, followed by miRNA. However, the majority of the differentially expressed (DE; p-Adj ≤ 0.1) sncRNA species were miRNA followed by piRNA (Fig. 2C,D). A principal component analysis (PCA) did not reveal global differences in the sncRNA populations of the 4 groups (Fig. 2E). Indeed, there were only 2 DE piRNAs in IAV infection vs. uninfected tissue and 1 DE piRNA in BCG infection vs. uninfected tissue, and their expression in uninfected tissue was already low, indicating that their functional significance is uncertain (Fig. 2F,G). On the other hand, more substantial expression changes were observed in P. aeruginosa infection. 20 sncRNA were DE: 13 miRNA, 6 piRNA, and 1 snRNA (Fig. 2H). Among the DE miRNAs, mir-410-5p and miR-665 have been linked to increased cell proliferation and are candidate biomarkers for cancer. Decreased expression of one of the downregulated miRNA, miR-7974, has been reported in hypoxia (22) and oxidative stress (23), and increased expression in breast tumor tissue (24). A gene set enrichment (GSEA) analysis of predicted mRNA targets of the miRNA DE in P. aeruginosa infection suggested activation of pathways related to peroxisomes, inflammation, and cell cycle (Fig. 2I). However, statistical significance was modest. Nonetheless, DE of about 10% of the predicted mRNA targets could be verified in the long RNAseq dataset (Figure S1B) and DE of additional predicted targets would be expected at later time points. piRNAs were the 2nd most significantly DE sncRNA class in P. aeruginosa infection and were notably overrepresented among the most strongly DE sncRNA (e.g., piR_010299). piRNAs have been implicated in tissue responses to infections in humans (3) and animals (25), suggesting that their DE in this HLTE model is pathophysiologically significant.
Marked reprogramming in long RNA (mRNA, long noncoding RNA, pseudogene RNA) populations reveals signatures for bacterial vs. viral infection. Sequencing of long RNA (> 100 nucleotides long) revealed deep coverage in all samples in that we obtained a mean (SD) 42.1 x 106 (±7.3 x 106) reads mapping to protein-coding RNA (mRNA), long noncoding RNA (lncRNA and long intervening noncoding RNA, lincRNA), and pseudogene RNA species. To examine the spread of RNAseq data across the four groups, we generated a graphical representation of Fragments Per Kilobase Million (FPM) for each group. The plot illustrates that the data are equally distributed across all groups (Figure S2A). With respect to both total and DE RNA, the majority of reads mapped to protein-coding RNA, followed by lncRNA and then pseudogene RNA in infection with all three pathogens (Figure S2B-D), whereby the percentages of the four captured RNA classes (protein-coding RNA, lncRNA, lincRNA, pseudogenes) in the DE RNAs were similar across infection with the three pathogens (Figure S2E). A PCA based on all species of long RNA suggested similarities between IAV infection and control samples on the one hand and between the two bacterial infections on the other hand (Fig. 3A). Still, no clear separation was seen even though expression changes were substantially greater than in the sncRNA populations. Nonetheless, it was evident that the extent of differential expression increased from IAV (n = 131 DE RNAs) over BCG (n = 964 DE RNAs) to P. aeruginosa (n = 2423 DE RNAs) (Fig. 3B).
High levels of mRNA corresponding to all 8 IAV genome segments were detected in the IAV-infected samples, but only background signals in the other three groups, thus evidencing successful IAV infection of the HLTEs used for the RNAseq analysis (Figure S3). Volcano plots were then used to visualize differences between the three pathogens in DE with respect to uninfected tissue (Fig. 3C-H). In IAV infection, 131 RNAs were DE with respect to uninfected tissue. Several mRNAs were highly upregulated which are known to be strongly associated with host responses to influenza virus infection (Fig. 3C). In particular, CCL8, CXCL9, GBP1, IFI44L, IFI6, IFNL2, IL27, ISG15, MX1, and MX2 are all contained in the Gene Ontology (GO) term BP 0009615 response to virus (Fig. 3D), thus providing further evidence that the HLTE mounted the expected anti-IAV innate immune response. Furthermore, we noted 19 lncRNAs and lincRNAs among the DE RNAs (Fig. 3E). Of note, lncRNAs ADAM2-2 (ENSG00000253939), SIMALR (suppressor of inflammatory macrophage apoptosis lncRNA, ENSG00000226004), LINC02865 (ENSG00000229922), LGALS17A (ENSG00000226025), and LGALS14 (ENSG00000269460) were among the most highly upregulated transcripts in IAV infection (Fig. 3E). lnc-ADAM2-2 is a novel gene which encodes antisense RNA to the genes encoding the IDO1 and IDO2 enzymes, both of which have immunomodulatory functions in influenza virus infection (26, 27). lncRNA SIMALR and lincRNA LINC02865 (a suppressor of inflammatory macrophage apoptosis) are transcribed from overlapping sequences on chromosome 6. They have emerged as novel regulators of macrophage biology by suppressing inflammatory macrophage apoptosis via Netrin-1 (NTN1) (28) and are located just downstream of the gene encoding the anti-inflammatory deubiquitinase TNFAIP3 (also known as A20). LGALS17A and its downstream lncRNA lnc-LGALS14 share a genomic location on Chromosome 19 with galectins, which possess antiviral functionality in IAV infection (29). LGALS17A and lnc-LGALS14 are derived from a galectin pseudogene, and LGALS17A has previously been identified as an IFN-stimulated transcript (30). In addition, they are located downstream of a gene cluster relating to type III IFN responses, of which INFL2 was significantly upregulated in the current dataset, whereas upregulation of IFNL1 and 3 was less significant (p Adj > 0.1). We then aimed to validate the expression of a subset of the DE lncRNA, lincRNA and pseudogene RNAs by RT-qPCR. Robust PCR assays could be established for RP11-202G18.1, GBP1P1, SIMALR, AC093063.1, and LGALS17A. Indeed, expression of all five increased significantly 24 h after IAV infection, but to lower levels than CXCL10 mRNA, which was measured for comparison (Figure S4).
Transcriptomic tissue responses to the two bacterial pathogens were much more extensive than to IAV (Fig. 3F,G), and there was considerable overlap in DE RNAs, even though P. aeruginosa and BCG belong to two different phyla and differ significantly in their interactions with host cells. 781 long RNA species were DE (with respect to uninfected tissue) in infection by both bacterial pathogens, but not in IAV infection (Fig. 3B), and a four-way plot based on Wald statistic suggested that all genes that were DE by both bacterial pathogens were regulated in the same direction (Figure S5). This plot identified GJB2 as the gene that was most robustly upregulated in both bacterial infections. GJB2 encodes gap junction binding protein 2, a connexin which plays roles in epithelial barrier integrity (30) and is known for its importance in lung adenocarcinoma. Little is known about its role in inflammation or infections, but it is listed in GO:BP gene set RESPONSE_TO_BACTERIUM (GO:0009617). VNN1 encodes vanin 1, an oxidative stress sensor that is part of the tissue response to oxidative stress and inflammation in a variety of organs (31). The most specifically and strongly regulated RNAs were found in P. aeruginosa infection. In addition to GJB2 and VNN1, the most significantly (Wald test statistic) upregulated mRNAs were DUSP4 (encoding dual specificity phosphatase 4, a protein tyrosine/serine/threonine phosphatase involved in mitogenic signal transduction), IL10 (encoding a predominantly anti-inflammatory cytokine (32)), and SERPINB7 (encoding a member of the serpin family of protease inhibitors, which have anti-inflammatory and anticoagulant properties and counteract proteolytic tissue destruction (33)). The RNAs that were most significantly downregulated by both pathogens were PKMYT1 (encoding a kinase that inhibits G2/M transition in the cell cycle (34)), GGTA1P (encoding an enzymatically inactive alpha-1,3-galactosyltransferase involved in immune regulation (e.g., (35)), and S100A4 (encoding the S100 Ca++-binding protein S4, which furthers cell motility, migration, and invasion (36)). All three genes were downregulated more significantly by P. aeruginosa than by BCG. In spite of the strong agreement between P. aeruginosa and BCG in terms of directionality of DE, some degree of DE between the two bacterial infections was apparent. This is evidenced in the four-way plot (see RNAs marked in green [significant only in BCG infection] and violet [significant only in P. aeruginosa infection] in Figure S5 as well as in the volcano plot shown in Fig. 3H.
In hierarchical clustering analysis based on the 75 most significantly DE RNA selected by ANOVA across all four experimental groups (p values: 3.6E-06–7.8E-04), there was a remarkably clear separation in that the bacterial infections clustered in one clade and within it BCG and P. aeruginosa into separate subclades (Fig. 4). Separation of IAV infection from uninfected tissue was also excellent in that only one IAV-infected and one control sample were classified in the same subclade. Even though there were much fewer DE RNAs in IAV infection than in the bacterial infections, this hierarchical clustering analysis revealed a distinct “viral” signature, which consisted of 5 upregulated RNAs (marked red in Fig. 3I). Three of these (lnc-DDX60-1 [ENSG00000248601,] lnc-THOC3-2 [ENSG00000248596)], and RP11-202G18.1 [ENSG00000227531] are transcripts of unknown function which may arise during transcription of neighboring genes related to cellular responses to infection and inflammation. SLC8A1 encodes the Na+/Ca++ exchange channel NCX1, which is expressed in a variety of cell types and plays roles in mediating airway smooth muscle responses to inflammation (37) and thrombin-induced increased endothelial permeability (38). SNX10 (sorting nexin 10) is an important differentiation factor of myeloid cells and supports a proinflammatory phenotype of macrophages (39).
GSEA analysis based on Hallmark pathways revealed pathogen-specific regulation of functional pathways, but also common features shared among the three types of infection (Fig. 5). Several pathways were commonly upregulated by all three pathogens. These related to innate immune system signaling (IFNG response, IFNA response, inflammatory response, IL6_JAK_Stat3 signaling, TNFA signaling via NFKB), other aspects of innate immune activation (allograft rejection, complement activation), and apoptosis. One commonly downregulated pathway was cell cycle progression: E2F targets, indicating inhibition of cell proliferation by all three pathogens. The two bacterial infections could be distinguished from IAV infection by weaker induction of type I and type II IFN signaling, but stronger enrichment of signaling pathways such as IL2_Stat5_signaling, IL6_STAT3 signaling, programmed cell death, and KRAS signaling (upregulated genes). In contrast, depletion of myc targets, reflecting decreased cellular processes and cell cycle, and enrichment of KRAS signaling (downregulated genes) were unique for IAV infection. The pathways uniquely regulated in P. aeruginosa infection reflected a more vigorous induction of compensatory tissue responses such as formation of tight junctions, protein secretion, and metabolism of xenobiotics. Adipocyte development and p53 pathway were the only two pathways that were uniquely differentially regulated in BCG-infected tissue only. Taken together, these results provide further evidence that this HLTE model recapitulates classic tissue responses to infectious agents and reveals pathophysiologically plausible differences between the responses to viral vs. bacterial pathogens.
Cell type-specific responses to IAV infection. We then tested whether the HLTE could be used to dissect interactions between IAV and lung tissue at the single cell level and applied the 10x genomics single-cell RNA sequencing (scRNAseq) platform to single cell suspensions from IAV infected and uninfected HLTE from the same donor (a patient with indication for lung transplantation due to COPD) 24 h after IAV or mock infection. To correct for differences among individual HLTE pieces, single cell suspensions from 12 infected and 12 uninfected HLTE pieces were pooled to yield one infected and one uninfected sample for single cell analysis. After excluding dead cells by FACS, we loaded cells (n = 4500) for a targeted recovery of 3500 cells from each pooled sample at a sequencing depth to obtain 50 k reads per cell. Using cell-specific markers to identify cell clusters allowed identification of 11 cell types (Fig. 6A). These included the expected functional cell types of lower airways (airway epithelial cells, multiciliated cells, endothelial cells) and the major immune cell types: macrophages, monocytes, mast cells, B cells, CD4 + and CD8 + T cells, and NK cells. Mito_high T cells correspond to apoptotic T cells with high mitochondrial DNA content, indicating cell death, and were excluded from the host cell gene expression analyses. We then aimed to identify cells containing viral RNAs (Figure S7). Since IAV is a negative sense RNA virus, viral mRNA is transcribed from the negative sense (-) strand, and identification of a viral transcript (+ strand) therefore reflects viral transcription. However, by scRNAseq it is not easily possible to distinguish between intracellular viral mRNA synthesis in infected cells and uptake or binding of free viral mRNA from the extracellular environment. Of the eight viral genome segments, RNA corresponding to the segment encoding M1 and M2 proteins could not be detected for unknown reasons. Of the remaining seven segments, NS and HA mRNAs were the most abundant, which is consistent with the viral life cycle. Strongest viral mRNA expression per infected cell was seen in epithelial cells, the main host cell supporting productive IAV infection, followed by macrophages, which are known to support a nonproductive IAV infection and contribute to antiviral immune responses. However, there was clear evidence of viral transcripts also in association with the other cell types, whereby levels were lowest in endothelial cells, mast cells, and monocytes. Comparing DE between the same cell types in single cell suspensions from infected and uninfected HLTE enabled us to assess transcriptomic responses in the 10 cell types (Fig. 6B-G, Figure S8). CD4 T cells (72 DE genes), CD8 T cells (43), macrophages (42), NK cells (26), and mast cells (25) mounted the briskest transcriptional responses, whereas only weak responses were seen in multiciliated cells (0), B cells (1), endothelial cells (2), epithelial cells (1), and monocytes (1). Thus, the most robust antiviral host responses did not originate from the main target cells supporting viral transcription (epithelial cells) but from other cell types with less or no evidence of viral mRNA.
A comparison of DE genes identified by bulk RNAseq and scRNAseq showed that scRNAseq identified 19% of the DE genes that were identified by bulk RNAseq, but it also demonstrated that 89 DE genes were revealed by scRNAseq only (Fig. 7A). Of note, six genes (IFI44L, MX1, MX2, IRF7, IFI6, and ISG15) were identified by scRNAseq as DE in all major immune cells as well as by bulk RNAseq and thus constituted the core of the anti-IAV transcriptomic tissue response in this early time window of IAV infection (Fig. 7B).
Thus, this analysis supported the notion that IAV infects the expected cell types in this HLTE model, that antiviral and inflammatory responses are to some extent due to responses in noninfected “bystander” cells, and that scRNAseq may identify DE transcripts that are not detected as differentially expressed by bulk RNAseq.
We then assessed differences among the cell types in terms of functional responses to IAV infection (Fig. 7C). GSEA (Hallmark) revealed that type 1 and 2 IFN responses were induced in all cell types except B cells. CD4 + T cells and macrophages were the most responsive in that the most pathways (six in each cell type) were regulated. Notably, TNFA signaling induced by NFKB was enriched only in macrophages and monocytes. We then compared GSEA in each of the cell types to the GSEA performed at the tissue level (bulk RNAseq). Regulation of hallmark pathways (40) related to interferon response, inflammation, and MYC targets was detected at the tissue and the single-cell level, whereas five pathways were detected only by scRNAseq.