Allergic inflammatory memory in human respiratory epithelial progenitor cells

Abstract

Barrier tissue dysfunction is a fundamental feature of chronic human inflammatory diseases¹. Specialized subsets of epithelial cells-including secretory and ciliated cells-differentiate from basal stem cells to collectively protect the upper airway^2-4. Allergic inflammation can develop from persistent activation⁵ of type 2 immunity⁶ in the upper airway, resulting in chronic rhinosinusitis, which ranges in severity from rhinitis to severe nasal polyps⁷. Basal cell hyperplasia is a hallmark of severe disease^7-9, but it is not known how these progenitor cells^2,10,11 contribute to clinical presentation and barrier tissue dysfunction in humans. Here we profile primary human surgical chronic rhinosinusitis samples (18,036 cells, n = 12) that span the disease spectrum using Seq-Well for massively parallel single-cell RNA sequencing¹², report transcriptomes for human respiratory epithelial, immune and stromal cell types and subsets from a type 2 inflammatory disease, and map key mediators. By comparison with nasal scrapings (18,704 cells, n = 9), we define signatures of core, healthy, inflamed and polyp secretory cells. We reveal marked differences between the epithelial compartments of the non-polyp and polyp cellular ecosystems, identifying and validating a global reduction in cellular diversity of polyps characterized by basal cell hyperplasia, concomitant decreases in glandular cells, and phenotypic shifts in secretory cell antimicrobial expression. We detect an aberrant basal progenitor differentiation trajectory in polyps, and propose cell-intrinsic¹³, epigenetic^14,15 and extrinsic factors^11,16,17 that lock polyp basal cells into this uncommitted state. Finally, we functionally demonstrate that ex vivo cultured basal cells retain intrinsic memory of IL-4/IL-13 exposure, and test the potential for clinical blockade of the IL-4 receptor α-subunit to modify basal and secretory cell states in vivo. Overall, we find that reduced epithelial diversity stemming from functional shifts in basal cells is a key characteristic of type 2 immune-mediated barrier tissue dysfunction. Our results demonstrate that epithelial stem cells may contribute to the persistence of human disease by serving as repositories for allergic memories.

Extended Data Figure 1 |. Consistency of cell capture and identification in surgical EthSin scRNA-seq patient cohort

(a) Number of unique molecular identifiers (nUMI) and genes identified, and fraction of reads mapping to mitochondrial or ribosomal genes across recovered cell types; 3,222 basal cells, 4,362 apical cells, 2,192 glandular cells, 498 ciliated cells, 835 T cells, 2,976 plasma cells, 1,724 fibroblasts, 1,143 endothelial cells, 811 myeloid cells, 273 mast cells. (b) nUMI and genes identified, and fraction of reads mapping to mitochondrial or ribosomal genes across patient samples; 789 Polyp 1 cells, 1,309 Polyp 2 cells, 1,153 Polyp 3 cells, 913 Polyp 4 cells, 1,219 Polyp 5 cells, 1,141 Polyp 6A cells, 1,334 Polyp 6B cells, 1,314 Polyp 7 cells, 1,286 Polyp 8 cells, 1,481 Polyp 9 cells, 2,988 Polyp 11 cells, 3,109 Polyp 12 cells. (c) tSNE plot as in Fig. 1b colored by cell types across all patients and then separated by sample; 18,036 single cells (n=12 samples). (d) The percentage of each cell type recovered within each sample. (e) Select marker gene overlays displaying binned count-based UMI-collapsed expression level (log(scaled UMI+1)) on a tSNE plot from Fig. 1b for key cell types identified (see Supplementary Table 3 for full gene lists); area under the curve (AUC) 0.998 to 0.7 for all markers displayed.

Extended Data Figure 2 |. Top marker genes for cell types by scRNA-seq and bulk tissue RNA-seq from EthSin recovers expected T2I and eosinophilic modules

(a) Row-normalized heatmap of the top-10 marker genes identified by ROC-test (AUC>0.73 for all) over all cell types (Fig. 1b,c) with select genes displayed on y-axis and cells on x-axis (see Supplementary Table 3 for full gene lists); maximum 500 cells/type. (b) An overlay of CLC displaying binned count-based expression level (log(scaled UMI+1)) amongst myeloid cells (a pathognomonic gene for eosinophils); 811 myeloid cells from n=12 samples. (c) A row-normalized and row-clustered heatmap over the top 100 positively and negatively differentially-expressed genes (50 in each direction) in bulk tissue RNA-seq of 27 samples from non-polyp (n=10) and polyp (n=17) tissue with select genes displayed; DESeq2 Wald Test, all p<9.03×10⁻⁵ for genes displayed, corrected for multiple comparisons by Benjamini procedure, samples ordered as in Fig. 3g, (see Supplementary Table 3 for full gene list and associated statistics). (d) The top differentially regulated pathways identified by Ingenuity Pathway Analysis (see Methods) over the top 1,000 differentially expressed genes, as determined by p<0.05 corrected for multiple comparisons by Benjamini procedure, across polyp and non-polyp tissue. (e) Predicted upstream regulators based on differentially expressed gene modules in polyp tissue relative to non-polyp determined using Ingenuity Pathway Analysis (see Methods).

