RNAseq analysis of olfactory neuroepithelium cytological samples in individuals with Down syndrome compared to euploid controls: a pilot study

Cognitive evaluation

The cognitive datum was an additional information of the DS group, already present at AGBD association, as in our previous work [19]. Cognitive evaluation was performed by an expert psychologist, by means of the Vineland II scale (Vineland Adaptive Behavior Scales-II-second Edition) [28] and the Leiter-R scale (Leiter International Performance Scale-Revised) [29], considering both verbal and non-verbal abilities. Recently, in DS, a high interindividual variability is reported and various cognitive profiles could emerge in the verbal and non-verbal domains [30, 31].

The Vineland-II scale is a valid and reliable method to measure a person’s adaptive level of functioning. It is helpful in diagnosis and in classifying intellectual and developmental disabilities and other disorders, such as developmental delays, and it is organized within a three-domain structure: communication, daily living, and socialization. The communication scale domain was used as a measure of verbal intelligence.

Although the Leiter-R scale is a test designed for children and adolescents (ages 2–18), it can yield an intelligence quotient (IQ) and a measure of logical ability for all ages. This test provides a non-verbal measure of general intelligence by sampling a wide variety of functions from memory to non-verbal reasoning. A remarkable feature of the Leiter scale is that it can be administered completely without the use of oral language, including instructions, and requires no verbal response from the participant. Because of the exclusion of language, it claims to be more accurate than other tests when testing subjects who cannot or will not provide a verbal response. Leiter contains 20 subtests organized into two domains: visualization and reasoning (VR) and attention and memory (AM). The VR domain is the only domain routinely used at the AGBD Association. Through the different subtests, it is possible to obtain a series of measures connected to intelligence (i.e., reasoning and problem solving).

Olfactory evaluation

Olfactory function was assessed by means of a standardized test battery, the “Sniffin’ Sticks Extended test” (Burghart Company, Wedel, Germany). One DS subject presented quite severe intellectual disability and reduced speech so that finally, 6 out of 7 DS individuals were able to undergo this assessment.

This validated procedure consists of three subtests, namely, threshold (the concentration at which the odor is reliably detected), discrimination (the subject’s ability to distinguish odors), and identification (the subject has to identify 16 different odors, choosing among different options of answer each time). In order to increase the reliability of the measurements, each subject must give an answer (forced-choice paradigm).

During each subtest, the experimenter removes the pen’s cap and the pen’s tip is held for around 3 s approximately 1 cm under both nostrils. All participants were tested blindfolded by a sleeping mask to prevent visual identification of the odorant-containing pens, for the threshold and discrimination test, as required by the procedure. During the identification assessment for DS people, subjects were asked to choose the answer option that they think to be correct after the odor had been presented, with the possibility to read a paper words list linked to pictures of the four choice options, as previously reported in DS [19].

Scores of the three subtests are presented as a composite “TDI score,” the sum of results obtained for threshold, discrimination, and identification measures. This global score represents a reliable measure to estimate the degree of olfactory function and allows for the detection of normosmia (TDI ≥ 30.3), hyposmia (30.3 > TDI > 16), and functional anosmia (TDI ≤ 16). Kobal et al. introduced the term “functional anosmia” in 2000. This definition means that subjects with a TDI score below 16 are considered completely anosmic or to have some olfactory function left, even if not useful in daily life [32–34]. Indeed, a subject with functional anosmia may still perceive a few odors, be able to discriminate between some of them, or even show olfactory event-related potentials [35]. However, this residual olfactory function does not contribute to the enjoyment of food/drink or to the detection of spoiled food or gas leaks.

Olfactory swabbing

After olfactory assessment, all participants underwent olfactory swabbing (DS = 7; euploid controls = 10). An expert otolaryngologist explored nasal cavities using a 30° rigid endoscope. Olfactory swabbing sampling was performed using a sterile disposable nasal swab (Copan Flock Technologies, Brescia, Italy) in both nasal cavities. The human olfactory neuroepithelium is located on the cribriform plate, the superior part of the nasal septum, and the superior and middle turbinate [36, 37]. To minimize discomfort in this kind of individuals, the swab was gently rolled on the mucosa surface at the level of the middle turbinate. We previously showed that, from a surface of 2 cm², it is possible to collect ~ 2 × 10⁶ cells for both nostrils, and among these cells, ~ 30% are ONs [20, 22]. This technique is a gentle approach to collect in vivo olfactory epithelial cells and no participant evaluated this procedure as painful. In particular, the most frequent DS individuals’ opinion after the procedure was of “no pain and no discomfort at all.” Swabs were then immersed in RNA stabilization solution for total RNA extraction. Additional swabs were immersed in a fixative solution (Diacyte, Diapath, Italy) for cytological quality control of the samples.

RNA extraction

Nasal swabs (DS = 7; euploid controls = 10), collected in TRIzol reagent (Invitrogen, Italy), were frozen and stored at − 80 °C. RNA was extracted from each sample using Direct-zol™ RNA Miniprep Plus kit (Zymo Research, CA, USA) following the manufacturer’s instructions. Briefly, the tubes were shaken for 8 s and the swab head removed. Then, each tube was centrifuged at 14,000 × g to remove particular debris and the supernatant was transferred into an RNase-free tube. An equal volume of ethanol (95–100%) was added to each sample and transferred into the columns, and then centrifuged. After two washes, the RNA was eluted. Purification of total RNA was performed using Agencourt RNAClean XP beads (2 × the volume per RNA volume; Beckman Coulter Genomics, Danvers, MA, USA). The concentration and purity of the total RNA samples were measured using the NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE). RNA integrity was assessed with an Agilent 2100 Bioanalyzer and the RNA 6000 LabChip kit (Agilent Technologies, Palo Alto, CA). Total RNA was then verified on Bioanalyzer 2100 (Agilent Technologies) to assess its quality and integrity, to a final RIN of 6.6 ± 1, and then quantified using Qubit RNA HS assay.

RNA preparation and sequencing

Samples were further processed with Lexogen RiboCop rRNA depletion kit to remove the ribosomal content and prepared for sequencing using Lexogen SENSE RNA-seq kit following the manufacturer’s protocol (Lexogen). The 17 samples were sequenced by using the Illumina NextSeq 500 applying the 75-single-end chemistry. The data were deposited with links to BioProject accession number PRJNA789170 in the NCBI BioProject and SRA databases.

RNAseq data analysis

