Characterisation of Neurospheres-Derived Cells from Human Olfactory Epithelium

A major problem in psychiatric research is a deficit of relevant cell material of neuronal origin, especially in large quantities from living individuals. One of the promising options is cells from the olfactory neuroepithelium, which contains neuronal progenitors that ensure the regeneration of olfactory receptors. These cells are easy to obtain with nasal biopsies and it is possible to grow and cultivate them in vitro. In this work, we used RNAseq expression profiling and immunofluorescence microscopy to characterise neurospheres-derived cells (NDC), that simply and reliably grow from neurospheres (NS) obtained from nasal biopsies. We utilized differential expression analysis to explore the molecular changes that occur during transition from NS to NDC. We found that processes associated with neuronal and vascular cells are downregulated in NDC. A comparison with public transcriptomes revealed a depletion of neuronal and glial components in NDC. We also discovered that NDC have several metabolic features specific to neuronal progenitors treated with the fungicide maneb. Thus, while NDC retain some neuronal/glial identity, additional protocol alterations are needed to use NDC for mass sample collection in psychiatric research.


Introduction
Compared to most other disease research, the arsenal of tools for psychiatric research is limited by the lack of reliable biomarkers, the small selection of animal models, and the poor availability of brain tissue. In recent years, studies of genetic and epigenetic factors of mental illness have revealed that most of them are associated with impaired embryonic development of the brain [1][2][3][4]. Since it is impossible to obtain foetal brain samples from future mentally ill people, the study of neuronal stem cells in vitro is the next best thing. This makes the use of neuronal cells in vitro a priority object for model research in psychiatry [5,6]. A common procedure for obtaining such cells is to grow neuronal progenitors and neurons from induced pluripotent cells. This model was used to show a decrease in the expression of protocadherins in schizophrenia [7], impaired β-catenin/BRN2 cascade for idiopathic autism [8], lithium-dependent hyperexcitability [9] and altered calcium signalling pathway [10] in bipolar disorder. Such studies are usually limited to small numbers, often less than 10 samples. The problem with this small sample

Sample Collection
The study included 11 inpatients (5 women, 6 men, 25-41 y.o.) admitted to the Sverzhevsky Institute to fix the curvature of the nasal septum or hypertrophic rhinitis. Sampling of biological material from the nasal cavity was performed under general anaesthesia in the area of the upper third of the nasal septum (superior turbinate), which is a convenient site for a biopsy that contains most of the olfactory epithelium according to our and other studies [33,34]. Inclusion criteria for the study were age (20-45 y.o.) and the absence of severe somatic and mental illnesses. Exclusion criteria were the absence of chronic rhinosinusitis or acute allergic rhinitis. The collected material was placed in a tube containing 1 mL of DMEM/HAM F12 medium supplemented with foetal bovine serum (10%) and Antibiotic-Antimycotic, Gibco (1%). Samples were delivered to the laboratory within 3 h of biopsy collection.

NS and NDC
NS were obtained according to the protocol described earlier [19][20][21][22]. Biopsies were incubated with Dispase II (2.4 IU/mL) for 1 h at 37 • C. Further, the olfactory epithelium was separated from the lamina propria. Then the lamina propria was divided into samples with a thickness of 200 to 500 µm and each sample was placed in a standard 6-well plate under a cover glass, where it was cultured for 18 days (5% CO 2 , 37 • C) in DMEM/HAM F12 supplemented with 10% foetal calf serum and antibiotics. The medium was changed every 2-3 days. After that, the cells were removed with trypsin and transplanted into flasks coated with poly-L-lysine with a density of 16,000 cells per cm2 to form NS (culture medium based on DMEM/F12 with the addition of GlutaMAX, 1% ITS, 50 ng/mL EGF, 25 ng/mL FGF2). After 15-20 days, the formed neurospheres were collected. Some of them were frozen for RNA isolation using the protective reagent RNALater, the remaining ones were dissociated with trypsin and transplanted onto Petri dishes, cultured in DMEM/F12 medium supplemented with 10% FBS, GlutaMAX, and an antibiotic. After three passages, the cells were harvested and frozen for RNA isolation.

