Transcriptomics unravels molecular changes associated with cilia and COVID-19 in chronic rhinosinusitis with nasal polyps

Chronic rhinosinusitis with nasal polyps (CRSwNP) is a common upper respiratory tract complication where the pathogenesis is largely unknown. Herein, we investigated the transcriptome profile in nasal mucosa biopsies of CRSwNP patients and healthy individuals. We further integrated the transcriptomics data with genes located in chromosomal regions containing genome-wide significant gene variants for COVID-19. Among the most significantly upregulated genes in polyp mucosa were CCL18, CLEC4G, CCL13 and SLC9A3. Pathways involving “Ciliated epithelial cells” were the most differentially expressed molecular pathways when polyp mucosa and non-polyp mucosa from the same patient was compared. Natural killer T-cell (NKT) and viral pathways were the most statistically significant pathways in the mucosa of CRSwNP patients compared with those of healthy control individuals. Upregulated genes in polyp mucosa, located within the genome-wide associated regions of COVID-19, included LZTFL1, CCR9, SLC6A20, IFNAR1, IFNAR2 and IL10RB. Interestingly, the second most over-expressed gene in our study, CLEC4G, has been shown to bind directly to SARS-CoV-2 spike's N-terminal domain and mediate its entry and infection. Our results on altered expression of genes related to cilia and viruses point to the de-regulation of viral defenses in CRSwNP patients, and may give clues to future intervention strategies.

RNA extraction and cDNA synthesis.The 50 tissue samples from CRSwNP patients and healthy controls were directly put in RNA later medium (Thermo Fisher Scientific Inc., CA, USA) and processed for extraction and purification of total RNA using the Kingfisher RNA kit together with the Kingfisher instrument (Thermo Fisher Scientific Inc., CA, USA) according to the manufacturer's instructions.The quantity and quality of isolated RNA was determined using the NanoDrop 2000 Spectrophotometer (Thermo Fischer Scientific) and 2200 TapeStation Automated Electrophoresis System (Agilent Technologies).Samples with an RNA integrity number of greater than 6.0 were chosen for sequencing, with the exception of one polyp sample at RIN 4 which was included.
Whole genome RNA sequencing.Preparation of libraries was carried out using the TruSeq Stranded Total RNA Sample Preparation Kit with Ribo-Zero Gold (Illumina, Inc., San Diego, CA), using 60-1000 ng of total RNA input.The Novaseq 6000 platform was used (Illumina, Inc., San Diego, CA) and 100 bp paired-end reads were generated by Clinical Genomics (Gothenburg, Sweden).
Adapters and low-quality tail were trimmed from reads prior to read alignment.Clean sequence reads aligned to the human genome were used to assemble transcripts, estimate the abundance of these transcripts and detect differential expression among samples.For mRNA and lncRNA analyses, the reference genome build GRCh38/ hg38 was chosen as the annotation references.Fragments per kilo-base of exon per million fragments mapped (FPKM) of both lncRNAs and mRNAs in each sample were calculated based on the length of the fragments and reads count mapped to this fragment.

Data availability statement.
Raw sequence reads in FASTQ format were evaluated in terms of read quality (per base sequence quality, per base G + C content, sequence length distribution, sequence duplication levels, Kmer content and low complexity sequences).The quality statistics were gathered using FastQC 21 (version 0.11.2) available at https:// www.bioin forma tics.babra ham.ac.uk/ proje cts/ fastq c/.The resulting reports were Differential expression analysis.For differential gene expression discovery, the DESeq2 R package was used 29 .A variance stabilizing transformation (vst) was applied to gene expression counts.To identify differentially expressed genes (DEGs), an adjusted p-value < 0.01 and 0.001 were set as thresholds to define the significance in the three comparisons: Control vs non polyp mucosa from CRSwNP patient, non-polyp CRSwNP patient mucosa vs polyp mucosa from the same patient with CRSwNP and finally control vs polyp mucosa from CRSwNP patient.Due to an uneven matching on age and gender between the CRSwNP patient and controls, gender and age were used as covariates in all the statistical analyses using DESeq2.Heatmaps were generated using "pheatmap" in R, R version 4.1.3.

