Differentiation of ciliated human midbrain-derived LUHMES neurons

Primary cilia are non-motile, hair-like cell surface protrusions present in a single copy on polarized cell types. The cilium consists of a microtubular core structure, the axoneme, surrounded by a specialized membrane. The axoneme is a continuous outgrowth from the centriole-derived basal body, which anchors the ciliary shaft within the cytoplasm of the cell. A transition zone at the base of the axoneme forms a ciliary gate where cargo is sorted for import and export to and from the ciliary compartment. The maintenance of cilia depends on selective transport of protein components mediated by the intraflagellar transport (IFT) machinery (Ishikawa and Marshall, 2011).

Primary cilia are not restricted to terminally differentiated cell types, but can also be found on proliferating cells, such as neural progenitors. The centriole functions both as basal body and axonemal template in ciliogenesis and as microtubule-organizing center (MTOC) of the centrosome. Cilia are thus tightly associated with the dynamics of the cell cycle (Malicki and Johnson, 2017). Ciliogenesis is linked to the G0 and G1 phases of the cell cycle. Prior to cell division the cilium is resorbed to free a pair of centrioles from the ciliary base. During S and G2 phases, centrioles are duplicated and mature into functional centrosomes, which organize the formation of the bipolar spindle for sister chromatid separation in M phase (Pletz et al., 2013; Ford et al., 2018; Joukov and De Nicolo, 2019).

Primary cilia are antennae that receive and transduce signals from the immediate environment of a cell. The primary cilium as a signaling hub integrates several aspects of the sonic hedgehog (SHH) pathway, but is also involved in other signal transduction pathways, such as WNT, Notch, platelet-derived growth factor (PDGF) or multiple G-protein-coupled receptor pathways (Goetz and Anderson, 2010; Lauter et al., 2018; Liu et al., 2018; Wheway et al., 2018). In addition to being a receiver, the primary cilium has also gained attention as a source of signals to its cellular environment (Garcia et al., 2018; Nachury and Mick, 2019).

Cilia are highly conserved within the eukaryotic kingdom (Piasecki et al., 2010), and in humans they can be found on many different cell types (Reiter and Leroux, 2017). Primary cilia are a common feature of polarized cells, and they are abundant in the brain, both during development and at mature stages, extending from progenitor cells as well as from terminally differentiated neurons and astrocytes (Bishop et al., 2007; Arellano et al., 2012; Sarkisian and Guadiana, 2015). Disruption in the structure, maintenance or function of primary cilia gives rise to a group of highly divergent disorders collectively termed ciliopathies.

Owing to the widespread distribution of cilia, ciliopathies display pleiotropic features, including brain phenotypes (Mitchison and Valente, 2017; Reiter and Leroux, 2017). A number of neurodevelopmental disruptions that result from SHH-dependent patterning defects and certain aspects of neuropsychiatric dysfunctions, as observed in Joubert syndrome and autism-spectrum disorders, have been connected to ciliary aberrations (Massinen et al., 2011; Lauter et al., 2018; Youn and Han, 2018; Thomas et al., 2019). The associations between changes in ciliary structure and function and the manifestation of a particular (brain) disorder vary in solidity and strength. The causality of the link is also a matter of debate. Thus, investigations of the underlying molecular and cellular mechanisms, possibly cilia-associated, are urgently needed.

Currently, multiple mammalian cell culture models exist for the investigation of cilia biology, including ciliary aspects of sensory perception and signal transduction. Despite the undisputed value of existing models, no easily amenable human neuronal cell model exists that allows study of ciliogenesis, cilia structure and function in (neural) progenitor cells during all stages of the cell cycle, as well as during development and neuronal differentiation. To fill this gap, we set up the Lund human mesencephalic (LUHMES) cell line as a model for studying primary cilia in human neurons. Ciliated LUHMES cells are easy to maintain and can be differentiated into fully functional neurons within about one week.

The LUHMES cell line is an easy-to-maintain human neuronal cell model and a candidate for studying cilia biology (Scholz et al., 2011). We used immunocytochemistry to identify neuronal cilia during both cell proliferation and after various days past the induction of differentiation into functional neurons. ADP-ribosylation factor-like GTPase 13B (ARL13B) is a widely used marker for cilia that is specifically and strongly enriched in the cilium (Caspary et al., 2007, 2016; Ferent et al., 2019). Immunofluorescence stainings revealed that ARL13B is expressed and specifically localized in LUHMES cells at the stage of cell proliferation (0 d of differentiation, d0) as well as during differentiation into neurons (1 d to 6 d of differentiation, d1 to d6) (Fig. 1A), suggesting that LUHMES cells and neurons are ciliated. To confirm ciliary localization we looked at double stainings of ARL13B together with other known ciliary markers. Pericentrin (PCNT) localizes to centrioles and to centriole-derived ciliary basal bodies (Pitaval et al., 2010; Ishikawa and Marshall, 2011; Pletz et al., 2013). ARL13B and PCNT colocalized at the ciliary shaft and basal body (Fig. 1B), corroborating the presence of cilia. Proper structure and function of cilia depend on IFT (Ishikawa and Marshall, 2011). As expected, IFT88 protein colocalized together with ARL13B along the ciliary shaft (Fig. 1C), further supporting the presence of neuronal cilia. Other ciliary markers, such as detyrosinated tubulin and RAB8A, also showed expression and specific localization to cilia (data not shown). Using immunofluorescence stainings we could thus demonstrate the presence of cilia on LUHMES cells and neurons.

To unambiguously confirm the presence and the nature of cilia on LUHMES cells and neurons, we performed transmission electron microscopy (TEM). Using TEM we could clearly verify the presence of extended cilia (Fig. 2). We detected all the defining parts of cilia in cross-sections as well as longitudinal sections through the structure. Basal bodies displayed the characteristic configuration of nine microtubule triplets (Fig. 2A,B), whereas more distal parts, including the ciliary axoneme, consisted of the typical microtubule doublet formation (Fig. 2C) (Pletz et al., 2013). Basal body appendages and transition zone structures (Fig. 2B–E) as well as a pronounced ciliary pocket were clearly distinguishable (Fig. 2C–G). Of note, the apparent lack of dynein arms between the microtubule doublets (Fig. 2C) and the arrangement of these microtubule doublets in a 9+0 configuration confirmed that LUHMES cilia are non-motile, primary cilia.

