Distinct Gene Set Enrichment Profiles in Eosinophilic and Non-Eosinophilic Chronic Rhinosinusitis with Nasal Polyps by Bulk RNA Barcoding and Sequencing

Chronic rhinosinusitis with nasal polyps (CRSwNP) is a chronic inflammatory disease with a high symptom burden, including nasal congestion and smell disorders. This study performed a detailed transcriptomic analysis in CRSwNP classified as eosinophilic CRS (ECRS), nonECRS according to the Japanese Epidemiological Survey of Refractory Eosinophilic Chronic Rhinosinusitis (JESREC) criteria, and a group of ECRS with comorbid aspirin intolerant asthma (Asp). Gene expression profiles of nasal polyps and the uncinate process in CRSwNP patients and normal subjects (controls) were generated by bulk RNA barcoding and sequencing (BRB-seq). A differentially expressed genes (DEGs) analysis was performed using DESeq2 software in iDEP to clarify any relationship between gene expression and disease backgrounds. A total of 3004 genes were identified by DEGs analysis to be associated with ECRS vs control, nonECRS vs control, and Asp vs control. A pathway analysis showed distinct profiles between the groups. A Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) showed distinct phenotype-specific pathways of expressed genes. In the specific pathway of “cytokine–cytokine receptor interaction”, the differentially expressed genes were widely distributed. This study indicates that transcriptome analysis using BRB-seq may be a valuable tool to explore the pathogenesis of type 2 inflammation in CRSwNP.


Introduction
Chronic rhinosinusitis (CRS) is a chronic inflammatory upper airway disease with a high symptom burden, including nasal congestion, smell disorders, nasal discharge, and facial pain/pressure [1,2]. Based on the presence or absence of nasal polyps, CRS can be classified as without nasal polyps (CRSsNP) or with nasal polyps (CRSwNP), which have distinct inflammatory profiles and pathogenesis. CRS can be subclassified into further categories by its underlying cellular and molecular pathophysiologic mechanisms [3], and CRSwNP can also be subclassified into ECRS and nonECRS in terms of the degree of eosinophilic infiltration of sinonasal tissue [4]. Compared with nonECRS, which is supposed to be characterized by Th1-dominant inflammation, ECRS is a phenotype of CRSwNP characterized by Th2-dominant inflammation, which is difficult-to-treat and leads to poorer surgical outcomes [4,5]. Additionally, ECRS patients with comorbid aspirin 2 of 22 intolerant asthma (AIA) are thought to have more difficult-to-treat CRS due to more severe type 2 inflammation [4,6]. Therefore, the type of inflammation in CRSwNP is thought to be a key factor in classifying CRSwNP as a difficult-to-treat CRS. Currently, JESREC criteria are used to classify ECRS into mild, moderate, and severe ECRS [4]. From a transcriptome analysis of nasal polyps, it was certified that nasal polyps can be segregated into two major endotypes showing a type 2 or non-type 2 endotype [7]. However, a previous study using a transcriptome analysis with three ECRS, three nonECRS, and three control subjects showed no distinct differentiation of expressed mRNAs between ECRS and nonECRS [8]. Given the discrepancy between the inflammatory differences of CRSwNP and differentiation of expressed mRNAs in transcriptome analysis, we presume that the severity of eosinophilic infiltration does not enhance the severity of type 2 inflammation or pathological endotype. This study performed a detailed transcriptomic analysis among ECRS and nonECRS based on the JESREC criteria and a group of ECRS with comorbid AIA (Asp) to reveal whether the diagnosis of ECRS and/or Asp correlates with the severity of type 2 inflammation or pathological endotype. Because RNA-seq, a method for transcriptomic analysis, is still too laborious and expensive, we used bulk RNA barcoding and sequencing (BRB-seq), which is a comparatively novel method for RNA sequencing and uses early multiplexing to produce 3 cDNA libraries for dozens of samples [9]. BRB-seq has a comparable performance to the standard TruSeq approach while showing greater tolerance for lower RNA quality and being up to 25 times cheaper [9]. In this study, we applied unsupervised clustering methods from comprehensive RNA gene expression data using BRB-seq and analyzed the difference between these three groups of CRSwNP to gain insights into the molecular profiles and endotype categorization of CRSwNP.

