In vitro atlas of dorsal spinal interneurons reveals Wnt signaling as a critical regulator of progenitor expansion.

Restoring sensation after injury or disease requires a reproducible method for generating large quantities of bona fide somatosensory interneurons. Toward this goal, we assess the mechanisms by which dorsal spinal interneurons (dIs; dI1-dI6) can be derived from mouse embryonic stem cells (mESCs). Using two developmentally relevant growth factors, retinoic acid (RA) and bone morphogenetic protein (BMP) 4, we recapitulate the complete in vivo program of dI differentiation through a neuromesodermal intermediate. Transcriptional profiling reveals that mESC-derived dIs strikingly resemble endogenous dIs, with the correct molecular and functional signatures. We further demonstrate that RA specifies dI4-dI6 fates through a default multipotential state, while the addition of BMP4 induces dI1-dI3 fates and activates Wnt signaling to enhance progenitor proliferation. Constitutively activating Wnt signaling permits the dramatic expansion of neural progenitor cultures. These cultures retain the capacity to differentiate into diverse populations of dIs, thereby providing a method of increasing neuronal yield.

INTRODUCTION

Spinal cord injuries (SCIs) can result in the loss of both coordinated movement and somatosensation (Ionta et al., 2016; Lenggenhager et al., 2012), as well as the ability to perceive the environment. Without somatosensation, patients can no longer communicate through touch, sense danger through pain, or control the functions of the gut and bladder (Furusawa et al., 2011). Although stem-cell-based cellular replacement therapies have shown promise for SCI (Antonic et al., 2013; Kong et al., 2021), the focus so far has been on the recovery of coordinated motor function (Faravelli et al., 2014; Trawczynski et al., 2019). Less attention has been paid to regaining sensory functions, despite their critical roles ensuring physical and emotional well-being, and providing modulatory feedback to the motor system through proprioceptive pathways. Somatosensory information is received peripherally and then relayed centrally by different populations of dorsal spinal interneurons (dIs) in the spinal cord. Each dI population is specialized for distinct somatosensory modalities and forms microcircuits that span different layers of the dorsal horn (Gupta and Butler, 2021; Lai et al., 2016). For example, dI1s (proprioceptors) (Lai et al., 2011), dI2s, dI3s (both mechanosensors) (Bui et al., 2013; Haimson et al., 2020), and dI6 (gait) (Perry et al., 2019) are located in the deep layers of the dorsal horn, while dI4s and dI5s regulate pain, heat, and itch in the superficial layers of the dorsal horn (Koch et al., 2018). To replace damaged spinal tissue, it will be critical to derive a complete complement of dorsal sensory interneurons that molecularly and functionally phenocopy endogenous dIs.

Toward this goal, we are working to identify the developmental mechanisms that specify sensory dIs in the spinal cord and then apply these insights to develop protocols to derive dIs from stem cells (Andrews et al., 2017; Gupta et al., 2018, 2021). Distinct classes of dIs arise during embryonic development, when growth factors, including retinoic acid (RA) and the bone morphogenetic proteins (BMPs), pattern six progenitor domains (dP1–dP6), which differentiate into dI1–dI6 (Andrews et al., 2019; Gupta and Butler, 2021). The addition of RA and BMP4 in embryoid body (EB) protocols (Wichterle et al., 2002) can induce some dI fates in both mouse and human stem cell cultures (Andrews et al., 2017; Duval et al., 2019; Gupta et al., 2018). However, it has remained unclear how the addition of RA and BMP4 specifically induces dI fates, and whether the dIs generated from these protocols mirror their endogenous counterparts. BMPs pattern multiple organ systems throughout development, including the non-neuronal cardiac mesoderm (Kattman et al., 2011; Ladd et al., 1998), as well as osteogenic tissues (Kawaguchi et al., 2005). The BMPs also have reiterative activities both patterning dPs and controlling their proliferation to ensure that the precise number of specific dIs are generated (Andrews et al., 2017; Ille et al., 2007). Although the Smad second messenger complex mediates dI patterning (Hazen et al., 2012; Le Dreau et al., 2012), the downstream signals regulating proliferation remain unresolved.

Here, we define RA ± BMP4 protocols that generate mouse embryonic stem cell (mESC)-derived dIs via a neuromesodermal progenitor (NMP) intermediate (Gouti et al., 2014) that faithfully mimic the normal developmental programs of the neural tube (Delile et al., 2019; Rodrigo Albors et al., 2018; Sagner et al., 2018). This approach has successfully generated a complete atlas of in-vitro-derived dIs. Transcriptional profiling of these dIs demonstrated that they possess sensory function-specific gene signatures that strikingly resemble endogenous dIs, indicating that our protocols generate bona fide dIs. Using these protocols as an in vitro model for dorsal spinal cord development, we identified the hierarchy of fate specification decisions mediated by RA and BMP signaling. We further identify Wnt/β-catenin signaling as an immediate downstream response to BMP signaling that maintains neural progenitors in a mitotic state. Remarkably, elevating Wnt signaling, using the CHIR99021 (CHIR) agonist, can dramatically extend the proliferative capacity of dPs, while preserving their ability to differentiate into specific sensory neurons on demand.

Taken together, this study provides a mechanistic understanding of the developmental trajectories needed to generate the full complement of different dI identities. These insights were used to develop a protocol to expand the propagation of multipotential dPs, thereby paving the way for the production of in-vitro-derived dIs needed for both drug screening and cellular transplantation therapies.

RESULTS

Establishment of 2-d protocols to direct mESCs toward dorsal spinal interneuron fates

To direct dI fates through an NMP intermediate, we adapted a mESC protocol previously used to derive ventral spinal neurons (Gouti et al., 2014; Sagner et al., 2018). In this protocol, treatment with basic fibroblast growth factor (bFGF) and the Wnt agonist CHIR (Figure 1B) results in ~90% of mESCs acquiring a bipotential NMP identity by day 3 in culture (Figure S1A) distinguished by the co-expression of both mesodermal (brachyury (T)) and neural progenitor cell (NPC; Sox2) markers (Figure S1A).

Figure 1.

Distinct temporal combinations of RA and BMP4 direct different mESC identities

(A) Schematic of the source of patterning signals in the neural tube.

(B) Schematic of three RA ± BMP4 protocols used to direct mESCs toward dorsal spinal cord identities through an NMP intermediate.

(C–E) Protocols 2 and 3, but not protocol 1, induce phosphorylated Smad1/5/8 (p-Smad1/5/8) by day 5 showing active BMP signaling.

(F–K) Rosettes expressing Sox2 (NPCs, red, I, K), Pax3 (dPs, red, F, H), Pax6 (dPs, green, F, H), and Tuj1 (neurites, green, I, K) are observed in protocols 1 and 3, but not in protocol 2 (G, J), by day 9.

(L) Timeline for bulk RNA-seq sample acquisition (n = 3).

(M) Principal-component analysis (PCA) identified similar trajectories for protocols 1 (red, R-branch) and 3 (blue, B-branch), which are distinct from the trajectory of protocol 2 (green, C-branch).

(N) Pearson correlation analysis showing the cells from protocols 1 (R9) and 3 (B9) are transcriptionally similar and distinct from those in protocol 2 (C9) at day 9.

(O) Weighted gene co-expression network analysis (WGCNA) identified the successive downregulation of a rRNA biogenesis module implicated in the loss of pluripotency (Watanabe-Susaki et al., 2014; Woolnough et al., 2016). Conversely, neural modules are sequentially upregulated in the R- and B-branches, while a cardiac module is sequentially upregulated in the C-branch.

See Figure S1E for other modules. Scale bars: 100 μm (C–K).

Using the developmentally relevant growth factors (Figure 1A) previously shown to direct ESCs toward dorsal spinal fates (Andrews et al., 2017; Duval et al., 2019; Gupta et al., 2018; Wilson et al., 2004), we then treated day 3 NMPs with 100 nM RA for 48 h either alone (protocol 1, Figure 1B) or in combination with 10 ng/mL BMP4 starting at either day 3 (protocol 2) or day 4 (protocol 3). Treatment with RA alone did not activate BMP signaling (distinguished by the absence of phosphorylated [p] Smad1/5/8 staining; Figure 1C) and resulted in dorsalized neural cultures, with Pax3+ Pax6+ Sox2+ neural rosettes extending Tuj1+ neurites (arrowheads, Figures 1F and 1I). When BMP4 was added in combination with RA, BMP signaling was now activated (arrowheads, Figures 1D and 1E). However, the timing of BMP4 addition critically affected NMP identity. When BMP4 was added concomitantly with RA at day 3, no Pax3+ (Figure 1G) or Sox2+ (Figure 1J) cells were observed. The aggregates exhibited Tuj1 staining but had no apparent neuronal morphology (Figure 1J). In contrast, when RA was added at day 3, followed by BMP4 at day 4, we again observed Pax3+ Pax6+ Sox2+ neural rosettes, which robustly extended Tuj1+ neurites (arrowheads, Figures 1H and 1K).

Unbiased transcriptomic analyses identify the developmental trajectories leading to either spinal cord or cardiac fates