Primary cilia are involved in the transduction of a variety of signaling pathways, including, in vertebrates, the coordination of the canonical SHH signaling pathway (Anvarian et al., 2019). In vertebrates, the cilium is required for activation of the canonical SHH signaling cascade (May et al., 2020, preprint), whereas the cilium is not essential for SHH signaling in select species such as the invertebrate Drosophila (Chen et al., 2009), where only certain types of sensory neurons are ciliated (Lauter et al., 2018). Here, we examined the output of canonical SHH signaling in LUHMES, our ciliated human neuronal cell model. We treated differentiating (d2) LUHMES neurons for 24 h with either an SHH pathway activator, smoothened agonist (SAG), or an antagonist, cyclopamine, that modulate the canonical SHH pathway through the GLI transcription factors (Stanton and Peng, 2010). After treatment (on d3), we evaluated by reverse transcription qPCR (RT-qPCR) the gene expression of four downstream SHH pathway target genes: GLI1, HHIP, PTCH1 and PTCH2. Samples treated with the SHH pathway activator SAG (500 nM) showed a significant upregulation of target gene expression, as compared to expression in the samples treated with the antagonist cyclopamine (10 µM) or the control samples (0.25% DMSO) (Fig. 3). We have obtained very similar results by treating LUHMES cells and neurons at different time points with either SAG (d0) or N-terminus recombinant human SHH protein (d2) (data not shown). These observations confirm that human LUHMES cells and neurons transmit molecular signals, exemplified by the canonical SHH signal transduction pathway.

To establish the LUHMES cell model as a useful tool for studying human, neuronal cilia we adopted a simplified protocol for differentiating LUHMES cells into functional neurons: within about one week (d0–d6) proliferating LUHMES cells can be differentiated into a state at which different neuronal differentiation and maturation markers are readily detectable and even solid electrophysiological assays are possible (Scholz et al., 2011). To better mimic in situ neuronal development, we used unsynchronized culture conditions for neuronal differentiation.

By following representative groups of several hundreds of cells through image capture and by manual counting, we determined that upon induction of neuronal differentiation at d0, LUHMES cells divide only once from d0 to d1 and then become postmitotic (Fig. 4A–D). The complete arrest of all cell proliferation after about one day post-induction of neuronal differentiation confirmed that the applied culturing protocol pushed all the cells into full differentiation mode. We did not observe any significant cell death (tested up to d6), and cell viability was consistently >95% (Fig. 4E). These features suggest that the adopted culture conditions (Scholz et al., 2011) were in the optimal range.

To confirm that LUHMES cells develop cilia as they become postmitotic and differentiate, we followed the expression of both ciliary (ARL13B) and neuronal differentiation markers. From loosely aggregated colonies of proliferating LUHMES cells at d0 (Fig. 5A), differentiating neurons dispersed and progressively grew neurite extensions to eventually form a well-developed network of neurites, easily discernible on d5 or d6 (Fig. 5). At any given time point during this developmental timecourse, LUHMES cells and neurons display primary cilia (d0 to d6) (Figs 1, 4 and 5). The percentage of ciliated cells and neurons was not uniform between stages. Proliferating LUHMES cells progressing through the cell cycle in our unsynchronized culture conditions displayed a much lower percentage of ciliation (∼20–25% at d0 and 30–35% at d1, after exit from the cell cycle) as compared to after the induction of (postmitotic) neuronal differentiation. Ciliation of these neuronal populations quickly rose to well above 50% at d2 to d5 of differentiation. Thus, ciliation slightly preceded and then coincided well with the differentiating neuron population attaining a complex network of neurites during the same time period (d2 to d5).

LUHMES differentiation into functional neurons has been described elsewhere, and our work faithfully replicated the findings, including the overall time frame (Scholz et al., 2011; Shah et al., 2016; Harischandra et al., 2020). Accordingly, we observed that, in addition to the ciliary marker ARL13B (Fig. 5A), differentiating LUHMES neurons expressed the general neurite marker β3-tubulin (TUBB3) and the synaptic marker postsynaptic density protein 95 (PSD95, also known as DLG4) (Fig. 5B), the neuronal transport marker kinesin 17 (KIF17) (Fig. 5C) as well as tyrosine hydroxylase (TH), the terminal differentiation marker for dopaminergic neurons (Fig. 5D) (Wulle and Schnitzer, 1989; Lee et al., 1990; Hunt et al., 1996; Franker et al., 2016). The expression of various markers for neuronal differentiation, maturation and function suggested that our culture conditions successfully induced neuronal differentiation, in accordance with previous findings.

Proliferating, progenitor-type LUHMES cells (d0) differ strongly from differentiating LUHMES neurons (d1 to d6) (Scholz et al., 2011) with regard to anatomy and morphology, physiology and function, and a number of molecular aspects. One main difference is the absence of developing or differentiated neurites (dendrites and axons) in proliferating LUHMES cells. Despite these differences, the position and localization of the primary cilium remained constant during both proliferation and differentiation (Fig. S1). Using regular monolayer-type culture conditions and immunofluorescence staining (ARL13B), we found by visual inspection and distance measurements (Fig. S1A,B) that cilia always located at the cell body of proliferating cells and of neurons with differentiating neurites. Using microscopy image stacks and 3D rendering along the full z-axis encompassing entire LUHMES cells or neurons, we found the cilium to always locate in the middle sections of the stacks (Fig. S1C; Movie 1). Primary cilia in LUHMES are, thus, ideally located for sensory tasks concerning neighboring cells or the culture medium. The length of LUHMES cilia can reach up to 4–5 µm (data not shown). On average, LUHMES cilia are slightly shorter than cilia from, for example, the human retinal pigment epithelial 1 (RPE1) cell line (Pitaval et al., 2010).

