IDO1⁺ fibroblasts and LY6D⁺ club cells contribute to neutrophil recruitment in CRSwNP

CRSwNP patients exhibit both eosinophilic and neutrophilic inflammation¹⁹. Increased neutrophilia and eosinophilia were detected in the mucosa of NP tissues from CRSwNP patients (Fig. 2a). We also used the xCell algorithm to quantify neutrophil infiltration in a bulk RNA-seq dataset of CRSwNP (GSE179265), and again detected the significantly elevated neutrophilia in CRSwNP samples as compared to healthy tissue samples (Fig. 2b). Seeking to identify epithelial and stromal cell subsets contributing to neutrophil infiltration, we generated an integrated dataset built from the data from this study and the data of neutrophils from another CRSwNP scRNA-seq dataset (HRA000772)^20–22(Supplementary Fig. ^6a), and subsequently used an algorithm combining Networkx, Community, and Pygraphviz to plot chemokine-chemokine receptor interaction networks and infer the strongly interacting cell subset pairs. Notable signals from the network included a superlatively strong interaction between IDO1⁺ fibroblasts and neutrophils (Fig. 2c), with MMP7⁺ GCs and LY6D⁺ club cells also interacting strongly with neutrophils (Fig. 2c).

Open in a separate window

Fig. 2

Identification of cell subsets involved in neutrophil recruitment in CRSwNP.

a Representative image and quantification of immunofluorescence staining for neutrophils (upper panel: MPO-red, the nuclear marker DAPI-blue) and eosinophils (lower panel: ECP-red, the nuclear marker DAPI-blue) in IT tissue from healthy individuals (HC-IT) and NP tissue from CRSwNP patients (CRS-NP) (n=10 per group). Neutrophils are indicated by white arrows. Eosinophils are indicated by orange arrows. Scale bar, 50μm (upper), 20μm (below). The data are presented as the means±SEM. Each point represents one sample. HPF, high-power field. The P-value is calculated and reported using the two-tailed Student’s t-test. b Neutrophil scores obtained using the xCell algorithm for CRSwNP (n=17) patients and healthy individuals (n=7). The data are presented as the means±95% CI. The P-value is calculated and reported using a two-tailed Student’s t-test. c Cell-cell interaction and network representation analysis based on chemokine-chemokine receptor interactions. The nodes with a degree of zero and a connection strength less than the average of all the edges were eliminated. The sizes of the nodes were defined as the log2 (counts+1) of the cell subsets. The thickness of the link reflects the degree of the interaction. Cell subsets that strongly interact with neutrophils are indicated in red boxes. d Bubble heatmap for chemokine-chemokine receptor interactions between immune cells and epithelial/stromal cells. Previously validated interactions are indicated by colored straight lines. Chemokines predominantly expressed in LY6D⁺ club cells and IDO1⁺ fibroblasts are indicated in red boxes. e Bar plot depicting the cell composition of epithelial cell subsets in the indicated anatomical regions of human nasal mucosa from CRSwNP patients and healthy individuals. f Immunofluorescence staining for LY6D (red), KRT13 (white), MPO (green), and the nuclear marker DAPI (blue) in CRS-NP tissue and HC-IT tissue. Staining is repeated on 6 participants per group. The epithelium is indicated with orange arrows. Neutrophils are indicated with white arrows. Scale bar, 20μm. g Immunofluorescence staining for COL1A2 (red), IDO1 (white), MPO (green), and the nuclear marker DAPI (blue) in HC-IT tissue and CRS-NP tissue. Staining is repeated on 6 participants per group. The epithelium is indicated with orange arrows. Neutrophils are indicated with white arrows. Scale bar, 20μm.

In particular, the chemokine receptors enriched in neutrophils (CCR1 and CXCR1/2/4) matched extensively with chemokines highly expressed in IDO1⁺ fibroblasts (such as CXCL1/2/3/5/6/8 and CCL5/7/8/11)^23,24(Fig. 2d). The interaction between neutrophils and MMP7⁺ GCs was characterized by high CXCR2 expression in neutrophils and strong CXCL2/3 expression in MMP7⁺ GCs. LY6D⁺ club cells interacted with neutrophils by expressing high levels of S100A8/A9, and their receptor TLR4 was expressed mainly on neutrophils (Fig. 2d). However, MMP7⁺ GCs did not show much difference in the proportion of total epithelial cells in different anatomical regions (Fig. 2e). They were probably a subset of cells with an intermediate state based on their low pseudotime ct values calculated by RNA velocity (Supplementary Fig. ^4a–d). The gene expression pattern of the MMP7⁺ GCs in each tissue was exhibited in four groups of tissues (Supplementary Fig. ^4e). Genes widely expressed in glands, such as PIP and STATH²⁵, were observed in MMP7⁺ GCs from healthy inferior turbinates. Differentially expressed genes in CRS-IT showed abnormalities in cell proliferation and differentiation (TUBA1C, LMNA)^26,27, ECM remodeling (ANXA2, TM4SF1)^28,29 and tissue repair (GADD45A, S100A10)^30,31. Characteristic genes in MMP7⁺ GCs in CRS-MT were mainly involved in immune defense (SOD2, LYZ, LTF)^32,33. MMP7⁺ GCs in CRS-NP expressed factors involved in cell proliferation (LY6E, FOS)³⁴ and antiviral activity (IFITM2 and IFITM3)³⁵. Therefore, MMP7⁺ GCs were not considered to be associated with neutrophil infiltration in CRSwNP.

In contrast to that of MMP7⁺ GCs, the proportion of LY6D⁺ club cells was greater in the IT and NP tissues of CRSwNP patients than in the IT tissue of healthy individuals (HC-IT), suggesting their potential role in CRSwNP development (Fig. 2e). Consistent with these findings, using the HRA000772 dataset, we also detected higher proportions of LY6D⁺ club cells and IDO1⁺ cells in the CRSwNP with higher neutrophil numbers (neCRSwNP) than those with lower neutrophil numbers (eCRSwNP) (Supplementary Fig. ^6b–d). We then conducted immunofluorescence analyses on NP tissues from 6 CRSwNP patients (CRS-NP) and IT tissues from 6 healthy controls (HC-IT), and the results revealed a preferential distribution of neutrophils (MPO⁺ cells) in the LY6D⁺ club cell-rich and IDO1⁺ fibroblast-rich regions (Fig. 2f, g), supporting their ability to recruit neutrophils in an inflammatory milieu in CRSwNP. These results collectively support that LY6D⁺ club cells and IDO1⁺ fibroblasts facilitate neutrophil recruitment in CRSwNP.