To assess the transcriptional changes resulting from the timed addition of RA ± BMP4, we conducted bulk RNA sequencing (RNA-seq) analyses along the timelines of the three protocols (Figure 1L). Principal-component (PC) analysis (PCA) (Son et al., 2018) revealed that PC1 and PC2 contribute 44% of the variation in the data and separate the least (D0, D3) from the most (D7, D9) differentiated samples (Figure S1B), while PC3 contributes 15% variance and separates the samples from the different protocols (Figure S1C). A three-dimensional PCA plot combining PC1-PC3 reveals three differentiation trajectories (Figure 1M), connected by bifurcation points at days 3 and 4. The C-branch (RA + BMP4, protocol 2) follows a distinct trajectory from the R-branch (RA, protocol 1) after bifurcating at day 3, while the B-branch (RA + BMP4, protocol 3) bifurcates from the R-branch at day 4 (Figure 1M). The R- and B-branches follow similar trajectories and terminate in adjacent PCA spaces, distinct from the C-branch, suggesting that their differentiation program is largely equivalent. Supporting this conclusion, a Pearson correlation demonstrated R- and B-branch transcriptomes were the most similar to each other (Figure 1N).

We next examined whether the three protocols direct distinct differentiation outcomes, by assessing the groups of genes upregulated over time using weighted gene co-expression and network analysis (WGCNA) (Langfelder and Horvath, 2008). Distinct eigengene modules (Farhadian et al., 2021; Panahi et al., 2020), i.e., a collection of genes with similar expression trends, were identified for each branch (Figure S1E). Of these, the yellow-green, dark gray, and green modules show increasing expression from day 0 to 9 for protocols 1, 3, and 2, respectively (Figure 1O). A functional enrichment analysis of these modules using Enrichr (Chen et al., 2013) revealed that both the yellow-green and dark-gray gene modules were highly enriched for nervous-system-specific Gene Ontology (GO) terms (Z score > 30) (Figure 1O), supporting the hypothesis that protocols 1 and 3 promote neural identities. In contrast, the top GO term in the green module is ‘‘outflow tract septum’’ (Figure 1O), a structure that regulates the blood flow in the heart (Webb et al., 2003). Other cardiac-related terms, such as actomyosin, contractile actin filament, and stress fibers terms, are represented in the cellular component categories (Figure S1D). We additionally observed that multiple cardiac-specific genes are upregulated by day 9 in cells derived from protocol 2 (Figure S1F). Taken together, these analyses reveal the timeline over which RA and BMPs act sequentially to direct NMPs toward the neural lineages, rather than endocardial fates.

Single-cell RNA-seq (scRNA-seq) of mESC-derived neurons identifies the full complement of spinal dIs

To assess whether protocols 1 and 3 direct distinct classes of dIs, we obtained single-cell transcriptomes from 9,704 cells using protocol 1 and 13,121 cells using protocol 3, taken at day 9 of the differentiation procedure (Figure 2A). We performed Louvain clustering followed by Uniform Manifold Approximation and Projection (UMAP) plots to visualize groups of transcriptionally distinct cell types (Becht et al., 2018). This pipeline identified 21 clusters of cells in both protocol 1 (Figure 2C) and protocol 3 (Figure 2D). In both conditions, ~30% of cells identify as being stressed; they cluster separately and display <2,000 RNA counts (Figures 2C, 2D, S2A, and S2B). Approximately 2%–4% of cells are pluripotent stem cells, marked by Sox2, Pou5f1 (Oct4), and Nanog expression, while 5%–8% of cells are non-neural, either of mesodermal or cardiac neural crest identity, expressing Twist1, Runx2, and Hand2 (Han et al., 2021; Soldatov et al., 2019; Vincentz et al., 2013) (Figures 2C, 2D, S2A, and S2B). The remaining 60% of cells have a neural identity, either Sox2 NPCs or Tubb3 (class III β-tubulin, Tuj1) differentiated neurons (feature plots, Figures 2C and 2D).

Figure 2.

Single-cell sequencing (sc-seq) identifies relevant spinal-function-specific modules in in-vitro-derived dIs

(A) Timeline for sc sequencing of protocols 1 and 3 (n = 1).

(B) Schematic showing the position of dPs and dIs in the spinal cord. dI1–dI6 mediate distinct functionalities and can be distinguished by distinct transcription factors.

(C and D) UMAP projection of day 9 cells derived from RA (C) and RA + BMP4 conditions (D). Feature plots show the Sox2+ NPCs and Tubb3 (Tuj1+) differentiated neurons.

(E) Protocol 1: sub-clustering of the Tubb3+ neuronal clusters (C) identified 14 clusters (3,394 cells). Feature plots show Pax2+Lhx1+ dI4/dI6, Lmx1b+ dI5, and Dmrt3+ dI6 populations. The pie chart depicts the percentage of dI1–dI6s obtained, only ~8.5% of neuronal cells of unknown identity (see Table S1). Concurring with the in vivo functionality (B), the dI4/dI6 clusters express Gad2 (inhibitory neurons), while the dI5 clusters express Slc17a6 (VGlut2, excitatory neurons).

(F) Protocol 3: sub-clustering of the Tuj1+ neuronal clusters shown in (D) identified 13 clusters (3,967 cells). Feature plots show Lhx2+Barhl1+ dI1, Foxd3+ dI2, and Isl1+ dI3 populations. The pie chart represents percentage of dI1–dI6 clusters; ~9.5% of cells are unknown neural identity (see Table S1, cluster 3). As in in vivo (B), the dI1–dI3 clusters express excitatory markers, Slc17a2 and Grin2b.

(G) Gene Ontology (GO) analysis of the genes upregulated (logFC > 0.2) in the indicated clusters reveal that the relevant sensory/motor signatures are enriched in the in-vitro-derived dIs.

To distinguish the neuronal subtypes, we subclustered the Tubb3ells, resulting in 14 (protocol 1) and 13 (protocol 2) transcriptionally distinct clusters (Figures 2E and 2F). We then used a well-characterized panel of transcription factors (Figure 2B) to assign dI identities to these clusters (Gupta and Butler, 2021; Lai et al., 2016). This analysis revealed that protocol 1 primarily generates dI4, dI5, and dI6 that mediate pain, itch, and heat perception, while protocol 3 most notably generates dI1, dI2, and dI3 that regulate proprioception, gait, and mechanosensation. Specifically, in protocol 1, the majority of subclustered Tubb3 cells divide into two major groups of closely aligned clusters (Figures 2B, 2E, and S2C): one group (clusters 7, 2, 0, and 1) expresses Lmx1b, Prrxl1, and Tlx1/3, which define dI5 identity (Lai et al., 2016); the second group (clusters 6, 12, 13, 4, 5, and 3) co-expresses Pax2 and Lhx1, which denote dI4/dI6 identity. Of these clusters, only cluster 6 also expresses Dmrt3, an established marker of dI6s (Andersson et al., 2012). The neurotransmitter profile supports these fate designations: phenocopying endogenous neurons (Figure 2B), the dI4/dI6 clusters express Gad2, a GABA-synthesizing enzyme (Pillai et al., 2007), while the dI5 clusters express Slc17a6 (VGlut2), a glutamate transporter specific to excitatory neurons (Figure 2E) (Cheng et al., 2004). Of the remaining four clusters, clusters 8 and 9 map to dI1 and dI2 identities (Figure S2C), while cluster 10 is enriched for neural- and neuroendocrine-specific genes, and cluster 11 is enriched for ribosomal genes (Table S1).

In protocol 3, the subclustered Tubb3 cells divide into three groups of closely aligned clusters (Figures 2B, 2F, and S2D): clusters 7, 1, 0, and 8 express Lhx2/9 and Barhl1/2, denoting dI1 identity; clusters 11, 5, and 6 express the dI2 markers Foxd3, Lhx1, and Lhx5, while clusters 12 and 10 express Isl1 and Tlx3, indicating dI3s. As in vivo, all three classes of in-vitro-derived neurons are excitatory, expressing either Slc17a6 or Grin2b (Figure 2F). Of the remaining four clusters, clusters 2 and 9 express Atoh1 and Neurog1, respectively, which mark the dP1 and dP2 state and suggest these cells are immature dI1s and dI2s (Figure S2D); cluster 4 expresses dI4-specific markers (Figure S2D), and cluster 3 expresses multiple ribosomal and mitochondrial genes indicative of dying neurons (Ilicic et al., 2016) (Table S1).

In summary, protocol 1 (RA) generates ~33% dI4s, 40.5% dI5s, 8% dI6s, 5% dI1s, 5% dI2s, and 8.5% unknown neural cell types (pie chart, Figure 2E). In contrast, protocol 3 (RA + BMP4) generates ~45.5% dI1s, 29% dI2s, 8% dI3s, 8% dI4s, and 9.5% unknown cell types (pie chart, Figure 2F). Both immunohistochemical and qRT-PCR analyses further support these designations (Figures S2E–S2H).

GO analyses of ESC-derived dIs identify sensory modality-specific ontology modules

We next assessed whether the functional identities of these in-vitro-derived neurons could be predicted from the transcriptomic data. In vivo, the dIs relay diverse somatosensory modalities either locally within the spinal cord (Fernandes et al., 2016; Todd, 2017) or by long-range afferent connections (Wercberger and Basbaum, 2019). Generally, the most dorsal populations regulate proprioception (dI1/dI2) (Sakai et al., 2012; Yuengert et al., 2015) and touch-induced mechanosensation (dI3s) (Bui et al., 2013), while the more intermediate classes form relay circuits for pain, itch, heat, and touch (dI4/dI5) (Koch et al., 2018; Todd, 2010). The dI6s regulate coordinated movement by forming inhibitory synapses with spinal motor neurons, while themselves receiving cholinergic and glutamatergic inputs (Perry et al., 2019).