Cilia localization at the cell body being ideal for sensory tasks, and the rise in ciliation during neuronal differentiation, prompted us to search for functional roles of cilia in differentiating LUHMES neurons. Our unsynchronized culture conditions allowed us to compare neurons with detectable cilia versus neurons without detectable cilia during their timecourse of differentiation. We observed that the presence of cilia strongly affected axonal connectivity patterns. LUHMES neurons with cilia display increased and earlier axon branching as compared to their counterparts without detectable cilia (Fig. S2). Our findings in human neurons therefore strongly support an important role of cilia during neuronal differentiation. Furthermore, our findings are corroborated by work in the mouse, where it was found that ciliary signaling promotes axon development (Guo et al., 2019) and coordinates neuronal migration (Matsumoto et al., 2019; Stoufflet et al., 2020).

To monitor gene expression changes during the differentiation process, from proliferating cells to differentiated neurons, we applied the single-cell tagged reverse transcription RNA-sequencing (STRT RNA-Seq) method, which captures the 5′ ends of poly-A transcripts (transcript far 5′ ends, TFEs) (Islam et al., 2011, 2014; Töhönen et al., 2015; Krjutškov et al., 2016). We sampled total RNA from d0 to d6 in two independent timecourse experiments with 5–6 replicates per time point; in total 81 samples were sequenced, split across three sample libraries (Table S1; Figs S3–S6). After having removed eight samples that failed quality checks and three technical replicates, a subtotal of 70 samples were subjected to principal component analyses (PCA) based on the expression levels of TFEs. We corrected biases between the three sample libraries using ComBat (Johnson et al., 2007) (Fig. S7A,B). In the resulting PCA plot, we confirmed that all the samples clustered based on the time points during the differentiation timecourse (Fig. 6A). Notably, d0 samples (proliferating LUHMES cells) cluster separately from all other time points, indicating that gene expression changed drastically between d0 and d1, reflecting the induction of neuronal differentiation. Samples from d5 and d6 cluster very closely, indicating that the gene expression profile of neuronal differentiation became stable around d5, reflecting the emergence of fully differentiated neurons. Correlations between samples also supported this notion, showing low correlation between d0 and samples from other time points, but high correlation around and after d5 (Fig. 6B).

Following previous reports about the LUHMES cell model (Delp et al., 2018; Pierce et al., 2018), we performed gene-based differential expression analysis between d0 and d6. We found 5800 genes with significantly upregulated expression, whereas expression of 3026 genes was downregulated (Fig. S7C). We performed gene ontology (GO) term enrichment analysis on these gene lists (Tables S2, S3), and confirmed that terms related to nervous system development, vesicular transport and axonogenesis were significantly enriched in the upregulated genes, appropriately reflecting that LUHMES cells differentiated into functional neurons. Interestingly, autophagy-, proteasome- and ubiquitylation-related terms were significantly enriched in the upregulated genes. It is known that these processes play essential roles in development and cell differentiation by recycling cellular components to meet biochemical demands (Mizushima and Levine, 2010). Autophagy- and proteasome-related features have also been reported for neurogenesis (Morgado et al., 2015), for neuronal differentiation (reviewed in Casares-Crespo et al., 2018), and also in connection to ciliogenesis (reviewed in Malicki and Johnson, 2017). Importantly, the ‘cilium assembly’ term itself was also significantly enriched in the upregulated genes (Table S2). In contrast, the downregulated genes were significantly enriched with terms broadly related to ribosomal functions, in agreement with previous reports (Delp et al., 2018; Pierce et al., 2018) (Table S3).

In a complementary approach, we classified 5596 genes with highly fluctuating expression into four clusters based on their expression pattern changes during neuronal differentiation, and identified classes of genes that were overrepresented using the core contributing genes in each cluster (shown as red lines in Fig. 7). Genes in cluster 1 were highly expressed until d1 and monotonically downregulated after d2. Cluster 1 is highly enriched with ‘ribosomal function’ genes (Fig. 7; Table S4), indicating that translation is intensely activated on d1, when the cells begin to convert from progenitors into neurons. Genes in cluster 2 showed the highest expression on d0 and their expression was then drastically downregulated, starting on d1. Cluster 2 is highly enriched with replication and cell cycle-related genes (Fig. 7; Table S5; compare to Fig. 4A–D). Expression of genes in cluster 3 was drastically upregulated between d0 and d1. Cluster 3 is highly enriched with genes associated with transcription, transport and catabolic processes (Fig. 7; Table S6), indicating that cells very dynamically changed their status during the early parts of neuronal differentiation. Interestingly, cluster 3 also contains sets of significantly enriched ciliary genes (Table S6). Finally, we found that genes in cluster 4, most strongly upregulated from d2 to d5, were highly enriched for the terms ‘axonogenesis’ (and similar) as well as ‘cilium assembly’ (Fig. 7; Table S7). In summary, we note that highly upregulated genes during LUHMES neuronal differentiation (clusters 3 and 4) are significantly enriched with ciliary genes, including 100 ciliary genes already described in the SYSCILIA gold standard list (van Dam et al., 2013).

To assess the similarity between LUHMES neuronal differentiation in culture (d0 to d6) with available in vivo data from human (neural) progenitor cells and differentiated neurons, we compared our gene expression differentiation timecourse results with single-cell RNA-seq data from human fetal midbrain (La Manno et al., 2016) (Fig. 6C), using the tool SingleR, which assigns the cellular identity of single-cell transcriptome data (Aran et al., 2019). We found that d0 samples (proliferating LUHMES cells) showed high similarity with progenitor-type cells from human midbrain, a similarity that was gradually lost during the neuronal differentiation timecourse. In contrast, LUHMES cells differentiating into neurons (d1 to d6) progressively acquired high similarity with differentiated human midbrain neurons including dopaminergic, serotonergic and GABAergic neurons. Furthermore, LUHMES cells and neurons showed strong dissimilarity with other types of cells such as pericytes and endothelial cells (Fig. 6C). These quantitative observations confirm that LUHMES neuronal differentiation in culture closely mimics the in vivo situation in human midbrain.