Extended Data Figure 3 |. Sub-clustering of myeloid, fibroblast and endothelial cell types from the EthSin T2I inflammatory ecosystem

(a) tSNE plot of 811 myeloid cells (n=6 non-polyp, n=6 polyp samples), colored by clusters identified through shared nearest neighbor (SNN) analysis (Supplementary Table 3; Methods), from CRS-EthSin; select marker gene overlays displaying count-based (unique molecular identifier (UMI) collapsed) expression level (log(scaled UMI+1)) on a tSNE plot (see Supplementary Table 3 for full gene lists; genes identified via ROC test with AUC 0.689 for S100A8, 0.763 for CD1C, 0.927 for C1QC); a clustered correlation matrix of marker genes identified in single-cell data from myeloid cells; and violin plots for the expression value (log(scaled UMI+1)) of selected markers of myeloid activation state. (b) tSNE plot of 1,724 fibroblasts (n=6 non-polyp, n=6 polyp samples), colored by clusters identified through shared nearest neighbor (SNN) analysis (Supplementary Table 3; Methods), from CRS-EthSin; select marker gene overlays displaying count-based (unique molecular identifier (UMI) collapsed) expression level (log(scaled UMI+1)) on a tSNE plot (see Supplementary Table 3 for full gene lists; genes identified via ROC test with AUC 0.691 for CTGF, 0.683 for CXCL12, 0.726 for MYH11); and a clustered correlation matrix of marker genes identified in single-cell data from fibroblasts. NB: Clusters 4 and 5 likely represent doublets with epithelial cells and endothelial cells, respectively. While we exclude these clusters from further formal analyses, we note that there may be interesting biology within pairs of cells found to interact more frequently than by chance. (c) tSNE plot of 1,143 endothelial cells (n=6 non-polyp, n=6 polyp samples), colored by clusters identified through shared nearest neighbor (SNN) analysis (Supplementary Table 3; Methods), from CRS-EthSin; select marker gene overlays displaying count-based (unique molecular identifier (UMI) collapsed) expression level (log(scaled UMI+1)) on a tSNE plot (see Supplementary Table 3 for full gene lists; genes identified via ROC test with AUC 0.742 for SELE, 0.706 for PODXL, 0.822 for PLAT); and a clustered correlation matrix of marker genes identified in single-cell data from endothelial cells.

Extended Data Figure 4 |. Mapping T2I mediators within EthSin non-polyp or polyp ecosystems and the identities of T cells

(a) Dot plots of chemokines and lipid mediators with known roles in T2I mapped onto cell types divided by non-polyp or polyp disease state. Dot size represents fraction of cells within that type expressing, and color intensity binned (log(scaled UMI+1)) gene expression amongst expressing cells (related to Figure 1d). (b) Dot plot of inducers and effectors of T2I mapped onto cell types divided by non-polyp or polyp disease state. Dot size represents fraction of cells within that type expressing, and color intensity binned (log(scaled UMI+1)) gene expression amongst expressing cells (related to Figure 1d). (c) tSNE plot of re-clustered T cells with select gene overlays displaying binned count-based expression level (log(scaled UMI+1)) for Th2A-specific genes (top row) and canonical T cell markers (bottom row); 835 T cells from n=6 non-polyp and n=6 polyp samples. (d) Violin plot of five identified T cell clusters scored for expression of T cell receptor complex genes (e.g. TRAC and CD3E, see Methods, Supplementary Table 4); dots represent individual cells, 835 total T cells. (e) Dot plot of inducers and effectors of Type 1 immunity across all cell types (NB: IL17F not detected). (f) Dot plot of select GWAS risk alleles for allergic disease, mapped onto cell types divided by non-polyp or polyp disease state. Dot size represents fraction of cells within that type expressing, and color intensity binned (log(scaled UMI+1)) gene expression amongst expressing cells (related to Figure 1d).

Extended Data Figure 5 |. Relationship of EthSin epithelial cell clusters and secretory/glandular distinctions