Gene Expression Profiles of NPs between ECRS, NonECRS, and Asp
PCA, Heatmap of the Correlation Matrix, and Hierarchical Cluster Analysis PCA was performed among all groups and revealed that the control group was clearly segregated from the other groups, but ECRS, nonECRS, and Asp were not segregated from each other ( Figure 1A). A heatmap of the correlation matrix also clarified the difficulty of the segregation among the three CRSwNP groups ( Figure 1B).

Segregation All CRSwNPs into Two Types of Clusters by Hierarchic Cluster Analysis
Due to the close association and similar clusters between ECRS, nonECRS, and Asp, we segregated these CRSwNPs into two clusters (cluster1 and 2) using the hierarchical cluster analysis with all CRSwNPs (Figure 6). The two clusters were not correlated with the comorbid of asthma. Additionally, height, weight, BMI, sex, age, comorbid of olfactory disturbance, otitis media, smoking, blood eosinophil, total IgE, nasal polyp score, CT score, tissue eosinophil, JESREC score, Oral FeNO, and Nasal FeNO were not correlated with the segregation (data not shown).
we segregated these CRSwNPs into two clusters (cluster1 and 2) using the hierarc cluster analysis with all CRSwNPs (Figure 6). The two clusters were not correlated the comorbid of asthma. Additionally, height, weight, BMI, sex, age, comorbid of olfa disturbance, otitis media, smoking, blood eosinophil, total IgE, nasal polyp score score, tissue eosinophil, JESREC score, Oral FeNO, and Nasal FeNO were not corre with the segregation (data not shown).  The cluster-specific upregulated and downregulated DEGs (Table S1) were ide between cluster1 vs cluster2, and the KEGG pathway analysis by DAVID with up lated DEGs revealed that the disease-associated terms "cytokine-cytokine receptor action" was commonly highly ranked between cluster1 vs control (cluster1-Ctrl) an ter2 vs control (cluster2-Ctrl) (Table S2).