GO term enrichment analyses (Tables S2–S7) found the term ‘cilium assembly’ (and similar) significantly enriched in the upregulated genes and in defined subsets of genes with highly fluctuating expression (clusters 3 and 4; Fig. 7), molecularly supporting the ciliogenesis observed during the LUHMES neuronal differentiation timecourse (d0 to d6). We confirmed neuronal differentiation and ciliary gene STRT RNA-seq data using an independent method: RT-qPCR normalized against the housekeeping gene GAPDH. Pearson correlation coefficient plots demonstrated effectively identical differential gene expression patterns for a replication group of 17 genes (Figs S3–S6; Table S8) (Tammimies et al., 2016; Reiter and Leroux, 2017; Lauter et al., 2018; Youn and Han, 2018): four neuronal differentiation markers (DCX, MAP2, TUBA1A and TUBB3) (Fig. S3); two cilia and ciliopathy genes (CCDC28B and IFT20) (Fig. S4A,B); two candidate ciliogenic transcription factor genes (RFX2 and RFX5) (Fig. S4C,D); six cilia and ciliopathy genes with brain or neuron phenotypes (BBS1, BBS2, CC2D2A, IFT81, PDE6D and TCTN2) (Fig. S5); a ciliary dyslexia candidate gene (DYX1C1, also known as DNAAF4) (Fig. S6A); and two neuronal Parkinson's disease-associated genes (MAPT and SNCA) (Fig. S6B,C). Ciliogenesis was found as a prominent aspect during LUHMES neuronal differentiation. Similar ‘ciliary’ observations were made for the differentiation from human neuroepithelial stem (NES) cells and induced pluripotent stem cells (iPSCs) into neurons (Bieder et al., 2020).

Using the SYSCILIA gold standard list of ciliary genes (van Dam et al., 2013), we conducted gene-set enrichment analysis (GSEA) by comparing each day of LUHMES neuronal differentiation with the undifferentiated state on d0. Enrichment scores (ES) were highest when we compared d0 with d5 (normalized ES=2.57, P<0.001) (Fig. 8A; Fig. S8), demonstrating that ciliary genes were most significantly enriched in the upregulated genes on d5. This outcome is in accordance with the gene expression patterns of highly fluctuating genes in cluster 4, which showed strong upregulation prominently including d5 of differentiation (Fig. 7). These results prompted us to explore which regulatory factors might control the expression of ciliary genes during LUHMES neuronal differentiation.

We applied motif activity response analysis (MARA) (FANTOM Consortium and Riken Omics Science Center, 2009; Balwierz et al., 2014) to find candidate transcription factors responsible for gene expression changes during LUHMES neuronal differentiation. MARA computationally predicts key transcription factors that drive gene expression changes by integrating gene expression data with transcription factor binding motif statistics in gene promoters. MARA output consists of a binding motif activity Z-value, which quantifies the contribution of a given binding motif to the gene expression variations of its target genes across individual time points of a timecourse. In this analysis ‘MXI1_MYC_MYCN’, ‘SOX3_SOX2’, ‘RFX3_RFX2’, and ‘E2F7_E2F1’ motifs showed exceptionally high Z-values in all the three libraries (Fig. 8B; Table S9). Conversely, the activity Z-values of ‘MXI1_MYC_MYCN’ and ‘E2F7_E2F1’ motifs decreased drastically on d1 (Fig. 8C). This observation is explained by the addition of tetracycline to the neuronal differentiation growth medium, which shuts down a v-myc expression cassette (see Materials and Methods). The E2F transcription factor family is known to play crucial roles in the control of the cell cycle, which is in accordance with the downregulation of cell cycle-related genes on d1 (cluster 2 in Fig. 7). In contrast, activity Z-values of ‘SOX3_SOX2’ and ‘RFX3_RFX2’ motifs increased during LUHMES neuronal differentiation (Fig. 8D). SOX2 is known to be required for the maintenance of neural stem cells (Pevny and Nicolis, 2010), and is downregulated during neuronal differentiation (Zhang and Cui, 2014). Here we also found that SOX2 expression was downregulated during LUHMES neuronal differentiation (Fig. S6D), indicating that SOX2 might act as a repressor for neuronal genes (Liu et al., 2014). The activity Z-value of ‘RFX3_RFX2’ motifs strongly increased during LUHMES neuronal differentiation until d5 (Fig. 8D), a signature very similar to the upregulated expression patterns of many ciliary genes (Figs S4A,B, S5 and S6A). RFX transcription factors are known to be key ciliogenic regulators in many different organisms (Senti and Swoboda, 2008; Piasecki et al., 2010; Choksi et al., 2014; Sugiaman-Trapman et al., 2018). As an example, the ciliary gene BBS1, commonly involved in a ciliopathy with brain phenotypes – Bardet–Biedl syndrome (Mykytyn et al., 2002), was also upregulated until d5 (Fig. S5A). Notably, the BBS1 gene promoter has a binding site for RFX2 (Yan et al., 2013) (Fig. 8E). For comparison, RFX3 showed a similarly upregulated expression pattern until d5, whereas RFX2 was most highly expressed on d1 but downregulated after d2 (Fig. S6D). These results suggest that during LUHMES neuronal differentiation, RFX2 might be responsible for the induction of ciliogenesis early on, whereas RFX3 is responsible for ciliary maintenance later on. We are aware of possible redundancies though, given that the human RFX transcription factor family consists of eight members (Sugiaman-Trapman et al., 2018) and that RFX binding motifs (X-boxes) from the different family members are very similar to each other (Fig. 8D).

We independently confirmed the above findings using motif enrichment analysis of the differentially expressed TFEs during LUHMES neuronal differentiation. E2F and MYC binding motifs were significantly enriched in downregulated TFEs between d0 and d1. RFX binding motifs were significantly enriched in upregulated TFEs between d3 and d4, as well as between d4 and d5 (Table S10). Taken together, these observations strongly suggest that RFX transcription factors play a crucial role in the regulation of ciliary genes during LUHMES neuronal differentiation.

We have thoroughly characterized the LUHMES ciliated human neuronal cell line. The LUHMES neuronal cell model is well established, and a large body of published work demonstrates its experimental and technical usefulness for basic neurobiology research. To an array of cell biological, molecular and biochemical studies, including high-throughput approaches, we add a series of experimental investigations of cilia. Below we highlight important advantages and possible limitations with a focus on cilia biology.

We have appropriated and successfully optimized standardized culture conditions (Scholz et al., 2011); work with LUHMES cells and neurons is technically straightforward, and the culture conditions used do not cause any relevant losses by cell death. LUHMES cells and their cilia can be experimented with under proliferative conditions, by maintaining progenitor-type cells going through the cell cycle. Alternatively, we have differentiated LUHMES cells into functional neurons with cilia within about one week by following published protocols (Lotharius et al., 2005; Scholz et al., 2011), whereby differentiated LUHMES neurons can be kept in culture for up to two weeks. Immunostainings with functional neuronal differentiation markers have been performed up to d9 (Shah et al., 2016) or d13 in our hands. Neuronal physiological parameters have been assessed using both metabolite uptake assays and electrophysiological recordings (Scholz et al., 2011). A particular strength of the LUHMES cell model is the fact that one can completely separate cell proliferation from neuronal differentiation within one day of changing the culture medium. Such strict experimental separation greatly facilitates comparative studies of ciliary structure and function (e.g. signal transduction) between two different biological situations, while keeping an identical genetic background. The LUHMES cell model is amenable to genetic modification. Concordant with Shah et al. (2016), lipid-based transfection methods were inefficient in our hands (<1%), whereas nucleofection-based methods with improved transfection efficiency (up to around 30%) are more practical, for example for CRISPR/Cas9-based genome engineering approaches (Shah et al., 2016; and our own unpublished work). One can thus examine, within a matter of weeks, the biological consequences of human brain disorder candidate mutations in a homologous system (a human neuronal cell line), although background genetic variation present in human patient material is obviously not attainable in the LUHMES cell model.

Karyotyping of the LUHMES cell model revealed a normal, female diploid set of chromosomes (Shah et al., 2016). Even though LUHMES cells and neurons are maintained in vitro in cell culture, our immunocytochemistry and detailed transcriptomics timecourse work demonstrated that both stages, proliferation and differentiation, showed exceptionally high similarities to the corresponding in vivo cell types, progenitors and differentiating neurons from human midbrain (La Manno et al., 2016), the LUHMES tissue of origin (Lotharius et al., 2002, 2005). Slight adjustments in LUHMES cell culture growth medium even allow for driving differentiation toward a specific neuron class, such as dopaminergic neurons (Scholz et al., 2011; Luk et al., 2015). Damage to and loss of dopaminergic neurons is a hallmark of Parkinson's disease patients (reviewed in Hegarty et al., 2013). Thus, the LUHMES cell model has commonly been used in Parkinson's disease studies and in neurotoxicity assays (Schildknecht et al., 2009; Stępkowski et al., 2015; Pierce et al., 2018; Ganjam et al., 2019; Kranaster et al., 2020). However, whether there is any relevant connection between cilia malfunction and Parkinson's disease is entirely unexplored. Our in-depth transcriptomics timecourse analyses of human neuronal differentiation complement and greatly extend similar transcriptomics work in the mouse (Guo et al., 2019). Using the LUHMES cell model, we provide step-by-step descriptions of the molecular flux from progenitor cell proliferation, to cell polarization and then neuronal differentiation, whereby we include a first description of the ciliary contributions during these neuro-developmental events. Importantly, we provide evidence, for the first time in human neurons, that during differentiation the presence of cilia exerts a strong impact on axon branching and connectivity patterns.

The length of LUHMES cilia typically reaches up to 4–5 µm. This may pose limitations on (fluorescence-based) live imaging of neuronal ciliary transport (IFT) or signal transduction, or render difficult the analysis of subtle ciliary structural defects. In all experiments, we have grown LUHMES cells and neurons in vitro in two-dimensional monolayer cultures. Although this approach is advantageous for microscopy, it does not entirely reflect the complex three-dimensional, in vivo situation in the human brain. Furthermore, the LUHMES cell model consists of a single cell type, either as progenitor cells or differentiating neurons. Again, this does not reflect the in vivo situation in the (developing) human brain, where a mixture of different cell types is present at any given time (La Manno et al., 2016). Although brain organoid approaches (Lehmann et al., 2019) would remedy this limitation, these would make impossible the detailed timecourse studies presented here. LUHMES neurons have successfully been co-cultured with astrocytes (Efremova et al., 2015) and in 3D cultures (Smirnova et al., 2016). However, brain organoids in connection to LUHMES neurons have not been reported yet.

In conclusion, we have successfully established a human neuronal cell model, LUHMES, where ciliogenesis, cilia structure and cilia function can easily be investigated, both in the context of neural progenitor cell type proliferation and the entire timecourse of differentiation into a mature, functional neuron. In particular, the speed and one-week time frame of LUHMES cell differentiation, in combination with an available neuronal and ciliary experimental toolbox, make for an ideal test bed to investigate (genetic) defects associated with neuronal cilia. LUHMES cells will thus be of great value for better understanding specific brain phenotypes of human ciliopathies including Joubert syndrome, certain neurodevelopmental disruptions as observed in dyslexia and aspects of neuropsychiatric dysfunctions as in autism-spectrum disorders, for all of which strong connections to ciliary dysfunction have been uncovered (Reiter and Leroux, 2017; Lauter et al., 2018; Thomas et al., 2019).

