Twenty-one subjects were enrolled in the Neurocovid project, conducted at the University of Padova, with the aim of investigating potential SARS-CoV-2 related transcriptional signatures in olfactory mucosa that might be linked to anosmia. Before entering the study, approved by the local Ethical Committee, all subjects signed a written informed consent form. The enrolled subjects had a mean age of 34 years old, both sexes and, in most cases, without significant lung involvement. In many of them, anosmia or hyposmia persisted at least 4-6 weeks after COVID-19 diagnosis. Following an anamnestic interview, all the patients were assessed using the Sniffin’ stick test (Burghart Medical Technology, Wedel, Germany). The Sniffin’s stick test is a psychophysical validated smell test consisting of three subtests, which allow the study of odor threshold, discrimination and identification, respectively. The sum of the three subtest results gives a composite score, known as TDI (threshold, discrimination, identification), which can indicate normal olfactory function, hyposmia, or anosmia (28) test scores range from 0 to 45, hyposmia (decreased olfactory function) is defined with a score 16<TDI<30, while a TDI score <16 is defined as anosmia (loss of olfactory function). Olfactory epithelium cells were sampled from each subject by using Copan 491CE.A swab (Copan Italy, Brescia, BS, Italy) under endoscopic guide at the University Hospital in Padova. Immediately after sampling, the swab was inserted in RNA/DNA protection buffer and placed on ice ( Supplementary Figure 1 ).

Total RNA, including miRNAs, was extracted from olfactory epithelium cells by using Monarch^® Total RNA Miniprep kit (New England Biolabs, Ipswich, MA, USA) according to the manufacturer’s instructions. Total RNA was quantified by NanoDrop 2000c (Thermo Fisher Scientific, Waltham, MA, USA). RNA integrity and content of miRNAs (%) in each sample were assessed by capillary electrophoresis with the RNA 6000 Nano LabChip and the Small RNA Nano LabChip, respectively, by using the Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA). Both nostrils were sampled independently, and the obtained cells were used for total RNA extraction ( Supplementary Figure 1 ). For each subject, only the best sample in terms of RNA quantity and quality was selected and used for further analysis. More detailed information about RNA samples is reported in Table 1 . Only samples with RNA Integrity Number (R.I.N.) values higher than six and a percentage of miRNA, calculated on the total of small RNA, between 40% and 65% were used for gene and miRNA expression analysis.

The analysis of miRNA expression profiles was performed using “Agilent SurePrint G3 human, 21^st version (8x60k)” microarray (Agilent Technologies, Santa Clara, CA, USA), which allows the detection of 2,549 human miRNAs (miRBase 21.0^th version) and 76 viral miRNAs (GEO platform N. GPL24741). Each miRNA was targeted by 16 to 20 array-probes of different sizes. Total RNA (200 ng) was labeled with pCp Cy3, according to the Agilent’s protocol, and unincorporated dyes were removed with MicroBioSpin6 columns (BioRad, Hercules, CA, USA) (29). Probes were hybridized at 55°C for 22 hours using the Agilent’s hybridization oven, which is suitable for bubble-mixing and microarray hybridization processes. Slides were examined using Agilent microarray scanner (model G2565CA) at 100% and 5% sensitivity settings. Agilent Feature Extraction software version 12.0.0.7 was used for image analysis of miRNA expression arrays. Raw miRNA data are available in the U.S National Centre for Biotechnology Information Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database with the Accession GSE209806.

The RNA-Seq of each individual sample was carried out from IGA Technology Services (Udine, Italy). cDNA libraries were constructed with 100 ng of total RNA by using “Universal Plus™ Total RNA-Seq with NuQuant kit” (Tecan Genomics, Redwood City, CA) following the manufacturer’s instructions. The workflow consists of fragmentation of total RNA and cDNA synthesis with a mixture of random and oligo (dT) primer, followed by end repair to generate blunt end, ligation of UDI adaptors, strand selection, AnyDeplete to remove unwanted transcript, such as ribosomal RNA, and PCR amplification to generate the final library. The libraries were quantified with the Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and quality tested by Agilent 2100 Bioanalyzer High Sensitivity DNA Assay. Sequencing was carried out in paired-end mode (150 bp) by using NovaSeq 6000 (Illumina, San Diego, CA) with a targeted sequencing depth of about 80 million reads per sample. Raw data were processed with the software CASAVA v1.8.2 (Illumina) for both format conversion and demultiplexing. Sequence reads are available on NCBI BioProject database with the accession number PRJNA806721. Raw reads were trimmed to remove adapter sequences using cutadapt (version 3.4). The abundances of all human transcripts annotated by GENCODE (release 38) were estimated using Salmon software (version 1.5.2) (30) and then summarized at the gene level using tximport (version 1.20.0) (31). Genes were filtered genes by their expression levels using the strategy described in Chen et al. (32), as implemented in the edgeR package with default parameters, 22.612 genes were retained. Sample P21 (Group 1) was removed from the analysis due to the low amount of total RNA available.