(a) A phylogenetic tree based on the average cell from each cluster of epithelial cell clusters in gene-space. (b) Violin plot of expression contribution to a cell’s transcriptome of basal cell genes (see Methods and Supplementary Table 4) across all epithelial cells; 794 cells cluster 12; 924 cells cluster 8; 1,504 cells cluster 2; 1,561 cells cluster 1; 1,600 cells cluster 0; 1,201 cells cluster 4; 725 cells cluster 13; 1,467 cells cluster 3; 498 cells cluster 16; *Mann-Whitney U-test, with Bonferroni correction, p<1.76×10⁻¹⁵, 12, 8 or 2 vs. the mean score of basal/apical epithelial cells; p= 0.5392, 1 vs. the mean score. (c) Canonical correlation analysis (CCA) displaying our cell type annotations for basal and apical cells derived through clustering and biological curation alongside CCA clusters in tSNE space; 7,584 basal and apical cells. (d) Violin plots for the count-based expression level (log(scaled UMI+1)) of selected marker genes for each identified epithelial cell subset; cell numbers as in (b). (e) Row-normalized heatmap of the top marker genes identified by ROC-test (AUC>0.6) within each cell type for each cell cluster with genes displayed on y-axis and cluster annotations on x-axis (see Supplementary Table 3 for full gene lists). (f) Select overlays on clusters 0 and 4 (differentiating/secretory) and 13 (glandular) displaying binned count-based expression level (log(scaled UMI+1)) in tSNE space for canonical goblet (MUC5B, MUC5AC, SPDEF, FOXA3) and secretory (SCGB1A1) genes; 3,526 cells. (g) A clustered correlation matrix of glandular, goblet, and secretory cell genes; Pearson’s abs(r)> 0.038 is p<0.05 significant based on asymptotic p-values.

Extended Data Figure 6 |. Glandular cell subsets, their relationship to apical secretory cells, and immune cells recovered through nasal scrapings

(a) tSNE plots of 5,928 single epithelial cells (n=6 non-polyp samples), and 4,346 single epithelial cells (n=6 polyp samples) colored by clusters identified through (left) shared nearest neighbor (SNN) analysis and (right) original biological curation of cell types (Supplementary Table 3; Methods) as illustrated in Figure 2a. NB: cluster colors in left panels of each disease are not comparable but curated clusters in right are, and glandular cells are highlighted for subsetting in next panel. (b) Violin plots for the count-based expression level (log(scaled UMI+1)) of selected marker genes identified through marker discovery (ROC test) for each subset of glandular cells; 2,114 total cells (791 cells cluster 1; 709 cells LCN2 cluster; 283 cells SERPINB3 cluster; 209 cells MUC5B cluster; 183 cells PRB1 cluster) with representation of every non-polyp patient in each cluster of cells (e.g. no cluster is unique to one patient) and AUC metric 0.800 for LCN2, 0.736 for SERPINB3, 0.985 for MUC5B, 0.973 for BPIFB2, and 0.908 for PRB1. (c) Samples were acquired through the two distinct methods of nasal scraping and ethmoid sinus surgical intervention. This allowed for sampling of left: healthy tissue from InfTurb (scraping), middle: of CRS-EthSin-non-polyp tissue (surgery), right: of CRS-EthSin-polyp tissue (surgery), of InfTurb of polyp-bearing individuals (scraping), and of CRS-EthSin-polyp tissue accessible for scraping (scraping). Anatomy of the nasal turbinates (healthy and CRS polyp) and ethmoid sinus (CRS non-polyp and CRS polyp) where samples were acquired is displayed below, highlighting the depth of cells recovered from each site related to Fig. 2. Healthy tissue is annotated with basal and apical cell types, including sub-mucosal glands. (d) Left: tSNE plot of 18,704 single cells from nasal scrapings (n=9 samples) colored by clusters identified through shared nearest neighbor (SNN) analysis (Supplementary Table 3; Methods); middle: tSNE plot colored by cell types identified through marker discovery (ROC test) and biological curation of identified clusters (Supplementary Table 3; Methods); right: colored by disease and tissue of origin from healthy InfTurb (7,603 cells; n=3 samples), polyp-bearing patient InfTurb (2,298 cells; n=4 samples), and polyp scraping directly from EthSin-polyp (8,803 cells; n=2 samples); with adjacent select marker gene overlays displaying count-based UMI-collapsed expression level (log(scaled UMI+1)) for apical epithelial (KRT8) and hematopoietic (PTPRC) genes. (e) Select marker gene overlays displaying count-based UMI-collapsed expression level (log(scaled UMI+1)) on a tSNE plot from (a) for key cell types identified (see Supplementary Table 3 for full gene lists); area under the curve (AUC) 0.946 to 0.705 for all markers displayed. (f) Violin plots for the count-based expression level (log(scaled UMI+1)) for key differentially expressed genes using ROC test within myeloid cells across disease states and tissues identified (Methods); 137 cells, n=3 healthy inferior turbinate; 157 cells, n=4 polyp inferior turbinate; 210 cells, n=2 polyp ethmoid sinus samples; AUC 0.67 for TXNRD1, 0.615 for RALA, 0.647 for TLR2, 0.619 for RIPK2, 0.747 for C1QA, 0.674 for FGL2.