To identify functionally relevant GOs, we selected the cluster for each dI population that highly expressed the differentiation markers (>0.2 log fold change [FC]) and subjected them to the DAVID pipeline (Huang et al., 2009a, 2009b), which identifies enriched biological processes and signaling pathways related to diseases or drugs (Figure 2G). Remarkably, each class of ESC-derived dIs displayed the relevant modality-specific signatures. For the dorsal-most neurons, we analyzed clusters 0 (dI1), 5 (dI2), and 10 (dI3) and specifically found enriched GO categories related to balance and walking behaviors, as well as excitatory synaptic signaling. Both dI1 and dI2 clusters additionally include proprioception-related terms (blue terms, Figure 2G; Table S2). For the more intermediate neurons, we analyzed clusters 4 (dI4), 1 (dI5), and 6 (dI6). The dI4/dI5 clusters include GO terms associated with pain and itch perception, specifically implicating both oxytocin, a neuropeptide that modulates pain processing (Boll et al., 2018), and β-alanine, an amino acid that induces itch (Liu et al., 2012). As in vivo, the dI4 cluster contains terms related to inhibitory synaptic signaling, while the dI5 cluster is enriched for excitatory synaptic terms (red terms, Figure 2G; Table S2). In contrast, the dI6 cluster contains terms associated with coordinated movement and synaptic signaling pathways, accurately reflecting their functional identity. Finally, we also found enriched terms in the dI clusters related to different drug addiction pathways, including amphetamine and endocannabinoids for dI1/dI2s, and amphetamine, morphine, and cocaine for dI4/dI5s (Figures 2G and S3A; Table S2). Such signatures suggest the dIs mediate the ability of psychoactive drugs to modulate pain and itch perception (Lipman and Yosipovitch, 2021) or proprioception (Downey et al., 2017). Other dI-subtype markers identified in this analysis include receptors, ion channels, and adhesion molecules (Figures S3A and S3B; Table S4).

ESC-derived dIs are transcriptionally indistinguishable from endogenous spinal interneurons

To assess whether the ESC-derived dIs are also transcriptionally similar to their in vivo counterparts, we used reciprocal PCA-based data integration and label transfer (Hao et al., 2021) to compare a single-cell dataset taken from embryonic stage (E) 9.5–E13.5 mouse spinal cords (Delile et al., 2019) with the in-vitro-derived cell types in protocols 1 and 3. Projecting the combined in vivo and in vitro datasets into UMAP plots revealed a remarkable degree of overlap (Figures 3A–3C). This overlap was observed only for the neural progenitor and neural identities; the blood, mesoderm, neural crest, and skin lineages present in the in vivo dataset were notably absent from the in vitro datasets (Figure 3A). Although cells from protocols 1 and 3 mapped directly on top of the spinal cord cells (Figure 3B), these three datasets were distinct from trachea, lung, and kidney cells (Tabula Muris Consortium et al., 2018) (Figure 3E).

Figure 3.

Stem-cell-derived sensory interneurons resemble their endogenous counterparts

(A) Integrated UMAP plot of the sc-seq datasets from the in vitro differentiations and the in vivo embryonic spinal cord (E9.5–E13.5) (Delile et al., 2019). Cell-type labels from the in vivo dataset were projected onto the in vitro dataset.