LY6D⁺ club cells drive IL-1 signaling-mediated neutrophilic inflammation in CRSwNP

We further compared the proportions of LY6D⁺ club cells across anatomical regions. An elevated proportion of LY6D⁺ club cells within the total epithelial cell population was noted in CRS-ITs compared to HC-ITs, and in NP tissues compared to adjacent MT tissues (Fig. 3a). We then performed immunofluorescence analyses to evaluate the distribution of LY6D⁺ club cells in different tissues. Immunostaining detected only a few LY6D⁺ cells in normal IT tissues, but more LY6D⁺ cells in NP tissues, reflecting the preferential induction of LY6D⁺ club cells in an inflammatory milieu (Fig. 3b). A higher expression of S100A8 was also detected in the epithelium of NP tissues, further suggesting the neutrophilic inflammation in CRSwNP³⁶ (Fig. 3c and Supplementary Fig. ^6e). LY6D⁺ club cells were highly conserved across three anatomical regions as indicated in the fate decision tree, and PITX1 ranked as the top differentially expressed transcription factor determining LY6D⁺ club cell differentiation (Fig. 3d, e).

Open in a separate window

Fig. 3

LY6D⁺ club cells regulate neutrophilic inflammation in CRSwNP.

a UMAP displaying the distribution of LY6D⁺ club cells in the indicated anatomical regions of human nasal mucosa from CRSwNP patients and healthy individuals. Enriched LY6D⁺ club cells are indicated with blue circles (left panel). The proportions of LY6D⁺ club cells in the indicated anatomical regions are shown on the right. The colors are consistent in the two panels. b Representative immunofluorescence staining for LY6D (green), S100A8 (red), KRT13 (white), and the nuclear marker DAPI (blue) in HC-IT tissue and CRS-NP tissue. White arrows indicate colocalization of LY6D, S100A8, and KRT13. Scale bar, 20μm. c Quantification of the data in (b) (n=5 per group) and Supplementary Fig. ⁶e (n=5 per group). Each point represents one sample. The data are presented as the means±SEM. HPF, high-power field. The P-value is calculated and reported using a two-tailed Student’s t-test. d A dendrogram of regulons for epithelial cell subsets in the indicated anatomical regions constructed from the fate decision tree analysis. The transcription factors at each branching point are representative regulons of subjacent groups. The brown triangles show the proximity between LY6D⁺ club cells (E11) and club cells (E03) during differentiation. The numbering of cell subsets is consistent with that in Fig. 1k. e TFs enriched in LY6D⁺ club cells aligned by TF specificity score. PITX1 (red) is the top transcription factor responsible for LY6D⁺ club cell differentiation. f Heatmap of gene expression analyzed by scRNA-seq displaying representative genes for 12 epithelial cell subsets. LY6D⁺ club cells are highlighted in cyan letters. g Volcano plot displaying the differentially expressed genes (DEGs) between LY6D⁺ club cells and other epithelial cell subsets. The P-value is calculated and reported using a two-tailed Student’s t-test. h UMAP with the epithelial cell compartment displaying the expression of four upregulated genes (IL1RN, S100A8, AREG, and EREG) in LY6D⁺ club cells. LY6D⁺ club cells are indicated with red circles. i Pathway enrichment analysis revealing the enriched signaling pathways in LY6D⁺ club cells when compared with those in other epithelial cells. j RNA velocity analysis based on RNA splicing information indicating that LY6D⁺ club cells are maturely differentiated club cells. k Heatmap displaying dynamic changes in the expression of functional genes and TFs during the maturation process of LY6D⁺ club cells.

Previous studies have shown that the expression of S100A8 and S100A9 is elevated in nasal polyps as compared to control tissues^37,38, and is associated with neutrophilic inflammation and CRS severity³⁹. By exhibiting epithelial cell subset marker gene expression via heatmap, we observed the upregulation of S100A8 in LY6D⁺ club cells (Fig. 3f). We next explored the differentially expressed genes (DEGs) in LY6D⁺ club cells as compared to other epithelial cells (Fig. 3g). In addition to LY6D and S100A8, S100A9 was also significantly upregulated in LY6D⁺ club cells (Fig. 3g). UMAP showed that LY6D⁺ club cells were the main cell source of S100A8 and S100A9 in the epithelium that may promote neutrophil chemotaxis in CRSwNP³⁶ (Fig. 3g, h). The high expression of EREG and AREG in LY6D⁺ club cells indicated their involvement in eosinophil reprogramming and goblet metaplasia in response to inflammation^40,41 (Fig. 3h).

Pathway enrichment analysis revealed that the transcriptome of LY6D⁺ club cells was enriched in genes induced by IL-1 signaling (Fig. 3i). The RNA velocity profile of total club cells indicated that LY6D⁺ club cells originated from resident club cells, suggesting that some club cells in the face of upregulated IL-1 signaling progressively acquired LY6D⁺ club cell identity in the mucosal epithelium in CRSwNP patients (Fig. 3j). The expression of several key functional genes and transcription factors upregulated during the maturation process of LY6D⁺ club cells was presented in the heatmap (Fig. 3k). IL1RN was inferred by RNA velocity, iteratively indicating that IL-1 signaling participates in the transdifferentiation of LY6D⁺ club cells (Fig. 3k). Pathway enrichment analysis revealed that genes involved in neutrophil degranulation were also enriched in LY6D⁺ club cells, reflecting the regulation of neutrophil recruitment by LY6D⁺ club cells (Fig. 3i). Taken together, these findings underscore the role of LY6D⁺ club cells in IL-1 signaling-mediated neutrophilic inflammation in CRSwNP.