Extended Data Figure 7 |. Changes in cellular composition between EthSin-non-polyp and EthSin-polyp tissue by scRNA-seq and flow cytometric gating and histological strategy for quantification and isolation of basal cells

(a) The frequency of each cell type recovered amongst all cells within each patient sample (n=6 non-polyp, n=6 polyp) grouped by disease state; *t-test, two-sided, p=0.0003 for apical, p<0.0001 for glandular, p=0.0047 for ciliated, p=0.00014 for plasma cell, p=0.0098 for myeloid and p=0.00018 for mast cell; all non-polyp vs. polyp with Holm-Sidak correction for multiple comparisons; mean±s.e.m. (b) The frequency of basal cells amongst epithelial cells captured in scRNA-seq data displayed for each sample and colored by non-polyp or polyp designation. (c) tSNE plots with each patient’s cells clustered independently over a common list of most variable genes identified from all epithelial cells and with clustering parameters set constant to 12 principal components and resolution set to 1.4; minimum 789 cells in each plot, Extended Data Fig. 1b and Supplementary Table 3 for specific cell numbers. (d) Simpson’s index of diversity over epithelial cell clusters identified in (c), an indication of the total richness present within an ecosystem, calculated for each patient; n=6 non-polyp and n=6 polyp samples; *t-test, two-tailed, p=0.0384, mean±s.e.m. (e) Correlation of Simpson’s index of diversity calculated over epithelial cells against the ranked order of samples based on clinical pathological evaluation; n=6 non-polyp and n=6 polyp samples; r=0.6824, p=0.009. (f) Simpson’s index of diversity over stromal and immune cell types and total cells, an indication of the total richness present within an ecosystem, calculated for each sample (n=6 non-polyp and n=6 polyp); points represent individual samples, *t-test, two-tailed, p=0.0015 stromal and immune, p=0.0145 total cells, non-polp vs. polyp; mean±s.e.m. (g) Reproduced from Figure 2a: tSNE plot of 10,274 epithelial cells, colored by clusters identified through SNN, with adjacent color bars representing related cell clusters, and overlays displaying binned count-based expression level (log(scaled UMI+1)) of selected genes used to negatively (CD45, EPCAM, THY1) and positively (NGFR, ITGA6, PDPN) identify basal cells. (h) Full flow cytometric gating strategy for quantification and isolation of basal cells from non-polyp and polyp tissue, (related to Fig. 3c). (i) Representative histology (5x magnification) of the glandular area detected in haematoxylin and eosin stained tissue sections from non-polyp or polyp patients; quantification in Fig. 3e. (j) Representative immunofluorescence of p63-staining cells (basal cell marker) relative to isotype control; quantification in Fig. 3d; scale bar 100μm. (k) Basal cell fraction of transcripts from bulk tissue RNA-seq data of our own data set (related to Fig 3g,h) and two GEO data sets containing healthy and healthy/polyp nasal mucosa biopsies; our data n=10 non-polyp samples, n=17 polyp samples; reference data n=6 healthy, n=6 polyp samples; *t-test, two-tailed p=0.0465 our data and p=0.0040 GEO data; mean±s.e.m. (l) Secretory cell fraction of transcripts from bulk tissue RNA-seq data of our own data set (related to Fig 3g,h) and two GEO data sets containing healthy and healthy/polyp nasal mucosa biopsies; our data n=10 non-polyp samples, n=17 polyp samples; reference data n=6 healthy, n=6 polyp samples; *t-test, two-tailed p=0.0465 our data and p=0.0040 GEO data; mean±s.e.m

Extended Data Figure 8 |. Epithelial cytokine signatures from CRS-EthSin tissue demonstrate T2I pattern, discovery of gene modules in the fibroblast niche which correlate with basal cell hyperplasia, and differential expression within myeloid and endothelial cells by polyp status