Discussion
CRS has been categorized using the pathological molecular evidence and syndromic phenotypes, and currently it is believed that strong pathological associations between endotype and phenotypes would be the key factors to resolve the clinical problems of CRS therapies [3]. ECRS involving Asp is known as the most difficult-to-treat type of CRSwNP because of the easy relapse of NPs after sinus surgery and resistance to medical therapy except for oral systemic corticosteroids; nonECRS is known as an easy-to-treat type of CRSwNP with good response to ESS and ordinal macrolide therapy [10]. The pathological differences between ECRS and nonECRS were recognized by the type of inflammation: whether type 2-dependent with eosinophil infiltration or type 1-dependent.
Currently, the criteria for categorizing ECRS have been readily accepted, especially in terms of the disease severity [11][12][13]. To detect clear differences in the type of inflammation and other factors, we performed the transcriptome analysis for ECRS, nonECRS, and Asp. The PCA and heatmap of the correlation matrix showed clear segregation between the control group and the others, but the others were not segregated. Additionally, hierarchical cluster analysis of the DEGs revealed no distinct segregation among the three CRSwNP groups, and some ECRS and nonECRS patients were segregated in a similar cluster as a control group. These results suggested that CRSwNP was widely distributed, but endotype separation using the current criteria for ECRS, which mainly use a count of eosinophil infiltration in NP and blood eosinophil count or positive of comorbid asthma and/or AIA, could not well-segregate CRSwNP into the pathological endotypes. Furthermore, these results also suggested that new criteria for the segregation of CRSwNPs to enhance the pathological endotypes should be proposed to analyze the effectiveness of various therapies.
In our study, a difference in significant enrichment pathways among the CRSwNP groups could not be detected using the DEGs because of no or a low number of DEGs. However, a difference in category-oriented significant enrichment pathways, which came from the comparison between each CRSwNP groups vs the control, could be identified, and "cytokine-cytokine receptor interaction" was mainly upregulated. Additionally, the DEGs associated with "cytokine-cytokine receptor interaction" among the three comparison groups were widely distributed, and the upregulated genes of CCL13 and CCL26 (eotaxin-3) in ECRS and Asp and CCL18 in ECRS were identified in either ECRS or Asp compared with nonECRS. In addition, CCL11(eotaxin) was upregulated in ECRS and nonECRS, and CCL24 (eotaxin-2) was not identified in any CRSwNPs. Previous reports mentioned that type 2 inflammatory transcript expression of CCL13, CCL18, and CCL26 were upregulated in the type 2 endotype [7] and ECRS compared with nonECRS [8]. It is known that CCL11, CCL13, and CCL26 recruit eosinophils, basophils, mast cells, and CD4(+) T helper 2 cells [14], and CCL18 recruits naive T cells, B cells, immature dendritic cells, and activate fibroblasts [15]. CCL18 is thought to be involved in T(H)2-related inflammatory diseases, including asthma and atopic dermatitis, and the number of CCL18(+) cells was significantly increased within NPs [15]. Production of CCL26 is induced by IL4/13 in epithelial cells [16], and CCL26 induces not only eosinophil infiltration but also NK cell migration, which contributes to trafficking NK cells to the allergic upper airway mediated by CX3CL1 and the tyrosine kinase pathway [17]. Higher expression of type 2 inflammatory molecules (CCL13, IGHE, CCL18, CCL23, CCR3, and CLC) along with lower levels of LACRT, PPDPFL, DES, C6, MUC5B, and SCGB3A1 are related to a stronger clinical response to glucocorticoids [18]. In terms of the above-mentioned evidence, these differences would explain the difference of each phenotype as to whether they mainly have type 2 inflammation or not. From the upregulation of these gene expressions, we concluded that ECRS have higher type 2 inflammation with the upregulation of CCL13, CCL18, and CCL26 than nonECRS, and eosinophil infiltration in ECRS would be mainly induced by CCL13 and CCL26 rather than CCL11 and CCL24.
Furthermore, the upregulated expression of TNFRSF18, IL1RL1, and IL2RA in ECRS and Asp, INHBB in ECRS, and TNFRSF10D, IL6, IL11, and IL31RA in Asp is also associated with "cytokine-cytokine receptor interaction". IL1RL1 is the receptor of IL33, which induces T helper 2-type inflammatory responses. IL1RL1-positive cells in the inflammatory cells of the subepithelial layer are significantly higher, and the expression of IL1RL1 in polyps is significantly increased in the ECRS compared with controls. Many IL1RA1-positive eosinophils are observed only in the mucosa of ECRS but not of nonECRS [19]. Therefore, overexpression of IL1RA1 in ECRS and Asp suggested that higher type 2 inflammation would be induced by IL33 and IL1RA1 with CD4(+) T helper 2 cells recruited by CCL11, CCL13, and CCL26. Furthermore, TNFRSF10D is expressed by macrophages/monocytes and promotes eosinophil survival via inhibiting the spontaneous apoptosis [20]. IL11 is selectively expressed in eosinophils and epithelial cells in patients with moderate and severe asthma, and correlates directly with disease severity. IL11 also causes nodular mononuclear infiltrates [21]. IL31 and IL31RA induce MUC5AC gene expression in human airway epithelial cells and play important roles in mucus overproduction [22]. These findings suggested that eosinophil accumulation caused by prolonged survival of eosinophil with TNFRSF10D, nodular mononuclear infiltration with IL11, and mucus overexpression with MUC5AC via IL31RA would be the characteristic features of Asp. In the same way as IL11 enhances the disease severity in asthma, overexpression of IL6 and IL11 in Asp may reflect the severity of IL6-mediated inflammation and type 2 inflammation enhancing the disease severity in CRS.
Overexpressions of TNFRSF18, IL2RA, and INHBB in either ECRS or Asp groups may also reflect the disease severity and pathology of the CRSwNP phenotype via interference of cytokine-cytokine receptor interaction despite poor reports of pathological analyses between the overexpressions of these genes and the phenotype of ECRS and/or Asp.
However, the KEGG pathway analysis for the downregulated DEGs of ECRS-Ctrl, nonECRS-Ctrl, and Asp-Ctrl also revealed that "metabolic pathways" would be another factor for the CRSwNP phenotype. The gene count differences of "metabolic pathways", namely ECRS-Ctrl of 55, nonECRS-Ctrl of 36, and Asp-Ctrl of 65, also pave the way to elucidate the differences in each phenotype associated with disease severity.
Hierarchical cluster analysis using the DEGs revealed no distinct segregation among the three CRSwNP groups, but it did suggest that these CRSwNPs segregate into two clusters in terms of the similarity of up and downregulated DEGs. Two clusters (cluster 1 and 2) among all CRSwNPs, after modifying a hierarchical cluster analysis, showed clear separation in the PCA analysis. The cluster-specific upregulated DEGs and the KEGG pathway analysis by DAVID of two clusters and the control group revealed that the diseaseassociated term "cytokine-cytokine receptor interaction" was commonly highly ranked. Compared to cluster1-Ctrl and cluster2-Ctrl in the DEGs of "cytokine-cytokine receptor interaction", cluster2 was more upregulated in the genes of the CC subfamily and CXC subfamily in "chemokines", gamma-chain utilizing, and IL-4 like and IL-6/12-like genes in "the class I helical cytokines", "IL1-like cytokines", and "TNF family". In the IL4-like genes in "the class I helical cytokines", a comparison between cluster1-Ctrl and cluster2-Ctrl revealed that cluster2 was more upregulated in the expression of CSF2RB and IL4R than cluster1. These genes are deeply associated with type 2 inflammation via IL4, IL5, and IL13 signaling [23,24]. Previously, dupilumab, which is a dual inhibitor of IL4 and IL13 signaling and controls the type 2 inflammation, was shown to decrease several symptoms of CRSwNP despite the diversity in the degree of eosinophilic infiltration of sinonasal tissue by classifying CRSwNP into the categories from nonECRS to severe ECRS [25]. Given the result of the effectiveness of dupilumab against all CRSwNP, all CRSwNP groups should have type 2 inflammation, which could differ in degree regardless of the criteria for ECRS. The upregulated expression of CSF2RB and IL4R in cluster2 is expected to be strongly associated with the effect of dupilumab by blocking this signaling. In addition, the segregation of DEGs into cluster1 and cluster2 will elucidate the broad effectiveness of dupilumab in all CRSwNPs in spite of the phenotypes of nonECRS, which are suggested to be deeply associated with type 1 inflammation. The segregation of DEGs into cluster1 and 2 will become the criteria for true endotype segregation, and further studies are needed to categorize CRSwNP enhancement by/of IL4, IL5, and IL13 signaling.
There are some limitations to our study. First, due to the small sample size in Asp, we could not demonstrate more clear differences in DEGs between ECRS and nonECRS, and the tendency as to whether Asp could be classified as more cluster1 or 2. Second, this study only included Japanese patients; hence, caution should be taken when extrapolating our results to other ethnic groups. Third, we could not measure expression protein levels. A previous report mentioned that the protein expression levels of CCL13 and CCL26 were minimal or unchanged compared with the control despite the overexpression of these mRNAs [26]. Further studies are needed to analyze the correlation among the CRS phenotype, protein expression levels, and transcriptional overexpression of CCL13 and CCL26. Additionally, some protein expression levels may also be minimal or unchanged despite transcriptomic overexpression. Further studies are necessary to explore the correlation between the DEGs detected in BRB-seq and protein expression levels.