Fluorescent Microscopy
NS obtained according to the protocol described above were collected and dissociated with 0.25% Trypsin-EDTA solution, after which they were planted in a 24-well culture plate treated with poly-L-lysine, using a nutrient medium based on DMEM/F12 supplemented with glutamine, ITS 1%, EGF (50 ng/mL), FGF2 (25 ng/mL), and incubated for a day with 5% CO2, 37 • C. NDCs were seeded in 24-well plates using the standard design described above. NS and NDC were washed three times with PBS and fixed in 4% paraformaldehyde in PBS (pH 7.4) for 15 min at room temperature, after which they were washed again with PBS. Further, the cells were supplemented with primary antibodies diluted (v/v 1:100) in a "block solution" (DPBS, 10% FBS, 0.1% Triton X-100, 0.01% Tween 20) and incubated for 12 h at 4 • C. After incubation, the cells were washed three times with PBS solution, after which they were supplemented with secondary antibodies diluted (v/v 1:1000) in a block solution and incubated for 1 h at 37 • C. Then, the cells were washed three times with PBS solution and DAPI solution was added for nuclear staining for 15 min at room temperature. Afterwards, the cells were once again washed with PBS solution. Control wells were processed in the same way, except for the addition of primary antibodies, as well as any antibodies for "pure" control.

RNA-seq and Expression Analysis
RNA was isolated from the cells using RNeasy column kits (Qiagen, Hilden, Germany), then treated with DNase I and purified using the CleanRNA Standard kit (Evrogen, Moscow, Russia). The quality and quantity of the RNA were assessed with the Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). cDNA libraries were constructed with 125 ng (NS) or 500 ng (NDC) of total RNA using NEBNext Ultra II Directional RNA Library Prep (NEB, Ispwich, MA, USA) with NEBNext Poly(A) mRNA Magnetic Isolation Module and barcoded with dual index primers for pooled sequencing. The sequencing was performed to yield 10-20 million reads (2 × 150 base pairs) per library. The gene-level estimated counts were obtained with the "salmon" tool [35]. The analysis for differentially expressed genes between 11 pairs of NS and NDC was performed with the "DESeq2" R package with sample ID used as a covariate (Love, Huber, and Anders 2014). Transcriptomic cell typing was performed using t-SNE with publicly available datasets; the SRA accession numbers are in Supplementary Table S1. Public data were reanalysed from raw FASTQ reads with salmon the same way as described above. For plots on Figures 2B and 4 we used combined TPM (transcripts per million) normalised data filtered by ensembl genes which have TPM > 5 in at least 90% of all datasets. The log-transformed TPM values were used for dimensional reduction. We utilised R packages "Rtsne" for t-SNE with default settings (Barnes-Hut implementation of t-SNE with initial PCA with 50 principal components and perplexity parameter specified in the text) and "umap" for UMAP with default settings ("naive" method). The same data were used for Spearman's rank correlation computations. Gene Set Enrichment Analysis was performed with "msigdbr" and "clusterProfiler" R packages with default parameters on a list of genes, ordered by -log(p-level,adjusted) * sign(FC), where p-values and FC (fold change of expression) were taken from Deseq2 results. Net plots visualisations were generated with the "enrichplot" R package [36] and the visualisation of biochemical pathways in Figure S5 was created with the help of "pathview" R package [37].