Sequenced reads were trimmed by using cutadapt v1.16 [38] to remove the first 9 nucleotides associated with the library preparation. Trimmed reads were mapped with Salmon v0.9.1 [39] to the Ensembl Homo sapiens GRCh38 cDNA (ftp.ensembl.org/pub/release-92/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz) using v92 of the Ensembl gene annotation (http://ftp.ensembl.org/pub/release-98/gtf/homo_sapiens/Homo_sapiens.GRCh38.98.gtf.gz). Automatic selection of library type (-l A) and aggregated gene-level abundance estimation (–geneMap) were added to standard salmon parameters. Salmon was executed into a quasi-mapping-based mode (salmon quant), and to improve the read mapping process, whose performance could be reduced by a short read length (~ 66 bp), a k-mers size of 21 was chosen to calculate salmon genome index.

The obtained tables were aggregated into a unique file of raw counts and further normalized to account for sequencing depth between samples, using the procedure implemented in the DESeq2 package [40]. Data analysis, statistical testing, and plotting were performed with Python3 and R, exploiting appropriate libraries and packages.

Differential expression analysis

Differential expression (DE) analysis was performed with DESeq2 version 1.22.1 [40] with standard parameters. The full DESeq2 pipeline was applied to raw gene counts to characterize DE genes (DEGs). Genes with an adjusted p-value lower than 0.1 were considered differentially expressed. No specific filter on the fold change was applied.

Olfactory evaluation and olfactory swabbing

In accordance with our previous work [19], all DS individuals (n = 6) showed a clear olfactory deficit in all the three assessed domains (Threshold, Discrimination, Identification). In particular, all DS were markedly hyposmic with one case at the limit of functional anosmia (TDI score: 16.5). All euploid controls (n = 10) were normosmic (Table ​(Table2).2). Olfactory swabbing was bilaterally performed in all the recruited individuals (7 DS individuals and 10 euploid controls) and all of the obtained samples were of good quality, showing both neuronal and non-neuronal cellular component in both DS and euploid controls at light microscopy check, as previously reported [22].

Sequencing results and DE analysis

To evaluate if significant differences in gene expression among DS individuals and euploid controls were detectable, we sequenced the RNA depleted from rRNA of 17 people (7 DS individuals and 10 euploid controls) by means of Illumina NextSeq 500, after the proper RNA quality check (RIN 6.6 ± 1). For each sample, 66.8 ± 18.3 million reads were produced, with a minimum of 37.9 and a maximum of 91, thus granting a high coverage of sequenced transcripts. Several alignment and feature association pipelines were tested (data not shown), finding the best choice in the transcript-level quantification of salmon, accounting for the percentage of the assigned reads to the features (61.2 ± 5.1).

Raw count tables were processed using DESeq2 in order to identify genes showing a differential expression (adj p-value < 0.1) between controls and DS individuals. A total of 52 differentially expressed genes (DEGs) were detected (Table S1), the majority of which has a |logFC|> 0.5, although no filter on fold change was applied. As expected, several DEGs are located in chromosome 21, and genes such as APP, DYRK1A, and DOPEY2 were significantly upregulated in DS individuals (Fig. 1). In addition, we also noticed that the misregulation is spread along the entire genome (Fig. S1): in fact, other interesting DEGs, including MUC16, S100PBP, CREB3L2 and CREB5, which could play a role in the DS related olfactory peripheral impairment, are located outside of the chromosome 21.

Open in a separate window

Fig. 1

Down syndrome DEG heatmap. Cluster heatmap shows samples in columns and genes in rows. The level of expression is represented by the background color, where blue means low and red means high expression. Experimental conditions are shown in green and purple, illustrating euploid controls and Down syndrome (DS) individuals respectively

Regarding the downstream steps, no significant result was showed. A possible explanation for the non-significant enrichment analysis results (related to Reactome pathways and GO) could be the relative low number of detected DEGs. Nevertheless, pathway analysis allowed us to better characterize data found in the previous step. At this regard, we noticed modifications of glycosylation processes, mainly O-linked and mucin associated, mediated by MUC16 and POFUT2. In addition, APP was found to be involved in different processes so that its upregulation could trigger an increase in inflammation and neuronal dysfunction (Fig. S2; Table S2; Table S3).

Olfactory test scores and gene correlation analysis

Pearson’s analysis showed different significant correlations (Fig. 2). In particular, a subset of interesting genes involved in neuronal function, cellular regeneration, and mucus physiology located in the chromosome 21 and also in chromosomes other than 21 was considered for correlations (i.e., APP, DYRK1A, DOPEY2, S100PBP, CREB3L2, CREB5, POFUT2, MUC16, PAX7, KLF9, SPARCL1, ITGA6, STATH). Among the DS upregulated genes (i.e., APP, DYRK1A, DOPEY2, S100PBP, CREB3L2, CREB5, POFUT2, MUC16), a strong significant negative correlation with the global olfactory TDI score emerged. Thus, when gene expression increases, olfactory performance decreases (APP ρ = − 0.87, p-value = 0; MUC16 ρ = − 0.67, p-value = 0; CREB5 ρ = − 0.52, p-value = 0.04; CREB3L2 ρ = − 0.79, p-value = 0; DYRK1A ρ = − 0.71, p-value = 0; DOPEY2 ρ = − 0.62, p-value = 0.01; POFUT2 ρ = − 0.81, p-value = 0; S100PBP ρ = − 0.66, p-value = 0.01). On the other hand, looking at the downregulated genes, a significant positive correlation with the TDI score emerged, namely when gene expression increases, olfactory performance increases. This is particularly clear for the PAX7 and ITGA6 genes (PAX7 ρ = 0.73, p-value = 0; ITGA6 ρ = 0.61, p-value = 0.01).

Open in a separate window

Fig. 2

DEG and olfactory scores correlation matrix. Pearson’s coefficient (lower triangle) and p-value (upper triangle) between the expression of interesting differentially expressed genes in Down syndrome (DS) individuals (upregulated: MUC16, CREB5, CREB3L2, DYRK1A, DOPEYE2, APP, POFUT2, S100PBP; downregulated: PAX7, KLF9, SPARCL1, ITGA6, STATH) and olfactory test scores (TDI, T, D, I) are shown in the matrix. Blue background indicates negative correlation while red indicates positive. p-value background ranges from pale yellow (high significance) to green (no significance). T, threshold; D, discrimination; I, identification; TDI, total smell score

Analyzing the single subtest olfactory scores, the correlation between the threshold score and MUC16 expression (ρ = − 0.52, p-value = 0.04) is particularly interesting. Regarding discrimination score, a significant correlation is showed both with the CREB3L2 (ρ = − 0.68, p-value = 0), DYRK1A (ρ = − 0.73, p-value = 0), POFUT2 (ρ = − 0.77, p-value = 0), and S100PBP (ρ = − 0.56, p-value = 0.02) genes. The identification (I) score shows a strong correlation with both the CREB3L2 (ρ = − 0.72, p-value = 0), DYRK1A (ρ = − 0.68, p-value = 0), and DOPEY2 (ρ = − 0.62, p-value = 0) genes.

Regarding the relationship between downregulated genes and olfactory scores, PAX7 expression was found to be strongly correlated with both threshold and discrimination (ρ = 0.63 and ρ = 0.62 respectively, p-value = 0.01), while ITGA6 showed a positive correlation with the discrimination score (ρ = 0.68, p-value = 0).

In addition, the aforementioned results were also retained by means of the Benjamini/Hochberg FDR method (p-value = 0.1) showing a robust significant correlation between olfactory scores and genes’ expression for the global TDI score (APP ρ = − 0.87, p-value = 0; MUC16 ρ = − 0.67, p-value = 0.01; CREB5 ρ = − 0.52, p-value = 0.06; CREB3L2 ρ = − 0.79, p-value = 0; DYRK1A ρ = − 0.71, p-value = 0.01; DOPEY2 ρ = − 0.62, p-value = 0.02; POFUT2 ρ = − 0.81, p-value = 0; S100PBP ρ = − 0.66, p-value = 0.02) and for the following single subtests: threshold score and MUC16 expression (ρ = − 0.52, p-value = 0.06); discrimination and CREB3L2 (ρ = − 0.68, p-value = 0.01); discrimination and DYRK1A (ρ = − 0.73, p-value = 0.01); discrimination and POFUT2 (ρ = − 0.77, p-value = 0); discrimination and S100PBP (ρ = − 0.56, p-value = 0.05); identification with CREB3L2 (ρ = − 0.72, p-value = 0.01), DYRK1A (ρ = − 0.68, p-value = 0.01), and DOPEY2 (ρ = − 0.69, p-value = 0.01) genes. The same thing occurs considering the DS downregulated gene PAX7 and the TDI score (ρ = 0.73, p-value = 0.01) and ITGA6 and the TDI score (ρ = 0.61, p-value = 0.03) as well as for the single subtest scores: both threshold, discrimination with PAX7 (ρ = 0.63 and ρ = 0.62 respectively, p-value = 0.03), and discrimination with ITGA6 (ρ = 0.68, p-value = 0.02).

The authors warmly thank the “AGBD: Associazione Genitori Bambini Down,” Marzana (Verona), Italy, and all the families who contributed to this research. The authors also thank the Copan Company, Brescia, Italy, for donating the flocked nasal swabs, the “Centro Piattaforme Tecnologiche” of Verona University for providing facilities and technical assistance, and Dr. Alessandro Marcon, Department of Diagnostics and Public Health, Unit of Epidemiology and Medical Statistics, University of Verona, Verona, Italy, for his useful advices. The authors would also like to thank Dr. Stefania Vendramin for proof-reading the paper and Dr. Maria Jimena Ricatti for the graphical support.