(a) Violin plots of IL-4- or IL-13-uniquely induced gene signatures in respiratory epithelial cell clusters or grouped by disease state presented as expression contribution to a cell’s transcriptome (see Methods, Figure 4b for shared genes, and Supplementary Table 4); 794 cells cluster 12; 924 cells cluster 8; 1,504 cells cluster 2; 1,561 cells cluster 1; 1,600 cells cluster 0; 1,201 cells cluster 4; 725 cells cluster 13; 1,467 cells cluster 3; 498 cells cluster 16; *Mann-Whitney U-test, p<2.2×10⁻¹⁶, 0.305 effect size IL-4 polyp vs. non-polyp and −0.448 effect size IL-13 polyp vs non-polyp. (b) Violin plots of IFNα- or IFNγ-induced gene signatures in respiratory epithelial cell clusters or grouped by disease state presented as expression contribution to a cell’s transcriptome (see Methods, and Supplementary Table 4); cell numbers as in (a); *Mann-Whitney U-test, p=4.98×10⁻⁶, −0.156 effect size IFNα polyp vs. non-polyp; *Mann-Whitney U-test, p<2.2×10⁻¹⁶, 0.161 effect size IFNγ polyp vs non-polyp. (c) Selected genes detected in fibroblasts from single-cell data which correlate with the samples ranked by basal cell frequency detected within each ecosystem; n=6 non-polyp, 6 polyp samples, all genes used: Spearman correlation, abs(r)>0.7651, p<0.0037. NB: to determine genes correlated in specific cell types (e.g. fibroblasts) with the frequency of basal cells present in a cellular ecosystem, we correlated the average log-normalized single-cell count data for each gene to the rank of samples determined by increasing frequency of basal cells in each ecosystem (8.2% to 19.1% for non-polyp and 27.9% to 70.1% for polyp samples, Extended Data Fig. 7b). (d) A clustered correlation matrix of genes identified as per (c) in single-cell data from fibroblasts; Pearson’s abs(r)>0.048 is p<0.05 significant based on asymptotic p-values. (e) Row-normalized heatmap for myeloid cells from ethmoid sinus with select genes displayed on y-axis including a core myeloid signature (ROC-test myeloid cells vs. rest of cells, AUC>0.8), and genes differentially expressed (bimodal test) by disease state, with disease state annotations on x-axis; bimodal test, all non-core genes p<0.0002 or less with Bonferroni correction for multiple hypothesis testing based on number of genes tested. (f) Row-normalized heatmap for endothelial cells from ethmoid sinus with select genes displayed on y-axis including a core basal signature (ROC-test endothelial cells vs. rest of cells, AUC>0.75), and genes differentially expressed (bimodal test) by disease state, with disease state annotations on x-axis; bimodal test, all non-core genes p<2.43×10⁻⁶ or less with Bonferroni correction for multiple hypothesis testing based on number of genes tested.

Extended Data Figure 9 |. Pseudotime analysis on basal and differentiating/secretory cell clusters from EthSin, transcriptional motif enrichments in non-polyp and polyp basal cells, and the identity of cell types in air-liquid interface cultures

(a) Pseudotime analysis using diffusion mapping (see Methods) of selected clusters of epithelial cells, here displaying diffusion pseudotime (related to Fig. 4d); 3,516 cells (clusters 8/1/4); 4,064 cells (clusters 12/2/0); and n=6 non-polyp, n=6 polyp samples; diffusion map and diffusion coefficients (DC) are calculated over the set of basal and apical marker genes identified in Fig. 1a, see Supplementary Table 3. (b) The top 60 negatively correlated genes expressed in non-polyp cells with pseudotime trajectory and Pearson correlation values for genes in polyp cells also displayed; differential correlation coefficient analysis using Fisher’s Z-statistic, accounting for number of cells in each group (specific genes highlighted all > 2 Z, full results including Bonferroni corrected p-values in Supplementary Table 3). (c) Correlation matrices (row and column clustered) of the normalized read counts per sample in motif associated-peaks for non-polyp or polyp samples; Pearson correlation, n=3 non-polyp, n=7 polyp. (d) A column-normalized heatmap (row and column clustered) for the fraction of peaks with a motif corresponding to accessibility of the respective transcription factor displayed by patient; n=3 non-polyp, n=7 polyp. (e) IGV tracks for ATF3 and KLF5 based on peaks detected and averaged by non-polyp and polyp samples from ATAC-seq profiling. (f) IGV tracks for S100A9 and MUC4 based on peaks detected and averaged by non-polyp and polyp samples from ATAC-seq profiling. (g) Violin plots for the count-based expression level (log(scaled UMI+1)) for key marker genes using ROC test across cell types identified in (Fig. 5a; Supplementary Table 3); 1,345 basal; 6,420 secretory; 6,381 hybrid; and 2,027 ciliated cells; from n=2 non-polyp and 2 polyp patients; AUC 0.943 for KRT5, 0.667 for TP63, 0.644 for LYPD2, <0.55 for SPDEF, <0.55 for KRT8, 0.602 for BPIFA1, 0.813 for PIFO, 0.73 for FOXJ1. (h) Row-normalized heatmap for ALI-secretory cells (subsampled to 300 cells per donor) as in Fig. 2f of the top in vivo secretory marker genes identified by ROC-test (AUC>0.662) with select genes displayed on y-axis including a core secretory signature (ROC-test secretory cells vs. rest of cells), and then within secretory cells a ROC-test used to identify marker genes within each disease/location category, and basal-cell derived annotations on x-axis (see Supplementary Table 3 for full gene lists, all AUC>0.65 for markers displayed in Fig. 2f). (i) Quantification of flow cytometry for the ratio of basal to Epcam^hi cells (gating as in Extended Data Fig. 7h) from ALI cultures at 21 days stimulated with IL-13 over the indicated doses; points represent individual biological replicates; n=6 non-polyp, n=5 polyp samples for each dose; *2-way ANOVA, n.s. between disease groups at any dose tested; *2-way ANOVA, p=0.0224 for IL-13 dose; mean±s.e.m. (j) Expression levels for CTGF (Log2 expression value of log-normalized count data) in basal cells from non-polyp or polyp individuals across doses of cytokines displayed; n=4 samples each dose; 2-way ANOVA p<0.0260 for CTGF; all conditions non-polyp vs. polyp except 0.1 ng/mL IL-4 dose for CTGF.