Inter-array normalization of miRNA expression levels was performed with cyclic Lowess for miRNA (33), the average of replicates being used. Feature Extraction software (Agilent Technologies) was employed to obtain spot quality measures for evaluating the quality and the reliability of the hybridization. In detail, the flag “glsFound” (set to 1 if the spot had an intensity value significantly different from that of the local background, 0 otherwise) was used to filter out unreliable probes: a flag equal to 0 was noted as “not available” (NA). In order to make a robust and unbiased statistical analysis, probes with a high proportion of NA values were removed from the dataset. NA (44%) was used as threshold in the filtering process, a total of 210 available human miRNAs being obtained. Differentially expressed miRNAs were identified with two class-Significance Analysis of Microarray (SAM) algorithm (34) with default settings. SAM, which uses a permutation-based multiple testing algorithm, associates a variable false discovery rate (FDR) with the significant genes. FDR, which refers to the percentage of error that can occur in the identification of the statistically significant differentially expressed miRNAs in multiple comparisons, can be manually adjusted (FDR < 0.05). Samples P16 (Group 1), P6 (Group 3), and C9 (Group 4) were excluded from the analysis due to issues in the miRNA labelling step; sample P16 was replaced with sample P12, belonging to the same Group but not selected for gene expression analysis.

Gene-level counts deriving from RNA-Seq were normalized using RUVseq (version 1.26.0; RUVg method, k=7 confounding factors) (35). Differential expression was tested with edgeR (version 3.34.1) (36), using the GLM model. Genes with an adjusted p-value (FDR) < 0.10 after correction for multiple testing (Benjiamini-Hochberg method) were considered differentially expressed. Finally, to analyze the functional relationship of these differentially expressed genes, a Gene Ontology (GO) functional enrichment analysis through the ShinyGO tool (FDR < 0.05) was performed (37). Sample C10 (Group 4) was excluded from the analysis due to his low sequencing depth.