Patient Recruitment
Patients with or without CRS, who underwent endoscopic sinus surgery in Hiroshima Medical University Hospital between October 2016 and October 2021, were enrolled in this study. The diagnosis of CRS was based on computed tomography (CT) scanning, patient history, clinical symptoms, and endoscopic findings. The inclusion criteria for CRSwNP were as follows: treatment without oral/nasal steroids within 4 weeks before surgery, no improvement in continuous nasal drip, post-nasal drip, and nasal congestion after medical treatment including low-dose macrolide therapy. Patients with CRSwNP were clinically diagnosed as ECRS or nonECRS based on the diagnostic criteria of the JESREC study [4]. The diagnostic criteria for ECRS (JESREC score) include (1) side affected: both sides, 3 points; (2) with nasal polyps, 2 points; (3) CT changes: ethmoid/maxillary ≥ 1, 2 points; (4) peripheral blood eosinophil count (%): 2 < and ≤ 5%, 4 points; 5 < and ≤ 10%, 8 points; 10%<, 10 points. Eosinophil infiltration into the NPs was diagnosed by calculating the mean cell count of the 3 most dense areas of eosinophils in the hematoxylin and eosin (HE) stained sections under ×400 magnification. A diagnosis of ECRS was made by both the total JESREC score of 11 points or higher and eosinophil infiltration into the NPs of 70 or more [4]. The group of ECRS with comorbid aspirin intolerant asthma was independently classified as Asp (N = 3), and the others were classified as ECRS (N = 9) and nonECRS (N = 8). Control patients (Ctrl) (N = 6) were also diagnosed based on the anatomical abnormalities but showed no inflammatory mucosal change or bacterial infection around the uncinate process. The exclusion criteria for CRSwNP were as follows: treatment with oral/nasal steroids within 4 weeks before surgery, fungal/allergic fungal rhinosinusitis, and primary ciliary dyskinesia. Demographics and clinical characteristics were obtained from the medical records. The detailed patient demography is shown in Table S3.