(B) Integrated UMAP plot highlighting the overlap between the in vitro and in vivo cells. An atlas of these combined datasets with dI annotation can be found here (https://www.dropbox.com/s/kau8oc7t6i5lchl/dImap_Briscoe_in vitro combined_3D_interactive.html?dl=0; file opens after being downloaded).

(C) Clustering of the integrated dataset yields 24 shared clusters, with 6 clusters unique to the in vivo dataset.

(D) Dot plot of the 24 shared clusters showing the expression levels of dI marker genes.

(E) Integration of the in vivo spinal cord and in vitro datasets with lung, kidney, and trachea datasets (Tabula Muris Consortium) shows the in vitro cells overlap only with spinal cord cells.

(F) UMAP plots of the neuronal portion of the integrated datasets. The in vivo dIs are colored using the in vivo dataset annotations. Each dI population was then overlaid onto the in vitro datasets (black) from either protocol 1 (RA) or 3 (RA + BMP4). Substantial overlap was seen for dI1–dI6, while no overlap was observed for motor neurons (pink).

(G) dI-specific cells were extracted from both the in vitro protocols and in vivo spinal cord dataset and then plotted in the integrated UMAP space with in vitro cells on top of in vivo cells.

To compare the distribution of cell types in the in vivo and in vitro datasets, we performed unsupervised clustering on an integrated dataset, yielding a total of 29 clusters (Figure 3C); 24 clusters (0–10 and 12–25) were common to all datasets, while 6 clusters were unique to the spinal cord dataset. A dot plot analysis of dI-specific markers within the shared 24 clusters showed striking similarities between the in-vivo- and in-vitro-derived datasets (Figure 3D). Although there are modest differences, there is overall concordance between both the number of cells expressing a given dI marker and the strength of marker expression in the in vivo and in vitro cell clusters. This concordance continues as the neurons mature, demonstrated by the expression of neuropeptides, which modulate the intensity of sensory inputs in the adult spinal cord (Lai et al., 2016; Russ et al., 2021). Specific neuropeptides, including prepronociceptin (pnoc; dI2, dI4), neuromedin S (nms; dI4), neuropeptide Y (npy; dI4), and gastrin-releasing peptide (grp; dI1, dI5), are enriched in the same populations of dIs in both the in vitro and in vivo datasets (Figures S3C and S3D). Thus, the ESC-derived dIs appear to proceed along the same maturation process as endogenous spinal neurons.

Finally, to assess whether the in-vitro-derived dI populations specifically mirror their endogenous counterparts, we compared the protocol 1 versus protocol 3 datasets with an atlas of the neuronal classes of the spinal cord (Figure 3F), generated using the annotation provided by Delile et al. (2019). The in vivo dI1–dI3s and dI4–dI6s map to identical positions in protocol 3 (RA + BMP4) and protocol 1 (RA), respectively (Figure 3F). In contrast, there was no overlap between the ventral spinal neurons, such as the motor neurons (MNs), with the in vitro datasets (Figure 3F and data not shown). We also computationally extracted the cells expressing specific dI subtype markers from the in vitro and in vivo datasets and projected them in the same UMAP space (Figure 3G). This analysis further demonstrated that all six classes of in-vitro-derived dIs are transcriptionally indistinguishable from their spinal cord-derived counterparts.

Taken together, these studies suggest that our in vitro protocols are accurately replicating the in vivo developmental program of the spinal cord permitting the generation of bona fide sensory interneurons.

RA activates a transcriptional network that specifies a dI4–dI6 default state

Using the directed differentiation protocols, we investigated the unresolved mechanisms by which RA and BMP directdI fate specification. We first assessed the genetic basis of RA-mediated fate decisions. During the 24-h pulse of RA at day 3, Pax3 levels gradually increase, whereas T expression concomitantly declines (Figure S5A). By day 4, ~90%–95% of cells are Pax3+ (Figure S5B), indicating NMPs have transitioned to a dorsal neural progenitor state (dP state1) (Figure S4A). To identify the transcriptional changes occurring during this first commitment step, we performed bulk RNA-seq at day 3 (NMP state), day 3.25 (early transcriptional response to RA), and day 4 (dP state 1) (Figure S4B).

The differential expression analysis identified 138 upregulated and 19 downregulated genes after 6 h of RA treatment (Figure S4C). Upregulated genes include known RA-regulated factors, such as Meis1, Meis2, Stra8, Cyp26a1, and Rarb, confirming that cells are directly responding to RA. By 24 h of RA exposure, 389 genes were significantly upregulated, including genes present in intermediate dPs and the dI4/dI5s (Figure S4D), while 590 genes were significantly downregulated, which include mesodermal-specific genes (Figures S4E–S4G). These findings support the hypothesis that RA induces NMPs toward intermediate spinal identities, while also directing them away from mesodermal fates, thereby ensuring the transition from an NMP to Pax3+ dP state 1.

We next assessed whether the dP state 1 can be defined as a function-specific interaction network, by subjecting the upregulated genes to a Metascape protein-protein interaction (PPI) analysis (Zhou et al., 2019). At 6 h, a nascent PPI network exists, with one functional module for embryonic patterning containing the RA-target genes Meis1 and Meis2 (Figure S4H, inset). However, by 24 h, the PPI network has expanded into 13 functional modules (Figures S4I and S5C). Four modules, insulin-like signaling, pattern specification, sensory organ morphogenesis, and peptide ligand receptor, contained genes expressed in the intermediate spinal cord (insets, Figures S4I and S5C). Notably, Meis1 and Meis2 interact with Pax6, Lmx1b, Msx1, and Msx2 in the pattern specification module. Using both qRT-PCR (Figure S5D) and transcriptional profiling of untreated day 4 cells (no RA control), or after a 24-h pulse of RA ± AGN193109, a pan-RAR inhibitor (Figure S4J), we independently validated that Meis1 and Meis2, but not Meis3, are upregulated by RA, including many patterning genes. Thus, Meis1/2 are good candidates for the key transcription factors that link RA signaling to the initiation of spinal cord patterning, i.e., dP state 1.

Wnt signaling is induced as an immediate consequence of BMP4 treatment

We hypothesized that the addition of BMP4 to NPCs in dP state 1 directs them to a more dorsal dP state, i.e., dP state 2 (Figure S4A). To define the transcriptional changes occurring after BMP4 addition, we performed bulk RNA-seq at day 4 (dP state 1), day 4.25 (early transcriptional response to RA + BMP4), and day 5 (dP state 2) (Figure 4A). We identified 106 genes upregulated by 6 h and 305 genes upregulated after 24 h of BMP4 exposure (Figure S6A). Seventy-eight genes, out of 411, were common between both time points (Figures 4B and S6B). The GO analysis of these 78 common genes included the BMP signaling pathway (Figure 4C; Table S3), confirming that the analysis identified BMP-responsive genes.

Figure 4.

Wnts are upregulated as an immediate response to BMP4 signaling in mESC-derived NPCs

(A) Timelines for bulk RNA-seq analysis of protocol 3 (n = 3).

(B) Venn diagram showing the overlap of genes significantly upregulated after 6 and 24 h of RA + BMP4 treatment (logFC > 2, false discovery rate [FDR] < 0.01). 78 genes were common to both comparisons.

(C) GO analysis of the 78 common genes shows the upregulation of the Wnt signaling pathway in all three GO categories: biological processes, molecular processes, and pathway categories.

(D) Heatmap showing the expression (FPKM, N = 3) of 19 Wnt ligands at day 4 (dP state1), day 4.25 (immediate response to BMP4), and day 5 (dP state2).

(E) Volcano plots showing upregulated genes at 6 and 24 h of RA + BMP4 treatment. Multiple Wnt ligands and Wnt pathway genes are highly upregulated in both conditions.

(F) qRT-PCR validation of selected Wnt genes under RA, RA + BMP4, and RA + BMP4 + noggin conditions at day 5. Wnt expression was significantly upregulated by BMP4 but reduced in the presence of noggin. n = 3 differentiations. The data are presented as the mean ± SEM. Significance is determined by two-way ANOVA (Tukey’s multiple comparison test), *p < 0.05, **p < 0 .005, ***p < 0.0005.

In addition, we also found enriched GO terms related to Wnt signaling across three functional categories (Figure 4C). Both Wnt ligands, including Wnt1, Wnt2, Wnt3, Wnt3a, and Wnt4, and Wnt receptors, frizzled (Fzd) 8, Fzd3, Fzd10, and Lrp6 (Figure S6C), were upregulated in response to BMP4 treatment, along with Wnt signaling regulators, such as Dkk, Lmx1a (Hoekstra et al., 2013), and Bambi (Zhao et al., 2020) (Figures 4D, 4E, and S6B). Selected Wnt-signaling genes were further assessed in a qPCR analysis of day 5 NMPs to validate whether they are BMP4 responsive. We found that expression of all tested genes significantly increased after the addition of BMP4, compared with the RA alone condition (Figure 4F). The upregulation of Wnt expression was suppressed if Noggin, a BMP inhibitor, was concomitantly added with BMP4 on day 4 (Figure 4F).

Many of the upregulated Wnts, Wnt1, Wnt3a, and Wnt4, are expressed in the developing dorsal spinal cord in vivo (Agalliu et al., 2009; Daneman et al., 2009) and signal through the canonical Wnt signaling pathway (MacDonald et al., 2009), suggesting that canonical Wnt signaling mediates the ability of BMP4 to direct the dP state 1-to-dP state 2 transition.

Inhibiting canonical Wnt signaling blocks the induction of the dorsal-most fates

To assess the role of canonical Wnt signaling in dI fate specification, we added endo-IWR1 (IWR1e), a small-molecule inhibitor that blocks β-catenin function (Chen et al., 2009) to dPs from day 4 to 6 (Figure 5A). By day 9, the RA + BMP4 condition robustly induces dI1s–dI3s (Figures 5F–5I and 5N–5P). In contrast, the addition of IWR1e resulted in the loss of dI1s and dI2s (Figures 5J, 5K, 5N, and 5O), while dI3 differentiation was less affected (Figures 5L and 5P).

Figure 5.

Wnt/β-catenin signaling is required for BMP4-mediated neuronal diversity

(A) Timeline to assess the requirement for Wnt signaling in the induction of dI fates.

(B–Q) Protocol 1 induces Pax2+ dI4/dI6s (E and Q) and Lmx1b+ dI5s (E), while protocol 3 induces the Lhx2+ dI1s (F and N), Foxd3+ dI2s (G and O), and Isl1+ dI3s (H and P). The addition of IWR1e to protocol 3 dramatically reduces the number of dI1s and dI2s (J, K, N, and O), more modestly decreases the dI3s (L and P), and concomitantly increases the dI4–dI6s (M and Q). n = 3 differentiations for qRT-PCR.

Significance tests: one-way ANOVA, *p < 0.05, **p < 0 .005, ***p < 0.0005. Scale bar: 100 μm.

BMP4 directs the dorsal-most fates (dI1–dI3), while concomitantly suppressing RA-dependent intermediate dorsal fates (dI4–dI6; Figures 5E, 5I, S6D, and S6E). Supporting a role for the Wnts mediating this latter activity, there is no suppression of Pax2+ dI4/dI6 and Lmx1b+ dI5 fates when IWR1e is added to the dP cultures, along with RA + BMP4 (Figures 5M and 5Q). The addition of IWR1e did not erode Pax3+ pan-dorsal identity or Pax7 expression, but did suppress Olig3 (Figure S6F), which may contribute to the loss of dI1/dI2 and perdurance of the dI4–dI6 fates.

Together, these results demonstrate that the BMP4-dependent dI fates require canonical Wnt/β-catenin signaling. In its absence, BMP4-mediated suppression is lost, and dPs retain their intermediate spinal cord identity.

Wnt signaling promotes BMP4-dependent fates by regulating proliferation, not patterning

Previous studies have suggested that Wnts act either as mitogens (Dickinson et al., 1994; Megason and McMahon, 2002) or patterning factors (Alvarez-Medina et al., 2008; Muroyama et al., 2002) in the specification of the dorsal spinal cord. We assessed these models, first by determining the effect Wnts (Wnt1, Wnt2, of canonical Wnt3a, and Wnt4), a non-canonical Wnt (Wnt9b), and a small-molecule Wnt agonist (CHIR99021) on dI identity in the RA ± BMP4 protocols. When Wnts were added together with RA, during the period when NPCs are competent to respond to patterning factors, they were unable to induce the dorsal-most dIs (Figure S7A). Similarly, when the Wnts were added concomitantly with BMP4, there was no synergistic effect on the specification of dorsal spinal identity (Figures 6A and 6B). Thus, Wnts are not sufficient to direct dI fate identity.

Figure 6.

Wnt/β-catenin signaling mediates proliferation, not patterning, of mESC-derived spinal NPCs

(A) Timeline to assess whether Wnt signaling modulates the BMP-mediated dI fates through patterning activities.

(B) qRT-PCR analyses for dI1–dI3 marker genes show that neither Wnts (100 ng/mL) nor CHIR (5 μM) affect (one-way ANOVA) the BMP4-mediated dI fates. Expression levels were normalized to the day 0 and RA condition (n = 2).

(C) Embryoid body (EB) differentiation timeline to assess the effect of activating Wnt signaling on the NPC proliferation. Samples were collected at day 9 for immunohistochemistry (IHC) and qRT-PCR.

(D–I) Both EdU incorporation (red, D–G; S-phase) and phosphorylated histone H3 staining (blue, D–G; pH3, M-phase) show that CHIR treatment increases cell proliferation by ~3-fold in RA ± BMP4-patterned EBs, compared with DMSO control, resulting in increased numbers of Sox2+ NPCs (green, D–G). Significance was determined using one-way ANOVA (Kruskal-Wallis test) (n = 15–25 EBs).

(J–P) IHC (J–M) and qRT-PCR (N–P) for Pax3 (red, J–M; all dPs), Olig3 (green, J–M; dP1–dP3), and Pax7 (blue, J; dP4–dP6) demonstrated that CHIR treatment increases the numbers of NPCs but does not significantly (unpaired t test) alter their dorsal-ventral identities (n = 3 differentiations for qRT-PCR).

Significance values: *p < 0.05, **p < 0 .005, ***p < 0.0005. Scale bar: 200 μm.

We next asked whether activating Wnt signaling affected the mitotic index of mESC-derived NPCs after the patterning period. To better visualize proliferation, we allowed mESCs to form EBs (Andrews et al., 2017; Duval et al., 2019; Gupta et al., 2018) and then treated them with RA ± BMP4 over the same timeline as protocols 1 and 3 (Figure 1B) to induce dorsal spinal cord patterning. EBs were then cultured with either DMSO or CHIR from day 5 to 9 (Figure 6C) and finally pulsed with 5-ethynyl-2’-deoxyuridine (EdU) to label S-phase cells. Strikingly, there was an ~3-fold increase in EdU+ cells in CHIR-treated EBs (Figures 6E, 6G, and 6H) compared with DMSO control (Figures 6D, 6F, and 6H) in both protocols. Similarly, the CHIR-treated EBs had an ~2- to 4-fold increase in the number of pHistone H3 (pH3)+ cells in M-phase, compared with control (Figures 6D–6G and 6I). Many of the dividing cells in the CHIR-treated cultures were Sox2+ NPCs (Figures 6E and 6G), while NPCs are largely depleted by day 9 in control cultures (Figures 6D and 6F). Thus, extending the window of Wnt signaling prolongs NPC proliferation without affecting their original identity, i.e., NPCs from protocol 1 expressed Pax3 and Pax7, while NPCs from protocol 3 expressed Pax3 and Olig3 (Figures 6J–6P).

Together, these studies support the hypothesis that Wnt/β-catenin signaling functions as a mitogen for dorsal spinal NPCs, rather than specifying distinct dI identities.

Wnt/β-catenin signaling can be modulated to extend the timeline of dP patterning

Finally, we assessed whether the mitogenic properties of Wnt signaling can be leveraged to expand the pool of mESC-derived NPCs. RA ± BMP4-patterned EBs were serially passaged at 1:3 dilution every 7 days in the presence of DMSO (control) or CHIR (Figure 7A). Although the number of EBs was comparable in the DMSO and CHIR conditions on the first passage, they declined sharply by the second passage in the controls (Figures 7B, 7C, 7E,7F, 7G, 7I, S7B, and S7C). By the third passage, few EBs were observed in the DMSO condition, while the number of CHIR-treated EBs had increased by ~3- to 6-fold (Figures 7E, 7I, S7B, and S7C). The ability of CHIR to expand the number of EBs was observed in both RA ± BMP4 conditions but was most robust for the RA + BMP4-patterned EBs (Figures 7E and 7I). CHIR-treated EBs remained Sox2+ after multiple passages, suggesting that they maintained their NPC identity (Figures 7D, 7H, S7D, and S7E).

Figure 7.

Activating canonical Wnt signaling expands mESC-derived NPC populations

(A) Timeline for the expansion protocol for mESC-derived spinal cord progenitors.

(B–I) RA-treated (B–E) and RA + BMP4-treated (F–I) EBs were passaged with DMSO (B and F) or CHIR (C, D, G, and H). Bright-field images (B, C, F, and G) show that passaging with CHIR, but not DMSO, causes >6-fold increase in the number of EBs (E and I; one-way ANOVA Kruskal-Wallis test, n = 2 independent differentiations). All CHIR-treated EBs maintain a Sox2+ NPC identity after passages 1 and 2 (D and H), but the expansion of RA-patterned EBs was less robust than that of RA + BMP4-patterned EBs (E and I).

(J) Schematic of the timeline and experimental approaches to assay the differentiation potential of CHIR-treated EBs. The left timeline shows conditions to induce neural differentiation in CHIR-expanded EBs following their dissociation. The right timeline determines the extent to which dI patterning can be restored in CHIR-treated EBs, with a 2-day pulse of RA ± BMP4.

(K–Q) IHC analysis for Sox2+ NPCs and Tuj1+ neurites on passage 6 EBs (48 days) revealed that removal of CHIR induces spontaneous neural differentiation in both RA-treated (L) and RA + BMP4-treated (O) EBs. The addition of DAPT had a modest effect on the number of Tuj1+ neurites (M and P). qRT-PCR analyses confirm the upregulation of Tubb3 and NeuN in CHIR-depleted conditions. (Q) However, the number of Sox2+ cells remains unchanged.

(R–X) Assessing the patterning profile of expanded EBs showed they can differentiate into dI5s (R and X; RA control) or dI3s (U; RA + BMP4 control), but not dI1s, dI2s, and dI4s (R, T, and X). A 2-day pulse of RA ± BMP4 from day 3 to 5 restored dI1/2 differentiation (V and X) with an ~5-fold increase in the number of Pax2+-dI4/dI6 neurons (S and X). The differentiation potential of passage (p) 1 and p6 (X) EBs did not show significant change over time, except for dI2s, where differentiation potential modestly improves over time (n = 6–10 EBs, Mann-Whitney test). Significance values: *p < 0.05, **p < 0.005, ***p < 0.0005. Scale bars: 2 mm (B, C, F, and G); 100 μm (Q–V); 2,000 μm (J–O).

To examine whether expanded NPCs can differentiate, we dissociated RA ± BMP4 CHIR-treated EBs at passages 6 (Figure 7J) or 7 (data not shown) onto Matrigel-coated plates and cultured them in the presence of CHIR, DMSO (i.e., CHIR withdrawal), or DAPT, a pro-differentiation agent (Crawford and Roelink, 2007) (Figure 7J). By day 10, while the CHIR-treated cultures remained Sox2+ (Figures 7K, 7N, and 7Q), the DMSO-treated cultures show robust induction of post-mitotic neural markers, neuronal nuclei (NeuN/Rbfox3), and Tubb3 (Figure 7Q) and contained many Tuj1+ neurites (Figures 7L and 7O). Thus, the removal of CHIR triggers spontaneous neural differentiation. The addition of DAPT did not significantly augment NeuN and Tubb3 expression compared with DMSO (Figure 7Q), although there was increased sprouting of Tuj+ neurites (Figures 7M and 7P).

We next asked whether specific dI fates are preserved in the expanded NPC pool. By day 10 (timeline in Figure 7J), passage 6/7 cultures derived from RA-patterned EBs contained many Lmx1b+ dI5s, but only a few Pax2+ dI4/dI6s, while cultures from the RA + BMP4-patterned EBs contained both Isl1+ dI3s and Lhx1b+ dI5s, but almost no Lhx2+ dI1s or Foxd3+ dI2s (Figures S7I, S7J, 7R, 7T, and 7U). Thus, although some patterning information is preserved over many passages, prolonged CHIR treatment erodes its fidelity, with RA + BMP4 dPs now expressing high levels of Pax7 and low levels of Olig3, as if they have been returned to an intermediate spinal identity (Figures S7G and S7H). However, dorsal patterning information can be restored by including a pulse of RA ± BMP4 between days 3 and 5 during the differentiation of the expanded EBs on Matrigel (Figure 7J). The addition of RA resulted in an ~6-fold increase in the number of dI4/dI6s (Figures 7S and 7X), while the addition of RA + BMP4 induced ~6-fold more dI1s and ~2-fold dI2s (Figures 7V and 7X), compared with control (Figures 7R, 7T, and 7X). Unexpectedly, the pulse of RA + BMP4 suppressed the dI3 fate (Figures 7U, 7W, and 7X). This strategy was successful at restoring almost all dI identities in cultures taken from EBs both at the earliest passage and at the latest passage assessed (passage 7; Figure 7X).

Thus, these data show that the chemical activation of the canonical Wnt signaling can be successfully used to expand the size of the mESC-derived NPC population, and thereby to derive large numbers of specific dI populations.

DISCUSSION

Stem-cell-based directed differentiation protocols result in both a source of in-vitro-derived cell types and a model system to dissect developmental mechanisms. We have built on previous studies that derive spinal motor neurons in vitro from an NMP intermediate (Gouti et al., 2014; Sagner et al., 2018), to develop protocols that produce a complete in vitro atlas of dorsal spinal interneurons. Our transcriptomic analyses suggest that these neurons are indistinguishable, both molecularly and functionally, from their endogenous counterparts. We then used these protocols to investigate the mechanisms that regulate dI fate choices. These studies identified the hierarchy of dI fate specification and resolved the role of Wnt signaling as a mitogen in this process. We leveraged this mechanistic understanding to expand the numbers of dPs as a means of increasing dI yield.

Hierarchy of dI fate specification

Previous studies have shown that NMPs contribute to spinal cord and paraxial mesodermal lineages both in vivo (Attardi et al., 2018; Tzouanacou et al., 2009) and in vitro (Gouti et al., 2014). Our transcriptomic analyses revealed that the dI differentiation program conforms to the canonical Waddington model (Rajagopal and Stanger, 2016; Waddington, 1957), where NMPs are at the top of a ‘‘hill’’ of possible fates, and the time and sequence in which they receive RA and BMP4 signals control their differentiation into dIs versus mesodermal fates. Thus, the addition of RA at day 3 channels NMPs toward the intermediate dI fates. However, if BMP4 is added concomitantly with RA, the differentiation trajectory alters to specify endocardial lineages. BMP4 can induce only the dorsal-most dI fates after 24 h of RA treatment, suggesting that NMPs must move through an RA-induced early-competency state (dP state 1) before BMP4 can relay dorsal spinal patterning information. These studies thereby inadvertently developed a directed differentiation protocol (protocol 2) that is a starting point to derive the cardiac cell types needed to regenerate the heart (Duelen and Sampaolesi, 2017). BMP4 may direct NMPs to form lateral plate mesoderm (Row et al., 2018), which contributes multiple cell types to the embryonic heart. This finding also suggests that NMPs may have more pluripotency for mesodermal fates than previously realized. In this case, the use of alternative growth factors in protocol 2 will direct NMPs toward additional mesodermal derivatives.

Our studies suggest that RA is required to establish dP state 1, before specifying the dI4–dI6 fates. The Meis1/2 genes, which are dynamically expressed in the early spinal cord neuroepithelium (Sánchez-Guardado et al., 2011), may be central to this competency state. Meis1/2 are among the earliest genes to be induced by RA (Gouti et al., 2017; Oulad-Abdelghani et al., 1997), and they become central to a gene interaction network that contains the genes that specify intermediate spinal identity. Thus, the upregulation of Meis1/2 may drive the multipotential progenitor state (dP state 1) that differentiates into dI4–dI6. The addition of BMP4 both suppresses dI4–dI6 fates and induces another multipotent state, dP state2, that differentiates into dI1–dI3. Interestingly, dP state1 may be a default state, which was first suggested by earlier studies in the intermediate spinal cord (Diez Del Corral et al., 2003; Novitch et al., 2003). When RA + BMP4-patterned progenitors are expanded, they appear to revert to the default dP state1 identity. This reversion can in turn be reversed: a pulse of RA + BMP4 remains sufficient to convert dP state 1 back to dP state 2 and elicit formation of the dorsal-most dIs.

Extent to which in-vitro-derived dIs phenocopy endogenous sensory interneurons

Our studies significantly extend previous studies seeking to generate dIs (Andrews et al., 2017; Duval et al., 2019; Gupta et al., 2018) by directly assessing the heterogeneity of the cultures and documenting the extent to which stem-cell-derived dIs resemble their endogenous counterparts. Using scRNA-seq analyses, we identified that our protocols generate the full complement of dIs, that there is no heterogeneity beyond dorsal spinal derivatives in the Tubb3 cells, and that the transcriptomes of in-vitro- and in-vivo-derived dIs are indistinguishable. In-vitro-derived dIs express the correct patterning factors and neurotransmitters, and display the relevant functional GO terms, suggesting they encode the correct sensory modalities. For example, balance-related GO terms were enriched in dI1, dI2, dI3, and dI6, which regulate proprioception and gait, while pain and itch GO terms were enriched for dI4 and dI5, which relay noxious stimuli in vivo. We also find GO terms related to serotonergic and dopaminergic synapses, which may permit communication between the spinal cord and enteric nervous system. This signature is associated with the synapses formed between spinal neurons in the dorsal horn and the descending afferents from the brainstem (Schwaller et al., 2017) and diencephalon (Skagerberg and Lindvall, 1985), which regulate both pain perception (Gautier et al., 2017) and sensation in gut (Travagli et al., 2006) and bladder (Hou et al., 2021). These circuits are of considerable therapeutic importance, given that spinal injuries can lead to both irritable bowel syndrome (Holmes and Blanke, 2019) and loss of bladder control (Hou et al., 2021; Taweel and Seyam, 2015). In addition, we identified gene signatures related to addictive psychoactive compounds like morphine, cocaine, and amphetamine, suggesting the dIs are the cellular targets for these drugs, elevating their importance for drug screening platforms.

Role and application of Wnt signaling as a mitogen

The BMPs have reiterative roles regulating dI fate specification and modulating the cell cycle (Andrews et al., 2017; Ille et al., 2007). Here, we demonstrate that the control of proliferation is mediated through the activation of Wnt signaling. Like the BMPs, multiple Wnt genes, especially those that activate the canonical Wnt/β-catenin pathway, are expressed in the developing dorsal spinal cord (Agalliu et al., 2009). The BMP and Wnt pathways are known to interact in other contexts, including during neural crest delamination (Burstyn-Cohen et al., 2004) and the assignment of hematopoietic fates (Lengerke et al., 2008). However, the role of the Wnts in the spinal cord has remained unresolved. Previous in vivo studies ascribed both dI fate specification (Muroyama et al., 2002; Zechner et al., 2007) and mitogenic roles (Ille et al., 2007; Megason and McMahon, 2002) to the dorsal Wnts. Our data unambiguously showed that the Wnts tested so far are not sufficient to pattern dI fates, rather they function as mitogens in this context, perhaps controlling the size of the dorsal-most dI populations. It remains unresolved as to why there are so many Wnts present in the spinal cord in vivo. One possibility is that they have signal-specific activities, like the BMPs (Andrews et al., 2017), sculpting the number of cells in each dI population. This mechanistic insight enabled us to expand mESC-derived spinal cultures by adding CHIR to constitutively activate Wnt/β-catenin signaling after the patterning stages. Remarkably, this manipulation sustained dorsal spinal neural identity while enabling the dP to proliferate for more than eight passages (the maximum number we attempted).

This protocol modification may represent a key step toward the long-sought goal of being able to supply specific neural populations in a limitless manner. Previous neurosphere studies have suggested that immortalized cells lose their neurogenic potential over time (Olivos-Cisneros et al., 2021; Shen et al., 2006). We rather find, in this expansion protocol, that even the progenitors from the latest passages we tried (passages 5–7) retain the ability to differentiate as neurons. However, while neurogenesis remains stable, the patterning information does erode over time. Cell-cycle duration has been shown to regulate cell fate by acting as a filter for long transcripts (Abou Chakra et al., 2021; Singh et al., 2013). Such mechanisms may erode cellular identity (i.e., dP state 2) when NPCs are held in an extended proliferative phase. Nonetheless, dI fates can be again restored with a 2-day pulse of RA ± BMP4.

Implications of stem-cell-derived sensory neurons for drug screening and regenerative therapies

Our studies have identified distinct psychoactive and analgesic drug-related signaling pathways in different dI populations. These findings raise the possibility of using stem-cell-derived dIs in drug screening studies to identify additional analgesics and to study the mechanistic action of drugs such as cocaine and amphetamine at the cellular level. Targeting a spinal cellular substrate for therapeutics that concomitantly spare cortical activation may result in pain relief that avoids the addictive side effects that have led to opioid abuse.

It is also now critical to assess the regenerative potential of these cells, by transplantation into an injured or diseased spinal cord. At present, it remains unresolved whether the pure neuronal populations (Fortin et al., 2016) will be more functionally beneficial than replacing populations of multipotential progenitors or mixtures of neurons (Kawai et al., 2021). The identification of dI-specific cell surface markers in this study (Table S4) enables the sorting of specific neural populations, which will make it possible to distinguish between these possibilities.

In summary, these studies uncover further mechanistic insights into the process of dI differentiation in vitro to derive bona fide spinal sensory neurons in the large numbers needed for clinical and drug discovery applications. They represent a critical step toward using ESC-derived dIs for studying the effects of various addictive compounds and in vitro disease modeling for pain disorders.

Limitations of the study

Although this study shows that directed differentiation protocols can generate bona fide sensory modality-specific dorsal spinal interneurons, some hurdles remain. In each protocol, the dIs arise as heterogeneous populations, and the yield of dI3 and dI6, the populations that regulate mechanosensation and gait, remain modest. Further transcriptomic analyses are needed to assess the mechanisms by which dP state1 resolves into dI4–dI6 and dP state2 resolves into dI1–dI3. The identification of any dI-specific regulators might permit the generation of large numbers of a specific dI population using the CHIR-mediated expansion protocol. Finally, in-vitro-derived dIs need to be electrophysiologically assessed to determine their maturity; such an analysis requires an understanding of the electrophysiological signatures of endogenous spinal sensory neurons, which are not yet well defined.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Samantha J. Butler, PhD. Email: butlersj@ucla.edu.

Materials availability

Further information regarding the resources and reagents should be directed to the lead contact- Samantha Butler, PhD (butlersj@ucla.edu).

Bulk and Single-cell RNA-seq data have been deposited at GEO and are publicly available. Accession numbers are listed in the key resources table. Microscopy data reported in this paper will be shared by the lead contact upon request.

KEY RESOURCES TABLE

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Antibodies
Goat polyclonal anti-Lhx2	Santa Cruz	Sc-19344; RRID:AB_2135660
Goat polyclonal anti-Isl1	R&D Systems	AF1837, RRID:AB_2126324
Rabbit polyclonal anti-Pax2	Invitrogen	71–6000, RRID:AB_2336046
Mouse monoclonal anti-Lhx1/5	Developmental Studies Hybridoma bank	4F2, RRID:AB_531784
Goat polyclonal anti-Pax3	R&D Systems	AF2457, RRID:AB_416599
Mouse monoclonal anti-Pax7	Developmental Studies Hybridoma bank	AB_528428, RRID:AB_528428
Rabbit monoclonal anti-phosphoSMAD1/5	Cell Signaling Technology	11971, RRID:AB_2797785
Rabbit polyclonal anti-phospho Histone H3	Cell signaling Technology	9701, RRID:AB_331535
Mouse monoclonal Tubulin III beta (clone-)	Covance	MMS-435P, RRID:AB_2315514
Goat polyclonal anti-Sox2	Santa Cruz	sc-17320, RRID:AB_2286684
Rabbit polyclonal anti-Pax6	MBL life science	PD-022, RRID:AB_1520876
Goat polyclonal anti-T (Brachyury)	R&D Systems	AF-2085, RRID:AB_2200235
Guinea pig polyclonal anti- FoxD3 sera	gift from Thomas Müller; (Müller et al., 2005)	N/A
Guinea pig polyclonal anti- Lmx1b sera	gift from Thomas Müller; (Müller et al., 2005)	N/A
Guinea pig polyclonal anti- Olig3 sera	gift from Thomas Müller; (Müller et al., 2005)	N/A
Chemicals, peptides, and recombinant proteins
CHIR99021	Tocris Bioscience	4423
IWR1e	Millipore	681669
XAV939	Stem Cell Technologies	72674
Matrigel membrane matrix	Fisher Scientific, Corning	CB-40234
DAPT	Sigma Aldrich	565784
Human recombinant BMP4	Thermo Fisher Scientific	PHC9534
Retinoic acid	Sigma Aldrich	R2625
Human recombinant Wnt1	Thermo Fisher Scientific	PHC1804
Human recombinant Wnt2	Millipore	SRP-6560
Human recombinant Wnt3a	R&D Systems	5036-WN-010
Mouse recombinant Wnt4	R&D Systems	475WN005
Human recombinant Wnt5b	R&D Systems	7347-WN
Mouse recombinant Wnt9b	R&D Systems	3669-WN-025
Human basic FGF (bFGF)	Thermo Fisher Scientific	PHG0023
DMEM:F12	Hyclone	SH30023
Neurobasal medium	Thermo Fisher Scientific	21103049
B27 supplement	Fisher Scientific	17–504-044
N2 supplement	Thermo Fisher Scientific	17502048
Mouse LIF	Millipore	ESG1107
Trypsin Inhibitor 1x solution	Sigma-Aldrich	T7659
Paraformaldehyde 16%	Sigma-Aldrich	28906
Phosphate buffer saline (PBS) 1X	Hyclone	SH30028.02
Prolong Gold	Thermo Fisher Scientific	P36930
Critical commercial assays
RNeasy plus mini kit	Qiagen	74106
Superscript IV First-Strand synthesis kit	Thermo Fisher Scientific	18091050
Dead cell removal kit	Miltenyi Biotec Inc.	130–090-101
EdU labelling assay	Invitrogen	C10338
Qia shredder	Qiagen	79654
Software and algorithms
Corel DRAW version 20	Corel DRAW	www.coreldraw.com
Zeiss Blue lite	Zeiss	https://www.zeiss.com/microscopy/us/products/microscope-software/zen-lite.html
Prism 9	Graphpad	https://www.graphpad.com/scientific-software/prism/
Image J	National Institute of Health	https://imagej.nih.gov/ij/
BioRender		https://biorender.com
R-Package V.3.2.5	The R project for Statistical computing	https://www.r-project.org/
Cell Ranger v 3.0.2	10x Genomics	https://support.10xgenomics.com/single-cell-gene-expression/software/
Seurat v. 3.1.1	Satija Lab	https://satijalab.org/seurat/
EnrichR	Ma’ayan Lab	https://maayanlab.cloud/Enrichr/
Experimental model: Cell lines
MM13 mouse embryonic stem cell line	Tom Jessell Lab, Columbia University (Wichterle et al., 2002)	N/A
Deposited data
Bulk and single-cell RNA-Seq	Gene Expression Omnibus (GEO)	GEO: GSE185891
Oligonucleotides	Mouse primers for qPCR analysis
	Forward primer	Reverse primer
GAPDH	GGCCTTCCGTGTTCCTAC	TGTCATCATACTTGGCAGGTT
Lhx2	CAGCTTGCGCAAAAGACC	TAAAAGGTTGCGCCTGAACT
Foxd3	CCCCAACACTGACCAACAG	GTTTGCTCCGCCAGCTTA
Isl1	ATCGAGTGTTTCCGCTGTGT	ATGGGTTCGGCTGCCATTT
Lhx5	CGGAGCAAAGTCTTCCACCT	GGACGACACCGAGTTGAGAC
Lbx1	GCCGAAGGCCGCGATG	TGCGTGATTTTCGCCGTTTC
Lmx1b	GCATCAAGATGGAGGAGCAC	CTCGTTGACTCGCATCAGG
Meis1-spliced variant A	ACCCTGATGGACAGCCAATG	GCATACTTTGCAGCCCTGG
Meis1-spliced variant B	CAATCGAAAGCCAGGGGGAA	CTCGGCTGTCCCATACTCAC
Meis2	GCACAGGATGCGTTTCAGTG	TGTCTAACCCATCGCCTCCA
Meis3	TGGGTTCAAAGGGAAGGAAGG	GAGGTGAAGTCCTGAAGCCC
Irx1	CGGCCGCGGTCACCT	CAGTTCATACTGAGAGCCCATCT
Irx3	ACAGGTGACAAGAAAAATATCGC	CAGCGTCCAGATGGTTCTGT
NeuN	ATCGTAGAGGGACGGAAAATTGA	GTTCCCAGGCTTCTTATTGGTC
Tubb3	ATGGACAGTGTTCGGTCTGG	AGCACCACTCTGACCAAAGAT
Sox2	CAAAAACCGTGATGCCGACT	CGCCCTCAGGTTTTCTCTGT
Pax3	AAACCCAAGCAGGTGACAAC	GCGTCCTTGAGCAATTTGTC
Pax7	TGGCCAAACTGCTGTTGATT	TAGGCTTGTCCCGTTTCCAC
Olig3	ACGTTAGGGCTGTTCCAAGG	GGCAAGGGGGTAGGATGAAC
Wnt1	ACTACGTTGCTACTGGCACT	GTTCACGATGCCCCACCATC
Wnt2	TGATGTAGACGCAAGGGGGT	CATGTACCACCATGAAGAGCTGA
REAGENT or RESOURCE	SOURCE	IDENTIFIER
Wnt3a	CTACCCGATCTGGTGGTCCT	ACAGAGAATGGGCTGAGTGC
Wnt4	CAGGAAGGCCATCTTGACACAC	GTCTTTACCTCGCAGGAGCC
Wnt5b	CCGGCCACCTTTCTCTCTTA	CGGAGCTCCCACGCTGAT
Wnt10b	ACGGTTTAAGCCCCAAGAGC	GAGTCAGTCAGCGCCTCC
Axin2	GCTGCGCTTTGATAAGGTCC	CGCGAACGGCTGCTTATTTT
Tcf7	GGACGGGAGCACACTTCG	GCCTTCAGGCCGTCCTCT

All original code is available from the lead contact upon request.

Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

ESC maintenance and two-dimensional (2D) differentiation

The mouse ESC line, MM13, was routinely maintained in ES cell media (DMEM+20% FBS) with 100u/ml of mouse leukocyte inducing factor (LIF). Cells were maintained on mitotically inactive irradiated mouse embryonic fibroblasts (MEFs) and were passaged at least twice before starting the differentiation. To prepare cells for the differentiation, ESC colonies were dissociated with 0.25% trypsin and plated on gelatin coated plates at 1:10 dilution to reduce MEF transfer. Cells cultured on gelatin coated plates and allowed to proliferate for 1–2 days. To initiate differentiation, cells were plated on 0.1% gelatin coated 24-well CellBIND dishes (Corning) at 20,000 cells/ml density in 0.5mL/well N2/B27 medium that contains 10ng/ml basic fibroblast growth factor (bFGF). The N2/B27 medium contains 1:1 portions of Advanced Dulbecco’s modified Eagle Media (DMEM) F12 (Hyclone) and Neurobasal media (Thermo Fisher Scientific) supplemented with 0.5x N2, 1x B27 supplement (with Vitamin A), 2mM L-glutamine, 0.1mM 2-mercaptoethanol (ME), and 1% BSA. On day 1, small colonies of cells can be observed attached to the bottom of the wells. On day 2, cells were supplemented with 10ng/ml bFGF and 5μM CHIR99021 for 24 h to induce a neuromesodermal identity (Gouti et al., 2014). On day 3, media was changed as follows: for protocol 1, cells were exposed to 100nM RA for 48hrs. For protocol 2, 10ng/ml human recombinant BMP4 (Thermo Fisher Scientific) was added concomitantly with 100nM RA (RA + BMP4) for 48hrs. For protocol 3, cells were exposed to 100nM RA for 24hrs, followed by a pulse of RA + BMP4 for another 24 h. To induce terminal differentiation in all three conditions, the medium containing growth factors was replaced with basic N2/B27 medium at day 5, and cell allowed to differentiate to day 9. At the end of the differentiation, the cultures were fixed for immunohistochemical analysis or lysed to obtain RNA for RNA-Seq or quantitative reverse transcriptase PCR analysis.

METHOD DETAILS

Bulk-RNA sequencing and data processing

Total cell lysate was obtained in buffer RLT (Qiagen) at different time points from three independent differentiations (biological replicates) conducted parallelly. RNA extraction was performed using RNeasy mini kit (Qiagen) and the quality was determined using Agilent Technologies 2100 Bioanalyzer and only samples with a RIN score >8.0 were sequenced. Stranded libraries were constructed using Universal plus mRNA sequencing kit (NuGEN) and sequenced onto 2 lanes of NovaseqS1 to generate 30 million reads/sample. Obtained reads were aligned to the latest mouse reference genome (mm10) using the STAR spliced read aligner. Pearson correlation and principal component analysis (PCA) was performed to determine the similarities among the samples and replicates. Differentially expressed genes were identified using DESeq and EdgeR package on R platform. The candidate genes were selected using cutoff false discovery rate (FDR <0.05) and log fold change (>2). Differentially expressed gene lists were further subjected to Enrichr, Metascape (Zhou et al., 2019) and DAVID functional annotation tools to identify enriched gene ontology (GO) categories.

WGCNA analysis

We used WGCNA package (Langfelder and Horvath, 2008) to identify the expression dynamics of gene modules associated with differentiations for protocol 1, 2, and 3. WGCNA was performed separately on RNA-Seq samples from protocol 1 (R-branch), 2 (C-branch), and 3 (B-branch) with day 0 and day 3 as common samples. First, differentially expressed genes were extracted for each group to determine the most variable genes. The Pearson correlation matrices were calculated for all RNA pairs in each group and were transformed into adjacency matrices using the power function. A dynamic tree-cut algorithm was used to identify gene co-expression modules where modules were defines as branches cut off from the tree and labeled in unique colors. The module eigengene (ME) represents the first principal component for each module. Gene ontologies associated with the modules were determined by supplying top 100 genes to the Enrichr platform.

Protein-protein interaction (PPI) network analysis

We used the metascape algorithm (www.metascape.org) to identify PPI networks induced by 6hrs and 24hrs of RA exposure. First, significant differentially expressed (DE) genes were extracted by pairwise comparison with the day3 dataset (day3 v/s 6hrs RA and day 3 v/s 24hr RA) using log fold change> 2 and false discovery rate (FDR) < 0.01 as cut off for the upregulated genes. The entire gene list obtained from the DE analysis was then submitted to Metascape platform with Mus musculus as both input and analysis species. Both gene ontologies (GO) and PPI network were then retrieved for each gene list using the default Metascape parameters applied under the express analysis tool (minimum network size = 3, maximum network size 500).

Single-cell RNA sequencing and analysis

Preparation of single cell suspension

Day 9 cultures from both protocol1 (RA) and protocol3 (RA + BMP4) were dissociated with 0.25% cold trypsin for 5 min. Trypsin activity was stopped with the addition of trypsin inhibitor (Sigma-Aldrich) in 1:1 ratio. Cells were then pelleted by centrifugation at 1000RPM and resuspended in 1x PBS. The ratio of viable cells was determined using Trypan blue staining. Dead cells were removed using MACS dead cell removal kit (Miltenyi Biotec Inc.) by following the manufacturer’s protocol. Eluted live cells were then suspended in 1X PBS containing 0.04% BSA solution (400μg/ml) for the library preparation.

Library preparation, and sequencing

~10,000 live cells/conditions were used to construct single-cell specific cDNA libraries using protocol described in 10x Genomics chromium single cell 3’ reagent kit (v3.1 Chemistry). Briefly, cells were partitioned into nanoliter-scale Gel-Beads-in-emulsion (GEM) using the 10x Chromium controller. Each GEM contains a unique barcode which is shared among the cDNA generated from a single cell. Cells were then lysed, and cDNA synthesis and feature barcoding were performed in the GEMs. The sequencing libraries were recovered using Magnetic separation and the quantity and quality of cDNA were assessed by Agilent 2100 expert High Sensitivity DNA Assay. cDNA samples were sequenced on 1 lane of NovaSeq 6000 S2 flow cell and reads were mapped to mouse mm10 genome using Cell Ranger v.3.0.2 to generate fasq files.

Filtering genes and cells

Seurat package v3.1.1 was used for the analysis. For each condition, genes expressing in less than 3 cells were removed and cells with >250 features were retained. PCA was used for dimension reduction (npcs = 50) with top 3000 most variable genes. Raw count data were normalized using regularized negative binomial regression with SCT transformation. Cell clustering is done using Shared Nearest Neighbor (SNN) Graph method and cluster specific markers were identified by Wilcox Rank-Sum test. Each cell cluster was annotated by a combination of the following methods 1) the expression of known canonical cell-type markers 2) Gene Set Enrichment Analysis (GSEA) of cluster specific markers determined by FindMarkers function in Seurat.