IDO1⁺ fibroblasts secrete chemokines that facilitate neutrophil recruitment in CRSwNP

To identify the stromal cell subsets responsible for inflammation in CRSwNP, we subclustered the stromal cell compartment into 27 clusters and annotated them into 15 distinct cell subpopulations (Supplementary Fig. ^{7a, b} and Fig. 1g). Endothelial cells, pericytes and smooth muscle cells did not show much variation in the proportions of cell subsets across different anatomical regions except for an increased proportion of arterial endothelial cells and decreased proportion of lymphatic endothelial cells in CRS-related tissues as compared to those in HC-ITs, suggesting weakened lymphatic infiltration but enhanced angiogenesis in inflammatory tissues (Fig. 4a and Supplementary Fig. ^7c–f). Fibroblast subsets exhibited substantial disparities in cellular proportions within the stromal cell compartment. The proportions of IDO1⁺ and OXTR⁺ fibroblasts were markedly higher in NPs than in other tissues (Fig. 4a).

Open in a separate window

Fig. 4

IDO1⁺ fibroblasts contribute to IL-1 signaling-mediated neutrophilic inflammation in CRSwNP.

a Bar plot depicting the cell composition of stromal cell subsets in the indicated anatomical regions of human nasal mucosa from CRSwNP patients and healthy individuals. b Transcription factor fate decision tree analysis of stromal cells displaying two distinguishable modules consisting of six fibroblast subsets with remarkable differences in TF patterns. Module 1 and module 2 are indicated by blue and green boxes, respectively. c UMAP displaying the expression of marker genes of the two modules in the stromal cell compartment. d UMAP of stromal cell compartment displaying fibroblasts in the indicated anatomical regions. Bar plot displaying differences in the proportions of seven fibroblast subsets in the indicated anatomical regions. e UMAP displaying the predominant expression of CXCL1 and CXCL8 in IDO1⁺ fibroblasts. f Bubble heatmap depicting the expression of representative genes of IDO1⁺ fibroblasts across different stromal cell subsets. g Radar plot displaying the pathway enrichment analysis results for the four fibroblast subsets in module 1. The colors in the circles reflect the P-values. The P-values are calculated using the Wilcoxon rank-sum test, applying a two-sided hypothesis, and are adjusted for multiple testing using the Benjamini–Hochberg correction. h Representative immunofluorescence staining for CXCL8 (white), IDO1 (green), COL1A2 (red, a marker gene of fibroblasts), and the nuclear marker DAPI (blue) in HC-IT tissue and CRS-NP tissue. Staining is repeated on 6 samples per group. The white arrows indicate the colocalization of IDO1, CXCL8, and COL1A2 in the NPs. Scale bar, 20μm.

To better characterize the functionality of the fibroblast subsets, we proceeded to construct a transcription factor fate decision tree for different cell subsets within the stromal cell compartment. Six out of the seven distinct fibroblast subsets were categorized into two main modules (Fig. 4b, c). We noticed that module 1 comprised fibroblasts that exerted pro-inflammatory effects by enhanced production of certain chemokines, such as CXCL1 and CXCL8, while module 2 encompassed fibroblasts that mainly reside in the adventitia and are essential for regulating the integrity and function of the vessel structure^42–44. The fibroblast clusters were displayed according to different anatomical regions (Fig. 4d). Within the cell clusters in module 1, IDO1⁺ and OXTR⁺ fibroblasts were enriched in inflammatory tissues, mostly in NPs, and they were barely detected in healthy tissues (Fig. 4d).

Considering the potent interaction detected between IDO1⁺ fibroblasts and neutrophils, we explored the gene expression patterns of different stromal cell subsets. The transcriptome of IDO1⁺ fibroblasts was enriched in chemokines (CXCL1/2/3/8) that are relevant to neutrophilic inflammation (Fig. 4e, f). We subsequently conducted pathway enrichment analysis on the fibroblast subsets within module 1, whose gene expression pattern was associated with inflammatory responses, to scrutinize the regulatory pathways in which IDO1⁺ fibroblasts are implicated (Fig. 4g). Interestingly, both the IL-1 signaling pathway and the NF-κB signaling were enriched in IDO1⁺ fibroblasts, indicating the upstream regulation of IDO1⁺ fibroblasts by IL-1 signaling in inflammation development. The high expression of MMP3 and LIF in IDO1⁺ fibroblasts also indicated the regulation of IDO1⁺ fibroblasts by IL-1 signaling^45,46 (Fig. 4f). Immunofluorescence staining revealed increased CXCL8 protein level and a greater number of IDO1⁺ fibroblasts (IDO1⁺ COL1A2⁺ cells) in CRS-NP samples as compared to HC-IT samples (Fig. 4h). Together, these findings elucidate IL-1 signaling as a common pathway inducing transdifferentiation of LY6D⁺ club cells and IDO1⁺ fibroblasts to facilitate neutrophil recruitment in CRSwNP.

IL-1β−induced transdifferentiation of LY6D⁺ club cells and IDO1⁺ fibroblasts promotes neutrophil recruitment

IL-1 signaling can be activated by the interaction between IL-1β and IL-1 receptor (IL-1R), leading to various immune responses including neutrophilic inflammation⁴⁷. Here we deployed recombinant IL-1β on air-liquid interface (ALI) cultures generated from primary human nasal epithelial cells (HNEs) (Supplementary Fig. ^8a). Bulk RNA-seq data revealed that the addition of IL-1β elicited the LY6D⁺ club cell state of ALI-cultured HNEs, as IL-1β-stimulated ALI-cultured HNEs highly expressed genes that were also upregulated in LY6D⁺ club cells detected by scRNA-seq, such as LY6D, SPRR2F, S100A9, and LYPD3 (Fig. 5a–c). Immunofluorescence staining of ALI-cultured HNEs revealed the colocalization and elevation of LY6D and S100A8 upon IL-1β stimulation (Fig. 5d). These results suggested that IL-1β induced transdifferentiation of LY6D⁺ cells in vitro. ELISA detected the increased secretion of S100A8/A9 protein from ALI-cultured HNEs upon IL-1β stimulation (Fig. 5e). Considering the ability of S100A8, S100A9, and S100A8/A9 to promote neutrophil activation, chemotaxis and adhesion, we performed a chemotaxis assay to determine whether the secretion of S100A8/A9 contributes to IL-1β-mediated neutrophil recruitment³⁶. The results showed that the media from ALI-cultured HNEs stimulated with IL-1β exhibited a stronger neutrophil chemotactic capacity compared to the control media. However, the application of a neutralizing antibody (anti-S100A8/A9 antibody) significantly reduced the neutrophil chemotactic capacity of the conditioned media. (Fig. 5f and Supplementary Fig. ^8b).