Results
We collected nasal biopsies from eleven people. We then obtained NS from the collected material and utilised a protocol described in the Methods section to grow a Cells 2021, 10, 1690 5 of 15 monolayer of NDC. NDC on serum-rich media produce an abundant monolayer of cells, which can be harvested and stored in nitrogen for later use.
NDC and NS were both stained with vascular markers, SMA and TAGLN ( Figure 1), but NDC cultures have a portion of cells with a much more prominent SMA signal ( Figures 1A and S1). The expression data support this observation: NDC have much higher ACTA2, SMA gene, (log2FC = 3.65, p = 4.46 × 10 −22 ), and TAGLN expression (log2FC = 4.87, p = 9.5 × 10 −43 ). It seems that while NDC still retain glial (GFAP) and neuronal (beta-tubulin and nestin) markers, they are shifted more towards other cell types.
From the expression analysis, it is clear that the transcriptomes of NS and NDC are different ( Figure 2A). To ascertain the cellular identity of NDCs, we compared the obtained transcriptomes with previously published bulk RNAseq data of different neural progenitor cells (NPC), glial, nasal, and other (listed in Supplementary Table S1, note that there are no available RNAseq data for OEC). The clustering picture of transcriptomes on t-SNE plot suggests that while NS have transcriptomes very similar to previously published expression data from nasal biopsies, NDC are closer to the transcriptomes of cells from smooth muscle tissue and fibroblasts and not NPC ( Figure 2B). Rank correlations between tested expression profiles reveal that NS to NDC transition leads to loss of similarity with glial cells (ρ falls from 0.74 to 0.58 on average) and NPC (0.51 to 0.38). NDC, like NS, remain similar with smooth muscle cells (0.87 to 0.78 with maximum similarity with fetal aorta cells, 0.82 to 0.84 [38]) and to less extend with fibroblasts (0.82 to 0.73) (Supplementary Figure S2A). Note that the same analysis for Matigian et al. array expression data suggests that expression profiles from ONS also resemble data from fibroblasts more than NPC, glial cells, and cells in this study (Supplementary Figure S2B).
The differential expression analysis (Deseq2) of NDC vs. NS revealed that there were multiple changes in gene expression ( Figure 2C, Supplementary Table S2). We had enough data for gene set enrichment analysis (GSEA) which we used to elaborate the processes involved in NS to NDC transformation. GSEA on gene ontology collection of cellular components (GO:CC) shows that energy and protein metabolism, followed by ontologies related to neuronal components, are affected the most during NS to NDC transition (Supplementary Figure S3 and Table S3). GSEA on gene ontology collection of biological processes (GO:BP) revealed a plethora of affected pathways (Supplementary Figure S4 and Table S4). There are three main themes among them: upregulated in NDC metabolic processes (mainly energetic), downregulated vascular transformation, and neuronal morphogenesis (Supplementary Figure S4). The same three themes are present in the GSEA analysis of pathway enrichment for KEGG [39] (Supplementary Table S5), WikiPathways [40] (Supplementary Table S6), and MSigDB Hallmark [41] (Supplementary Table S7) pathways collections. Some representative GSEA enrichment score curves are presented in Figure 3A. The interaction of affected genes in the "vascular" and "neuronal" categories are presented on the net plot in Figure 3B. Note that both groups of genes contain considerable intersection, involved both in neuro-and vasculogenesis, namely MAPK1, VEGFA, APOE, NRP1, and others.  correlations between tested expression profiles reveal that NS to NDC transition leads to loss of similarity with glial cells (ρ falls from 0.74 to 0.58 on average) and NPC (0.51 to 0.38). NDC, like NS, remain similar with smooth muscle cells (0.87 to 0.78 with maximum similarity with fetal aorta cells, 0.82 to 0.84 [38]) and to less extend with fibroblasts (0.82 to 0.73) (Supplementary Figure S2A). Note that the same analysis for Matigian et al. array expression data suggests that expression profiles from ONS also resemble data from fibroblasts more than NPC, glial cells, and cells in this study (Supplementary Figure S2B).  Table S1 in the Supplementary File. (C) Volcano plot of NDC vs. NS results of differential expression analysis with Deseq2.
The differential expression analysis (Deseq2) of NDC vs. NS revealed that there were multiple changes in gene expression ( Figure 2C, Supplementary Table S2). We had enough data for gene set enrichment analysis (GSEA) which we used to elaborate the processes involved in NS to NDC transformation. GSEA on gene ontology collection of cellular components (GO:CC) shows that energy and protein metabolism, followed by ontologies related to neuronal components, are affected the most during NS to NDC transition (Supplementary Figure S3, Supplementary Table S3). GSEA on gene ontology collection of biological processes (GO:BP) revealed a plethora of affected pathways (Supplementary Figure S4, Supplementary Table S4). There are three main themes among them: upregulated in NDC metabolic processes (mainly energetic), downregulated vascular To understand why NDC are further from neuronal cell identity than NS, we compared our cells with a large publicly available cytotoxicity RNAseq study of in vitro 2D-cultures of neuronal progenitors [42]. Various dimensional reduction plots demonstrate that the transcriptomic position of NPC treated with the fungicide maneb (manganese ethylenebis-dithiocarbamate) is very robustly found colocalized with positions of data from NDC ( Figure 4). Another notable feature (less prominent by a wide margin though) is cadmiumtreated NPC; the effect is best seen on PCA plot ( Figure 4B). Spearman's rank correlation of expression profiles of NDC and Kuusisto et al. data were not the highest for NPC treated with maneb (the highest correlation was for trans-retinoic acid, ρ = 0.73), but for maneb we observed the highest and the only major increase in similarity from NS to NDC (0.59 to 0.7) for all of the tested datasets.
themes are present in the GSEA analysis of pathway enrichment for KEGG [39] (Supplementary Table S5), WikiPathways [40] (Supplementary Table S6), and MSigDB Hallmark [41] (Supplementary Table S7) pathways collections. Some representative GSEA enrichment score curves are presented in Figure 3A. The interaction of affected genes in the "vascular" and "neuronal" categories are presented on the net plot in Figure 3B. Note that both groups of genes contain considerable intersection, involved both in neuro-and vasculogenesis, namely MAPK1, VEGFA, APOE, NRP1, and others. To understand why NDC are further from neuronal cell identity than NS, we compared our cells with a large publicly available cytotoxicity RNAseq study of in vitro 2Dcultures of neuronal progenitors [42]. Various dimensional reduction plots demonstrate that the transcriptomic position of NPC treated with the fungicide maneb (manganese ethylene-bis-dithiocarbamate) is very robustly found colocalized with positions of data from NDC (Figure 4). Another notable feature (less prominent by a wide margin though) is cadmium-treated NPC; the effect is best seen on PCA plot ( Figure 4B). Spearman's rank correlation of expression profiles of NDC and Kuusisto et al. data were not the highest for NPC treated with maneb (the highest correlation was for trans-retinoic acid, ρ = 0.73), but  for maneb we observed the highest and the only major increase in similarity from NS to NDC (0.59 to 0.7) for all of the tested datasets.