Data Processing of BRB-Seq
Read1 (barcode read) was extracted by using UMI-tools (ver.1.1.1) with the command "umi_tools extract -I read1.fastq -read2-in=read2.fastq -bc-pattern=NNNNNNNNNNCCC CCCCCC -read2-stdout". Adaptor sequences and low-quality sequences were removed and read length below 20 bp were discarded using Trim Galore (ver.0.6.7). Reads were mapped to the GRCh38 reference using HISAT2 (ver.2.2.1). Read counts for each gene were obtained by feature Counts (ver.2.0.1), and UMI duplication was removed by UMI-tools with the command "umi_tools count -method=unique -per-gene -per-cell -gene-tag=XT". Normalized read count value was obtained by using DESeq2 (ver.1.34.0) (Table S4). Genelevel expression data (read counts) were processed by using the web portal for integrated differential expression and pathway analysis (iDEP.95; http://bioinformatics.sdstate.edu/ idep, accessed on 25 April 2022). iDEP.95 was used for principal component analysis (PCA), hierarchical clustering with a heatmap, identification of differentially expressed genes (DEGs), and pathway analysis. DEGs were extracted with an FDR cutoff of 0.1 and min-fold change of 2. Pathway analysis was performed using GAGE with the genesets of GO Molecular Function. Geneset size was set at minimum 5 and maximum 2000, and pathway significance cutoff (FDR) was set at 0.2. To ascertain the functions of the differential expressed mRNAs, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was used with the Database for Annotation, Visualization, and Integrated Discovery (DAVID; https://david.ncifcrf.gov/, accessed on 25 April 2022).

Conclusions
This study revealed and shed light on the future gene expression targets for endotype categorization. Additionally, this study also revealed that BRB-seq was useful to elucidate the differences in CRSwNP. BRB-seq allowed us to perform a transcriptome analysis for as many as 26 patients because of the low cost of the sequencing. A greater number of patients could help us to segregate into new clusters and elucidate the differences between them. This method may provide a valuable tool for further studies to analyze the differences in CRSwNPs in terms of the phenotype and endotype categorization.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki, and was approved by the Ethics Committee of Hiroshima University Hospital, Hiroshima University, Hiroshima, Japan (Hi-136; August.11.2014 as date of approval).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.
Data Availability Statement: Not applicable.