Open in a separate window

Fig. 5

IL-1β induces transdifferentiation of LY6D⁺ club cells and IDO1⁺ fibroblasts to promote neutrophil recruitment.

a Volcano plot of DEGs between IL-1β -stimulated and PBS-treated air-liquid interface (ALI) -cultured primary human nasal epithelial cells (HNEs) identified with the cut-off criterion P<0.05 and |log2FC|≥ 1. The P-values are calculated and reported using two-tailed Student’s t-tests. Blue dots: significantly downregulated genes; red dots: significantly upregulated genes. b UMAP with epithelial cell compartment displaying the expression of four genes that are upregulated in LY6D⁺ club cells identified by scRNA-seq. The red circles indicate LY6D⁺ club cells in the epithelial cell compartment according to the scRNA-seq data. c Expression of four genes in panel b of ALI-cultured primary HNEs treated with the indicated conditions (n=5 per group, as determined by bulk RNA-seq). The length of the error bars is a 95% confidence interval for the mean in Fig. 5c. The P-values are calculated and reported using a two-tailed Student’s t-test. d Representative immunofluorescence staining for LY6D (green), S100A8 (red), and the nuclear marker DAPI (blue) in ALI-cultured primary HNEs upon the indicated stimulations. Scale bar, 50μm. e S100A8/A9 protein levels in the media (measured by ELISA) upon the indicated stimulations of ALI-cultured HNEs (n=6 per group). Each point represents one sample. The data are presented as the means±SEM. The P-value is calculated and reported using a two-tailed Student’s t-test. f Number of neutrophils passing through the membrane of a transwell insert by ALI-cultured HNEs-exposed media in the presence or absence of IL-1β stimulation and anti-S100A8/A9 antibody (n=8 per group). ALI-cultured HNEs-exposed media (control media) and fresh media containing IL-1β (control media+IL-1β) are used as a negative control. Data are presented as means±SEM. The P-values are calculated and reported using one-way ANOVA. g Heatmap displaying the expression of the indicated chemokines in primary fibroblasts derived from HC-IT tissue and CRS-NP tissue upon the indicated stimulations. h Representative immunofluorescence staining for IDO1 (green), CXCL8 (red), and the nuclear marker DAPI (blue) displaying the protein levels of IDO1 and CXCL8 upon the indicated stimulations. Scale bar, 100μm. i Quantification of the data in (h) (n=3 per group). Each point represents one sample. The data are presented as the means±SEM. The P-values are calculated and reported using one-way ANOVA. j CXCL8 protein levels measured by ELISA after indicated stimulations in primary fibroblasts derived from HC-IT tissues and CRS-NP tissues (n=4 per group). Each point represents one sample. The data are presented as the means±SEM. The P-values are calculated and reported using two-way ANOVA. k Number of neutrophils passing through the membrane of a transwell insert by fibroblasts-exposed media in the presence or absence of IL-1β stimulation and anti-CXCL8 antibody (n=4 per group). Each point represents one sample. Fibroblasts-exposed media (control media) and fresh media containing IL-1β (control media+IL-1β) are used as negative controls. Data are presented as means±SEM. The P-values are calculated and reported using two-way ANOVA.

To explore the effect of IL-1β on the induction of fibroblasts, we treated cultured primary fibroblasts isolated from IT tissues and NPs with IL-1β and performed bulk RNA sequencing. Bulk RNA sequencing results revealed high expression of genes encoding neutrophil chemoattractants (CXCL1, CXCL2, CXCL3, CXCL5, CXCL6, and CXCL8) in fibroblasts upon IL-1β stimulation (Fig. 5g). These genes were also upregulated in IDO1⁺ fibroblasts according to scRNA-seq analysis, suggesting that the primary fibroblasts acquire the identity of IDO1⁺ fibroblasts upon IL-1β stimulation. Immunofluorescence staining demonstrated that IL-1β was capable of activating fibroblasts and inducing the expression of CXCL8 and IDO1 in both IT-derived and NP-derived fibroblasts (Fig. 5h, i). In line with the results of bulk RNA sequencing, ELISA showed an increase in CXCL8 secretion in fibroblasts treated with IL-1β compared to those treated with PBS (Fig. 5j). As expected, culture media from IL-1β-exposed human nasal primary fibroblasts resulted in an increase in the transmigration of purified blood neutrophils compared with media of normal fibroblasts or fresh media mixed with IL-1β, whereas the application of a neutralizing antibody (anti-CXCL8 antibody) significantly reduced the neutrophil chemotactic capacity of the conditioned media (Fig. 5k). Therefore, we reason that IL-1β induces both epithelial cells and fibroblasts to promote the recruitment of neutrophils to the sites of inflammation in CRSwNP.

Preparation of single-cell suspensions