Extended Data Figure 10 |. In vivo blockade with an anti-IL-4Rα monoclonal antibody shifts secretory cell state towards healthy-associated genes

(a) left: tSNE plot of 8,764 single cells (related to Fig. 6g) from the nasal polyps of an anti-IL-4Rα (dupilumab) treated individual (1 patient, sampled at n=3 timepoints) colored by clusters identified through SNN analysis (Supplementary Table 3; Methods); middle: tSNE plot colored by timepoint and tissue of origin from polyp pre-dupilumab scraping (5,731 cells), from polyp post-dupilumab scraping (647 cells), and polyp post-dupilumab surgical sample (2,386 cells); and right: tSNE plot colored by cell types identified through marker discovery (ROC test) and biological curation of identified clusters (Supplementary Table 3; Methods)). (b) Select cell-type specific score overlays for cell types indicated in original core data set (see Supplementary Table 3 for full gene list). (c) Row-normalized heatmap for myeloid cells of the top marker genes identified by ROC-test (AUC>0.8) with select genes displayed on y-axis including a core myeloid signature (ROC-test myeloid cells vs. rest of cells), and then genes found to be differentially expressed from (Fig. 5f) in basal cells, and treatment annotations on x-axis; bimodal test, * denotes differential genes in both basal cells and myeloid cells pre- vs post-treatment p<0.003 or less with Bonferroni correction for multiple hypothesis testing based on number of genes tested. (d) Violin plots for basal cells (200 cells pre-dupilumab and 151 cells post-dupilumab, noted in (a)) for the count-based expression level (log(scaled UMI+1)), except where indicated for gene scores, fraction of transcriptome and z-score (see Methods, Supplementary Table 4 for gene set used) for key basal cell genes for selected biological processes, or from the baseline upregulated gene set from polyp basal cells in vitro (Fig. 5c); differential expression testing for decreased expression post-treatment using bimodal test n.s. unless denoted by * for p<0.00087 or less with Bonferroni correction for multiple hypothesis testing based on number of genes tested; see Supplementary Table 3 for full list; Basal in vitro score Pre vs Post: *t-test, two-tailed, p<3.897×10⁻¹⁵, effect size 0.822. (e) tSNE plot of 4,486 single cells (related to Fig. 2e, and Fig. 5e) from the inferior turbinate or nasal polyps of an anti-IL-4Rα (dupilumab) treated individual (n=4 samples) colored by timepoint and tissue of origin from inferior turbinate pre-dupilumab scraping (643 cells), from inferior turbinate post-dupilumab scraping (1,596 cells), polyp pre-dupilumab scraping (1,600 cells), and polyp post-dupilumab scraping (647 cells); and tSNE plot colored by cell types identified through marker discovery (ROC test) and biological curation of identified clusters (Supplementary Table 3; Methods)); black outline indicates cells considered in (g). (f) Select deconvolution score overlays for cell types indicated in original core data set (see Supplementary Table 3 for full gene list). (g) Violin plot for the gene set score over Wnt pathway (z-score) and expression contribution to a cell’s transcriptome over IFNα- and IL-4/IL-13-commonly induced gene signature in secretory cells grouped as in (e) and sub-sampled to a maximum of 150 cells from each disease/location category from inferior turbinate pre-dupilumab scraping (150 cells), from inferior turbinate post-dupilumab scraping (23 cells), polyp pre-dupilumab scraping (150 cells), and polyp post-dupilumab scraping (38 cells); see Methods, Supplementary Table 3, Supplementary Table 4 for gene lists used; *t-test, two-tailed, Wnt score Pre vs Post Polyp Tissue: effect size 1.02, p=1.091×10⁻¹⁴; Wnt score Pre vs Post Inferior Turbinate Tissue: effect size −0.17, p=0.3706; IL-4/IL-13 score Pre vs Post Polyp Tissue: effect size 1.17, p<2.2×10⁻¹⁶; IL-4/IL-13 score Pre vs Post Inferior Turbinate Tissue: effect size −0.17, p=0.163; IFNα score Pre vs Post Polyp Tissue: effect size −1.25, p=4.254×10⁻⁰⁵; IFNα score Pre vs Post Inferior Turbinate Tissue: effect size −0.304, p=0.2766; differential expression testing for decreased expression post-treatment using bimodal test denoted by * and p<7.81×10⁻⁰⁶ or less between pre- and post-treated polyp. (h) Violin plots of secretory cells grouped as in (e) and sub-sampled to a maximum of 150 cells from each disease/location category from inferior turbinate pre-dupilumab scraping (150 cells), inferior turbinate post-dupilumab scraping (23 cells), polyp pre-dupilumab scraping (150 cells), and polyp post-dupilumab scraping (38 cells) for the count-based expression level (log(scaled UMI+1)) and for secretory cell genes from the gene set used in Fig. 2f affected by treatment within anatomical regions indicated by heading; differential expression testing for decreased expression post-treatment using bimodal test n.s. unless denoted by *, all p<6.36×10⁻⁵ or less except KLF5 (p=0.0033) and FOSB (p=0.0053) with Bonferroni correction for multiple hypothesis testing based on number of genes tested, see Supplementary Table 3 for all genes tested.