Pathway enrichment analysis.
Lists of the most differentially expressed genes (DEGs) in each comparison were imported into GeneTrail 3.0 30 (https:// genet rail.bioinf.uni-sb.de/ start.html) for pathway enrichment analysis, and the top significantly enriched KEGG 31 , Wiki, Panther and Reactome pathways as well as tissue type/cell marker overrepresentation, provided by GeneTrail 3.0, were determined.

Results
Profiling DEGs in the polyp mucosa of CRSwNP patients and adjacent non-polyp mucosa from the same CRSwNP patient.Between polyp mucosa and non-polyp mucosa from the same CRSwNP patient a total of 1968 genes were significantly different using an adjusted p-value of < 0.001.The most significant DEGs between polyp mucosa from CRSwNP patients and adjacent non-polyp mucosa from the same CRSwNP patient are shown in Table 1 and Fig. 1a.Results from the over-represented pathway analyses between these two comparisons are shown in Table 2. Top significant DEGs up-regulated in the polyp mucosa include genes involved in T-helper 2 (Th2) type response and eosinophil migration, such as CCL13, CCL18, CCL8 12 .C-Type Lectin Domain Family 4 Member G (CLEC4G) gene was also up-regulated and represents the second most significant DEG.Comparative analysis of the polyp and non-polyp mucosa from the same patients revealed Solute carrier family 9 member A3 (SLC9A3) and its antisense gene SLC9A3-AS1 among the top 5 DEGs (Table 1 and Fig. 1a).Additionally, the HISLA gene (alias LINC01146) that is considered as a main "expression regulatory hub" 17 were significantly up-regulated in polyp mucosa from CRSwNP patients, albeit also up-regulated in nonpolyp mucosa of the same patient compared with healthy mucosa of control individuals.The molecular pathway which was the most statistically significant between polyp mucosa and adjacent non-polyp mucosa of the same patient was "ciliated epithelial cell" (p = 8.5 × 10 -78 ).Other pathways, which included significant DEGs were "Metabolic pathways" (108 DEGs, adjusted p = 0.001), "Microglia Pathogen Phagocytosis Pathway" (16 DEGs, adjusted p = 9.97 × 10 -07 , "inflammation mediated by chemokine and cytokine signaling pathway" (24 DEGs, adjusted p = 0.003 and "Neutrophil degranulation" (54 DEGs, adjusted p = 1.32 × 10 -05 ).A heatmap with the top 40 DEGs between polyp mucosa and non-polyp mucosa from the same CRSwNP patient is shown in Fig. 2a.A volcano plot showing the result from all expressed genes is shown in Fig. 3a.

Profiling DEGs between mucosa from healthy control vs non-polyp mucosa from patients with
CRSwNP.Overall, 1494 RNAs were differentially expressed between healthy mucosa of control vs non-polyp mucosa of patients with CRSwNP (Genes p-adjusted < 0.01).The most significant DEGs are shown in Fig. 1b and Table 3. Results from the over-represented pathway analyses between this comparison are shown in Table 4.The most significant pathways were the ones involving natural killer T (NKT) cell (p = 2.4 × 10 -45 ) and herpes simplex virus 1 (HSV-1) infection (4.3 × 10 -19 ).Altered expression of genes related to NKT-cells and HSV-1, may point to the deregulation of viral defense in CRSwNP.
A heatmap with the top 40 DEGs between healthy mucosa of control vs non-polyp mucosa of patients with CRSwNP is shown in Fig. 2b.A volcano plot showing the result of all expressed genes is shown in Fig. 3b.Top significant up-regulated DEGs were PTP4A1, GP2, ELF2, FOSB and FOS.Additionally, AC096642.1,TCEANC2, KLLN and PTCD2 were among the most significantly down-regulated genes in polyp mucosa.
Profiling DEGs between healthy mucosa from control vs polyp mucosa from patients with CRSwNP.Between polyp mucosa and healthy mucosa from control patients, a total of 1733 genes were significantly differentially expressed (adjusted p-value of < 0.001).In Supplementary Table S1 and supplementary Fig. S1, the most significant DEGs between polyp mucosa from CRSwNP patients and healthy mucosa from the same CRSwNP patient are shown.The most significant pathways from the over-represented analyses were, considering the comparison within the same CRSwNP patient, the ones involving cilium processes and cilium Vol:.( 1234567890 -32 ).Out of the ten most associated genes (Supplementary Table S1 and Fig. S1), seven were also among the top associated genes when comparing healthy mucosa and polyp tissue from the same patient (Table 1, Fig. 1a).The significant DEGs in nasal mucosa within these gene-regions were the genes IFNAR1, IFNAR2 and IL10RB located on chromosome 21, LZTFL1, CCR9, RN7SL145P, SLC6A20 and XCR1 located on chromosome 3, MYDGF SEMA6B and TNFAIP8L1 located on chromosome 19 and CFAP73, OAS1, OAS2, OAS3 and RASAL1 on chromosome 12 (Fig. 5a-d).
The impact of gender and age on the transcriptomics profile of nasal mucosa in healthy controls.A total of 1683 genes were significantly differently expressed by age using gender as a covariate, while gender had a significant influence on the expression of 1313 genes (p < 0.05) when correcting for age.
Since age and sex influence the expression of a substantial proportion of genes, we use them as covariates in order to identify differential expressed genes regardless of the influence of age or sex.

Discussion
Using whole transcriptomics analyses with a relatively large sample size of 50 nasal mucosa biopsies, this study has expanded our current knowledge of the transcriptome profile of polyp mucosa of CRSwNP patients integrating sex and age as covariates in the analysis.Our results confirm, and expand, previous studies suggesting the induction of genes involved in inflammation and the dysregulation of genes involved in cilia generation in polyp mucosa [12][13][14][15][16][17][18][19] .We found that the expression of the cilia inhibitory gene LZTF1, located in the most associated gene region of COVID-19, is up-regulated in polyp mucosa, supporting a potential common mechanism in the pathogenesis of CRSwNP and viral infections like COVID-19.Previous studies have shown that loss of cilia related proteins cause anosmia 32 a common symptom of COVID-19 33 , and coronaviruses selectively target ciliated cells, and can lead to the withdrawal of the cilia into the cell body, suggesting that the loss of cilia is likely to cause rhinorrhoea. 34.
The second most up-regulated gene in polyp mucosa was a member of C-type lectin protein family called CLEC4G (alias: LSECtin).This protein interacts with surface glycoproteins from several viruses including the spike protein of SARS coronaviruses 35 .CLEC4G has been shown to bind directly to spike's N-terminal domain of SARS-CoV-2 and allow entry in an ACE2-independent fashion 36 .
In non-polyp mucosa of CRSwNP patients, we have additionally identified DEGs associated with NKT-cell and HSV-1 infection pathways.The most differentially expressed gene was however the protein tyrosine phosphatase 4a1 (PTP4A1).PTP4A1 is believed to play a role in the development and maintenance of differentiating epithelial tissues and enhances cell proliferation, cell motility and invasive activity, and promotes cancer metastasis (UniProt Q93096) 37 .PTP4A1 is higher expressed in non-polyp mucosa from CRSwNP patients compared to healthy controls, as well as in tumours and is strongly down-regulated upon tetrodotoxin treatment 38 .
In the pathway enrichment analysis, HSV-1 infection was the top Kegg pathway identified when controls were compared with non-polyp mucosa from CRSwNP patients.Feng Lan et al. has suggested that the inadequate response of CRSwNP may be associated with a deeper intrusion of viruses 39 and, epithelial damage.Intrusion of HSV-1 into nasal mucosa has been shown to be more extensive in nasal polyp tissue compared to that of non-CRS controls 40 .Nasal polyp tissue may thereby be more susceptible to virus invasion and studies have reported a higher prevalence of respiratory viruses in the nasal mucosa of CRS patients compared to that of controls 41,42 .Additionally, a study by Zaravinos et al. 43 has shown an increase in the prevalence of human papilloma virus and www.nature.com/scientificreports/human herpes virus types 1-7 in human nasal polyposis and that the presence of these viruses likely influences the pathogenesis of the benign nasal tumors.In line with this notion, the transcriptomics study by Peng et al. 19 found that the top GO sets enriched in non-polyp mucosa from patients with CRSwNP compared with healthy controls were the interferon signalling pathway and pathways involved in viral responses.

Sex
The upper airways are specialized in removing airborne pathogens and allergens via effective mucociliary clearance, which is essential to protect the airways from pathogenic insults and prevent pulmonary injury.Ciliary dysfunction may also contribute to the pathogenesis of chronic airway inflammatory diseases such as asthma, allergic rhinitis and CRS, possibly due to the negative impacts of chronic inflammation on mucociliary clearance 44 .Our analysis suggest that gene pathways related to ciliated epithelial cells are significantly altered in nasal polyp tissue, in line with previous observations which show aberrant ciliary marker expression and mislocalization of ciliogenesis markers in CRSwNP patients.This includes the study by Peng et al., 2018  45 that showed the regulator of motile cilia formation FOXJ1 is mis-localized in CRSwNP patients, which correlated with disease severity and the co-existence of allergic rhinitis or asthma.Further, the ciliary ultrastructural marker DNAH5 is reported to be mis-localized in patients with nasal polyps and that the negative modulator of ciliogenesis cp110 is upregulated in CRSwNP patients 46,47 .
Importantly, our study revealed that the gene LZTFL1, which is involved in cilia inhibition and located in the most associated gene region of COVID-19 48 , shows an upregulated expression in polyp tissue from CRSwNP patients.LZTFL1 interacts with a Bardet-Biedl syndrome (BBS) protein complex known as the BBSome, and negatively regulate ciliary trafficking of this complex 49 , and similarly to cp110, negatively influence ciliogenesis.BBS proteins are vital to maintain ciliary function by mediating protein trafficking to cilia, and ciliary dysfunction has been implicated in at least 35 different diseases which collectively affect nearly all organ systems 50 .Higher expression of LZTFL1 in polyp tissue thereby suggests down-regulation of cilia formation in this tissue.However, the mucosa next to the polyps in CRSwNP patients have no different expression of this gene compared to controls.
In vertebrates, cilium-dependent signalling orchestrates important developmental pathways, such as limb development and kidney morphogenesis, and is required for vision, hearing, and smell 51 .Cilium dependent signalling also appears to play a role in obesity 52 and type II diabetes, which have a high incidence amongst patients with compromised cilia function 53 .Endothelial cells, which are ciliated, have been implicated in flow-sensing and vascular hypertension, intracranial blood vessel formation, and atherosclerosis prevention [54][55][56] .Recent studies have further unveiled a novel role of primary cilium in preventing vascular regression 57 where there seems to be a mechanosensory role of primary cilia in vascular hypertension and atherosclerosis.Obesity, type II diabetes and vascular hypertension are all associated with an increased risk of developing a severe illness from COVID-19 58,59 .
Loss of function of ciliary proteins in mice has been shown to affect the olfactory epithelium, causing severe reduction of the ciliated border 31 and as previously mentioned, clinical manifestation suggested to arise from ciliary defects includes anosmia 31 .Secretory and ciliated cells also contain the highest fraction of www.nature.com/scientificreports/SARS-CoV-2-infected cells 60 , and a recent study reported that SARS-CoV-2 preferentially replicates in ciliated cells, damages the ciliary layer, which may result in an impaired mucociliary clearance 61 .Ciliary dysfunction is an emerging theme in COVID-19 pathogenesis and given that ciliated epithelial cells is the most differentially regulated pathway in CRSwNP patients, this could suggest that CRSwNP patients may have an altered susceptibility to be infected with SARS-CoV-2.Taken together, cilium function is altered in nasal polyp tissue and could potentially be an important component of the ability of virus to infect cells after initial contact.Withdrawal of cilia at the cell surface could be a defense mechanism leading the loss of smell as well as fewer infected cells upon viral contact.Polyp formation could therefore perhaps be part of a faulty response or perhaps over exaggerated response to a perceived or real viral threat.
A strength of this study is using a relatively large number of samples from both nasal polyp and non-polyp tissues obtained from the same patients as well as being able to compare with healthy control mucosa.However, this study has a few limitations.Firstly, due to lack of blood samples from the patients we only evaluated upper airway inflammation, but not systemic inflammation (i.e.eosinophilic versus non-eosinophilic) via blood sampling.Secondly, all CRSwNP patients participating in our study were treated with corticosteroids locally in their nasal cavity, which could serve as a confounding factor.A study by Benson et al. 62 investigated the expression of over 20,000 genes in nasal polyps before and after glucocorticoid treatment and found that the largest functional non−polyp mucosa from the same CRSwNP patient.DESeq2 R package analysis, including sex and age as covariates.P−values are unadjusted and shown for the comparisons between non−polyp (unaffected) mucosa versus polyp mucosa and control patient mucosa.

Figure 1 .
Figure 1.Top 20 differentially expressed genes (DEGs) between non-polyp CRSwNP patient and control patient.(a) Polyp mucosa versus non-polyp (unaffected) mucosa from the same CRSwNP patient and (b) Healthy control patient mucosa versus non-polyp (unaffected) mucosa from CRSwNP patient.p-values are unadjusted and shown for the comparisons between non-polyp mucosa versus polyp mucosa and control patient mucosa.DESeq2 R package analysis, including sex and age as covariates.

Figure 2 .
Figure 2. Heatmap of the 40 most differentially expressed genes (DEGs).(a) Polyp mucosa versus non-polyp mucosa from the same CRSwNP patient and (b) Healthy control patient mucosa versus non-polyp mucosa from CRSwNP patient.DESeq2 R package analysis, including sex and age as covariates.Heatmaps were generated using "pheatmap" in R, R version 4.1.3.

Figure 3 .
Figure 3. Gene expression of all expressed genes.Vulcano plot showing the negative logarithm of the p-value on the y-axis and the fold change of all expressed genes from the comparison of (a) Polyp mucosa versus non-polyp mucosa from the same CRSwNP patient and (b) Healthy control patient mucosa versus non-polyp mucosa from CRSwNP patient.DESeq2 R package analysis, including sex and age as covariates.

DEGs associated with COVID-19. We
The distribution of samples on the PC1-PC2 planes visualize the difference between healthy control and CRSwNP samples and if gene expression can be used to correctly classify the patients with CRSwNP.A PCA plot showing the overall relationship of the three tissue types investigated in this study, CRSwNP-polyp, CRSwNP-adjacent non-polyp tissue and healthy control tissue is shown in Fig.4.PCA plots of sex and age are included in the supplementary information, Supplementary Figs.S3 and S4.
also analyzed genes specifically located in four chromosomal regions on chromosome 3, 12, 19 and 21, which show genome wide significant association with COVID-19.

Table 1 .
The 40 top differentially expressed genes (DEGs) in non-polyp nasal epithelial mucosa from patients with CRSwNP versus polyp mucosa in the same patient.The mRNA levels (baseMean) of expressed genes and p-values adjusted using the Benjamini-Hochberg method.BaseMean mean RNA count, lfcSE log2 Fold Change Standard Error, stat Wald statistic Z-score, padj p-value adjusted.

Table 2 .
Enriched pathways among the top 1698 differentially expressed genes (DEGs) from the DESeq2 analysis of polyp versus adjacent non-polyp mucosa in the same CRSwNP patient.Analyzed using Genetrail 3.029 p-values adjusted using the Benjamini-yekutieli method.(Genes p-value adjusted < 0.001, n = 1698).Significant values are in bold. KEGG-

Table 3 .
The 40 top differentially expressed genes (DEGs) in non-polyp nasal epithelial mucosa from patients with CRSwNP versus mucosal tissue from controls.The mRNA levels (baseMean) of expressed genes and p-values adjusted using the Benjamini-Hochberg method.BaseMean mean RNA count, lfcSE log2 Fold Change Standard Error, stat Wald statistic Z-score, padj p-value adjusted.