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.

Experimental design and functional validation of the human lung tissue explants (HLTE) model. A Outline of the experimental procedure. After removal from the explanted lung, tissue was divided into pieces of approx. 30 mg and then incubated in medium overnight. HLTE pieces were infected with 2.0 × 105 FFU/ml of influenza virus strain A/Giessen /6/2009 H1N1 (in short IAV), 5 × 106 CFU/ml M. bovis strain H37Rv (BCG), or 1 × 108/ml P. aeruginosa strain PA14. All analyses were performed 24 h post infection (p.i.) unless indicated otherwise. B Comparison of two methods of tissue preservation before RNA extraction. HLTE pieces were either snap frozen in liquid nitrogen (n = 20) or incubated in RNAlater at 4 °C overnight (n = 38) before storage at − 80 °C, followed by RNA extraction. Y-axis = RNA integrity number (RIN). C Time course of LDH release from uninfected or IAV-infected HLTE pieces (n = 3 per group). LDH was measured in tissue culture supernatants and is expressed as % of LDH extracted from lysed tissue control. D Time course of IAV hemagglutinin (HA) mRNA levels (RT-qPCR), indicating transcription of viral RNA (n = 6). E Differences in cytokine/chemokine induction by infection with IAV, BCG, and P. aeruginosa. Upper row: protein concentrations measured by EIA. Bottom row: mRNA levels measured by RT-qPCR relative to mock-infected HLTE, using GAPDH as internal reference. n = 9. *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001; ****p ≤ 0.0001; B–D = unpaired parametric T test; E = Mann–Whitney U test. Data represent means ± SEM

Bulk RNA sequencing. Using HLTEs from five separate donors, we performed 5 independent infection experiments, each featuring infection with all three pathogens in parallel. We then applied RNAseq to measure expression of sncRNA, long noncoding, pseudogene, and protein-coding RNA in uninfected HLTEs and in HLTEs infected with IAV, BCG, or P. aeruginosa. To reduce the impact of intra-donor variability observed in the IAV infection time course shown in Fig. 1D, we pooled RNA from 2 to 4 HLTE pieces from each donor so that five pooled RNA samples per condition were available for RNAseq (final n = 20).

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.

Modest reprogramming in sncRNA expression in HLTE in viral and bacterial infections. sncRNA populations were determined by sncRNAseq 24 h after infection with IAV, BCG, and P. aeruginosa, or after mock infection (n = 5 per condition). Results of the DE analysis with DEseq2 are found in Table S1. A Detection efficiency (no. of reads) of all sncRNAs vs. sncRNAs mapping to human genome hg38. B Abundance of the major sncRNA subpopulations miRNA, piRNA, snoRNA, and snRNA. C Abundance of sncRNA subtypes detected at a total of ≥ 20 reads in all samples. D Abundance of differentially expressed sncRNA subtypes (pAdj. ≤ 0.1). E PCA showing poor separation of the four groups based on sncRNA. F–H Volcano plots showing DE sncRNAs (pAdj ≤ 0.1, fold change [FC], ≥|2|). F IAV vs. uninfected tissue. G BCG vs. uninfected tissue. H P. aeruginosa vs. uninfected tissue. I, Hallmark GSEA analysis of predicted targets of DE miRNA in P. aeruginosa infection

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 × 106 (± 7.3 × 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).

Differential reprogramming of RNA populations by viral and bacterial infection. Using RNA from the same total RNA samples as used for small RNA sequencing (Fig. 2), long RNA populations were determined by RNAseq 24 h after infection or mock treatment (n = 5 per condition). A PCA indicating somewhat better separation than with sncRNAseq (compare Fig. 2E). B Venn diagram showing the number of DE mRNAs (pAdj ≤ 0.1, FC ≥|2| with respect to uninfected tissue) unique to each pathogen, shared between any two pathogens, or common to all three. C–H Volcano plots identifying DE mRNA and lncRNA species with respect to uninfected tissue in infection with IAV (C–E), BCG (F), and P. aeruginosa (G), and comparing P. aeruginosa vs. BCG infected HLTEs (H). Cutoffs were set at FC =|2| (i.e. log2 =|1|) and p-Adj 0.1 (i.e. − log10 = 1)

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. 760 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 to 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. 4). 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].

Viral and bacterial infections leave different RNA signatures in HLTEs. Unsupervised hierarchical biclustering analysis based on the 75 most significantly up- or downregulated mRNA and lncRNA (p < 7.8E−04, ANOVA across all four groups). Red bracket: RNAs comprising a 5-transcript “viral signature”. The color scheme in the 3 columns on the left shows p-Adj values for DE due to infection with each pathogen compared with uninfected tissue

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.

Gene set enrichment analysis (GSEA) based on Hallmark pathways identifies common and distinct signaling pathways in bacterial and viral infection of explanted human lung tissue. mRNA expression data were extracted from the bulk RNAseq data set used for Fig. 3, and a GSEA was performed using the prerank gene list from Deseq2 analysis against the Hallmark gene set collection. Significantly enriched or depleted Hallmark pathways are indicated by the triangles as detailed in the legend next to the figure. A GSEA based on Gene Ontology Biological Process (GO:BP) is shown in Figure S6

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 10 × 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. Extensive expression of markers for tissue resident lymphocytes, CD69 and CD103 [40], indicated that most of the captured lymphocytes were indeed tissue resident and not remnants of peripheral blood cells that had not been removed by washing the tissue (Fig. 6B, Figure S7). We then aimed to identify cells containing viral RNAs (Figure S8). 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. 6C, Figure S9). CD4 T cells (72 DE genes), CD8 T cells (43 DE genes), macrophages (42 DE genes), NK cells (26 DE genes), and mast cells (25 DE genes) 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.

Single-cell RNA sequencing identifies strongest induction of antiviral host responses in CD4, CD8 + T, NK cells, mast cells, and macrophages. Single cell transcriptomes were determined in pooled single cell suspensions derived from HLTEs from one donor with emphysema due to COPD 24 h after IAV infection or mock infection. A Identification of 11 cell subpopulations by uniform manifold approximation and projection (UMAP). B CD69 mRNA expression (see Figure S7 for CD103 expression). C Volcano plots identifying DE genes in the major immune cell types and vascular ECs. The dotted lines identify cutoffs FC =|2| (i.e. log2 =|1|) and p-Adj 0.1 (i.e. − log10 = 1). DE in the other cell types is shown in Figure S9

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

Comparison of DE Genes and functional pathways in IAV-infected HLTEs detected at the tissue or the single cell level 9. A, B Venn diagrams illustrating the relationships between DE genes identified by bulk RNAseq and by scRNAseq in different cell types. A Bulk RNAseq vs. scRNAseq (pool of all cell types); B bulk RNAseq vs. scRNAseq (individual cell types indicated in the diagram). The arrow points to 6 mRNAs making up the core tissue response identified by bulk and scRNAseq. C GSEA of IAV infection based on the bulk RNAseq data or scRNAseq data. Hallmark pathways [41] are listed on the vertical axis in reverse alphabetical order (top to bottom). Significantly (FDR ≤ 0.05) enriched or depleted pathways are marked with triangles as indicated in the legend on the right

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 via 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 [41] 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.