Figure 1 |. Mapping the T2I inflamed human sinus cellular ecosystem by scRNA-seq

a, Clinical disease spectrum (n=12 samples) and experimental workflow leading to a tSNE plot displaying 18,036 single cells, colored by shared nearest neighbor (SNN) clusters and (b) cell types (ROC-test; Supplementary Table 3; Methods) from respiratory tissue. c, Heatmap of top-10 marker genes by ROC-test (AUC>0.73) for indicated cell types; maximum 500 cells/type (Extended Data Fig. 2a annotated; Supplementary Table 3 full gene list). d, Dot plot of T2I mediators mapped onto cell types across all samples (Extended Data Fig. 4 by disease state).

Figure 2 |. Single-cell transcriptomes of epithelial cells in T2I highlight shifts in secretory cell states across health and disease

a, tSNE plot of 10,274 epithelial cells (n=12 samples), colored by SNN-clusters (Fig. 1; Extended Data Fig. 6 re-clustered) with blue color bars representing cell types determined per Extended Data Fig. 5, and (b) heatmap of marker genes by ROC (AUC>0.65; Supplementary Table 3, full list). c,d, tSNE plot (a) colored by disease (n=6 non-polyp, n=6 polyp samples) and (d) violin plots (Methods: all violins generated using standard Seurat implementation with default smoothing, density generated at >25% positive values, widest aspect centre of positive measures, minima/maxima within scale representing all points) for differentially expressed genes across disease state in differentiating/secretory cells; 2,566 cells, n=6 non-polyp; and 1,796 cells, n=6 polyp samples; *bimodal test, all p<2.03×10⁻⁵⁵ or less with Bonferroni-correction (Supplementary Table 3, exact values). e, tSNE plot of 18,325 re-clustered single cells from merged nasal scrapings (n=9) and surgical samples (n=12) by (left) location (healthy InfTurb (3,681 cells, n=3 samples), polyp-bearing patient InfTurb (1,370 cells, n=4 samples), non-polyp EthSin surgical samples (5,928 cells, n=6 samples), and polyp surgical and scraping samples directly from polyp in EthSin (7,346 cells, n=8 samples)) and (right) cell type (3,152 basal, 3,089 differentiating, 8,840 secretory, 1,105 ciliated, and 2,139 glandular cells). f, Heatmap of secretory cells (1,000 cells/location) displaying select genes (AUC>0.65; Supplementary Table 3). g, Violin plots of IFNα, IFNγ, and IL-4/IL-13 gene signatures for secretory cells; healthy InfTurb (3,414 cells), polyp-bearing patient InfTurb (1,239 cells), non-polyp EthSin surgical samples (1,048 cells), polyp surgical and scraping samples directly from polyp in EthSin (3,139 cells); effect size −1.16, −1.05, 1.32, respectively polyp EthSin vs healthy; *Mann-Whitney U-test, p<2.2×10^-16.

Figure 3 |. Reduced epithelial ecological diversity and basal cell hyperplasia in nasal polyps