All the heat maps were obtained by Morpheus software (https://software.broadinstitute.org/morpheus, Broad Institute, USA) using an unsupervised two-dimensional hierarchical clustering approach with the average linkage method and Euclidean correlation.

First-strand cDNA synthesis was performed with “GoScript™ Reverse Transcriptase kit” (Promega, Madison, WI, USA) starting from 500 ng of total RNA in a final volume of 20 μL according to the manufacturer’s instructions. qRT-PCRs were performed in triplicate using Bio-Rad CFX 384 Touch System (Bio-Rad, Hercules, CA, USA) and “GoTaq^® qPCR Master Mix” chemistry (Promega, Madison, WI, USA). The qPCR cycling conditions were 95°C for 2 min, 39 cycles (95°C for 25 s and 60°C for 1 min), and a final step at 72°C for 3 min. The 2^-ΔΔCt (RQ, relative quantification) and ΔCt (Ct_(GOI) – Ct_{(end, ctl)}) method were used to calculate the relative expression ratio and individual expression level, respectively. The oligonucleotides employed are shown in Supplementary Table 1 , Beta-2-microglobulin (B2M) was used as endogenous control. Samples P4 (Group 2) and P22 (Group 3) were not included because of insufficient remaining amount of total RNA.

For the quantification of miRNA expression levels, cDNA was synthesized using the miRCURY LNA RT Kit (Qiagen, Hilden, Germany) starting from 10 ng of total RNA with the addition of 0.5 μL of UniSp6 as exogenous miRNA spiked-in control. PCR was performed in a 10 μL volume containing 5 μL 2x miRCURY SYBR Green Master Mix (Qiagen, Hilden, Germany), 3 μL cDNA of a dilution 1:60, 1 μL RNase-free water, and 1 μL of one of the following miRCURY LNA PCR primer sets (Qiagen, Hilden, Germany): hsa-miR-16-5p (ID YP00205702), hsa-miR29a-3p (ID YP00204698), hsa-miR-21-5p (ID YP00204230) and UniSp6 (ID ZP00004674). The qPCR reactions were performed in a Bio-Rad CFX 96 Touch System (BioRad, Hercules, CA, USA). The qPCR cycling conditions were 95°C for 2 min and 40 cycles (95°C for 10 s and 56°C for 1 min). Three replicates of each sample were amplified for each real-time PCR reaction. The relative expression levels between samples were calculated using the comparative delta Ct (threshold cycle number) method (2^-ΔΔCt). hsa-miR-16 and C6 patient were used as endogenous control and as calibrator sample respectively.

For this study, 21 subjects were enrolled after at least 4 months from the infection. All the subjects were tested for smell sensitivity by adopting the Sniffin’ sticks and sorted into four Groups according to their olfactory symptoms ( Supplementary Figure 1 ). Group 1 was composed by 6 patients recovered from COVID-19 infection by at least 3 months but still manifesting olfactory symptoms. Group 2 was composed by 5 patients who experienced olfactory symptoms during the infection, that resolved 3-4 weeks after recovery. Group 3 was composed by 5 patients who did not experience olfactory symptoms during and after the infection. Group 4 was composed of 5 healthy subjects who were never been diagnosed with SARS-CoV-2 infection and never complained of smell loss. Clinical details of the enrolled subjects relevant for this study are reported in Table 1 .

In order to check the enrichment in olfactory epithelium cells, the expression levels of Keratin 5 (krt5), marker of horizontal basal cells, and cytochromes cyp2j2 and cyp2a13, markers of sustentacular cells, were quantified and compared to those of respiratory epithelium samples; the ΔCts were used for a principal component analysis (PCA) (38–40). As shown in Figure 1 , PCA discriminates all the olfactory and respiratory epithelium samples in two different Groups, demonstrating a clear enrichment in olfactory epithelium cells in the samples of interest, with the only exception of P7 sample who has been removed from the following analysis due to its poor enrichment in olfactory epithelial cells.

A microarray analysis was performed to identify specific microRNA (miRNA) expression signatures in patients with persistent olfactory symptoms (Group 1) and patient who never had olfactory symptoms (Group 3), using healthy subjects (Group 4) as controls. An unsupervised hierarchical clustering analysis, by using the expression levels of 210 detected miRNAs, was able to clearly separate SARS-CoV-2 patients from healthy controls, but not the patients with different olfactory symptoms ( Figure 2A ). 57 differentially expressed miRNAs were found between Group 1 and healthy controls (Group 4) ( Supplementary Table 2 ), 21 miRNAs were differentially expressed between Group 3 and Group 4 (12 of which were also found in the comparison between Groups 1 and 4, Figure 2B ) ( Supplementary Table 3 ) and 4 miRNAs were differentially expressed between Group 1 and Group 3 with FDR < 0.05 ( Supplementary Table 4 ). The large majority of differentially expressed miRNAs in these comparisons were under-expressed in patients with respect to healthy controls, even months after negative PCR tests, suggesting that the expression levels of these miRNAs is affected by SARS-CoV-2 infection well beyond its acute phase ( Figures 2C, D ). Among the overexpressed miRNAs, we identified miRNAs involved in neuronal development (let-7 family, including let-7a-5p, let-7f-5p, and let-7g-5p, Supplementary Figure 2 ) (41) and miRNAs with pro-inflammatory activity (miR-34 family, including miR-34a-5p and miR-34b-5p, Supplementary Figure 2 ) (42). In addition, miR-21-5p and miR-29a-3p, that take part in immune cells differentiation, were also overexpressed in patients with respect to healthy controls ( Figure 3A ). The expression levels of miR-21-5p and miR-29a-3p were also confirmed by using qRT-PCR ( Figure 3B ).

To assess whether the SARS-CoV-2 infection was accompanied by differential expression of specific set of genes, we performed an RNA-Seq analysis of total RNA extracted from 18 subjects ( Table 1 ). A principal component analysis (PCA) clearly revealed that patients belonging to different Groups are characterized by specific gene expression signatures. The PCA was able to clearly stratify patients by attributing them to the Groups in which they were classified from a clinical point of view. Gene expression signatures clearly separate SARS-CoV-2 patients from healthy controls and also patients with different olfactory symptoms ( Figure 4 ). Interestingly, patient P5 (Group 2) is localized closer to the patients of Group 1 and this could be due to the fact that this patient at the moment of olfactory cell sampling reported to the clinicians to have yet some olfactory symptoms even if its TDI (35.25) indicated normosmia. We decided to exclude this sample from the subsequent analysis because it might have been improperly included in the Group 2. Paired comparisons revealed 206 (G1 vs. G2, 118 up and 88 down), 637 (G1 vs. G3, 459 up and 178 down), and 307 (G1 vs. G4, 195 up and 112 down) (LFCT=1, BH adjusted p-value < 0.05) differentially expressed genes (DEGs). Further, 412 (316 up and 96 down) and 810 (379 up and 431 down) DEGs were found between patients with a full recovery of olfactory perception (Group 2) and those with no olfactory symptoms (Group 3) and the healthy ones (Group 4) respectively. Finally, 676 DEGs (147 up and 529 down) were found between healthy controls and patients with no olfactory symptoms (Group 3) ( Supplementary Table 5 ). Interestingly, the large majority of differentially expressed genes were downregulated in SARS-CoV-2 patients (Group 2 and 3) with respect to healthy controls. On the contrary, we observed a clear up-regulation of DEGs in patients with persistent olfactory symptoms (Group 1) with respect to patients with transient olfactory symptoms (Group 2), the ones with no olfactory symptoms (Group 3) and healthy controls (Group 4).

Both the comparisons between patients without olfactory symptoms at the time of the cells sampling (Group 2 and Group 3) versus healthy controls (Group 4) are characterized by DEGs that showed an enrichment of biological processes involved in vasculature and endothelium development as well as in immune activities ( Figure 5A ) known to be triggered by bacterial/Lipopolysaccharide (LPS) stimuli ( Supplementary Figures 3, 4 , Supplementary Tables 6–8 ). In addition, the comparison between patient recovered from the olfactory symptoms and healthy controls (Group 2 vs. Group 4) also highlighted an enrichment in interleukin-17 (IL-17) and tumor necrosis factor (TNF) pathways, suggesting an immune response characterized by Th17 cells in the olfactory epithelium ( Figure 5B , Supplementary Figure 5 , Supplementary Table 9 ) of those individuals. Finally, the comparisons between patients with persistent olfactory symptoms (Group 1) and healthy controls (Group 4) showed an enrichment in biological processes involved in neutrophils, myeloid and granulocytes activities ( Figure 5C , Supplementary Figure 6 , Supplementary Table 10 ). Nine differentially expressed genes obtained from these comparisons were successfully validated by qRT-PCR, confirming the RNA-Seq data ( Supplementary Figure 7 ).

Interestingly, the lists of differentially expressed genes obtained in the following comparisons G1 vs. G4, G2 vs. G4 and G3 vs. G4 highlighted common DEGs ( Figure 5D ) that could be related to infection with SARS-CoV-2. Only two genes (ITGAX and FCGR2A) of the list of common DEGs (73 genes, and Supplementary Table 11 ) were found among G1 vs. G4 belonging to “neutrophil and granulocyte activation” biological process. If we consider the common DEGs between the G1 vs. G4 and G3 vs. G4 comparisons (107 genes, Supplementary Table 11 ), in addition to ITGAX and FCGR2A genes we find three other genes (HBB, OLR1, SLC11A1) that are involved in the inflammatory process. Overall, this analysis highlights an inflammatory/immune response triggered by SARS-CoV-2 infection. While this finding was not unexpected, it is surprising that this molecular signature is still present after months from COVID-19 diagnosis.

Notably, the GO analysis performed on DEGs between patients with persistent olfactory symptoms (Group 1) and patients who never experienced these symptoms (Group 3) after the infection revealed a statistically significant enrichment of genes involved in detoxification of inorganic compound (MT1E, MT1F, MT1G, MT1H and MT2A, Figure 5E , Supplementary Figure 8 , Supplementary Table 12 ). These genes encode for different metallothioneins, small proteins which bind and sequester zinc to reduce its toxic effects in the cell. Accordingly, their expression levels were shown to rapidly increase after zinc exposure or inflammatory processes (43). Interestingly, we observed a relationship between expression of metallothioneins and the olfactory symptoms reported by patients of different Groups. The level of expression of metallothioneins, obtained by RNA-Seq ( Figure 6A ) and qRT-PCR ( Figure 6B ), were the lowest in patients who never experienced olfactory disorders, intermediate in patients fully recovered from olfactory disorders, and the highest in patient with persistent olfactory symptoms.