Reclustering of neuronal clusters

From the entire single-cell dataset, neuronal cells were first identified by the expression of pan-neuronal markers-Tubb3 and Elavl1/2. Roughly 30% cells express these markers in both protocol 1 and protocol 3 and were thus annotated as the mature neurons. The neuronal populations were then extracted using the Seurat package in R, by clustering the cells expressing Tubb3 >2 (median normalized log expression levels) using the Shared Nearest Neighbor (SNN) Graph method. This manipulation resulted in 12 clusters in protocol 1 (RA) and 13 clusters in protocol3 (RA + BMP4). Clusters were visualized with Uniform Manifold Approximation Projection (UMAP) method and the unique cluster-specific markers were identified using FindMarker function in Seurat.

GO analysis for neuronal clusters

We used DAVID functional annotation tool (Huang et al., 2009b) (https://david.ncifcrf.gov/tools.jsp) for identifying enriched biological processes and pathways associated with different neuronal clusters. For this, cluster-specific gene lists were sorted according to log fold change (FC) values and all genes with log FC = >0.2 were selected and submitted to the DAVID platform with default settings and Mus musculus as the input species. Ontological terms were then surveyed to determine if terms reflect the known functions of the dIs that the given cluster belongs to. Selected GO terms were arranged in the decreasing order of −log10 (p value) and represented as bar plots.

Integration of single-cell datasets from different origins

The in vivo spinal cord single-cell dataset (Delile et al., 2019) was downloaded from https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-7320/and loaded into R. Gene names were standardized to match our dataset using the biomaRt function (Durinck et al., 2009). This data was then loaded into Seurat (version 4.0) and processed using the SCTransform pipeline (Hafemeister and Satija, 2019). In vivo cell type labels were simplified and transferred to the in vitro datasets. The in vivo dataset was then integrated with the in vitro protocol datasets using Seurat’s reciprocal PCA integration pipeline, and projected into a three-dimensional UMAP space. Tabula Muris reference datasets of kidney (2781 cells), lung (5449 cells) and trachea (11,269 cells) were downloaded from https://figshare.com/articles/dataset/Robject_files_for_tissues_processed_by_Seurat/5821263, loaded into R, and then processed through Seurat, as above. The three Tabula Muris datasets were then integrated with the in vitro and in vivo spinal cord datasets using Seurat’s reciprocal PCA integration pipeline and embedded into a three-dimensional UMAP space. Code for this analysis is available upon request.

dI cell type extraction based on the marker gene expression

Cutoffs for gene expression levels when extracting cells were initially determined through the cluster-specific violin plots from the integrated in vivo and in vitro spinal cord datasets for the SCT assay. Cutoffs were then refined by plotting the extracted cells on a UMAP plot adjusting the values to exclude non-specific cells with background levels of the genes. The cutoffs values for each gene are as follows: dI1 = Barhl1 > 1 (in vivo-417 cells, in vitro-1120 cells), dI2 = Foxd3 > 1 (in vivo-552 cells, in vitro-395 cells), dI3 = Isl1 > 1 Tlx3 >0.5 (in vivo-1669 cells, in vitro-89 cells), dI4 = Pax2 > 1 and Lhx1 > 1 (in vivo-1943 cells, in vitro-810 cells), dI5 = Lmx1b > 1 and Pou4f1 >1 (in vivo-433 cells, in vitro- 216 cells), dI6 = Dmrt3 > 1 (in vivo-73 cells, in vitro-60 cells). To extract Tubb3+ cells in the in vivo and in vitro combined dataset, expression by cluster was plotted, and the clusters expressing high levels of Tubb3 were extracted as a whole. Code is available upon request.

Embryoid body (EB) culture, and expansion

To initiate EB formation, dissociated embryonic stem cells were seeded in low attachment plates (Corning) at 105 cells/ml density in N2/B27 medium supplemented with 10ng/ml bFGF. To induce posterior spinal cord and dorsal spinal patterning, RA and RA + BMP4 were added using the same timeline as described for the 2D differentiation, and the media was changed at every other day after day 5.

For expansion cultures

EBs were dissociated in single cell suspension at day 6 by treating with 0.25% Trypsin, followed by trituration to break up the larger cell clumps. Cells were pelleted and passaged in N2/B27 medium containing either DMSO or 5μM CHIR at 1:3 dilution. Media was changed at every third day during the expansion culture.

Embryoid body (EB) differentiation

To induce differentiation, EBs were dissociated into smaller clumps by a 5–7-min treatment with cold 0.25% Trypsin (Gibco), followed by trituration. Trypsin was neutralized with 1:1 addition of Trypsin inhibitor (Sigma-Aldrich) and cells were pelleted by centrifugation at 1000RPM. Cell pellets were resuspended in N2/B27 medium (without growth factors) and plated onto a Matrigel precoated plate. In our experience, the single cell suspension of 100–150 EBs provided enough cells to seed one 24-well plate. Spontaneous neural differentiation is induced when dissociated EBs are cultured in N2/B27 medium for 10 days. At the end of 10 days, multiple neural processes can be visualized by phase contrast microscopy and cultures were either fixed with 4% PFA for immunohistochemistry or lysed to extract RNA in RLT buffer (Qiagen).

Immunohistochemistry

For 2d differentiations, adherent cultures were first washed with 1xPBS and fixed with fresh cold 4% PFA for 10 min in the well. Following fixation, cultures were washed twice with 1xPBS to remove any remaining PFA. Cells were blocked with 1xPBS with 1% heat inactivated horse serum for 1 h and primary antibodies are added in the blocking solution for an overnight incubation at 4C. Following washes, species appropriate secondary antibodies (Jackson Immuno Research Laboratories) were added in 1x PBT (PBS +0.1% Triton 20) for 1hr. The cultures were then washed with 1x PBT to remove traces of secondary antibodies and counterstained with DAPI to stain nuclei. Plates were then imaged on Zeiss LSM800 inverted confocal microscope at 10× magnification.

For EBs, EBs were collected and fixed with fresh 4% PFA for 20 min at 4°C. EBs were then equilibrated with 30% sucrose for overnight at 4C and embedded into the OCT compound for cryo-sectioning. 12um sections were collected onto positive charged slides (VWR), incubated with the antibody blocking solution for 1hr, and then stained as above. Sections were counterstained with DAPI, mounted in Prolong Gold and imaged using a Zeiss LSM800 inverted confocal microscope.

Quantitative (q) RT-PCR analysis

RNA was isolated with RNeasy Mini Kit (Qiagen, cat no. 74104), and converted to cDNA using Superscript IV First Strand Synthesis kit (Thermo Fisher Scientific). qRT–PCR was performed using SYBR Green PCR Master Mix (Applied Biosystems) on Roche 480 Light-cycler (Roche). The relative fold expression was calculated using 2-∆∆Ct method of comparing the expression of the target gene with that of the housekeeping gene glyceraldehyde 3-phosphate dehydrogenase (Gapdh). All experiments were repeated at least 2–3 times (technical replicates) from at least two independent differentiations (biological replicates), and data is represented as mean ± SEM. The primer sequences for all target genes are listed in the Table S5.

QUANTIFICATION AND STATISTICAL ANALYSIS

Image quantification

To count the number of nuclei, images were converted to 8bit image in ImageJ, their threshold intensity was adjusted to capture only the florescent nuclei, and the number of nuclei counted using the analyze particle tool in ImageJ. For each image, the area under the DAPI staining was determined using ImageJ and target cells were represented as % of cells in 100μm2 DAPI+ area.

Statistics

Data are represented as mean ± SEM (standard error of the mean). Tests for statistical significance were performed using Prism software (version 9). Values of p < 0.05 were considered significant in all cases.