Nasal mucosa was freshly sampled from the middle turbinates (n=7), inferior turbinates (n=9), nasal polyps (n=15) of CRSwNP patients, and inferior turbinates (n=2) of patients with cerebral spinal fluid leak. The nasal biopsies were washed in phosphate-buffered saline (PBS, 10010023, Thermo Fisher) to remove mucus and blood cells. Then, the nasal tissues were cut into approximately 0.5-mm³ pieces in RPMI-1640 medium supplemented with 1% penicillin/streptomycin, and then enzymatically digested with the Multi Tissue Dissociation Kit 2 (MACS# 130-110-203) at 37°C for 30min with agitation, according to the manufacturer’s instructions. Following cell dissociation, the resultant cell suspension was sequentially filtered through cell strainers with pore sizes of 70μm and 40µm (BD). Subsequently, the samples were centrifuged at 300×g for 10min. Subsequent to the removal of the supernatant, the cells forming the pellet were reconstituted in red blood cell lysis buffer (Thermo Fisher) and subjected to a 2-min incubation on ice to lyse the red blood cells. Following dual washes with PBS, the cellular pellets were resuspended in PBS supplemented with 0.04% bovine serum albumin (A7906, Sigma–Aldrich).

Single-cell RNA library construction and sequencing

DNBelab C Series High-throughput Single-cell System (BGI-research) was utilized for scRNA-seq library preparation. Briefly, the single-cell suspensions underwent a series of processes to generate barcoded scRNA-seq libraries. These steps encompassed droplet encapsulation, emulsion breakage, collection of beads containing the captured mRNA, reverse transcription cDNA amplification, and subsequent purification. The cDNA was subjected to fragmentation into shorter segments spanning 250–400 base pairs. Following this, the construction of indexed sequencing libraries was achieved in accordance with the protocol provided by the manufacturer. Qualification was performed using the Qubit ssDNA Assay Kit (Thermo Fisher Scientific) and the Agilent Bioanalyzer 2100. Subsequent to library preparation, all the constructs underwent sequencing using the DIPSEQ T1 sequencing platform in the China National GeneBank via pair-end sequencing methodology. The sequencing reads contained 30-bp read 1 (including the 10-bp cell barcode 1, 10-bp cell barcode 2, and 10-bp unique molecular identifiers [UMI]), 100-bp read 2 for gene sequences, and 10-bp barcodes read for sample index. Next, processed reads were aligned to the GRCh38 reference genome using STAR (v2.5.3). The identification of valid cells was achieved through an automated process utilizing the “barcodeRanks” function from the DropletUtils tool. This function was employed to eliminate background beads and those with UMI counts falling below a predetermined threshold, using the UMI number distribution characteristic of each cell. Finally, we computed the gene expression profiles of individual cells and subsequently generated a matrix of genes by cells for each library by means of PISA.

Alignment, quantification, and quality control of single-cell RNA sequencing data

The droplet-based sequencing data were subjected to alignment and quantification through the utilization of CellRanger software (version 3.0.2, designed for 3′ chemistry), employing the GRCh38.p13 human reference genome. The Python package Scanpy (version 1.7.1)⁶⁸ was employed to load the matrix containing cell-gene counts and to execute quality control procedures for both the newly generated dataset and the acquired datasets. For each sample, genes associated with mitochondria (indicated by gene symbols commencing with “MT-“) and ribosomal proteins (initiated by gene symbols commencing with “RP”) were eliminated from consideration. After that, cells possessing less than 2000 UMI counts and 250 detected genes were identified as empty droplets and subsequently excluded from the datasets. Finally, genes demonstrating expression in fewer than three cells were excluded from further analysis.

Doublet detection

In order to rule out doublets, we implemented the Scrublet software (version 0.2.3)⁶⁹, which facilitated the identification of artifactual libraries originating from two or more cells within each scRNA-seq sample, comprising both the newly generated dataset and the compiled datasets. The doublet score for each individual single cell, along with the threshold determined from the bimodal distribution, was computed using the default parameters (sim_doublet_ratio=2.0; n_neighbors=None; expected_doublet_rate=0.1, stdev_doublet_rate=0.02). After that, a comprehensive assessment was conducted on the remaining cells and cell subsets to identify potential false negatives from the scrublet analysis. This evaluation was guided by the following sets of criteria: (1) cells with more than 8000 detected genes, (2) subsets that expressed marker genes from two distinct cell types, which are unlikely according to prior knowledge (i.e., CD3D for T cells and EPCAM for epithelial cells). Any cells or subsets identified as doublets were excluded from subsequent downstream analyses.

Graph subsetting and partitioning cells into distinct compartments

The downstream analysis included normalization (scanpy.pp.normalize_total method, target_sum=1e4), log-transformation (scanpy.pp.log1p method, default parameters), cell cycle score (scanpy.tl.score_genes_cell_cycle method), cell cycle genes defined in Tirosh et al.⁷⁰, feature regress out (scanpy.pp.regress_out method, UMI counts, percentage of mitochondrial genes and cell cycle score were considered to be the source of unwanted variability and were regressed), feature scaling (scanpy.pp.scale method, max_value=10, zero_center=False), PCA (scanpy.tl.pca method, svd_solver=’arpack’), batch-balanced neighborhood graph building (scanpy.external.pp.bbknn method, n_pcs=20)⁷¹, leiden graph-based subseting (scanpy.tl.leiden method, Resolution=1.0)⁷², and UMAP visualization⁷³ (scanpy.tl.umap method) performed using scanpy. The initial categorization of the subsets encompassed a division into six distinct compartments, achieved through the utilization of marker genes established in the existing literature in conjunction with genes exhibiting differential expression. (scanpy.tl.rank_gene_groups method, method=’Wilcoxon test’). Specifically, the epithelial compartment was annotated using a gene list (EPCAM, KRT8, KRT18, KRT19, PIGR), T and ILCs compartment (CD2, CD3D, CD3E, CD3G, TRAC, IL7R), B cell compartment (JCHAIN, CD79A, IGHA1, IGHA2, MZB1, SSR4), MNPs compartment (HLA-DRA, CST3, HLA-DPB1, CD74, HLA-DPA1, AIF1), Mast cell compartment (TPSAB1, CPA3, TPSB2, CD9, HPGDS, KIT), and Stromal cell compartment (IGFBP7, IFITM3, TCF7L1, COL1A2, COL3A1, GSN). Subsequently, the epithelial compartment was subjected to sorting for subsequent downstream analysis. Detailed methods and marker genes were included in Supplementary Data ¹.

Transcription factor module analysis

The Python package pySCENIC workflow (version 0.11.0) with default settings was used to infer active TFs and their target genes in all human cells^74,75. Specifically, the pipeline was executed in three steps. Initially, the single-cell gene expression matrix was filtered to eliminate genes whose expression was detected in fewer than ten total cells. The retained genes were subsequently employed to construct a gene-gene correlation matrix, which facilitated the identification of co-expression modules through the application of a regression per-target approach utilizing the GRNBoost2 algorithm. Subsequent to the initial step, each identified module was systematically refined based on a regulatory motif in close proximity to a transcription start site (TSS). The acquisition of cis-regulatory footprints was facilitated through the utilization of positional sequencing methodologies. The binding motifs of the TFs were then used to build an RCisTarget database. Modules were retained based on the enrichment of transcription factor (TF)-binding motifs among their respective target genes. In cases where target genes lacked direct TF-binding motifs, they were excluded from consideration. In the third phase, we assessed the influence of each regulon on individual single-cell transcriptomes through the utilization of the area under the curve (AUC) score, employing the AUCell algorithm as the evaluative metric. The scores pertaining to transcription factor motifs within gene promoters and regions surrounding the transcription start sites, specific to the hg38 human reference genome, were acquired from the RcisTarget database. Concurrently, the list of transcription factor-associated genes was obtained from the Humantfs database⁷⁶.

Fate decision tree construction (regulon-based)

Dendrogram plots were generated for epithelial cells using the sc.pl.dendrogram method from the Scanpy package. These plots were generated based on the AUCell matrix comprising 608 regulons, aiming to visualize more nuanced alterations. We deciphered the diverging composite rules of a regulon-based dendrogram by testing each branching node for differential regulon importance. Thereafter, differential analysis of regulon expression was conducted for each node using the Wilcoxon test (implemented through the sc.tl.rank_gene_groups method with method=‘Wilcoxon test’), with the aim of deducing the sequence of regulon-driven propagation events. The hierarchical clustering (dendrogram tree) was generated at a single-cell level⁷⁷.

Datasets integration

In this study, we utilized a previously published scRNA-seq dataset of CRSwNP²⁰ (GSA: HRA000772), and integrated the neutrophils with the data of this study to investigate the expression of inflammatory factors in neutrophils in human nasal mucosal tissues. Specifically, we compared the downloaded fastq files with the barcodes-genes matrix utilizing Alevin-fry⁷⁸. The matrix underwent initial quality control, doublet removal, and normalization, applied in accordance with the dataset from the previous section. The gene expression and cell annotation of the dataset were modeled using CellTypist⁷⁹. Subsequently, the trained model was used to perform Label Transfer on the HRA000772 dataset. In particular, myeloid cells annotated by Label Transfer were manually reannotated based on marker genes, thereby identifying the neutrophils subset (FCGR3B⁺CXCR1⁺CXCR2⁺). The study also integrated the scRNA-seq data with the data of normal ethmoid and sphenoid sinuses (n=5) and nasal polyps (n=11) from the HRA000772 dataset for statistical analyses.

RNA velocity

Cells that met the quality control criteria were used to filter the loom file generated by the Velocyto python package based on the cell barcodes⁸⁰. This package was used to conduct splicing analysis on the bam file in preparation for subsequent RNA velocity analysis. The filtered loom file served as an input within the Scanpy pipeline, implemented as part of the CellRank pipeline⁸¹. The loom file derived from Velocyto was harnessed to compute RNA velocities for each cell according to standard parameters for the software. CellRank generates both stochastic and dynamic models of RNA velocity, which were compared via the computation of a consistency score for each cell, employing each modeling approach, in accordance with the guidance provided by the authors. Pseudotime was subsequently calculated based on the outcomes of RNA velocity analysis, while latent time was deduced from the dynamic velocity results.

Gene set scoring and identification of significant changes

We scored the gene sets of all cells and subsets using the Scanpy python package (sc.tl.score_genes method, ctrl_size=len(genesets), gene_pool=None, n_bins=25, use_raw=None). The score was the average expression of a set of genes subtracted from the average expression of a reference set of genes. The reference set was randomly sampled from the gene_pool for each binned expression value. To prevent highly expressed genes from dominating a gene set score, we scaled each gene of the log2 (TP10K+1) expression matrix by its root mean squared expression across all cells. After obtaining the score-cell matrix of the signatures, differential signature analysis (sc.tl.rank_gene_groups method, method=’Wilcoxon test’) was implemented to identify significant changes among different nasal anatomical regions. All pathways included in gene set enrichment analysis (Fig. 3i, Fig. 4g, and Supplementary Fig. ^{7e, f}) were obtained from Reactome⁸².

Cell-cell interaction and network representation analysis

To plot chemokine-chemokine receptor interaction networks, we employed the Networkx (version 2.5) (https://github.com/networkx/networkx), Community (version 1.0.0b1) and Pygraphviz (version 1.6) (https://github.com/pygraphviz/pygraphviz) python packages to construct a network defined using the count of interactions between cell subsets. The pipeline was implemented in three steps. First, the nodes with a degree of zero were eliminated. Second, any edges with a connection strength less than the average of all the edges were removed. Third, the sizes of the nodes were defined as the log2 (counts+1) of the cell subsets, and the network with the Kamada Kawai layout algorithm (networkx.kamada_kawai_layout method) was utilized to visualize the network. The thickness of the line connecting the two-cell subsets was directly proportional to the degree of interaction strength between them. The chemokines-chemokines receptor interaction data were obtained from IMEx Consortium⁸³, IntAct⁸⁴, InnateDB-All⁸⁵, MINT⁸⁶, and I2D⁸⁷ database.

Estimation of neutrophil infiltration in CRSwNP

In this study, we applied the xCell algorithm to determine the immune cell subsets in the RNA-seq dataset (GSE179265). The xCell algorithm represents a gene signature-based approach derived from learning from numerous pure cell types originating from diverse sources. This method adeptly enables a cell type enumeration analysis using gene expression data, providing a comprehensive assessment of 64 immune and stromal cell types. This attribute endows it with a commendable capability to accurately depict the intricate landscape of cellular heterogeneity within tissue expression profiles⁸⁸.

Animals

6-8 weeks-old female C57BL/6 mice (n=36) purchased from SPF Biotech were used in these experiments. The mice were maintained in individually ventilated cages in a specific pathogen-free facility under 12h light–dark cycles at 22–24°C and 50–60% humidity. The protocol for the animal studies was approved by the Laboratory Animal Ethical and Welfare Committee of Shandong University Cheeloo College of Medicine (23086).

Neutrophilic CRSwNP mouse model and treatment with an IL-1R antagonist (Anakinra)

Mice were randomly divided into three groups consisting of 6 individuals each. The construction of the mouse model of CRSwNP with neutrophilia was carried out following a previously described protocol⁴⁹. For the control group, 20µl of PBS was dropped into the nasal cavities three times a week for 3.5 consecutive months. Mice in the model groups received 10µg of LPS (from Escherichia coli; Sigma–Aldrich, Merck Millipore, Germany) in 20µl of PBS three times a week for 3.5 consecutive months. For the anakinra-treated group, starting on the 77th day, the mice were given 100µg of Anakinra (MedChemExpress, HY-108841, USA) in 20µl of PBS by intranasal instillation and 100µg of Anakinra in 200µl of PBS by intraperitoneal injection 30min after LPS stimulation for 2 weeks. For the following 2 weeks, only 100µg of Anakinra was intranasally administered in 20µl of PBS within 30min each after LPS stimulation. The animals were sacrificed 24h after the last nasal challenge. The graphic protocol is depicted in Supplementary Fig. ^8c. NLF was collected immediately from the sacrificed mice by washing the nasal cavity with 1mL of ice-cold PBS three times.

Immunofluorescence staining

We removed the skin on the heads of the mice and then excised the mandibles. The heads of the mice were fixed in 4% paraformaldehyde at room temperature for at least 24h and decalcified for 7 days. For human nasal tissues, biopsy samples were soaked in 4% paraformaldehyde for 24h. For both the murine and human studies, after dehydration and paraffin embedding, the tissue samples were cut into 4-µm-thick paraffin sections. The slides were incubated at 65°C for 1h, dewaxed, hydrated, and subsequently heated in antigen retrieval liquid for 15min in a microwave oven. After cooling to room temperature, the slides were permeated with PBS containing 1% Triton X-100 for 20min. The slides were washed in PBS 3 times and blocked with 5% bovine serum albumin at room temperature for 1h. After that, the slides were incubated with the primary antibody (see Supplementary Table ² for a complete list and dilutions) overnight at 4°C in a humidified chamber. The slides were gently washed with PBS 3 times, and incubated with a fluorescent secondary antibody at room temperature for 1h. After washing with PBS, the slides were stained with 4′, 6-diamidino-2-phenylindole (DAPI) (Solarbio, C006, China) for 10min. After another washing step with PBS, the slides were cover-slipped with an anti-fade mounting medium (Solarbio, S2100, China)⁸⁹. Image acquisition was performed using a fluorescence microscope (Olympus, VS200, Japan).

Multiplexed immunohistochemistry

Multiplexed immunohistochemistry (mIHC) assay was performed using the Opal 6-Plex Detection Kit (AKOYA #811001, USA) as described previously⁹⁰. Briefly, after dewaxing and hydration, the slides were boiled in AR6 buffer in a microwave oven for 15min. The tissue sections on the slides were incubated with blocking buffer for 30min and then with primary antibody (see Supplementary Table ² for a complete list and dilutions) for 2h at room temperature in a humidified chamber. Then the slides were washed with TBST twice and incubated with Opal polymer anti-rabbit/mouse horseradish peroxidase (HRP) for 10min at room temperature. Then, 100–300µl of Opal Fluorophore working solution was added to each slide. After washing with TBST twice, the slides were incubated at room temperature for 10min. The previous steps were repeated as needed. DAPI working solution was applied on the slides for 10min at room temperature. As a final step, the slides were washed and cover-slipped with an anti-fade mounting medium. Image acquisition was performed using the TissueFAXS imaging system (TissueGnostics, Germany).

Isolation and culture of primary human nasal epithelial cells (HNEs)

Human nasal epithelial cells were scraped from patients’ nasal mucosa during endoscopic sinus surgery. The cells were placed in an Eppendorf tube containing 1ml of bronchial epithelial cell medium (BEpiCM) (ScienCell, 3211, USA) supplemented with 1% penicillin/streptomycin and 1% bronchial epithelial cell growth supplement immediately upon acquisition. Cells were seeded within 6h in six-well plates pre-coated with Collagen Type I (Corning, 354236, USA) and maintained in a humidified incubator at 37°C containing 5% CO₂. The media was changed every two days. When cells reached 90% confluence in the well, they were transferred to the upper chamber of polyester Transwell inserts (0.4µm, 0.33cm², BIOFIL, TCS016012, China) pre-coated with Collagen Type I. After that, 1ml of BEpiCM was added into the lower chamber, and media was replaced every two days. At confluence, the media was replaced with differential media (BEpiCM: DMEM/F12=1:1) in the basal chamber and the apical surface was exposed to provide an air-liquid interface (ALI). Monolayers were grown at the ALI for an additional 21 days to promote differentiation into a nasal epithelium with basal, multiciliated, and secretory cells. On day 22, media containing PBS or recombinant IL-1β (10ng/ml) (Abbkine, PRP100051, USA) was added to the basal chambers for 3 days.

Isolation and culture of primary human nasal fibroblasts (HNFs)

The inferior turbinate or nasal polyp tissues were soaked in penicillin-streptomycin solution (Solarbio, P1400, China) for 3min and cut into small pieces. After digestion in Trypsin-EDTA solution (Macgene, CC017-500) for 10min, the tissues were put into cell culture flasks with DMEM media supplemented with 10% FBS. The cells were cultured in a humidified incubator at 37°C containing 5% CO₂, and the media was replaced every 2 days. The migrated cells were nasal mucosa-derived fibroblasts. When cells reached 90% confluence in the well, PBS or IL-1β (10ng/ml) was added into the wells, and the cells were cultured for 1 day.

Isolation of human peripheral blood neutrophils

Neutrophils were enriched from peripheral blood by means of Polymorphprep (Serumwerk Bernburg AG, 1895) density centrifugation. We carefully layered 5.0ml of anti-coagulated whole blood over 5.0ml of PolymorphPrep in a 15ml tube. The tubes were centrifuged at 500×g for 30min at 20°C. After centrifugation, two bands were visible, and the neutrophils were enriched in the lower band. The cells were aspirated to another clean tube and an equal volume of sterile normal saline solution was added. After incubating at room temperature for 10min, the tubes were put on a centrifuge at 500×g for 30min. The supernatant was discarded, and the cell pellet was resuspended in Roswell Park Memorial Institute (RPMI) 1640 media supplemented with 1% FBS.

Neutrophil chemotaxis assay

For the cell migration assay, after resuspension in RPMI-1640 media supplemented with 1% FBS, the neutrophils were seeded in the upper compartment of 24-transwell plates with 3-μm pores (Costar, 3415). Conditioned media from ALI-cultured HNEs or fibroblasts, in the presence or absence of IL-1β stimulation and neutralizing antibody blocking (see Supplementary Table ² for a full list of antibodies and dilutions), was added into the lower chamber of 24-transwell plates. After 3h of incubation at 37°C in 5% CO₂, the number of the migrated cells in the lower chamber was counted.

Enzyme-linked immunosorbent assay (ELISA)

ELISAs were performed using multiple ELISA kits (4A Biotech, CHE0011, CME0008, CME0004, China) according to the manufacturer’s instructions. In brief, the standards and samples were added to the antibody pre-coated 96-well ELISA plate, which was subsequently incubated at 37C for 2h. The liquid was removed, and the plate was washed 4 times with wash buffer. Then, an enzyme-linked antibody was applied to the plate, which was incubated at 37°C for 60min. After a washing step, the avidin-biotin-peroxidase complex was applied to each well, and the plate was incubated at 37°C for 30min. The plate was washed 4 times with wash buffer and the color developing reagent was added to each well of the plate and the plate was incubated at 37°C in darkness for 10–20min. The reaction was terminated by adding a stop solution and the optical density (OD) at 450nm was measured immediately using a microplate reader (Thermo Fisher, Varioskan Flash, USA). Analysis was performed using GraphPad Prism version 9.

Cytospin and Wright–Giemsa staining

NLF was collected immediately from the sacrificed mice by washing the nasal cavity three times with 1mL of ice-cold PBS. The collected lavage fluid was centrifuged at 1200×g for 5min at 4°C. The supernatant was collected and the precipitated cells were treated with red blood cell lysis buffer (Servicebio, G2015, China) for 5min. The precipitated cells were resuspended in 0.2mL of PBS and were counted using a cell counter (JIMBIO, China). NLF cell smears were prepared on cytospin slides and stained using the Wright–Giemsa staining kit according to the manufacturer’s instructions (Beyotime, C0131, China). Briefly, Wright–Giemsa staining solution was added to the NLF cell smears and stained for 45min. Then, the smears were washed with distilled water and dried at room temperature. The images were acquired by a Nikon Eclipse Ni-U epi-fluorescent microscope.

Immunohistochemistry

Paraffin-embedded sections were incubated at 65°C for 1h. Dewaxing, hydration, and antigen repair were performed sequentially as previously described⁹⁰. The endogenous peroxidase blocker was applied to the slides after they had cooled to room temperature. The slides were incubated for 20min at room temperature. The slides were then washed with PBS 3 times and incubated with the primary antibody (see Supplementary Table ² for a complete list and dilutions) in a humidified chamber at 4°C overnight. After washing with PBS, the sections were incubated with a reaction-enhanced solution. Following another wash, the sections were incubated with the secondary antibody for 10min, and the color reaction was developed using 3,30-diaminobenzidine tetrahydrochloride (DAB) (ZSGB-Bio, PV-9000, China). The slides were counterstained with hematoxylin. Finally, the slides were dehydrated and mounted. The images were acquired using a fluorescence microscope (Olympus VS200, Japan).

Hematoxylin and eosin staining (HE staining)

HE staining was performed using the HE staining kit (Beyotime, C0105S, China) according to the manufacturer’s instruction. Sections were dewaxed, hydrated and then washed with PBS. Then, the sections were incubated with hematoxylin for 10s and washed with distilled water for 10min. After that, the sections were differentiated with 1% hydrochloric ethanol for 20s. After a washing step with distilled water for 10min, the slides were stained with eosin for 1min. Following dehydration, clearing and mounting, the slides were ready for image acquisition under a microscope (Olympus VS200, Japan).

Statistics and reproducibility

No statistical analysis was performed to predetermine the sample size. The numbers of samples included in the analyses are listed throughout the figures. For the scRNA-seq data, statistical analyses and graphic production were performed using Python version 3.7.10. The experimental data are presented as mean±SEM or mean with 95% CI, as shown in the corresponding figure legends. Data distribution was assumed to be normal. In the quantitative statistical graph for all stained images, each data point represents the average result from at least three randomly selected fields of view within the same sample. Wilcoxon rank-sum test was used for the Radar plot, applying a two-sided hypothesis. One-way ANOVA, two-way ANOVA, and Welch ANOVA were used to compare multiple sets. Two-tailed Student’s t-tests were used for the comparisons between the two sets. Statistical analyses and graphic production were performed with GraphPad Prism version 9 (GraphPad Software Inc., San Diego, CA, USA). P<0.05 was considered statistically significant.

This research is supported by the National Natural Science Foundation of China (82171106, 82371120, and 81700890 to X.F.), Taishan Scholar Program of Shandong Province (tsqn202103166 to X.F.), and Natural Science Foundation of Shandong Province (ZR2022MH313 to P.W.).