a, scRNA-seq cell frequency (Fig. 2a; Extended Data Fig. 7a individual points) calculated for each sample; basal cell p=0.00023, glandular p<0.0001, ciliated p=0.0387, non-polyp vs. polyp; and (b) Simpson’s index (Methods, p=0.0021); n=6 non-polyp, n=6 polyp samples; *t-test, two-sided, mean±s.e.m. c, Flow cytometry quantification (Extended Data Fig. 7, full gating), n=6 non-polyp, n=7 polyp samples; *t-test, two-sided, p=0.0005; mean±s.e.m. d, Immunofluorescence and quantification for basal cells normalized to 1,000μm² of epithelium; n=5 non-polyp patients, 10 sections; n=8 polyp patients, 41 sections *Mann-Whitney U-test, p=0.0282, mean±s.d.; scale bar 100μm, and (e) quantification of glandular area, n=6 non-polyp, n=6 polyp patients; *t-test, two-sided, p=0.0022; (Extended Data Fig. 7, isotype and representative). f,g, Bulk-tissue RNA-seq deconvolution by (f) PCA and (g) heatmap over epithelial subset-specific genes (rows) with KNN-clusters (n=4; Methods), from n=10 non-polyp, n=17 polyp samples (columns). h, Violin plot of basal cell gene fraction in scRNA-seq epithelium (Methods, Supplementary Table 4); 5,928 cells, n=6 non-polyp; 4,346 cells, n=6 polyp; effect size 0.457, *Mann-Whitney U-test, p<2.2×10^-16.

Figure 4 |. T2I cytokines and developmental pathways converge at the epigenetic level in basal cells to intrinsically impair differentiation in vivo

a, Heatmap of select genes over basal cell clusters 8 and 12; 860 cells, n=6 non-polyp, 858 cells, n=6 polyp samples; *bimodal test, all displayed genes p<1.97×10⁻³⁹ or less with Bonferroni correction; (Supplementary Table 3). b, Violin plot (Methods; Extended Data Fig. 8; Supplementary Table 3 cell numbers, Supplementary Table 4 gene lists) for IL-4/IL-13 commonly-induced gene signature; *Mann-Whitney U-test, p<1.76×10⁻¹⁵, relative to mean score, with Bonferroni correction. c, Violins of shared IL-4/IL-13 signature (*Mann-Whitney U-test, p<2.2×10⁻¹⁶, effect size 1.305) and Wnt:Notch target gene proportion (*t-test, two-sided, p<2.2×10⁻¹⁶, effect size 0.334, NB: axis truncated, zero indicates equal scores) over 5,928 cells, n=6 non-polyp vs. 4,346 cells, n=6 polyp. d,e, Diffusion pseudotime (Methods) over epithelial cells with unified gene list (Supplementary Table 3 cell numbers, gene list) and (e) violin plot of pseudotime component with green (n=6 non-polyp) and purple (n=6 polyp) underlying distribution. f, Omni-ATAC-seq and HOMER motif enrichment over background peaks (Methods); all q-value<0.0002 Benjamini corrected, and (g) transcription factors from low-input RNA-seq on sorted basal cell populations; *t-test, two-sided, p<0.05 or less Holm-Sidak correction; mean±s.e.m; n=3 non-polyp (4 RNA-seq), n=7 polyp.

Figure 5 |. Transcriptional memory of IL-4/IL-13 exposure retained by basal cells ex vivo and in vivo IL-4Rα blockade partially resets state in polyps

a, tSNE plot of cell types from ALI-cultures (Extended Data Fig. 9g) over 16,173 single cells (8,483 non-polyp; 7,690 polyp) and (b) violin plots for secretory genes (from Fig. 2f) on ALI secretory cells (3,277 non-polyp; 3,143 polyp); not significant (n.s.) except MUC4, PSCA and SCGB1A1 greater in polyp secretory cells; bimodal p<4.04×10⁻¹⁶, Bonferroni correction. c, Basal cell inflammatory memory to IL-4/IL-13 performed (Methods) and displayed as PCA over variable genes and Venn diagram overlaps of differential expression; *t-test, two-tailed, Bonferroni corrected p<0.05, Supplementary Table 3 full lists. d, CTNNB1 expression and Wnt pathway z-score (Methods) in basal cells from (c); n=4 samples/dose; 2-way ANOVA; p<0.0001 for CTNNB1; p<0.0282 for Wnt pathway, (a-d) n=2 basal cell donors each non-polyp and polyp. e-g, scRNA-seq on anti-IL-4Rα-treated individual and (e) tSNE plot of 8,764 single cells from nasal polyps colored by pre- (5,731 cells) and post- (3,033 cells) treatment samples (Extended Data Fig. 10 for more), (f) heatmap of select genes (AUC>0.68 core, p<2.46×10⁻⁵ or less with Bonferroni correction) over basal cells (200 pre, 151 post), and (g) violin plots for select genes; n.s. save *bimodal p<0.00087 with Bonferroni correction Supplementary Table 3 for full list; Wnt score Pre vs. Post: *t-test, two-tailed, p<2.2×10⁻¹⁶, effect size 0.942.