Discussion
Our results show that the key problem of growing cells of neuronal origin from NS obtained from the olfactory epithelium is the balance between neural cells and non-neural cells (vascular or stromal cells, for example). The cells, originating from NS, often end up as a mixture of stem cells [43]. It is possible to further optimise cell development towards desired types. For example, in an article by Begum et al., the production of neurons from neurospheres was achieved by increasing the concentration of CO 2 with a serum-free

Discussion
Our results show that the key problem of growing cells of neuronal origin from NS obtained from the olfactory epithelium is the balance between neural cells and non-neural cells (vascular or stromal cells, for example). The cells, originating from NS, often end up as a mixture of stem cells [43]. It is possible to further optimise cell development towards desired types. For example, in an article by Begum et al., the production of neurons from neurospheres was achieved by increasing the concentration of CO 2 with a serum-free medium at the stage of neurospheres [44]. Possibly, under conditions of greater hypoxia with less reactive oxygen species (ROS) production, there are fewer stimuli for cells in the medium to develop towards vascularisation and angiogenesis [45][46][47].
Neuronal/glial and vascular markers are not mutually exclusive. For example, NPC express the smooth muscle marker SMA in mice which control their migration [48]. It was shown that OEC also express SMA and that the proportion of OEC simultaneously expressing SMA and GFAP grows with time of cultivation on the cell culture of model animals: rats, mice, hamsters, and monkeys [49,50]. Similar cells with glial, vascular, and partially neuronal markers are present with the production of spheres from other sensory organs in animals: mice hair follicle [51], porcine retina [52], and mice tongue fungiform papilla [53]. The multipotent state of those cells seems to be essential for sensory organ regeneration. The property is exploited in regenerative medicine. These kinds of cells were successfully employed for spinal injury therapy in animal models [54,55] and even in some human cases [25,26,56]. The particular combination of glial (GFAP) and vascular (SMA) markers was also observed for gliosarcomas, a subtype of gliomas [57][58][59]. It was found that the vascular component of these gliomas bears the same somatic p53 mutations as glial cells, pointing to the monoclonal origin of the glial and vascular components [60,61]. We observed that major cluster of ontologies of downregulated genes in NDC consists of both vascular and neuronal processes, with many shared genes ( Figure 3B, Supplementary Figure S4), suggesting that in the cells of this study, the vascular markers were acquired from or even present in the original NPC and/or OEC in the olfactory biopsies.
The apparent difference between expression in NDC and Matigian et al. ONS is worth consideration. There are small but important discrepancies in protocols that could contribute to a difference in ONS and NDC. Probably the most important of these is that in the Matigian et al. protocol, primary cultures were grown in suspensions from enzymatically and mechanically dissociated biopsies. It is a suitable method for targeting the numerous neuroblasts of the olfactory epithelium of rats and mice, but it may not be very effective in humans. We used the explant method, which is supposed to be more gentle. It is also possible that both protocols work well, but stem cells from NS are unstable in their undifferentiated state, and small changes in their handling can lead to completely different outcomes.
The reason that transcriptomes of NDC resemble NPC treated with maneb could be elevated ROS production (note, that the most enriched KEGG ontology in GSEA was "oxydative phoshorilation"). Indeed, the cytotoxicity of both maneb and cadmium is linked to elevated levels of ROS [64][65][66]. NDC have increased expression levels of genes of glutathione peroxidase, GPX1 (log2FC = 1.73, p = 2.99 × 10 −14 ) and superoxide dismutase 1, SOD1 (log2FC = 0.92, p = 1.02 × 10 −10 ). As we mentioned earlier, transient ROS production is known to promote vascular growth, which is what we see in the case of NS to NDC progression. The counterpoint is that many other substances in the Kuusisto et al. screen are known to cause oxidative stress, like lead [67], permethrin [68], and isotretinoin [69]. Furthermore, Roede et al. reported that in SH-SY5Y neuroblastoma cell line, it is paraquat rather than maneb that is responsible for ROS production [70]. Anderson et al. studied the biochemistry of acute treatment of neuroblastoma cell line SK-N-AS with maneb [71,72]. They found that while energy metabolism is affected by maneb, the effect is more specific. Maneb seems to act on thiol groups of mitochondrial proteins and alter glucose metabolism by promoting gluconeogenesis at the expense of ATP production. In our NDC, we also observed multiple changes in the expression of genes involved in energetic and biochemical pathways, namely oxidative phosphorylation, nucleotide metabolism, and the citrate cycle (Supplementary Tables S4-S7). Still, the origin of the similarity of expression profiles in NDC and NPC treated with maneb is unclear and worthy of additional investigation.
There are some limitations to this study. First, the presented RNAseq results are acquired on bulk RNA and are not aware of the heterogeneity of NDC and NS. The immunofluorescent pictures do confirm that most of the cells bear the same markers but to a different extent, like in the case of SMA. Different methods, such as scRNA-seq, are needed to study the contribution of different types of cells in the observed effects. Second, we used dimensionality reduction of transcriptome data to contrast obtained expression profiles against publicly available data. While data from similar cell types tend to cluster together, there is no way to correct for differences in technicalities used to produce data in different laboratories and interpretation of such results should be conducted with caution. On the other hand, our second method of expression profiles comparison with rank correlations could be too conservative and misleading when correlations are too small. Most importantly, the results concerning the similarity of NDC with NPC treated with maneb are based on a single observation, which should be independently replicated to draw a definitive conclusion. We plan to address this in future research.
In conclusion, we observed that gene expression profiles of NDC with present cultivation routine tend to alter their energetic metabolism, which reflects increased vascular and stromal processes associated with a partial loss of their neuronal identity. For now, it seems more practical to use original NS for studying ex vivo neuronal cells. Still, it is possible that by growing NDC in the presence of ROS, scavengers and/or different sources of energy cell cultures with more pronounced neuronal signatures can be obtained to render the material more useful in mass collection for psychiatric research.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/cells10071690/s1. Figure S1. Heterogeneity of SMA in NDC; Figure S2. Matrix of Spearman's correlations between analysed expression profiles; Figure S3: Map of mutually overlapping gene sets, enriched in GSEA analysis of gene ontologies of cellular components; Figure S4: Map of mutually overlapping gene sets, enriched in GSEA analysis of gene ontologies of biological processes; Figure S5: Pathview visualisation of fold change (NDC vs. NS) of genes involved in Parkinson's disease; Table S1: List of data used on dimensionality reduction plots; Table S2: Deseq2 results of differential expression of NDC vs. NS; Table S3: GSEA results for cell component gene ontologies; Table S4: GSEA results for biological process ontologies;   Informed Consent Statement: All study participants signed an informed consent form to participate in the research. Data Availability Statement: Raw RNA-seq sequencing reads generated in this study were deposited to the SRA database, accession number PRJNA732358 (https://www.ncbi.nlm.nih.gov/sra/PRJNA7 32358, accessed on 3 July 2021).

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analysis or interpretation of data, in the writing of the manuscript or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: