Transcriptomic analyses of NeuroD1-mediated astrocyte-to-neuron conversion.

Ectopic expression of a single neural transcription factor NeuroD1 can reprogram reactive glial cells into functional neurons both in vitro and in vivo, but the underlying mechanisms are not well understood yet. Here, we used RNA-sequencing technology to capture the transcriptomic changes at different time points during the reprogramming process. We found that following NeuroD1 overexpression, astroglial genes (ACTG1, ALDH1A3, EMP1, CLDN6, SOX21) were significantly downregulated, whereas neuronal genes (DCX, RBFOX3/NeuN, CUX2, RELN, SNAP25) were significantly upregulated. NeuroD family members (NeuroD1/2/6) and signaling pathways (Wnt, MAPK, cAMP) as well as neurotransmitter receptors (acetylcholine, somatostatin, dopamine) were also significantly upregulated. Gene co-expression analysis identified many central genes among the NeuroD1-interacting network, including CABP7, KIAA1456, SSTR2, GADD45G, LRRTM2, and INSM1. Compared to chemical conversion, we found that NeuroD1 acted as a strong driving force and triggered fast transcriptomic changes during astrocyte-to-neuron conversion process. Together, this study reveals many important downstream targets of NeuroD1 such as HES6, BHLHE22, INSM1, CHRNA1/3, CABP7, and SSTR2, which may play critical roles during the transcriptomic landscape shift from a glial profile to a neuronal profile.

INTRODUCTION

Mammalian central nervous system (CNS) is the control panel for the whole body, yet prone to heterogeneous traumas as well as pathological neurodegeneration. Once injured, adult CNS has very limited regeneration capability. Therefore, it is pivotal to develop remedies for neuronal replenishment and functional recovery after neuronal loss. Previous studies, including our own, have demonstrated that ectopic expression of transcription factors (TFs) can reprogram glial cells into neurons (Chen et al., 2020; Gascon et al., 2016; Ge et al., 2020; Guo et al., 2014; Heinrich et al., 2014; Li & Chen, 2016; Liu et al., 2020; Liu et al., 2015; Puls et al., 2020; Rao et al., 2021; Su et al., 2014; Torper et al., 2015; Wu et al., 2020; Xiang et al., 2021; Zhang et al., 2020). With this powerful neuronal conversion technology, researchers have attempted to repair damaged neural tissue and restore lost neuronal connectivity in multiple disease models. For example, we have first demonstrated that expressing a single neural TF NeuroD1 in 14‐month‐old mouse model for Alzheimer's disease can directly convert reactive astrocytes into functional neurons (Guo et al., 2014). Later on, this NeuroD1‐mediated astrocyte‐to‐neuron (AtN) conversion technology has been successfully applied in stab injury model to reverse glial scar tissue back to neural tissue (Zhang et al., 2020), and in ischemic stroke model to regenerate neural tissue and promote functional recovery (Chen et al., 2020; Tang et al., 2021). Recently, we have further demonstrated that NeuroD1‐based neuroregenerative gene therapy can successfully convert reactive astrocytes into neurons in nonhuman primate model (Ge et al., 2020), making this in vivo neural conversion approach one step closer to clinical trials. Besides NeuroD1, many other TFs including Neurogenin2 (Ngn2), Ascl1, Dlx2, and their combinations have been reported to convert glial cells into neurons both in vitro and in vivo (Barker et al., 2018; di Val Cervo et al., 2017; Gascon et al., 2016; Heinrich et al., 2014; Herrero‐Navarro et al., 2021; Liu et al., 2015; Su et al., 2014; Torper et al., 2015). These TF‐based transdifferentiation from glial cells into neurons provides a new path toward neural regeneration and repair.

Many studies have investigated the molecular mechanisms underlying intercellular transdifferentiation. For example, in the BAM (Brn2, Ascl1, Myt1l)‐mediated conversion of mouse embryonic fibroblast into neurons, Ascl1 acts as a pioneer factor and occupies many genomic sites on chromatin, while Myt1l can safeguard neuronal identity by repressing nonneuronal activators (Mall et al., 2017; Wapinski et al., 2013). Single‐cell RNA‐sequencing technology further found that Ascl1 triggered a relatively homogeneous initiation, forcing cells to exit cell cycle and adopt the neuronal maturation program (Treutlein et al., 2016). In the mouse astrocytes, Ngn2 overexpression can trigger a hierarchical sequence of gene activation and compete with REST for NeuroD4 promoter binding, which leads to the acquisition of neuronal identity (Masserdotti et al., 2015). A recent study revealed that Ascl1 targets Klf10, Myt1, NeuroD4, and Chd7 to convert astrocytes into neurons (Rao et al., 2021). Therefore, it seems that TFs operate in a concerted way in driving glial cell conversion.

NeuroD1 is a bHLH family TF originally identified by Weintraube and colleagues through injecting cDNAs into Xenopus oocytes that induced neural differentiation (Lee et al., 1995). Subsequent studies found that NeuroD1 not only promotes neural differentiation from neural stem cells but also promotes neuronal maturation and survival (Miyata et al., 1999). Besides a role in early brain development, NeuroD1 has also been found in adult neural stem cells in the mouse hippocampus (Gao et al., 2009; Kuwabara et al., 2009). Our lab initially tried to use Ngn2 to convert astrocytes into neurons in the mouse brain in vivo but found the efficiency was quite low. We then turned to NeuroD1 and achieved high AtN conversion efficiency in a series of neurological disease models in rodents as well as in nonhuman primates (Chen et al., 2020; Ge et al., 2020; Guo et al., 2014; Liu et al., 2020; Puls et al., 2020; Tang et al., 2021). Here, we used RNA‐sequencing technology to investigate the molecular mechanisms behind the NeuroD1‐mediated AtN conversion. Starting from cultured human astrocytes (HA), we analyzed the changes of transcriptome profile at different time points following NeuroD1 expression. Dissection of the differentially expressed genes (DEGs) within the first 2 weeks of the AtN conversion process revealed early response genes to NeuroD1 overexpression, including neural TFs and regulators of multiple signaling pathways such as Wnt, Notch, hedgehog, and MAPK. Our transcriptome analyses depicted a molecular roadmap of sequential upregulation of neuronal genes and downregulation of glial genes during AtN conversion. These findings will enhance our understanding of molecular programs instructing cell fate conversion, and facilitate the development of neuroregenerative therapies to treat neurological disorders.

RESULTS

Global transcriptomic changes after overexpressing NeuroD1 in HA

We have previously demonstrated that overexpressing a single neural TF NeuroD1 can efficiently convert astrocytes into functional neurons both in vitro and in vivo (Guo et al., 2014). To investigate the molecular mechanisms underlying such direct AtN conversion, we infected HA (HA1800, ScienCell) in culture with NeuroD1 retroviruses and performed RNA‐sequencing at day 1, 3, 5, and 14 following viral infection to interrogate the transcriptomic changes in the cultured cells (see Figure 1a for experimental design). Retroviruses expressing GFP alone served as controls. Consistent with our published work (Guo et al., 2014), more than 80% of NeuroD1‐infected HA were converted into neurons after 2 weeks of infection (Figure S1a). To avoid bias toward early or high expression cells, we collected RNA samples from the entire pool of cultured cells for high‐throughput sequencing analysis, similar to that performed following chemical reprogramming (Ma et al., 2019). The raw data count was examined by boxplot and density plot, and no outlier sample was identified (Figure S1b,c).

FIGURE 1

Experimental design and overall global transcriptomic analysis. (a) Illustration of the experimental design. Each time point has three replicates. (b) Hierarchical clustering of sample relationship based on global expression profile. (c) Principal component analysis of 18 RNA samples. PC1 = 43.9%, PC2 = 28.4%. (d) Bar plot of the numbers of up‐ or downregulated DEGs in pairwise comparisons with untreated astrocytes. (e) Venn diagram of comparisons shows that GFP control and NeuroD1‐GFP samples at D1 infection share the majority of DEGs caused by viral infection.

We performed sample comparison using hierarchical clustering (Figure 1b) and principal component analysis (PCA) (Figure 1c). On the whole transcriptome scale, the six samples at 1 day post viral infection clustered together, regardless of GFP or NeuroD1 retrovirus, suggesting that virus infection itself had a dominant effect at day 1 (Figure 1b). Day 3 samples were more closely related to day 5 than nontreated HA, and day 14 samples further shifted away from the rest of the samples, likely due to transdifferentiation into neuronal transcriptomes after 2 weeks of NeuroD1 expression (Figure 1b). This clustering relationship is also illustrated in the pattern of PCA analysis (Figure 1c). When plotted without the untreated HA samples, a trajectory from day 1 to day 14 can be observed (Figure S1d). We next performed the analysis of DEGs (fold change >2, expression base mean >50, adjusted p‐value < .01) (Figures 1d,e and S1e). Compared to untreated astrocytes, day 1 samples infected by GFP or NeuroD1 virus exhibited large number and similar DEGs, whereas day 3 samples showed less DEGs, with many from day 1 samples (Figure 1d,e).

From all pairwise comparisons throughout the 2‐week conversion process, we identified 2994 DEGs and divided them into three groups using hierarchical clustering based on the gene changing trend (Figure 2a). Gene ontology (GO) analysis identified three distinct categories with different set of genes highly upregulated during the neuronal reprogramming process. The first group was strongly related to virus infection and interferon signaling and highly expressed in day 1 samples infected by GFP and NeuroD1 virus (Figure 2a, red box). The second cluster had 827 genes involved in cell cycle and telencephalon development and was turned on during day 3–5 after NeuroD1 infection (Figure 2a, blue box). The third group contained 957 genes, which were expressed at late stage of reprogramming after 2 weeks of NeuroD1 expression and the GO categories showed neuronal gene clusters including “neuron projection” and “synapse assembly” (Figure 2a, green box).

FIGURE 2

DEG analysis revealed early inflammatory responses and transcriptomic shift toward neurons at late stage. (a) Heatmap of RNA‐seq data of 2994 DEGs, which were grouped into three clusters with most significant functional ontologies annotated. (b) Volcano plot of up‐ and downregulated genes in the control virus D1 group compared to the HA group (untreated astrocytes). The cutoff values are −1 and 1 for log2(fold change), and 2 for ‐log10(p‐adj). (c and d) Top up‐ and downregulated DEGs from panel (b)

To better understand the virus effect, we compared control GFP virus at D1 with untreated HA. The volcano plot showed 683 upregulated genes (URGs, red) and 510 downregulated genes (DRGs, blue; Figure 2b). Top URGs were mostly interferon‐induced genes (Figure 2c), such as 2′−5′‐oligoadenylate synthetase (OAS1), viperin (RSAD2), and viral RNA degradation gene IFI27. They were hardly expressed in untreated astrocytes, suggesting that they were viral‐specific response genes. The DRGs were involved in biological processes such as cell cycle phase transition and DNA damage check (Figure 2d), including E3 ubiquitin ligase producer UHRF1, cell proliferation regulator MYBL2, and checkpoint arrestor CLSPN. Together, these data suggest that viral infection of HA triggers strong immune responses within 24 h, which decay quickly when NeuroD1 expression starts to show its effects.

Early response genes induced by NeuroD1

We next tried to identify the immediate downstream targets of NeuroD1 by examining the DEGs that were up‐ or downregulated at day 1 after NeuroD1 expression (Figures 3 and S2). Compared to the control GFP virus at day 1, DEG analysis revealed 95 URGs and 28 DRGs induced by NeuroD1 at day 1 (Figure 3a). Among these DEGs, 13 TFs were identified (Figure 3b), including five members from bHLH family (NEUROD1, HES6, BHLHE22, MYCL1, and NHLH1) and four zinc finger C2H2 factors (INSM1, PRDM8, KLF4, and EGR1). INSM1 has been reported to cooperate with NEUROD1 and FOXA2 in pancreatic β‐cells, or with ASCL1 to regulate neurotransmitter synthesis in vertebrate hindbrains (Jacob et al., 2009; Jia et al., 2015). PRDM8 protein can form a repressor complex with BHLHE22 (also known as BHLHB5) to control axon targeting and circuit assembly (Ross et al., 2012). Thus, these 13 TFs may act synergistically to mediate the immediate response following NeuroD1 overexpression.

FIGURE 3

NeuroD1 overexpression rapidly activated multiple transcription factors and neuronal receptors within the first 24 h. (a) A total of 123 DEGs were picked from the comparison between day 1 NeuroD1 samples versus the GFP control samples. (b) Representative transcription factors regulated by NeuroD1 at day 1. (c) Heatmap showing the distribution of downregulated DEGs (DRGs) and two classes of upregulated DEGs (URGs) among different samples. (d–f) Distinct patterns of some of the representative genes among DRGs and URGs

We further clustered the DEGs at day 1 into DRG, URG class I, and URG class II, based on their decrease or increase pattern in the following days during reprogramming process (Figure 3c). Compared to the GFP control group, most DRGs showed a continuous reduction from day 1 to day 5, such as CPNE7 (copine 7, a Ca2+‐dependent phospholipid‐binding protein) and SCNN1B and SCNN1G (epithelial sodium channel subunits for Na+ absorption) (Figure 3d). The continuous downregulation of these DRGs during NeuroD1‐mediated cell conversion suggested that these genes may be specific to astrocytic functions and thus downregulated during neuronal conversion. The URGs were categorized into two classes: the class I was fast‐response genes that quickly peaked on day 1 but started to decrease at day 3 (Figures 3e and S2a), whereas the class II was slow‐response genes that started to increase at day 1 and further increased at day 3 (Figure 3f). The fast‐response class I URGs included NEUROD1 itself, which peaked at day 1 and then gradually decreased at day 3 and day 5 (Figure 3e). A surprising finding is that some neuronal receptor genes including acetylcholine receptor subunit CHRNA1 and CHRNA3 as well as dopamine receptor subunit DRD2 also showed a transient increase at day 1 followed by a decrease at day 3 and day 5 (Figure 3e). A transient increase of these receptor genes immediately following NeuroD1 expression might suggest a new function of these genes, which certainly warrants further investigation. Another immediate response gene is WISP1, which encodes a Wnt inducible signaling pathway protein and functions as a connective tissue growth factor to cause collagen linearization (also known as CCN4) (Figures 3e and S2b). WISP1 has also been reported to participate in apoptosis and stem cell proliferation through interacting with PI3K/AKT/mTOR/MAPK/JNK/GSK‐3β multiple signaling pathways (Maiese, 2014). In general, genes in this fast response cluster experienced the most significant increase in the first 24 h before starting to decrease at day 3, indicating that they may play an important role in the initiation stage and start the cascades of cell reprogramming.

The class II slow‐response URGs not only showed an instant increase in day 1 but also showed a sustained increase at day 3 and sometimes even day 5 or day 14 (Figure 3f). One interesting gene among this class II URGs is INSM1 (insulinoma associated 1) (Figure 3f), a zinc finger protein playing an important role in neurogenesis and neuroendocrine cell differentiation. Importantly, INSM1 is also acting as a transcriptional repressor of NEUROD1, which might explain the transient increase of NEUROD1 at day 1 followed by a decrease at day 3 and day 5 (Figure 3e). Another two slow‐response genes known to be enriched in neuroprogenitor cells are MFAP4, an extracellular matrix protein, and SSTR2 (somatostatin receptor 2), a G‐protein coupled receptor functioning as a negative regulator of proliferation (Buscail et al., 1995; Yuzwa et al., 2017). IGFBPL1 (IGF‐1 binding protein like 1) and PTCHD2 (patched domain‐containing protein 2) also showed continuous increase from day 1 to day 5, and both were reported in neural tissue to regulate neurogenesis and synaptogenesis (Guo et al., 2018; Konirova et al., 2017). In addition, our KEGG pathway analysis revealed that cholinergic synapse, circadian entrainment, gap junction, MAPK/cAMP/Rap1 signaling pathways, and neuroactive ligand–receptor interaction pathways were also significantly modulated on day 1 (Figure S2). These results reveal critical genes that showed rapid response to NeuroD1 expression within 24 h, which include TFs, neural receptors, and Hedgehog/Wnt/MAPK pathways.

NeuroD1 activation of endogenous neural TFs

After understanding the NeuroD1‐induced transcriptomic alterations with 123 DEGs on day 1, we next investigated further dynamic changes in the intermediate stage on day 3. Comparing NeuroD1 samples at day 3 versus day 1, volcano plot identified 774 DEGs, including 437 DRGs and 337 URGs (Figure 4a). These DEGs are enriched in various cellular functions and signaling pathways, such as cell proliferation, nucleotide catabolism, and immune system‐related genes (defense response to virus, NF‐κB signaling, JAK/STAT3, interferon‐γ response) (Figures 4b and S3). The downregulation of immune response genes at day 3 compared to day 1 suggested that the viral response had peaked at day 1 and now faded at day 3.

FIGURE 4

Activation of neurogenic factors along with the receding of viral responses after 3 days of NeuroD1 expression. (a) Volcano plot of genes at D3 versus D1 in the NeuroD1 group. (b) Enriched GO categories visualized in REVIGO. The p‐values are color coded. The size indicates the generalness of the GO term. Similar GO terms remained close together in the plot. (c) Bar plot of 21 transcription factors with fold change >3 from panel (a). Colors correspond to log2 expression base mean. (d) Activated and suppressed transcription factors with top fold changes

Among all the DEGs, we identified 52 TFs and half of them showed more than threefold change (Figure 4c). Interestingly, NeuroD1 itself showed a significant decrease, whereas NeuroD6 and NeuroD2 showed a significant increase at day 3 (Figure 4c), suggesting that NeuroD1 had passed its effects to the downstream targets including other NeuroD family members. Besides NeuroD family, the upregulated TFs induced by NeuroD1 were mostly involved in neurogenesis (Figure 4c,d). For example, SCRT2 (Scratch 2) is a transcriptional repressor that is important for neurogenesis and neuronal migration during embryonic development (Rodriguez‐Aznar et al., 2013). MYT1 (myelin TF 1) and NHLH1 are both targets of NeuroD1 (Seo et al., 2007), and play functional roles in early neural developmental stage (Figure 4d, top panel). The downregulated TFs include interferon regulatory factor IRF7 and immune system regulator BATF2, which indicated that the viral infection response has largely been reduced on day 3 (Figure 4c,d). Together, on day 3, neural TFs have been highly activated, while the virus‐induced immune response genes have started to fade.

Identification of NeuroD1‐related coexpression network

Many genes function together through close interactions and display similar expression profiles. Therefore, we constructed a weighted co‐expression network to detect gene functional modules. First, genes with coefficient of variation of expression levels less than 0.01 were removed, and a total of 15,215 genes were retained for coexpression analysis. A soft‐thresholding power of 10 was picked for scale‐free network construction (Figure S4a). Based on topological overlap matrix (TOM), hierarchical clustering generated a dendrogram of genes that was divided into 15 modules (Figure 5a). To examine the relationship between gene clusters and stage‐specific functions, we assigned three arbitrary traits, including time, infection, and initiation to each sample (Figure S4b). Time was set as day(s) from the beginning of the conversion and monotonically increases. Infection was assigned to day 1 and day 3 samples to mimic viral infection responses. Initiation was assigned to capture the early response following NeuroD1 expression. To assess the similarities among the 15 modules assigned by the three traits (Figure S4b), eigengenes of each module were calculated and clustered. Module eigengene (ME) is the first principal component of the expression matrix and can best represent the module. From the eigengene dendrogram (Figure 5b), five modules were close to trait initiation, among which purple module had the highest gene significance (0.726 ± 0.025; Figure S4c). These genes (see Figure 5d, purple color genes) were either activated or suppressed right after NeuroD1 overexpression, suggesting that they may directly interact with NeuroD1 and play important roles during the initiation stage. Figure 5c listed the number of genes and DEGs among each module. The turquoise module is the largest one with 3527 genes, and covers many members from STAT, ISG, and IFI families. These gene families are known to be induced by interferon and involved in JAK/STAT signaling (Morales & Lenschow, 2013; Schindler et al., 2007). For the trait time, module brown has the highest gene significance (0.757 ± 0.024) (Figure S4d). There are 765 DEGs out of 2979 genes in the brown module, including those related to neuronal function and maturation, such as DCX, GRIA, and RBFOX3 (NeuN).

FIGURE 5

Modular associations with conversion traits and network analysis of NeuroD1 correlations. (a) An example of gene dendrogram. The color block below the dendrogram indicates the corresponding module color. (b) Dendrogram of eigengenes from 15 modules together with three traits (initiation, time, and infection). (c) The number of genes and DEGs associated with each colored module. (d) Network of NeuroD1 with its 49 interacting genes. Line width: weight (min = 0.0022, max = 0.7175). Circle size: degree, that is, connected nodes (min = 1, max = 17). Filled color: betweenness centrality (min = 0.0, max = 0.2835).

Next, we constructed a NeuroD1‐interacting gene network based on interactions with a similarity score >0.15, which resulted in a network containing 50 nodes and 119 edges (Figure 5d). With gene names labeled in corresponding module colors, we found that most of the NeuroD1‐interacting genes were in purple and black modules. One of the central gene among the NeuroD1‐interacting network is CABP7, a calcium binding protein regulating lysosome clustering. Other central genes include SSTR2 (somatostatin receptor type 2, affecting neurotransmission and hormone secretion), KIAA1456 (regulating DNA methylation and cell cycle), GADD45G (a stress responding gene and a negative regulator of the Jak‐Stat3 pathway), INSM1 (a TF regulating embryonic neurogenesis), Sox4 (an oncogenic TF regulating epithelial–mesenchymal transition and cellular proliferation and differentiation), MFAP4 (an extracellular matrix protein regulating cell adhesion), LRRTM2 (a neurexin‐binding cell adhesion molecule regulating synapse formation), and many more (Figure 5d). The NeuroD1‐interacting genes identified here during AtN conversion process are partly consistent with previous report on NeuroD1 targets in cancer cell lines (Borromeo et al., 2016), suggesting that some of the interacting gene network is universal after NeuroD1 overexpression.

Since NeuroD1 is a bHLH family TF, which often activates a neurogenic TF network during neural development, we analyzed a total of 98 differentially expressed TFs and removed those related to virus infection (from module turquoise). Their correlation relationship is displayed in Figure 6a with node colors corresponding to the modules in Figure 5. Many of the well‐connected nodes, including NeuroD family members NeuroD2/NeuroD6, came from the brown module, which is closely related to trait time (Figure S4d). We then looked into the temporal expression pattern of cell type‐specific markers, focusing on astrocytes, neurons, and neuroprogenitors (Figure 6b). Astroglial genes were universally downregulated, whereas neuronal genes were universally upregulated, as expected (Figure 6b). The progenitor genes displayed a mixed pattern, with various factors turned on at different time points during conversion process, but then all downregulated at day 14 after NeuroD1 infection, suggesting that neuronal conversion is completed in 2 weeks. Consistently, enrichment score analysis also found that gene sets of Notch signaling and EMT (epithelial mesenchymal transition) were highly activated on day 14 (Figure 6c). Together, these results suggest that overexpression of NeuroD1 triggers a landscape change in many TFs, which together downregulate astroglial genes and simultaneously activate neuronal genes to change astrocytes into neuronal fate.

FIGURE 6

Regulatory network of transcription factors and cell type‐specific gene expression. (a) Network of 98 differentially expressed transcription factors from modules brown, blue, black, magenta, yellow, red, green, and pink. Circle size: degree (min = 1, max = 15). (b) Heatmap of progenitor‐, neuron‐, and astrocyte‐related genes over the reprogramming time course. (c) Gene set enrichment analysis of D14 compared to human astrocytes, showing Notch and EMT signaling pathways activated at day 14 after NeuroD1 expression.

Comparison between transcriptome changes induced by NeuroD1 and small molecules

In previous work, we have demonstrated that cultured HA can be converted into neurons using a combination of four chemicals (together called core drugs) and investigated the molecular processes using RNA‐seq (Ma et al., 2019; Yin et al., 2019). Comparing the RNA‐seq datasets, we found that the small‐molecule approach and the NeuroD1 reprogramming methods showed unique features as well as shared mechanisms. For example, in both conditions the astroglial genes were suppressed, while neurogenic factors were transiently activated, followed by sustained expression of typical neuronal genes. Figure 7a illustrates the heatmap of 1104 overlapping DEGs between the core drug‐induced and NeuroD1‐induced AtN conversion process. It is interesting to note that while many DEGs showed up in both conditions, their temporal expression patterns were quite different. For example, the first cluster (upper half of heatmap in Figure 7a) was mainly expressed in preconverted astrocytes, as indicated by the high level in HA, D0, and GFP groups. Representative glial genes including GFAP, ACTN1, SMAD9, and TGFBI were almost totally suppressed following core drug treatment and partially decreased following NeuroD1 treatment (Figure 7a), likely because NeuroD1 retrovirus only partially infected the cultured astrocytes but in core drug condition, all astrocytes were affected. The associated GO terms in the downregulated DEGs were extracellular matrix organization, cell adhesion, and TGF‐β response genes (Figure 7a). The second cluster (lower half of heatmap in Figure 7a) includes most upregulated DEGs and was also stronger in core drug condition than that in NeuroD1‐treated group. This cluster is related to GO terms of neurogenesis, neuronal maturation, synaptic signaling, and axonogenesis, with representative genes including TBR1, SEMA5A, and NLGN3.

FIGURE 7

Comparison between transcriptomic changes induced by NeuroD1 versus chemical reprogramming. (a) Heatmap of 1104 overlapping DEGs between NeuroD1 treatment and small molecule treatment (core drugs). Yellow indicates high expression; purple indicates low. The biological processes enriched in the top cluster and bottom cluster were annotated by GO analysis on the right. (b–e) Line plots of log‐read counts of representative up‐ and downregulated genes in the NeuroD1 group versus the core drug group. Solid lines for core drug group, and dashed lines for NeuroD1 group.

Next, we selected several signature DEGs among both groups and compared their temporal dynamics during conversion process (Figure 7b). In NeuroD1 group (dashed line), NeuroD1 and its target gene CHRNA3 reached peak on day 1 and then began to drop, whereas these two genes under core drug treatment showed a gradual increase from day 1 and peaked at day 5 (Figure 7b; NeuroD1, red; CHRNA3, greenyellow). Interestingly, core drug‐activated genes, such as NNAT (neuronatin, an important gene in brain development and neural differentiation) and PENK (proenkephalin), were also upregulated in NeuroD1 group but they reached peak 2 days in advance compared to that in the core drug group (Figure 7b; PENK, blue; NNAT, green). These patterns indicated that NeuroD1 conversion features a more expedited progress in some transcriptional changes compared to chemical reprogramming. On the other hand, some signaling genes including CNTN2 (contactin‐2, a GPI‐anchored neuronal membrane protein, red), IGFBPL1 (an IGF‐binding protein regulating mTOR phosphorylation, blue), NHLH1 (nescient Helix–Loop–Helix 1 protein playing a wide variety functions in development, greenyellow), and MAPK4 (mitogen‐activated protein kinase 4 regulating AKT/mTOR signaling pathways, green) showed continuous high upregulation from day 1 to day 5 under core drug treatment, but they were modestly changed in NeuroD1 group (Figure 7c). We also observed significant differences between these two reprogramming approaches. For example, the “starter” genes for chemical conversion, such as hedgehog genes ARHGAP36 and PTCH1 as well as GPR98 (a calcium‐binding G protein‐coupled receptor) and SOX21 (a regulator of neuronal gene expression), were highly upregulated upon core drug treatment at day 1 (Figure 7d, solid lines), but they were not turned on at all or even downregulated in NeuroD1 group (Figure 7d, dashed lines). The opposite direction was also observed, such as that in core drug treatment group, the ID1 (inhibitor of differentiation/DNA binding, a Helix‐Loop‐Helix protein), COL3A1 (collagen), DCN (decorin, an extracellular matrix protein), and PRPH (cytoskeletal peripherin) were downregulated by core drugs but slightly upregulated in NeuroD1 group (Figure 7e). Therefore, while both NeuroD1 and small molecules reprogram astrocytes into neurons, they alter the astroglial transcriptome profile by targeting different biological signaling pathways, yet in 2 weeks both achieved successful conversion into neuronal transcriptome profile.

DISCUSSION

We have previously demonstrated that ectopic expression of a neural TF NeuroD1 can reprogram astrocytes into neurons both in vitro and in vivo. In this work, we analyzed the transcriptomic changes at different time points following NeuroD1 expression and decoded the sequential changes of gene network during AtN conversion process. We found that NeuroD1 expression activated a cluster of neural TFs that strongly inhibited astroglial genes and simultaneously activated neuronal genes within the first 5 days. Expression of NeuroD1 in astrocytes also inhibited genes involved in cell proliferation within the first 3 days, and by 14 days following NeuroD1 expression, AtN conversion process was largely completed. These results revealed a molecular cascade that triggered a transcriptomic shift from astrocytic profile to neuronal profile after overexpression of a single neural TF NeuroD1.

Downregulation of glial genes by NeuroD1

After overexpressing NeuroD1 in astrocytes, we observed rapid downregulation of glial genes within 24–72 h of viral infection. Such rapid transcriptomic changes after overexpressing neural TFs are consistent with previous reports on Ascl1 or Ngn2 (Masserdotti et al., 2015; Rao et al., 2021; Wapinski et al., 2013). It is interesting to note that among the earliest downregulated genes, non‐voltage‐gated epithelial sodium channel subunit genes SCNN1B and SCNN1G stand out. SCNN1B and SCNN1G function in keeping the electrolyte homeostasis in epithelial cells (Voilley et al., 1995; Zhong et al., 2016). Since neuronal sodium channels are highly voltage dependent, these non‐voltage‐dependent sodium channels in astrocytes must be replaced by voltage‐dependent sodium channels in order to carry out the conversion process. Another interesting gene downregulated quickly by NeuroD1 in astrocytes is CPNE7, a Ca2+‐dependent phospholipid‐binding protein that regulates Ca2+‐dependent intracellular processes (Caudell et al., 2000; Seo et al., 2017). Ca2+ is an important second messenger and regulates many critical signaling pathways including Ca2+‐dependent protein kinases and phosphatases. Neuronal Ca2+ signaling is largely dependent on voltage‐gated Ca2+ channels and receptors such as NMDA receptors, whereas astrocytes rarely use voltage‐gated channels. Therefore, downregulation of CPNE7 may be an important step to prepare the converting cells adopting a neuronal Ca2+ signaling system. In addition, many other glial genes including ALDH1A3, CLDN6, EMP1, SLC7A1, SOX21, and TPST2 are also significantly downregulated within 3 days of NeuroD1 expression. Downregulation of these glial genes will likely pave the way for the NeuroD1‐infected cells to adopt neuronal structural and signaling compositions.

Neural TFs and neuronal genes activated by NeuroD1 in astrocytes

NeuroD1 is a bHLH family TF that has been originally reported to induce neural differentiation in Xenopus oocytes (Lee et al., 1995). Later studies found that NeuroD1 is critically involved in neuronal differentiation not only during embryonic brain development (Miyata et al., 1999) but also in adult neurogenesis (Gao et al., 2009; Kuwabara et al., 2009). We have previously reported in a series of studies that overexpression of NeuroD1 in astrocytes both in vitro and in vivo, and regardless the species of human, monkey, or rodents, can convert astrocytes into neurons (Chen et al., 2020; Ge et al., 2020; Guo et al., 2014; Liu et al., 2020; Puls et al., 2020; Tang et al., 2021; Wu et al., 2020; Xiang et al., 2021; Zhang et al., 2020). In this study, we employed RNA‐sequencing technology to investigate molecular mechanisms underlying such AtN conversion. We found that immediately following NeuroD1 expression in HA, a network of TFs was quickly activated within 24 h that might act synergistically to trigger wider transcriptomic changes in the infected astrocytes. Interestingly, the TFs activated by NeuroD1 expression in astrocytes are largely overlapping with those reported by NeuroD1 expression in embryonic stem cells, including Hes6, Insm1, Prdm8, Nhlh1, and Rcor2 (Pataskar et al., 2016). This similarity suggests that the NeuroD1 downstream effectors are rather conserved within different cells, and NeuroD1 can act as a master regulator to synchronize many other downstream TFs to promote neuronal fate determination.

Besides downstream targets, NeuroD1 itself appeared to be rapidly downregulated after reaching high level at 1 day post retroviral infection. This suggests that NeuroD1 transcription can be rapidly activated in cultured astrocytes using retroviruses. This high transcription level of NeuroD1 might include both the transcription from the retroviruses and from the activation of endogenous NeuroD1 in the astrocytes. A continuous downregulation of NeuroD1 expression at day 3 and day 5 suggests that there may be a negative feedback within astrocytes to prevent overactivation of NeuroD1. Similar to NeuroD1, previous studies have reported using other bHLH family TFs such as Ngn2 and Ascl1 to convert astrocytes or fibroblasts into neurons (Masserdotti et al., 2015; Rao et al., 2021; Wapinski et al., 2013). A common scene emerged from these different studies using different TFs is a rapid transcriptome change, typically within 1–3 days following TF expression, and consistent upregulation of neuronal genes within 1–2 weeks (Masserdotti et al., 2015; Rao et al., 2021; Wapinski et al., 2013). On the other hand, each TF appears to activate a unique cluster of downstream factors that eventually converge onto the same upregulation of neuronal genes such as NeuN, MAP2, and synapsin. Masserdotti et al. (2015) compared the transcriptome change induced by Ngn2 versus Ascl1 and found very little commonality between these two factors, except that NeuroD4 was identified as the shared downstream target. Recent study on Ascl1‐mediated astrocyte conversion confirmed the importance of NeuroD4 (Rao et al., 2021). However, in our NeuroD1‐mediated astrocyte conversion, we detected a significant upregulation of NeuroD6 and NeuroD2 but not NeuroD4, confirming that each pioneering factor (such as NeuroD1, Ngn2, Ascl1) may activate a distinct of downstream effectors to trigger AtN conversion. Interestingly, we detect a significant upregulation of INSM1 at day 1 following NeuroD1 overexpression, which is consistent with that reported for Ngn2 (Masserdotti et al., 2015), suggesting that INSM1 may be a common factor for reprogramming glutamatergic neurons.

While NeuroD1 activation of other TFs may not be surprising, one interesting finding is that some neurotransmitter receptors are rapidly upregulated by NeuroD1 within 24 h and then quickly downregulated. These include receptors for acetylcholine (CHRNA1, CHRNA3), dopamine (DRD2), and somatostatin (SSTR2). Since we do not expect that astrocytes would turn into neurons instantaneously after NeuroD1 expression, we wonder why these neurotransmitter receptors would be transiently activated so early and then decreased by day 3? Is it because astrocytes might have these receptors already? Surprisingly, such transient activation of neurotransmitter receptors is not a unique feature of astrocytes, but also has been reported as the target of NeuroD1 in human‐derived small‐cell lung cancer cell lines (Borromeo et al., 2016). Why these neurotransmitter receptor genes would be activated in cancer cells? What could be their function in these cancer cells? We speculate that such transient increase of receptor genes immediately following NeuroD1 expression might suggest a novel function of these genes that is beyond the classical role of receptors for binding with neurotransmitters.

In addition to upregulation of TFs and receptor genes, we also identified several well‐connected genes in the NeuroD1‐interacting network, such as CABP7, LRRTM2, and KIAA1456. For example, CABP7 is a calcium‐binding protein regulating cytokinesis and lysosome clustering in mammalian cells (Rajamanoharan et al., 2015). LRRTM2 is a neuronal‐specific leucine‐rich repeat transmembrane protein, and regulates the AMPA receptor expression as well as excitatory synapse formation (de Wit et al., 2009). KIAA1456 is also called TRMT9B, a tRNA methyltransferase that inhibits cell proliferation (Wang et al., 2017). Other significantly upregulated genes also include Wnt/IGF/MAPK/Hedgehog pathway regulators, such as WISP1, IGFBPL1, PTCHD2, and INSM1. The activation of neural TFs, neural receptor genes, and key signaling pathways by NeuroD1 forms a tightly interwoven network that drives astrocytes changing their glial properties toward neuronal properties.

Comparison between reprogramming induced by TFs and chemical compounds

Besides TF NeuroD1‐mediated AtN conversion, we have previously reported that combinations of small molecules can also convert HA into neurons with high efficiency (Yin et al., 2019; Zhang et al., 2015). We have also performed RNA‐seq to investigate the molecular mechanisms underlying the chemical reprogramming of astrocytes into neurons (Ma et al., 2019). Since both TF and small molecules can convert astrocytes into neurons, we wonder whether these two approaches share certain similarities in terms of transcriptomic changes. Interestingly, comparison of the NeuroD1‐induced transcriptomic changes with that induced by small molecules revealed quite different patterns. The most striking one is that while quite some genes are both upregulated or downregulated by NeuroD1 and small‐molecule strategies, they are regulated in very different paces during the conversion process. Not surprisingly, some genes are rapidly upregulated by NeuroD1 within 24 h, but rising rather slowly during small‐molecule treatment. On the other hand, some genes are up‐ or downregulated by small molecules, but not showing dramatic change following NeuroD1 expression. Therefore, NeuroD1 and small molecules may employ quite different mechanisms to achieve the same result of converting astrocytes into neurons.

CONCLUSION

NeuroD1 expression in astrocytes rapidly inhibits glial genes and activates other neural TFs simultaneously to orchestrate a landscape change of astroglial transcriptome profile toward neuronal transcriptome profile. Our transcriptomic analysis of NeuroD1‐induced AtN conversion is in sharp contrast to a recent work claiming that NeuroD1 cannot convert astrocytes into neurons (Wang et al., 2021). In fact, we have already reported successful conversion of the lineage‐traced astrocytes into neurons using the same Aldh1l1‐CreERT2 mice (Xiang et al., 2021). Another recent study also reported limited conversion of lineage‐traced astrocytes into neurons using the Aldh1l1‐CreERT2 mice (Leib et al., 2022). One potential explanation for the different conversion efficiency using the same TF NeuroD1 and the same transgenic mice may be due to different Adeno‐associated virus (AAV) vectors used in different labs that result in different NeuroD1 expression level in the lineage‐traced astrocytes. To solve the discrepancies among different labs, we recommend future work using a safe range of AAV titer (1011 to 1012 GC/ml) to investigate the expression level of NeuroD1 reached by different AAV vectors with different promoter and enhancer in the lineage‐traced astrocytes (Chen, 2021). We predict that lineage‐traced astrocytes will be converted into neurons in a NeuroD1 dose‐dependent manner.

EXPERIMENTAL PROCEDURES

Human cortical astrocyte culture

Human astrocytes HA1800 were purchased from ScienCell and cultured as previously described (Guo et al., 2014; Zhang et al., 2015). Briefly, cells were cultured on poly‐D‐lysine‐coated coverslips in 24‐well plates. The HA medium contains DMEM/F12, 10% fetal bovine serum (FBS), B27 supplement, 3.5 mM glucose, 10 ng/ml FGF2, 10 ng/ml EGF, and penicillin–streptomycin. Astrocytes were cultured to 90% confluent before split.

Retrovirus production and infection

The construct of the retroviral vectors pCAG‐NeuroD1‐IRES‐GFP and pCAG‐GFP‐IRES‐GFP was described in previous study (Guo et al., 2014). Virus particles pseudotyped with vesicular stomatitis virus G (VSV‐G) were packaged in HEK cells, and the titration was measured to be about 108 particles/μl.

When the astrocytes culture reached 80% confluent, 0.5 μl virus was added into each well. Twenty‐four hours later, the culture medium was completely changed into differentiation medium containing DMEM/F12, 0.5% FBS, N2, B27, 5 mg/ml vitamin C, and penicillin–streptomycin. Neurotrophic factors including 10 ng/ml BDNF, 10 ng/ml NT3, and 20 ng/ml IGF‐1 were also added to the culture.

Sample collection and RNA extraction

Control samples (HA) were collected from untreated HA. Twenty‐four hours after transfection, astrocytes infected with NeuroD1 or GFP control virus were collected, respectively (D1‐ND1, D1‐GFP). On days 3, 5, and 14, only astrocytes expressing NeuroD1 were collected (D3, D5, and D14). Each time point contained three biological replicates. RNA extractions were performed according to Macherey‐Nagel NucleoSpin® RNA kit protocol. It uses spin‐column‐based technology and the purified RNA was eluted with 40 μl RNase‐free water. The concentration was measured on NanoDrop™ 2000 spectrophotometers (Thermo Fisher).

RNA‐sequencing and analysis

All the RNA samples were sent to the UCLA Technology Center for Genomics & Bioinformatics for RNA quality check, mRNA enrichment, library construction, and single‐end 50 bp sequencing with HiSeq 3000. The raw data were checked by FastQC (v. 0.11.3) before aligning them against human reference genome hg38 by HISAT2 (v. 2.0.1) (Kim et al., 2015); they were summarized using featureCounts (v. 1.5.0) (Liao et al., 2014). Pairwise differential expression analysis was done using DESeq2 (v. 1.16.1) (Love et al., 2014). GO analysis was performed using Gene Ontology Consortium (Mi et al., 2017), Gorilla (Eden et al., 2009), Enrichr (Kuleshov et al., 2016), and g:Profiler (Reimand et al., 2007), and results were visualized using REVIGO (Supek et al., 2011). Enrichment plots were generated using GSEA (v.3.0) (Subramanian et al., 2005). Gene expression changes in significant KEGG pathways were shown by IPA (QIAGEN Inc.) and PathView (Luo & Brouwer, 2013). The correlation network was plotted using WGCNA (Langfelder & Horvath, 2008) and igraph (Csardi et al., 2010) and visualized in Cytoscape (Shannon et al., 2003).

CONFLICTS OF INTEREST

G.C. is a co‐founder of NeuExcell Therapeutics Inc. The other authors declare no conflict of interest.

Supporting Information

Click here for additional data file.