Lund human mesencephalic (LUHMES) cells constitute a human midbrain-derived neuronal cell model. The tissue source, the ventral midbrain A9 cluster of an 8-week-old female fetus, contains many dopaminergic neurons highly relevant for the development of the substantia nigra pars compacta (SNc). LUHMES cells were originally subcloned from MESC2.10 cells. Both the MESC2.10 and LUHMES cell lines have been developed at the Wallenberg Neuroscience Center at Lund University in Sweden (Lotharius et al., 2002, 2005). We have obtained LUHMES cells from the ATCC (https://www.lgcstandards-atcc.org; CRL-2927). Prior to all experimentation, LUHMES cells were collected, and the absence of mycoplasma contamination was confirmed. We used the same culture conditions as previously described (Scholz et al., 2011) with minor modifications. LUHMES cells were grown at 37°C in a humidified incubator (5% CO₂) and maintained in vessels (typically T75 flasks) pre-coated with poly-L-ornithine hydrobromide (Sigma P3655; 50 µg/ml) and human fibronectin (Sigma F1056; 1 µg/ml). The growth medium was DMEM/F-12 Ham (Sigma D6421-6) supplemented with L-glutamine (Sigma G7513; 2.5 mM), N-2 Supplement (Gibco 17502-048; 1×) and human basic Fibroblast Growth Factor (bFGF) (Thermo Fisher Scientific 68-8785-82; 20 ng/ml). LUHMES cells were grown to 80% confluency and passaged 1:6 using diluted trypsin-EDTA (0.25%; Thermo Fisher Scientific 25200056) with phosphate-buffered saline (PBS) at 1:1 volume. For all experiments except the routine maintenance of cells, cell culture vessels (typically 6-well plates) were pre-coated with higher concentrations of ornithine (100 µg/ml) and fibronectin (10 µg/ml) to improve cell attachment to the substrate. In the growth medium for LUHMES cell proliferation (d0), 40 ng/ml bFGF was added instead of 20 ng/ml, to accelerate cell propagation. For differentiating LUHMES cells into postmitotic neurons (d0 to d6), cells were pre-seeded at lower densities to reach a final cell density of 1.5×10⁵ cells/cm² during the differentiation process. 24 h after pre-seeding, on d0 of LUHMES differentiation into neurons, bFGF in the growth medium was replaced by tetracycline hydrochloride (Sigma T7660; 1 µg/ml), which shuts down a v-myc expression cassette used for cell line immortalization, as was previously described for −/− growth medium conditions (no addition of cAMP and GDNF; Scholz et al., 2011).

We point out that in the growth medium for cell proliferation (d0), LUHMES cells are proliferating in a progenitor-type cell state. In accordance with work from others (Lotharius et al., 2005; Scholz et al., 2011; Pierce et al., 2018) this statement is strongly supported by our STRT RNA-seq transcriptomics timecourse data (d0 to d6) (Fig. S6D), which demonstrate effectively absent or very low expression of several different classical neural stem cell markers. Of note, PAX6 is upregulated, consistent with its role in mediating axon outgrowth in postmitotic neurons (Sebastián-Serrano et al., 2012). SOX9, a marker for neurogenesis, is upregulated on d1, likely in a suppressive role during the transition from cell proliferation to neuronal differentiation (Vong et al., 2015) (Fig. S6D).

LUHMES cells and neurons were grown on coverslips (number 1.5 round-shaped borosilicate glass; VWR 631-0150) pretreated with concentrated hydrogen chloride (37%), washed with dH₂O and absolute ethanol, and then stored in 70% ethanol. Coverslips were placed in each well before the vessels were coated with ornithine and fibronectin. On the day of immunostaining, the wells with coverslips were fixed with 2% formaldehyde solution (Thermo Fisher Scientific, FB002) in PBS for 15 min at room temperature. After fixation, cells and neurons were permeabilized with 0.2% TritonX-100 in PBS for 12 min. The coverslips were then transferred to a humid chamber and incubated with the blocking buffer (2% bovine serum albumin and 0.1% Tween20 in PBS) for 30 min, followed by an overnight primary antibody incubation at 4°C and a 60 min secondary antibody incubation at room temperature. Cell nuclei were stained with DAPI solution (Hoechst, Invitrogen H3570; 50 µg/ml in PBS) for a few seconds and the plasma membrane with Cell Mask (Invitrogen HCS Deep Red Stain H32721; 2 ng/µl in PBS) for 30 min at room temperature, before the coverslips were mounted on a microscope slide with Mowiol mounting medium.

The following primary antibodies were used. Our main ciliary marker was ARL13B (rabbit; Proteintech 1711-1-AP; 1:5000). Other ciliary markers included PCNT (rabbit; Human Protein Atlas HPA019887; 1:1000), ACTUB (mouse; Sigma T7451; 1:5000) and IFT88 (rabbit; ProteinTech 13967-1-AP; 1:500). Neuronal differentiation markers were TUBB3 (Tuj1) (mouse; Covance, MMS-435P; 1:3000), PSD95 (rabbit; Cell Signaling Technology 2507; 1:100), KIF17 (rabbit; Abcam ab11261; 1:2000) and TH (rabbit; R&D Systems 779427; 1:50). Phalloidin was used to stain F-actin (Alexa Fluor 647-conjugated; Thermo Fisher Scientific A22287; 1×). Secondary antibodies (Alexa Fluor dye conjugated) were obtained from Thermo Fisher Scientific and used at a 1:600 dilution. All primary and secondary antibodies were diluted in blocking buffer during the incubation.

All images of samples used for immunocytochemistry were acquired at the Live Cell Imaging facility of the Karolinska Institute, Department of Biosciences and Nutrition (https://ki.se/en/bionut/welcome-to-the-lci-facility). A Nikon Ti-E inverted point scanning confocal microscope A1+ (Nikon Instruments) was used with a Plan Apo λ 20× NA 0.75 dry objective or a Plan Apo λ 60× NA 1.4 oil objective for higher magnifications. Image analysis was carried out using ImageJ software (Schneider et al., 2012).

Cell densities of cultures in six-well-plates were used to determine changes in population growth. Several different defined positions within a well were imaged on a bright-field microscope (Carl Zeiss, Axiovert 200M) on consecutive days of neuronal differentiation. The number of cells and neurons counted in these defined areas was extrapolated to the entire well area to derive the cell density of one well at a given time point during the differentiation timecourse. Three independent quantifications were conducted.

Bright-field pictures of cells and neurons were taken on a microscope equipped with epifluorescence (Carl Zeiss, Axiovert 200M) and dead cells were counted manually. Using an independent method, the same cell and neuron samples were stained by adding 2 drops/ml of NucGreen Dead 488 Ready Probes reagent (Thermo Fisher Scientific, R37109), a fluorescent dye that stains cells lacking full plasma membrane integrity. After 10 min of incubation, fluorescent dead cells were imaged by microscopy and manually counted.

Mitotic nuclei were counted on images of fixed DAPI-stained LUHMES cultures from d0 until d3 of neuronal differentiation. Immunofluorescence antibody detection (acetylated tubulin, ACTUB) was included as a differentiation marker. At least 150 nuclei were counted for each individual time point, whereby mitotic nuclei were visually identified.

To determine the exact position of cilia on LUHMES cells and neurons, these were stained and imaged as described above, using ARL13B as ciliary marker, DAPI for nuclei and Cell Mask to identify the cell-body boundaries. To discriminate the Cell Mask fluorescence signal between the cell body and developing neurites, a circle was drawn around the cell nucleus (Fig. S1A). This circle was then considered as the edge of the cell body for subsequent quantifications. The cilium–nucleus distance (d2, d3, d4) and the nucleus–cell-body edge distance (d0, d2, d3, d4) were measured using the ImageJ line tool, starting from the center point of the nucleus as visualized by DAPI staining. For all time points, a minimum of 50 ciliated cells or neurons were counted. For 3D rendering, snapshots of a series of 3D z-axis stacks were compiled using the 3D Viewer ImageJ plugin.

TEM work was carried out at the electron microscopy facility of the Karolinska University Hospital Huddinge (https://ki.se/en/research/the-electron-microscopy-unit-emil). LUHMES cells were differentiated into neurons for 3 d and subsequently harvested in a cell pellet. Differentiating neurons were fixed in 2.5% glutaraldehyde in 0.1 M phosphate buffer, pH 7.4, and stored at 4°C until further use. Following fixation, the cells were rinsed in 0.1 M phosphate buffer prior to post-fixation in 2% osmium tetroxide in 0.1 M phosphate buffer, pH 7.4, at 4°C for 2 h. The cells were then stepwise dehydrated in ethanol followed by acetone and embedded in LX-112 (Ladd; 21310). Ultrathin sections (between 50 and 60 nm) were cut using a Leica EM UC 7 microtome and contrasted with uranyl acetate followed by lead citrate. The sections were examined in a Tecnai 10 transmission electron microscope at 100 kV. Digital images were acquired using a Veleta CCD camera.

LUHMES cells were cultured as described above and pre-seeded in 12-well plates in proliferation growth medium for 24 h or overnight. At d0, the proliferation medium was changed to differentiation medium, and LUHMES cells were incubated for an additional 48 h (until d2 of neuronal differentiation). On d2, an SHH pathway activator, smoothened agonist (SAG; Merck 566660), or an SHH pathway antagonist, cyclopamine hydrate (Sigma C4116), was added to the wells for 24 h of treatment, prior to the extraction of total RNA on d3. SAG and cyclopamine (Stanton and Peng, 2010) were both dissolved in DMSO and were used at final concentrations of 500 nM and 10 µM, respectively. Control samples consisted of LUHMES cells and neurons cultured in differentiation medium for the same lengths of time, whereby only DMSO at 0.25% was added to the wells for 24 h of treatment. In parallel experiments, differentiating LUHMES neurons (d2) were treated with N-Terminus recombinant human SHH protein (R&D Systems, 1845-SH-025) at 400 ng/ml for 48 h, for SHH pathway activation.

Total RNA was isolated from LUHMES cells and neurons using a QIAGEN RNeasy Mini Kit and DNase Set, following the RNeasy Mini Handbook instructions. Initial RNA concentrations were measured with a Nanodrop ND-100 (Thermo Fisher Scientific). RNA quality was checked for an RNA integrity number>8 (Agilent Tech 2200 Tape Station; KI Bioinformatics and Expression Analysis core facility; https://ki.se/en/bionut/bea-core-facility). Prior to library preparation for transcriptomics analyses, RNA concentrations were measured with a Qubit 3.0 Fluorometer (Thermo Fisher Scientific), and RNA samples to be used further were diluted to 10 ng/µl (±0.5 ng/µl) in nuclease-free water.

RNA (500 ng) was reverse transcribed into cDNA using a RevertAid H Minus First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) with random hexamer primers, according to the manufacturer’s protocol. cDNA was diluted 1:10, and 2 µl of diluted cDNA was used per reaction. The mRNA expression was analyzed by RT-qPCR using a 7500 Fast Real-Time PCR instrument (Applied Biosystems) and FastStart Universal SYBR Green Master reagents (4913850001, Sigma-Aldrich). All primer sequences are listed in Table S8. Standard RT-qPCR runs consisted of a first passage at 95°C (10 min), followed by 40 cycles at 95°C (15 s) and 60°C (1 min). Melting curve analyses were performed to assess whether only single, specific PCR products had been obtained. The relative fold-increase of gene expression was calculated, using the 2^−ΔΔCt method, by subtracting the Ct values of the gene of interest from those of the housekeeping gene GAPDH (ΔCt), and then normalizing the resulting ΔCt values against the control samples at d0 (Livak and Schmittgen, 2001).

We used 20 ng of RNA to generate three 48-plex RNA-seq libraries employing a modified STRT method with unique molecular identifiers (UMIs) (Islam et al., 2011, 2014). Briefly, RNA samples were placed in a 48-well plate. A universal primer, template-switching oligonucleotides and a well-specific 6-bp barcode sequence (for sample identification) were added to each well of the plate (Krjutškov et al., 2016). The synthesized cDNAs from these samples were then pooled into one library and amplified by single-primer PCR using a universal primer sequence. STRT library sequencing was carried out with an Illumina NextSeq 500 System, High Output (75 cycles), at the Biomedicum Functional Genomics Unit (FuGU), University of Helsinki, Finland (https://www.helsinki.fi/en/infrastructures/genome-analysis/infrastructures/biomedicum-functional-genomics-unit).

The raw base-call (BCL) files were demultiplexed using Picard (v2.10.10; http://broadinstitute.github.io/picard/) ExtractIlluminaBarcodes and IlluminaBasecallsToSam to create unaligned BAM files. These BAM files were converted to FASTQ files using Picard SamToFastq, and aligned to the human reference genome hg19, human ribosomal DNA unit (GenBank: U13369) and ERCC spike-ins (SRM 2374) using the GENCODE (v28) transcript annotation by HISAT2 (v2.1.0) (Kim et al., 2019) with the option ‘--dta'. Aligned BAM files were then merged with the original unaligned BAM files to create UMI-annotated BAM files using Picard MergeBamAlignment. Subsequently, these BAM files corresponding to each sample derived from four lanes were merged using Picard MergeSamFiles. Finally, potential PCR duplicates were marked with Picard MarkDuplicates. The resulting BAM files were processed using featureCounts (v1.6.2) (Liao et al., 2014) to assign the reads to 5′ ends of genes or TFEs with options ‘-s 1 --largestOverlap --ignoreDup --primary'. For gene-based analysis, uniquely mapped reads that were within the 5′ UTR or 500 bp upstream of the NCBI RefSeq protein-coding genes, or within the first 50 bp of spike-in sequences, were counted. For TFE-based analysis, the mapped reads were assembled using StringTie (v1.3.3) (Pertea et al., 2015) with options ‘--fr -m 74 -c 10’, and those mapped reads within the first exons of the assembled transcripts were counted as previously described (Töhönen et al., 2015). FASTQ files after exclusion of duplicated reads were deposited in the ArrayExpress database at EMBL-EBI (https://www.ebi.ac.uk/arrayexpress) under the accession number E-MTAB-8908. The numbers of total and mapped reads for each sample are summarized in Table S1.

Among the 81 samples from LUHMES cells and neurons evenly distributed over three separate libraries, eight samples were excluded as outliers due to extremely high or low numbers of mapped reads. The bias in read counts between these three libraries was corrected by the ComBat function (Johnson et al., 2007) of the R (v3.6.0) package sva (v3.32.1) (Leek et al., 2012), using the remaining 73 LUHMES samples and 59 reference samples, also evenly distributed over the same three separate libraries, all of which were processed and sequenced in the exact same way. After removing three duplicates of technical replicates, 70 LUHMES samples were used for follow-up analyses. Corrected read counts were normalized by the sum of spike-in reads for each sample (Katayama et al., 2013). The significance of variation (fluctuation) of normalized expression levels was evaluated by comparing them with the technical variation of spike-ins, as described in Supplementary Text S1 of Krjutškov et al. (2016). In the line plots, normalized expression values were divided by the minimum non-zero value (5.78×10^-6) for visualization.

Pearson correlation coefficients between STRT RNA-seq and RT-qPCR data (Figs S3–S6) were calculated using the log₂ fold-change values from day 0 for each method (normalized read counts for STRT RNA-seq and 2^−ΔΔCt for RT-qPCR) per independent experiment.

Principal component analysis (PCA) was performed on significantly fluctuating TFEs [false discovery rate (FDR)-adjusted P<0.05] across all 70 LUHMES samples, using the pca function of the R package mixOmics (v6.8.0) (Rohart et al., 2017). Pearson correlation coefficients between samples were calculated using the log₂-scaled normalized expression values of significantly fluctuating TFEs. Differential expression analysis was carried out using the R package SAMstrt (v0.99.0) (Katayama et al., 2013), where we defined differentially expressed TFEs as those showing a differential expression with Benjamini–Hochberg (BH)-adjusted P<0.05 and fluctuation with FDR-adjusted P<0.05. Motif enrichment on differentially expressed TFEs was analyzed using the command findMotifsGenome.pl from HOMER (v4.10.3) (Heinz et al., 2010) with the option ‘-size -300,100 -mask’, using all the detected TFEs as background.

Clustering analysis by time-course expression pattern was carried out using the R package Mfuzz (v2.44.0) (Kumar and Futschik, 2007). A total of 5596 genes that showed highly significant fluctuation (FDR-adjusted P<1×10⁻¹⁰) for each experiment and high correlation between two experiments (Pearson correlation>0.7) were classified into four clusters, which were determined based on the minimum centroid distance. Gene ontology (GO) term enrichment analysis was performed using the R package enrichR (v2.1) (Kuleshov et al., 2016). Core genes in each cluster for the enrichment analysis were extracted based on high membership scores (>0.75) using the acore function of the Mfuzz package. Gene set enrichment analysis (GSEA) was performed with GSEA (v3.0; http://software.broadinstitute.org/gsea/) using the GSEAPreranked tool (Subramanian et al., 2005). Here, genes were ranked based on the combination of differential expression score and fold changes of their expression levels, and were then compared with the gene set of 302 ciliary genes from the SYSCILIA gold standard (SCGSv1) (van Dam et al., 2013). GSEA plots were made using ReplotGSEA.R (https://github.com/PeeperLab/Rtoolbox/blob/master/R/ReplotGSEA.R) with minor modifications. Integrative analysis with fetal midbrain single-cell RNA-seq data (GSE76381) (La Manno et al., 2016) was carried out with the CreateSinglerSeuratObject function of the R package SingleR (v1.0.1) (Aran et al., 2019) and Seurat (v3.0.2) (Satija et al., 2015), using the gene-based normalized expression data of LUHMES samples as the reference data set.

After removing duplicated, non-primary and unmapped reads from the BAM files, followed by sorting with SAMtools (v1.9), these processed BAM files from each library were subjected to the ISMARA (The Integrated System for Motif Activity Response Analysis) online tool (https://ismara.unibas.ch/mara/) (Balwierz et al., 2014) for the motif activity response analysis. Sequence logos were generated using the R package ggseqlogo (v0.1) (Wagih, 2017).

Parts of this study were performed at the Karolinska Institute Live Cell Imaging facility supported by grants from the Knut and Alice Wallenberg Foundation, the Swedish Research Council, the Centre for Innovative Medicine and the Jonasson Center at the Swedish Royal Institute of Technology. We acknowledge technical and personnel support from the electron microscopy facility of the Karolinska University Hospital, Huddinge. We acknowledge BEA, the Bioinformatics and Expression Analysis core facility, which is supported by the board of research at the Karolinska Institute and the research committee at the Karolinska Hospital. The Biomedicum Functional Genomics Unit at the University of Helsinki provided sequencing services. All computation for this work was performed on resources provided through the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under project contract SNIC 2017/7-317.

The peer review history is available online at https://jcs.biologists.org/lookup/doi/10.1242/jcs.249789.reviewer-comments.pdf