A controllable human spinal cord model with full dorsoventral patterning

Domain: physics
Published: 2026-03-28
Editorial outcome: major_revision

Abstract
Dorsoventral (DV) patterning of the spinal cord (SC) is orchestrated by morphogen gradients that specify distinct neural progenitor domains, the dorsal roof plate (RP), ventral floor plate (FP), and delaminating neural crest cells (NCCs). While foundational insights into SC patterning have been gained from animal models, key aspects - such as the role of retinoic acid (RA), dynamics of NCC lineage development, and human-specific features - remain poorly understood due to limitations in existing in vivo and in vitro models. Here, we present a human pluripotent stem cell (hPSC)-derived SC model, termed microfluidic SC-like structures (µSCLSs), generated by applying bioengineered, antiparallel morphogen gradients via a microfluidic platform over micropatterned hPSC colonies. The µSCLS robustly recapitulates complete DV patterning with human specific transcriptional signatures. Using this platform, we uncover a previously unrecognized RA-BMP signaling crosstalk that could explain conflicting reports on the role of RA in SC DV patterning. We further demonstrate lineage-specific ventral migration of NCCs in response to chemoattractant cues, enabling direct visualization and mechanistic dissection of sensory vs. other fate trajectories. This controllable, reproducible, and human-relevant model provides a powerful system for probing human SC development, neural crest biology, and disease modeling with unprecedented resolution.

Introduction
During vertebrate neural tube (NT) development, the caudal spinal cord (SC) region is patterned along the dorsal-ventral (DV) axis into 11 distinct neural progenitor domains and two key signaling centers - the dorsal roof plate (RP) and the ventral floor plate (FP) (Fig. 1a)1,2. Each SC progenitor domain gives rise to specific neuronal subtypes that contribute to the formation of the central nervous system (CNS)3. Concurrently, neural crest cells (NCCs) delaminate from the dorsal NT and migrate ventrally (Fig. 1a), forming components of the peripheral nervous system (PNS), including the dorsal root ganglia (DRG) and the sympathetic ganglia4.

Classic animal studies have demonstrated that DV patterning of the SC is largely orchestrated by opposing gradients of morphogens secreted from adjacent tissues: dorsally derived BMPs and WNTs, and ventrally derived SHH (Fig. 1a)5–10. These morphogens induce the expression of specific transcription factors (TFs) within SC progenitor domains in a dose-dependent manner, activating lineage-specific differentiation programs (Fig. 1a). Some TFs, such as OLIG2 and IRX3, also engage in mutually repressive interactions that help stabilize domain identities, sharpen boundaries, and ensure the segregation of distinct progenitor pools.

Retinoic acid (RA), secreted from the somites franking the SC, also plays a role in DV patterning, although its precise function remains unclear. Explant studies have shown that RA can upregulate both dorsal and ventral TFs in the SC11–14; notably, both OLIG2 and IRX3 are induced by RA treatment12,14. Disruption of RA signaling in vivo results in SC ventralization, suggesting a potential dorsalizing effect of RA - although this interpretation is confounded by the concurrent impairment of RP development11,12,14.

Shortly after NCCs delamination, these cells undergo lineage bifurcation into sensory or sympathetic fates15,16. Sensory neuron precursors halt their migration near the NT to form the DRG, while sympathetic neuron precursors continue migrating ventrally to the dorsal aorta17,18. Although several chemoattractants guiding sympathetic precursor migration have been identified19,20, the mechanisms governing sensory precursor arrest remain poorly understood. In vivo studies are limited by significant technical challenges, including imaging constraints, physiological complexity, and difficulties in controlling gene function or signaling dynamics. Therefore, an in vitro model that faithfully recapitulates NCC development – and is amenable to live imaging and genetic or chemical perturbation – would provide a powerful experimental tool for dissecting the cellular mechanisms underlying NCC migration and lineage development.

Recent studies have also identified human-specific features of SC development, such as unique gene expression profiles and the early emergence of oligodendrocyte progenitors21,22. These findings underscore the need for human-relevant SC models that are controllable, reproducible, ethically accessible, and compatible with hypothesis-driven research and therapeutic development23–25.

Human pluripotent stem cell (hPSC)-based developmental models have emerged as powerful, human-relevant platforms for studying development, physiology, and disease. However, most current models rely on uniform culture conditions without externally imposed signaling gradients, making regional patterning and tissue organization heavily dependent on self-organization. As a result, hPSC-derived SC-like tissues often exhibit inconsistent, incomplete patterning and low efficiency26–29. More robust models have recently been developed by introducing bioengineered gradients - such as those generated by microfluidics or synthetic organizers - to hPSC-derived neural tissues30–33. Notably, a recent model using synthetic organizers successfully recapitulated full DV patterning of the SC33. However, this approach has some limitations: the use of synthetic organizers restricts experimental control, the free-floating culture format hinders high-resolution imaging, and ventral migration of NCCs was not observed33.

In this study, we present a fully patterned and controllable human SC model, generated by applying microfluidic gradients of exogeneous morphogens to hPSC-derived SC-like tissues. We refer to these as microfluidic SC-like structures (µSCLSs). The µSCLS exhibits complete DV patterning - including all 11 progenitor domains as well as the dorsal RP and ventral FP - and recapitulates human-specific gene expression profiles. Using this controllable platform, we investigate the previously ambiguous role of RA in SC DV patterning and uncover a crosstalk between RA and BMP signaling. Furthermore, we demonstrate that NCCs within the µSCLS undergo lineage-specific ventral migration driven by chemoattraction, faithfully modeling in vivo behavior.

Results
DV patterned µSCLS
To generate DV patterned µSCLSs, we designed a microfluidic device composed of a central channel connected to two medium reservoirs, into which dorsalizing and ventralizing factors could be introduced, respectively (Fig. 1b and Supplementary Fig. 1a; see Methods). The reservoir containing dorsalizing factors was designated “dorsal reservoir,” while the one containing ventralizing factors was termed “ventral reservoir” (Fig. 1b).

The µSCLS protocol began with microcontact printing an array of rectangular Geltrex-coated adhesive islands onto a coverslip, which guided the formation of hPSC colonies with defined rectangular geometries. After colony seeding, the microfluidic structural layer was aligned and affixed to the coverslip (Supplementary Fig. 1b). Geltrex was chosen for its enrichment in basement membrane proteins characteristic of the NT environment17,18. The positions of hPSC colonies were precisely aligned with the microfluidic channel using positional markers engraved on the device (Fig. 1b and Supplementary Fig. 1b).

Two days after initial hPSC colony formation in mTeSR medium, the microfluidic channel was filled with 100% Geltrex to embed the colonies within a three-dimensional matrix, thereby supporting epithelialization and lumen formation (Fig. 1c and Supplementary Fig. 1b). At this point, the culture medium was switched to neural induction medium (NIM; see Methods). Because early activation of WNT signaling promotes SC fate during neural induction30,32,34, we supplemented NIM with the WNT agonist CHIR99021 (3 µM) from day 3 to day 5 to promote SC identity specification in the hPSC colonies (Fig. 1c ).

To induce DV patterning in µSCLS, opposing gradients of BMP4 (150 ng mL-1) and smoothened agonist (SAG; 250 nM), a SHH signaling agonist, were established by adding these factors to the dorsal and ventral reservoirs, respectively, starting on day 5 (Fig. 1c and Supplementary Fig. 1b). Passive diffusion of fluorescent dextran was used as a proxy to assess morphogen gradient formation, confirming that stable chemical gradients were established within approximately 48 hr (Supplementary Fig. 1c-e). Our µSCLS protocol supported continuous tissue growth over time (Fig. 1d&e). By day 12, µSCLSs had acquired a thoracic SC identity, as indicated by expression of the thoracic SC marker HOXC9 and the absence of the more caudal marker HOXC10 (Fig. 1f).

We also performed temporal immunofluorescence analysis and live imaging to examine tissue dynamics during µSCLS development. Between day 2 and day 6, we observed gradual cellular reorganization, characterized by progressive alignment of the long axes of cell nuclei along the thickness of µSCLS tissues, indicative of increasing epithelialization (Supplementary Fig. 1f). During this period, µSCLS acquired neural identity, as evidenced by gradual loss of OCT4 and NANOG expression and sustained expression of SOX2 (Supplementary Fig. 1f). Development of the µSCLS progressed through a transient neuromesodermal progenitor (NMP) state, with a significant portion of cells on day 4 co-expressing SOX2 and Brachyury (BRA) (Supplementary Fig. 1g). The SOX2+BRA+ putative NMP cells were restricted to the dorsal and ventral poles of the µSCLS on day 6 and became undetectable by day 10 (Supplementary Fig. 1g).

Live imaging using a Lifeact-GFP reporter hPSC line revealed the formation of a tubular µSCLS structure enclosing a single, continuous central lumen within 2 days (Supplementary Fig. 1h and Supplementary Movie 1). Immunostaining for ZO-1, a tight junction protein and apical polarity marker, confirmed that the apical surfaces of µSCLS faced the central lumen (Supplementary Fig. 1i). In vivo, neural progenitor cells in the pseudostratified neuroepithelium of the NT undergo interkinetic nuclear migration, where mitotic cells localize at the apical surface and nuclei in S phase are displaced toward the basal side. To assess whether µSCLSs recapitulate this organization, we performed EdU incorporation assay (marking S-phase nuclei) and immunostaining for phospho-histone H3 (pH3, a mitosis marker) on day 5 and day 12. EdU-labeled nuclei were preferentially located at the basal side, whereas pH3+ mitotic cells were concentrated near the apical surface (Supplementary Fig. 1j). To further confirm their pseudostratified cytoarchitecture, we generated µSCLSs by spiking single cytoplasmic EGFP-expressing hPSCs into wild-type cells. On both day 5 and day 12, EGFP-labeled cytoplasm of individual cells was observed spanning from the apical to the basal surfaces, consistent with the formation of a pseudostratified epithelium (Supplementary Fig. 1j).

Immunofluorescence analysis of day 12 µSCLSs revealed successful DV patterning, as evidenced by the spatially organized expression of canonical DV markers PAX3, PAX6, and OLIG2 along the DV axis (Fig. 1g). PAX6 and OLIG2 were expressed in mutually exclusive domains, whereas PAX3 and PAX6 exhibited slight overlap, consistent with the dorsal-intermediate regions of the SC21,35 (Fig. 1g and Supplementary Fig. 2a).

The dorsal-intermediate regions within the µSCLS were further identified by the co-expression of PAX6 and PAX7, which marks the dp4 - dp6 domains21,36 (Supplementary Fig. 2b). The ventral intermediate domain p2 was identified by co-expression of PAX6 and NKX6-1, along with the absence of OLIG221,37 (Supplementary Fig. 2c). Additional intermediate domains in the µSCLS were delineated through analysis of DBX1 and DBX2 expression using RNA fluorescence in situ hybridization (RNA-FISH) (Supplementary Fig 2d&e). DBX1 was localized to a distinct band dorsal to the OLIG2 expression domain, indicative of the intermediate domains p0 and dp621,38 (Supplementary Fig. 2d). DBX2, a broader intermediate domain marker, largely overlapped with DBX1 expression (Supplementary Fig. 2e). Notably, a DBX2+DBX1- region was observed dorsal to the DBX1+ domain, suggesting the presence of the dorsal intermediate domain dp521,38 (Supplementary Fig. 2e).

The presence of RP- and FP-like regions in day 12 µSCLSs was confirmed by LMX1a expression at the dorsal pole and FOXA2 expression at the ventral pole, respectively 2,39 (Fig. 1h and Supplementary Fig. 2f). Additionally, SOX10+ NCCs were observed near the dorsal pole, delaminating from the µSCLS (Fig. 1h and Supplementary Fig. 2f). The presence of FOXA2+ FP-like region at the ventral pole of the µSCLS prompted us to examine its ventral patterning. Ventral SC progenitor markers FOXA2, NKX2-2, and OLIG2, indicative of FP-, p3-, and pMN-like identities, respectively, showed a ventral-to-dorsal organization in the µSCLS (Supplementary Fig. 2g). While expression of these lineage markers was mostly exclusive, some cells co-expressed NKX2-2 and OLIG2, early emergence of which is a characteristic of human spinal cord development21,22 (Supplementary Fig. 3g). Additionally, some cells in the ventral µSCLS region co-expressed NKX2-2 and FOXA2; such cells have also been identified in e9.5 mice, corresponding to CS12 human embryos40 (Supplementary Fig. 2g).

µSCLSs exhibited high patterning efficiency across two human embryonic stem cell (hESC) lines and one human induced pluripotent stem cell (hiPSC) line, with an average efficiency of 97.9 ± 6.2% (mean ± s.d.). Patterning efficiency was defined by the spatial segregation of PAX3, OLIG2, and FOXA2 expression domains along the DV axis (Fig. 1i and Supplementary Fig. 3a-c; see Methods). In contrast, imposing a single gradient of either the dorsalizing signal BMP4 (150 ng mL-1) or the ventralizing cue SAG (250 nM) resulted in lower and more variable efficiencies of 77.6 ± 30.5% and 36.1 ± 37.6% (mean ± s.d.), respectively (Supplementary Fig. 3d-f).

DV patterning in µSCLSs was tunable through modulation of morphogen concentrations. Specifically, reducing BMP4 from 150 ng mL-1 to 50 ng mL-1 led to a contraction of the dorsal PAX3 expression domain (Supplementary Fig. 3g). Similarly, reducing SAG from 250 nM to 125 nM resulted in shortened expression domains of the ventral markers OLIG2 and FOXA2 (Supplementary Fig. 3h&i). Interestingly, while ventral OLIG2 expression remained unaffected by BMP4 reduction, FOXA2 expression at the ventral pole of µSCLSs showed an unexpected shrinkage, suggesting cross-regulatory effects between dorsalizing and ventralizing signals (Supplementary Fig. 3i). These findings underscore the complexity of SC DV patterning under the influence of opposing morphogen gradients.

DV patterning also promoted region-specific neurogenesis in µSCLSs. By day 12, neurite-like processes were observed extending from µSCLSs (Fig. 1j), and immunostaining for βIII-TUBULIN confirmed their neuronal identity (Fig. 1k). ISL1⁺HB9⁻ sensory neurons were identified among the dorsal NCC population, while ISL1⁺HB9⁺ motor neurons emerged at the periphery of the ventral OLIG2⁺ motor neuron progenitor domain, resembling the developing mantle layer 41 (Fig. 1l). Some motor neuron precursors retained OLIG2 expression, indicating an early stage of neuronal differentiation42 (Fig. 1l).

Transcriptome analysis of µSCLS
Next, we performed single-cell RNA sequencing (scRNA-seq) on day 12 µSCLSs to assess their cellular composition. Unbiased clustering using Uniform Manifold Approximation and Projection (UMAP) identified five distinct cell clusters, annotated as ‘Neuroepithelium’ (NE), ‘Neural Crest’ (NC), ‘Dorsal Root Ganglion’ (DRG), ‘Neuron’, and ‘Mesoderm’, based on lineage-specific marker expression (Fig. 2a&b, Supplementary Fig. 4a-c). The small population of Mesoderm cluster likely arose as a byproduct of early WNT activation during neural induction, which is known to transiently induce a neuromesodermal progenitor state34. The Mesoderm cluster contained only a small proportion of cells and showed high expression of sclerotome-related markers PAX1, PAX9, and NKX3-2 (Supplementary Fig. 4b)43. Immunostaining for the mesenchymal marker TWIST1 revealed the presence of these mesodermal cells mostly located near the ventral pole of µSCLSs (Supplementary Fig. 4d). Sclerotome specification requires the ventralizing signal SHH44,45, potentially accountable for the ventral localization of TWIST1+ mesodermal cells in µSCLSs. All cell populations in µSCLSs expressed thoracic markers HOXA9, HOXB9, HOXC9 and more anterior HOX genes whereas more caudal HOX genes were not expressed (Supplementary Fig. 4e), confirming a thoracic identity of µSCLS.

To determine the developmental stage represented by day 12 µSCLSs, we conducted principal component analysis (PCA) comparing the averaged gene expression profiles of µSCLSs with those of human embryonic SC tissues across developmental stages (Fig. 2c). This analysis revealed that the transcriptomic profile of day 12 µSCLSs most closely resembled that of the Carnegie stage (CS) 12 human SC (Fig. 2c).

To validate the authenticity of cell types within µSCLSs, we integrated the µSCLS transcriptomic dataset with that of human CS12 SC. The mesodermal population from the µSCLS dataset was excluded from this analysis due to its low abundance and limited relevance. Additionally, given the scarcity of DRG cells in CS12 human samples, DRG cells from CS14 human embryos were included to improve coverage in the integrated analysis. Upon integration, day 12 µSCLS clusters showed strong concordance with corresponding human SC cell types, with substantial overlap between µSCLS and CS12 human clusters (Fig. 2d). Consistently, Pearson correlation coefficients calculated for matched clusters - based on averaged gene expression - were highest between corresponding µSCLS and human SC cell types (Fig. 2e).

The NE cluster was further analyzed using a ‘knowledge matrix’ approach, which assigns cell identities based on the expression of established molecular markers22,46,47. This analysis revealed the presence of all 13 embryonic SC NE subtypes (Fig. 2f&g and Supplementary Fig. 4f&g). Upon integration with NE cells from human CS12 SC tissues, µSCLS-derived NE cells exhibited a coherent DV organization, with each subtype aligning closely with its in vivo counterpart in the integrated UMAP plot (Fig. 2f and Supplementary Fig. 4f). Likewise, Pearson correlation coefficients calculated from average gene expression profiles were highest between most corresponding NE subtypes from µSCLS and human SC (Fig. 2g). Some discrepancies were noted among dorsal intermediate progenitors, potentially due to their low abundance in the datasets (Fig. 2g and Supplementary Fig. 4f). Marker genes used to classify NE subtypes also showed highly similar expression patterns between NE cells from µSCLS and human SC tissues (Supplementary Fig. 4f&g). When NE progenitors were grouped into RP, dorsal (dp1-3), dorsal intermediate (dp4-6), ventral intermediate (p0-2), pMN, p3, and FP domains, the resulting cellular compositions were comparable between µSCLS and human SC NE populations (Fig. 2h).

The gene expression profiles of NE subtypes within µSCLS were consistent with our immunofluorescence and RNA-FISH results (Supplementary Fig. 4g). Specifically, the RP population exhibited upregulated expression of LMX1A, while the FP population strongly expressed FOXA2 (Fig. 1h and Supplementary Fig. 4g). PAX3 was broadly expressed in dorsal progenitors, whereas PAX6 was upregulated in intermediate progenitors but absent in the pMN domain and more ventral regions (Fig. 1g and Supplementary Fig. 4g). The pMN population showed elevated OLIG2 expression but lacked PAX 6 expression (Fig. 1g and Supplementary Fig. 4g). Additionally, WNT1 and SHH were expressed in the RP and FP cell populations, respectively, consistent with their established roles as signaling centers that secrete these morphogens2 (Supplementary Fig. 4g). Co-expression of PAX7 and PAX6 was observed throughout the dorsal intermediate domains (Supplementary Fig. 2b&4g). Within the p2 domain, cells co-expressed NKX6-1 and PAX6 but did not show upregulation of OLIG2 (Supplementary Fig. 2c&4g). DBX1 and DBX2 were co-expressed in the p0 and dp6 domains, with DBX2 additionally extending into the dp5 domain (Supplementary Fig. 2d&2e&4g).

‘Knowledge matrix’-based cell identification was also applied to the Neuron population. Upon integration of the transcriptomic data from µSCLS neurons with that of human CS12 SC neurons and projection onto a UMAP plot, µSCLS neurons exhibited an overall DV organization and colocalized with their in vivo counterparts (Supplementary Fig. 5a&b). Most neuronal subtypes were identified; however, dI1 neuron was undetected in both µSCLS and human CS12 SC (Supplementary Fig. 5b). Considering that dI1 neuron has been identified in CS14 human SC, it is possible that the developmental stage of day 12 µSCLS was too early for the appearance of dI1 neurons22. Similar to human CS12 SC, µSCLS contained a higher proportion of ventral neurons (V0-V3, MN) compared to dorsal interneurons (dI1–dI6) (Supplementary Fig. 5b&c). When interneurons were grouped into dorsal (dI1-dI3), dorsal intermediate (dI4-dI6), and ventral intermediate (V0-V2) categories, Pearson correlation coefficients of averaged gene expression levels were highest between the corresponding neuron subtypes from µSCLS and human CS12 SC (Supplementary Fig. 5d). Additionally, marker genes used to classify neuron subtypes displayed similar expression patterns across matched neuron subtypes from µSCLS and human CS12 SC (Supplementary Fig. 5e).

Human-specific gene expression signatures were faithfully recapitulated in hPSC-derived µSCLS. Of 51 primate-enriched protein-coding genes expressed in the human CNS - including 10 human-specific genes - µSCLS expressed 34, with expression patterns resembling those in human CS12 SC22,48,49 (Fig. 2i). Another example is NKX6-2, which is broadly expressed across ventral domains in human CS12 SC, in contrast to its more restricted p1 domain expression in stage-matched mouse E9.5 SC22,50. This broader expression pattern of NKX6-2 was also observed in the ventral regions of µSCLS (Fig. 2j). Previous studies have shown that NKX2-2 and OLIG2 co-expressing cells arise earlier in human SC than in mouse SC21,22. Consistent with observations made using immunofluorescence, cells co-expressing NKX2-2 and OLIG2 were identified in the µSCLS transcriptomic data (Fig 2k and Supplementary Fig. 2g). Importantly, among µSCLS cells expressing either NKX2-2 or OLIG2, 8% co-expressed both genes - higher than in mouse E9.5 SC (2%) and comparable to human CS12 SC (10%) (Fig. 2k). These co-expressing cells in µSCLS also exhibited expression patterns of oligodendrocyte markers similar to those in human CS12 SC and showed weak expression of neurogenic genes, in line with their identity as oligodendrocyte precursors21,22 (Fig. 2l).

Together, transcriptomic analyses demonstrate that µSCLS comprises all major cell types, including all 13 NE subtypes, and faithfully recapitulates the transcriptomic landscape of the human embryonic SC.

Effect of RA in SC DV patterning
RA has been reported to promote the expression of both dorsal and ventral markers, making its overall effect on SC DV patterning ambiguous 14,26,27,51,52 (Fig. 3a). To investigate the effect of RA on SC DV patterning, we compared DV patterning of the µSCLS cultured with or without RA (500 nM; Fig. 3b). To minimize endogenous RA production, we replaced the basal medium with N2B27 lacking vitamin A and reduced BMP4 concentration to 50 ng mL-1. RA treatment produced an overall dorsalizing effect: the expression domain of the dorsal marker PAX3 expanded ventrally, while expression of the ventral marker OLIG2 was markedly reduced and FOXA2 expression was completely abolished (Fig. 3c&d).

RA signaling has been shown to influence rostrocaudal axial identity13,34. To exclude this possibility in RA-treated µSCLSs, we stained RA-treated or untreated µSCLSs for the thoracic marker HOXC9 and the lumbar marker HOXC10 (Supplementary Fig. 6a). Irrespective of RA treatment, µSCLSs were HOXC9+HOXC10-, supporting a stable thoracic identity. Thus, the RA-induced alterations in DV patterning of µSCLSs are unlikely to result from changes in rostrocaudal axial identity.

The dorsalizing effect of RA contrasts with existing in vitro SC DV patterning models and motor neuron differentiation protocols, where RA has been reported to promote ventral marker expression, including FOXA2 and OLIG227,51,53. Importantly, those studies were conducted under uniform culture conditions containing only ventralizing cues. To test the hypothesis that RA and BMP4 may interact synergistically, we repeated the experiment in the absence of a dorsal-ventral BMP4 gradient. Under these conditions, RA did not affect the proportion of FOXA2⁺ or OLIG2⁺ cells, suggesting a crosstalk between RA and BMP signaling pathways (Fig. 3e&f). Further supporting this interaction, µSCLS patterned with antiparallel BMP4 and SAG gradients exhibited a 67.0% increase in SMAD1/5 phosphorylation upon RA treatment (Fig. 3g). To examine whether SHH signaling is also required for RA-mediated dorsalization of µSCLSs, we repeated the experiment in the absence of the ventral-dorsal SAG gradient. Even without SAG, RA successfully dorsalized µSCLSs, as the expression domain of the dorsal marker PAX3 expanded, whereas the expression domain of the intermediate marker PAX6 was reduced (Supplementary Fig. 6b&c). Together, these data suggest that SHH signaling is not required for RA to exert a dorsalizing effect on DV patterning of µSCLSs.

To investigate the mechanism underlying RA-mediated enhancement of BMP signaling, we re-analyzed data from a mouse ESC-based motor neuron differentiation study that included retinoic acid receptor (RAR) chromatin immunoprecipitation sequencing (ChIP-seq) for RA-treated and untreated samples54. This analysis revealed 16 RAR binding sites located near genes involved in BMP signaling regulation, with 5 sites specific to RA-untreated cells and 8 unique to RA-treated cells (Supplementary Table 1). These findings suggest that RA may enhance SMAD1/5 phosphorylation by modulating RAR binding to regulatory elements of BMP pathway-associated genes.

To identify which of the 16 genes were differentially expressed in µSCLSs in response to RA, we quantified mRNA levels of each gene in RA-treated and untreated µSCLSs using quantitative polymerase chain reaction (qPCR). Four of the 16 genes - GDF3, SKI, SKIL, and SFRP4 - were significantly upregulated following RA treatment (Supplementary Fig. 6d). To test whether these genes mediated the dorsalizing effect of RA, we individually knocked down each gene using small interfering RNA (siRNA) and compared DV patterning in RA-treated and untreated µSCLSs. Among the siRNAs tested, only knockdown (KD) of GDF3 resulted in a significant reduction of target mRNA levels (Supplementary Table 7). We therefore focused on examining DV patterning in RA-treated and untreated µSCLSs under GDF3-KD. Notably, in GDF3-KD µSCLSs, RA treatment no longer altered the expression domain of the ventral marker OLIG2 (Supplementary Fig. 6e&f). In contrast, RA treatment still reduced FOXA2 expression, suggesting that RA-mediated DV patterning involves domain-specific downstream mediators (Supplementary Fig. 6e&g). To further investigate the epistatic relationship among RA signaling, BMP signaling, and GDF3, we performed qPCR under conditions with or without RA and with or without a dorsal-ventral BMP4 gradient. GDF3 expression was significantly upregulated by RA treatment but was unaffected by BMP4, supporting a model in which RA directly induces expression of GDF3, which in turn directly or indirectly modulates BMP signaling to regulate DV patterning.

NCC migration
To track the delamination and migration of NCCs, we developed a SOX10::T2A-Cre lineage-tracing hESC line. A T2A-Cre cassette was knocked in downstream of the endogenous SOX10 locus in a hESC line harboring a Cre-inducible ZsGreen cassette at the AAVS1 locus (Fig. 4a). Live imaging from day 6 to day 8 of µSCLS derived from this lineage-tracing line revealed gradual delamination of NCCs from the dorsal side, with ZsGreen expression initiating after delamination (Supplementary Movie 2). This observation aligns with in vivo murine studies reporting that SOX10 expression in NCCs begins post-delamination15,19. Prolonged culture showed some NCCs detached from dorsal NCC clusters, became completely embedded in the surrounding matrix, and migrated ventrally (Fig. 4b&c, Supplementary Fig. 7a, and Supplementary Movie 3). Tracking individual NCCs demonstrated directional migration toward the ventral side (Fig. 4d), supported by a polar histogram of migratory trajectories skewed ventrally (Fig. 4e). Compared to a random walk simulation, migrating NCCs exhibited reduced circular standard deviation and increased directionality, consistent with directed migration (Fig. 4f). Their mean square displacement exhibited a parabolic trend, further characteristic of directional migration (Fig. 4g).

The earliest fate bifurcation of trunk NCCs occurs between the sensory and sympathetic lineages. Sensory neuron precursors halt their migration adjacent to the NT to form the DRG, whereas sympathetic neuron precursors migrate further ventrally toward the dorsal aorta15–18. Immunofluorescence analysis revealed that expression of SIX1, a DRG marker, was mutually exclusive with SOX10 (Fig. 4h). SIX1+ DRG cells were intermingled with SOX10+ NCCs in the dorsal region, whereas no DRG cells were observed among the more ventrally located, sparsely distributed NCCs, presumed to be migrating ventrally (Fig. 4h). For quantification, we defined cells located above the dorsal tip of the µSCLS as "dorsally located," and those below the midpoint as "ventrally located" (Fig. 4i). On average, the number of ventrally located SOX10+ non-DRG NCCs was 6.79 ± 9.74 per device, while DRG cells were rarely found ventrally (0.14 ± 0.36) (Fig. 4j). These ventrally migrating non-DRG NCCs, suspected to be sympathetic neuron precursors, were examined for ASCL1 expression by immunostaining but were negative (Fig. 4k), possibly due to a lack of maturation cues from the dorsal aorta17.

To determine the timing of NCC lineage specification, µSCLSs derived from the SOX10::T2A-Cre lineage tracer were immunostained for SOX10 and SIX1 at multiple time points (Supplementary Fig. 7b). SOX10 expression was detected on day 9 in delaminated NCCs, whereas ZsGreen expression appeared with a one-day delay, becoming detectable on day 10. SIX1 expression was observed in day 11 µSCLS, suggesting that NCC specification toward the DRG lineage occurs after delamination. ZsGreen expression in SIX1+ DRG cells further indicates that these cells were originated from a multipotent SOX10+ NCC pool. This observation aligns with previous studies showing that delaminating NCCs are multipotent progenitors and that fate bifurcation occurs during migration15.

Both sympathetic and sensory NCC lineages are reportedly responsive to chemoattraction by CXCL1219,20. Consistent with this, transcriptomic analysis of day 12 µSCLS revealed expression of CXCR4, the receptor for CXCL12, in both NC and DRG clusters (Fig. 4l). Additionally, CXCL12 was expressed in the FP, potentially driving ventral NCC migration (Fig. 4l). Chemical inhibition of CXCR4 with AMD3100 (200 nM), a CXCR4 antagonist, resulted in a dorsal shift in the localization of both non-DRG NCCs and DRG cells (Fig. 4m). Specifically, the number of ventrally located non-DRG NCCs was significantly reduced from 6.79 ± 9.74 to 1.54 ± 4.96 per device (Fig. 4n). Distribution analyses confirmed a dorsal shift for both populations (Fig. 4o-r). The proportion of dorsally located non-DRG NCCs increased from 40.70 ± 24.92% to 64.71 ± 29.95% (Fig. 4q), and that of DRG cells increased from 38.16 ± 27.14% to 69.24 ± 32.96% (Fig. 4r). These results are consistent with in vivo findings demonstrating CXCL12 as a chemoattractant for both sympathetic and sensory precursors19,20. Together, these findings show that µSCLS faithfully recapitulates the chemoattraction-based, fate-dependent migration behavior of NCCs.

Methods
Ethical statement
All experimental procedures used in this work with hPSCs were approved by the Human Pluripotent Stem Cell Research Oversight Committee at the University of Michigan. The research conforms to the 2021 Guidelines for Stem Cell Research and Clinical Translation recommended by the International Society for Stem Cell Research.

Cell line
hPSC lines used in this study were H9 hESC (WA09, WiCell; NIH registration number: 0062), H1 hESC (WA01, WiCell; NIH registration number : 0043), and 1196a (human induced pluripotent stem cell) lines from the University of Michigan Pluripotent Stem Cell Core. Additionally, a Lifeact-GFP H2B-mCherry WIBR3 hESC line62, a gift from O. Reiner, a cytoplasmic EGFP expressing H9 hESC line, and a SOX10::T2A-Cre lineage tracer H9 hESC line were used. hPSC lines used in this study were all validated by original sources and in-house by immunofluorescence staining for pluripotency markers along with successful differentiation to three germ layers. All hPSC lines were maintained in a feeder-free culture system for at least ten passages and were authenticated as karyotypically normal by Cell Line Genetics. All hPSC lines were not contamination by mycoplasma (LookOut Mycoplasma PCR Detection kit; Sigma-Aldrich).

Cell culture
All hPSC lines were maintained under standard feeder-free conditions using mTeSR medium (STEMCELL Technologies). Cells were cultured on tissue culture plates coated with dehydrogenase-elevating virus (LDEV)-free, hESC-qualified, reduced growth factor basement membrane matrix Geltrex (Thermo Fisher Scientific). Cultures were visually inspected at each passage for the presence of mesenchymal-like cells, which indicate spontaneous differentiation. All hPSCs were used between passages 60 and 80.

Cytoplasmic EGFP expressing hESC line
To generate cytoplasmic EGFP expressing hESC line, the lentiviral vector containing EGFP63 was co-transfected with pMDLg/pRRE, pRSV-Rev, and pMD2.G (Addgene IDs: 12251, 12253, and 12259) into 293 T cells (Takara Bio) using CalFectin (SignaGen Laboratories). 10 μM Forskolin (Selleck Chemicals) was added after 24 h of transfection. On the third day following transfection, the supernatant was collected, concentrated using Lenti-X Concentrator (Takara Bio), and used for the transduction of an H9 ES cell line. Then, EGFP-positive cells were sorted by flow cytometry with Bigfoot spectral cell sorter (Thermo Fisher Scientific).

SOX10::T2A-Cre lineage tracer
The Sox10::T2A-Cre lineage tracer hESC line was generated in the background of the previously published ZsGreen conditional reporter ES cell line32. Briefly, T2A-Cre was inserted before the stop codon of the SOX10 locus by using CRISPR-Cas9-mediated genome editing. The donor plasmid containing the left and right arms in the SOX10 locus, T2A-Cre, PGK-Hydromycin as well as gRNA plasmids targeting the SOX10 locus was generated as shown in our previous publication32. The gRNA sequence was: 5’ cactgtcccggccctaaagg 3’.
The donor and gRNA plasmids were transfected into the ZsGreen conditional reporter H9 hESCs and positive clones were selected in the presence of hygromycin (50 µg mL−1) for 7 days. Individual ES cell clones were genotyped and checked for a successful knock-in.

Device fabrication
The microfluidic device consisted of a polydimethylsiloxane (PDMS) structural layer attached to a coverslip. To fabricate PDMS structural layer, PDMS curing agent and base polymer (Sylgard 184; Dow Corning) were mixed at a ratio of 1:10. The mixed PDMS prepolymer was cast onto a microfabricated silicon mold and baked at 110 °C for 1 h. Harris Uni-Core punch tool (6 mm diameter, Ted Pella) was used to punch medium reservoirs in the PDMS structural layer. In parallel, PDMS stamps with rectangular micropatterns were fabricated by casting PDMS prepolymer (with the ratio of curing agent to base polymer of 1:20) onto a microfabricated silicon mold. This was baked at 110 °C for 1 h and 30 min. PDMS stamps were then peeled off, sterilized using 100% ethanol, and coated with 1% (v/v) Geltrex solution at 4 °C overnight. The next day, glass coverslips (Thermo Fisher Scientific) were cleaned by first sonicating in 2% (v/v) Hellmanex III (Hellma USA, 9-307-011-4-507) for 30 min, rinsing 3 times in deionized water, and then sonicating in 100% ethanol for 30 min. The cleaned coverslips were blow-dried and treated with ultraviolet ozone (Ozone cleaner; Jelight) for 7 min. The Geltrex-coated PDMS stamps were blow-dried and placed in conformal contact with ultraviolet ozone-treated coverslips to transfer Geltrex adhesive patterns onto the coverslips. Stamps were peeled off and PDMS structural layers were then attached to coverslips. Alignment of microcontact-printed Geltrex islands in the patterning region was guided by position markers engraved on the structural layer and was done using a desktop aligner64.

Generation of µSCLS
On day 0, colonies of human PS cells in tissue culture plates were dissociated using Accutase (Sigma-Aldrich) at 37 °C for 8 min before being suspended in DMEM/F12 (Gibco) as single cells. Cells were centrifuged and resuspended at a concentration of 20 × 106 cells per mL in mTeSR containing the ROCK inhibitor Y27632 (10 μM, Tocris) to prevent dissociation-induced apoptosis of hPSCs65. 7 μL of the hPSC suspension was injected into the microfluidic channel through openings to the reservoirs. The microfluidic device was incubated for 45 min to allow hPSCs to attach to the Geltrex adhesive islands. Unattached cells were flushed out by adding fresh mTeSR containing 10 μM Y27632 to one of the reservoirs. The resulting monolayer hPSC colonies were cultured one day more in mTeSR without Y27632. On day 2, mTeSR medium was aspirated from the microfluidic channel and 7 µL of 100% Geltrex was injected into the channel. The device was incubated for 10 minutes for gelation, establishing a 3D culture environment. Neural induction medium (NIM) was added to both reservoirs after gelation. NIM consisted of basal medium and two SMAD inhibitors for neural induction: TGF-β inhibitor SB431542 (10 μM, StemCell Technologies) and BMP inhibitor LDN193189 (500 nM, StemCell Technologies). The basal medium, N2B27, comprised 1:1 mixture of DMEM/F12 (Gibco) and neurobasal medium (Gibco), 1% N2 supplement (Gibco), 2% B-27 supplement (with vitamin A, Gibco), Glutamax (2 mM, Gibco), 1% non-essential amino acids (Gibco), 100 µM 2-mercaptoethanol (Gibco), and 1% antibiotic–antimycotic (Gibco). From day 3 to day 5, CHIR99021 (CHIR, 3 μM, StemCell Technologies) was further added to the NIM. From day 5 to day 12, antiparallel signaling gradients were imposed by adding basal medium supplemented with BMP4 (150 ng mL-1, R&D Systems) to the dorsal reservoir and basal medium supplemented with SAG (250 nM, StemCell Technologies) to the ventral reservoir. For prolonged culture from day 12 to day 18, both reservoirs were filled with basal media supplemented with neurotrophic factors: BDNF (10 ng mL-1; StemCell Technologies), GDNF (10 ng mL-1; StemCell Technologies), CNTF (10 ng mL-1; StemCell Technologies), IGF-1 (10 ng mL-1; PeproTech), cAMP (1 μM; Sigma-Aldrich), and AA (0.2 μg mL-1; Sigma-Aldrich). In all cases, culture media were replenished daily.

Confirmation of neural induction
To confirm neural induction in Supplementary Fig. 1f, µSCLS were stained for OCT4, NANOG, and SOX2 at different time points. A µSCLS was categorized as expressing OCT4 if at least one cell in the tissue exhibited nuclear localization of OCT4. Similar criterion was applied for NANOG and SOX2.

Confirmation of pseudostratified epithelium
To assess the structure of the neuroepithelium, cytoplasmic EGFP expressing H9 hESCs were mixed into wild type H9 hESCs with a ratio of 1 to 200 before injecting into the device. Then, all steps were carried out identical to the standard µSCLS generation protocol.

µSCLS patterning efficiency
Quantification of µSCLS patterning efficiency was based on marker expressions. A µSCLS was deemed patterned if dorsal marker PAX3 was expressed in the dorsal side without a separate expression domain further ventral to the main expression domain. Similarly, for samples stained for ventral markers OLIG2 or FOXA2, similar criterion was applied to assess patterning. The success rate was quantified as the ratio of number of patterned colonies to the number of total colonies.

Microscopy
Whole-mount immunostained µSCLS were imaged using an Olympus DSUIX81 spinning-disc confocal microscope, Nikon X1 Yokogawa spinning-disc confocal microscope, or a Nikon A1R confocal microscope. An array of partially overlapping images (15% overlap) were taken to cover the entire µSCLS area. Images taken with Olympus DSUIX81 were first corrected for background and shade using the Fiji (ver. 2.16.0.) plugin BaSiC66 and then stitched using the Fiji plugin MIST67. Images taken with Nikon A1R were first denoised using the ‘denoise AI’ functionality in Nikon NIS-Elements platform. Then, the denoised images were stitched using the ‘Stitch Image’ functionality in Nikon NIS-Elements platform. For 3D reconstruction, z stack images were acquired with a slice thickness of 1 μm. Low-magnification brightfield images were acquired using Olympus CKX53 inverted microscope. Live imaging was conducted using an inverted epifluorescence microscope (Zeiss Axio Observer Z1; Carl Zeiss MicroImaging) enclosed in an environmental incubator (XL S1 incubator, Carl Zeiss MicroImaging) that maintained the cell culture conditions at 37 °C and 5% CO2. Images taken using Zeiss Axio Observer Z1 were stitched using the Grid/Collection stitching plugin in Fiji.

Dextran diffusion assay
Morphogen diffusion in the microfluidic channel was characterized by using Texas Red-labelled dextran (70 kDa, with a hydrodynamic radius of 5.8 nm, Invitrogen) as proxy. From day 5 to day 8 of µSCLS culture, 10 µM Texas Red-labelled dextran was supplemented to dorsal culture medium (basal medium supplemented with BMP4). Fluorescence images covering the whole microfluidic channel were taken at different time points using an inverted epifluorescence microscope (Zeiss Axio Observer Z1; Carl Zeiss MicroImaging). Images were stitched using the Grid/Collection stitching plugin in Fiji 2.16.0. Average fluorescence intensities across the width of each patterning region were obtained using Plot Profile tool in Fiji 2.16.0.

Patterning through single gradient
To assess the patterning efficiencies of µSCLS when patterned through BMP4 gradient alone in Supplementary Fig. 3d&e, 150 ng mL-1 BMP4 was supplemented to the dorsal reservoir from day 5 to 12. The ventral reservoir was filled with basal medium without SAG supplementation. Similarly, to assess the patterning efficiency through SAG gradient alone in Supplementary Fig. 3d&f, 250 nM SAG was supplemented to the ventral reservoir from day 5 to 12 while the dorsal reservoir was filled with basal medium without BMP4 supplementation.

Signaling molecule titration assay
To demonstrate controllability of patterning through titration of signaling molecule concentrations in Supplementary Fig. 3g-i, three conditions were used: Standard, Less SAG, and Less BMP4 condition. For standard condition, from day 5 to 12, 150 ng mL-1 BMP4 was supplemented to the dorsal reservoir and 250 nM SAG was supplemented to the ventral reservoir. For Less SAG condition, the SAG concentration was decreased to 125 nM. For Less BMP4 condition, BMP4 concentration was decreased to 50 ng mL-1.

RA treatment assay
To observe the effect of RA in µSCLS DV patterning in Fig. 3c, d, g, RA (500 nM, STEMCELL Technologies) was added to both dorsal and ventral reservoirs from day 5 to 12, along with 50 ng mL-1 BMP4 in the dorsal reservoir and 250 nM SAG in the ventral reservoir. To compare DV patterning in the absence of BMP4 gradient, from day 5 to 12 the dorsal reservoir was filled with basal media, and the ventral reservoir was filled with basal media supplemented with 250 nM SAG. 500 nM RA was further supplemented to both dorsal and ventral reservoirs from day 5 to 12. To compare DV patterning in the absence of SAG gradient in Supplementary Fig. 6b&c, from day 5 to 12 the ventral reservoir was filled with basal media and the dorsal reservoir was filled with basal media supplemented with 50 ng mL-1 BMP4. Throughout the µSCLS process, the basal medium was switched to N2B27 without vitamin A, which is identical to the standard N2B27 except for the use of B27 without vitamin A (Gibco).

CXCR4 inhibition assay
To inhibit CXCR4, AMD3100 (200 nM, Selleckchem) was supplemented into both reservoirs together with BDNF (10 ng mL-1; R&D Systems), GDNF (10 ng mL-1; PeproTech), CNTF (10 ng mL-1; PeproTech), IGF-1 (10 ng mL-1; PeproTech) , cAMP (1 μM; Sigma-Aldrich), and AA (0.2 μg mL-1; Sigma-Aldrich) from day 12 to day 18.

Whole-mount immunofluorescence staining of µSCLS
First, µSCLS were fixed in 4% paraformaldehyde (PFA; buffered in 1× PBS) at room temperature for 1 h 30 min. Then, µSCLS were exposed by manually detaching PDMS structural layers from glass coverslips. Thinly sliced PDMS strips were placed between the PDMS structural layer and glass coverslip to provide space for solutions. The µSCLSs were washed 3 times in PBS for 5 min each. µSCLSs were permeabilized in 0.1% SDS solution (SDS dissolved in PBS) at room temperature for 3 h. After washing 3 times for 5 min in PBS, µSCLS were blocked in 4% donkey serum (Sigma-Aldrich) in PBS at 4 °C for 24 h followed by incubation with primary antibody solutions at 4 °C for another 24 h. µSCLS were washed in PBS 3 times for 5 min and labelled with donkey-raised secondary antibodies (1:400 dilution, Fisher Scientific) at 4 °C for 24 h. Both primary and secondary antibodies were dissolved in 4% donkey serum with 0.2% NaN3. 4',6-diamidino-2-phenylindole (DAPI, Thermo Fisher Scientific) was used to counterstaining cell nuclei. Antibodies used are listed in Supplementary Table 4.

Whole-mount RNA fluorescence in situ hybridization (RNA-FISH)
Similar to the whole-mount immunofluorescence staining protocol, the µSCLS was first fixed in 4% paraformaldehyde (PFA; buffered in 1× PBS) at room temperature for 1 h 30 min. Then, µSCLS was exposed by manually detaching PDMS structural layers from glass coverslips. Thinly sliced PDMS strips were placed between the PDMS structural layer and glass coverslip to provide space for solutions. The µSCLS was washed 3 times in PBS for 5 min each. µSCLS was permeabilized in methanol in -20 ℃ for 30 min and in room temperature for 5 min. The µSCLS was washed once with PBSTRG and once with PBSTR. PBSTR consists of PBS with 0.1% Tween-20 and 0.1 U μL-1 Protector RNase Inhibitor (Cat. No. 3335399001; MilliporeSigma), and PBSTRG consists of PBSTR with 0.1 M Glycine. Then, the µSCLS was incubated in probes (Supplementary Table 6; Integrated DNA Technologies) diluted to 4 µM in HCRTM Probe Hybridization Buffer (v3.0; Molecular Instruments) overnight in 37℃. The µSCLS was washed 4 times, 5 min each, with HCRTM Probe Wash Buffer (v3.0; Molecular Instruments) preheated to 37 ℃. The µSCLS was further washed 2 times for 5 min each with 5X SSCT, which consists of 5X sodium chloride sodium citrate buffer with 0.1% Tween-20. Then, the µSCLS was incubated in HCRTM Amplifier Buffer (v3.0; Molecular Instruments) in room temperature for 30 min. In the meantime, HCRTM Amplifiers (v3.0; Molecular Instruments) were snap-cooled for formation of hairpin structure by heating them up to 95 ℃ for 90 sec and cooling down to 12 ℃ for 30 min. The snap-cooled HCRTM Amplifiers were mixed in HCRTM Amplifier Buffer at 1/100 (v/v) dilution right before use. The µSCLS was incubated in the HCRTM Amplifiers mixed in the HCRTM Amplifier Buffer for 90 min in room temperature, avoiding light. Then, the µSCLS was washed with 5XSSCT 3 times, 5 min each. Nuclei of µSCLS was stained with 4',6-diamidino-2-phenylindole (DAPI, Thermo Fisher Scientific) by incubating it in DAPI diluted in 5XSSCT at a ratio of 1/400 (v/v) overnight at 4℃, avoiding light. The µSCLS was finally washed with 5XSSCT 3 times, 5 min each, before proceeding to tissue clearing.

Tissue clearing
For downstream microscopy after whole-mount immunofluorescence staining or whole-mount RNA-FISH, µSCLS were optically cleared by incubating in a refractive index (RI)-matching solution. The solution comprised of 6.3 mL ddH2O , 9.2 ml OptiPrep Density Gradient medium (MilliporeSigma), 4 g N-methyl-d-glucamine (MilliporeSigma) and 5 g diatrizoic acid (MilliporeSigma)68. 10 µL of RI-matching solution was applied onto the µSCLS, and PDMS structure layer was attached back to the glass coverslip. The re-assembled sample was incubated in RT for at least 30 min before imaging.

EdU cell-proliferacion assay
Cell proliferation was assessed using the Click-iT EdU Alexa Fluor 488 Imaging Kit (Thermo Fisher Scientific), following the manufacturer’s protocol. On days 5 and 12, PDMS structural layers were manually detached from the glass coverslips, exposing the µSCLSs that remained on the PDMS surfaces. The µSCLSs were then incubated for 1 hour in basal medium containing 20 μM EdU, followed by fixation, permeabilization, and incubation with the Click-iT reaction cocktail for 30 minutes. Afterward, samples were blocked with 4% donkey serum in PBS and incubated with primary antibodies. DAPI was used to counterstain cell nuclei. EdU-labeled proliferating cells were visualized using confocal microscopy.

Fluorescence intensity maps
To generate heatmaps in Fig. 1g and Supplementary Fig. 3g-i, a custom MATLAB (R2024b) script was developed to process fluorescence images of lineage markers and generate fluorescence intensity maps. First, stitched images of µSCLS were manually rotated to align the dorsoventral (DV) axis vertically. Binary masks were then generated from DAPI images by selecting pixels with intensities above 20% of the maximum signal, effectively outlining the µSCLS contours. Pixel intensities outside the DAPI mask were automatically set to zero. The remaining pixel intensities were normalized to the maximum intensity within each image to produce normalized fluorescence intensity maps. Normalized signal values were then averaged across the thickness of the tissue to generate a one-dimensional intensity profile along its length. PAX6 staining images had high cytoplasmic signals at the ventral tip. This led to misleading expression profiles with high expression at the ventral end, where the thickness of the tissues taper off. Thus, normalized PAX6 staining signals were summed up, instead of being averaged, across the thickness of the tissue. Because the tissue lengths varied across samples, each profile was interpolated to 1,000 equidistant points using the MATLAB function interp1. Each interpolated intensity profile was subsequently normalized by its minimum and maximum values. Finally, normalized profiles from different tissues were averaged to produce an average normalized intensity heatmap.
To generate heatmaps in Supplementary Fig. 1j showing cell cycle markers along the apical-basal axis, confocal micrograph of individual µSCLS was cut in half along the lumen. The lumen side was marked as apical and the exterior as basal. Fluorescence intensity maps were generated similarly and was averaged over the length of the tissue. The resultant expression profile along the apical-basal axis was interpolated to 1,000 equidistant points using the MATLAB function interp1. Each interpolated intensity profiles were normalized by its minimum and maximum values and averaged across different tissues.

Boundary determination
To determine the expression domain boundaries, the intensity profile of each marker obtained as described in ‘Fluorescence intensity maps’ section was first filtered using a moving average filter with a window size of 50 to smooth out local fluctuations. Then, the point where the expression value becomes 50% of the maximum was determined as the boundary of expression domain.

scRNA-seq
µSCLS was dissociated into single cells by manually detaching the PDMS structural layers and incubating the released µSCLS in Accutase for 2 hours. The resulting cell suspension was collected into PBS containing 0.5% BSA and centrifuged at 200 × g for 5 minutes. Cell pellet was resuspended in PBS with 0.5% BSA and filtered through a 40 μm cell strainer (pluriSelect USA, 43-50040-51) to remove debris and aggregates. Within 1 hour of dissociation, cells were loaded into a 10X Genomics Chromium system (10X Genomics). 10X Genomics v.3 libraries were prepared following manufacturer’s instructions. Then, libraries were sequenced using paired-end sequencing with a minimum coverage of 20,000 raw reads per cell using an Illumina NovaSeqXPlus. scRNA-seq data were aligned and quantified using Cell Ranger Single-Cell software suite (v.7.1.0, 10X Genomics) against Homo sapiens (human) genome assembly GRCh38 (version refdata-gex-GRCh38-2020-A, 10X Genomics).

Analysis of scRNA-seq data
scRNA-seq data analysis were conducted using the R (ver. 4.2.2.) package Seurat (v.5.1.0, https://satijalab.org/seurat/)69. Unless noted otherwise, default setups in Seurat were used. Briefly, scRNA-seq dataset was filtered based on the total number of genes detected and the percentage of mitochondrial genes. Raw counts were normalized by the total counts, multiplied by 10,000, and log-transformed. To regress out cell cycle information, cell cycle scores were first calculated using CellCycle Scoring function. Then, the cell cycle scores were regressed out during the data-scaling process through ScaleData. For dimensional reduction, top 2,000 variable genes were first identified using FindVariableFeatures. Identified variable genes were used for PCA analysis using RunPCA. The principal components were used for embedding the data into low-dimensional space with uniform manifold approximation and projection algorithm using RunUMAP. Neighboring cells were identified in the principal component space using FindNeighbors and clusters were identified by performing a shared nearest neighbor modularity optimization-based clustering algorithm using FindClusters. Differentially expressed genes (DEGs) in each cluster were identified using FindAllMarkers. Dot plots in Fig. 2b, Supplementary Fig. 4g, and Supplementary Fig. 5e were generated using DotPlot. Heatmap in Supplementary Fig. 4a was generated using DoHeatmap. Feature plots in Fig. 4l and Supplementary Fig. 4c were generated using FeaturePlot. Heatmap in Supplementary Fig. 4b was generated by first averaging the scaled expression levels of each gene in each cluster and using ggplot in conjunction with geom_tile of the R package ggplot2. Heatmap in Supplementary Fig. 4e was generated by first averaging the mRNA counts of each HOX gene in each cluster and using pheatmap function of the R package pheatmap.

Comparison with human embryo data
To deduce the developmental stage of µSCLS, PCA analysis was done on average counts of each transcript. Then, Euclidean distances were calculated in the principal component space between µSCLS and human SCs from different developmental stages22. To compare the transcriptome of µSCLS and stage matched human SC in Fig. 2d, the development systems ‘Neural crest progenitor’, ‘Neuron’, and ‘Progenitor’ were isolated from CS12 SC dataset. ‘Peripheral neuron’, corresponding to the DRG population, was obtained from CS14 SC dataset because of its scarcity in the CS12 SC dataset. For comparison of NE cluster in Fig. 2f&g and Supplementary Fig. 4f&g, clusters ‘RP’, ‘dP1’, ‘dP2’, ‘dP3’, ‘dP4’, ‘dP5’, ‘dP6’, ‘p0’, ‘p1’, ‘p2’, ‘pMN’, ‘p3’, ‘FP’ were isolated from human CS12 SC dataset. For comparison of neuron clusters in Supplementary Fig. 5, ‘Neuron’ cluster was isolated from human CS12 SC dataset. The datasets were integrated using a reciprocal PCA approach based on 30 dimensions and k.anchor = 20 through IntegrateLayers. Annotations were based on their original annotation before integration. Pearson’s correlations between cell clusters in Fig. 2e&g and Supplementary Fig 5d were calculated using the cor function in the R package Stats based on their integrated principal components averaged over each gene in each cluster. Red and gray feature plots marking subtype localization in UMAP plots in Supplementary Fig. 4f and Supplementary Fig. 5b were generated using FeaturePlot. Color-graded feature plots showing average expression of marker genes in Supplementary Fig. 4f and Supplementary Fig. 5b were generated by first averaging the scaled expression of marker genes used to identify each subtype (Supplementary Table 2) and using FeaturePlot. Heatmap in Fig. 2i was generated by averaging the scaled expression levels of each gene and using ggplot in conjunction with geom_tile of the R package ggplot2. Violin plots in Fig. 2l were generated by first isolating NE cells co-expressing OLIG2 and NKX2.2, determined by a threshold of one UMI count and using VlnPlot.

NE and neuron subtype partitioning
Partitioning of NE and neuron population into subtypes was done as previously reported22,46,47. First, expression levels of marker genes were binarized with a threshold of one UMI count. Then, Euclidean distance from a cell to each ideal subtype expressing all the marker genes was calculated in the binary space. The cell was labeled as the closest subtype. Markers used are listed in Supplementary Table 2.

Re-analysis of mouse data
Raw count matrices and metadata of scRNA-seq data from mouse embryos at E9.5 were re-analyzed for Fig. 2j&k. NE subtypes ‘RP’, ‘dP1’, ‘dP2’, ‘dP3’, ‘dP4’, ‘dP5’, ‘dP6’, ‘p0’, ‘p1’, ‘p2’, ‘pMN’, ‘p3’, ‘FP’ were isolated from mouse E9.5 SC dataset 46. In Fig. 2j, raw counts were averaged within each NE subtypes and normalized by maximum mean expression. In Fig. 2k, expression of NKX2.2 and OLIG2 were determined by a threshold of one UMI count.

Re-analysis of Chip seq
Table of aligned peaks from retinoic acid receptor chromatin immunoprecipitation sequencing experiment done on motor neurons differentiated from mouse embryonic stem cells with or without retinoic acid was reanalyzed54. Among the list of genes closest to peaks, BMP signaling pathway regulatory genes were identified. List of BMP signaling pathway regulatory genes was found from the Molecular Signature Database (MSigDB)70–72 (Supplementary Table 3).

Tissue area and length
To measure projected area and length of tissues in Fig. 1e, brightfield images of µSCLS at different days were first rotated so that the DV axis was aligned with the vertical axis. Then, tissues were manually contoured using freehand selection tool in Fiji 2.16.0. The height and area in pixels from Fiji measurements were converted to length and projected area in µm and µm2.

Cell number and position measurement
To specify position and number of NCCs in Fig. 4, or that of TWIST1+ mesodermal cells in Supplementary Fig. 4d, the fluorescent images were first binarized using a threshold. Then, adjoining cells were separated using watershed process in Fiji 2.16.0. The locations of separated cells were identified using analyze particle functionality in Fiji 2.16.0. Due to the high cell density in the NE, proportion of cells in Fig. 3f was estimated using proportion of pixel numbers. First, DAPI image was binarized by using Otsu’s thresholding. The positive pixels were counted and further used as a mask for fluorescence images of OLIG2 and FOXA2. After DAPI masking, OLIG2 and FOXA2 images were binarized by using Otsu’s thresholding and positive pixels were counted. Proportion of cells positive for OLIG2 was estimated by the proportion of OLIG2-positive pixels with to the number of DAPI-positive pixels.

Nuclei long axis angle quantification
To measure the angle of long axes of nuclei from DAPI images in Supplementary Fig. 1f, the DAPI images were first binarized by setting a threshold. Then, adjoining nuclei were separated using watershed process in Fiji 2.16.0. Each nucleus was fitted to an ellipse using analyze particle functionality in Fiji 2.16.0. The long axis of the ellipse was used to estimate the long axis of the nucleus.

Western blot
First, µSCLSs were dissociated into single cells. To do so, PDMS structural layer was manually detached from glass coverslips to release µSCLS, and released µSCLSs were incubated with Accutase for 20 min. Then, pellets were centrifuged and lysed for 10 minutes using ice cold radioimmunoprecipitation assay (RIPA) buffer supplemented with 1x protease and phosphatase inhibitor cocktail (Thermo Scientific). The lysates were vortexed every 2 min during lysis. The lysates were then centrifuged at maximum speed for 15 min. Total protein concentrations were measured using BCA Protein Assay Kit following manufacturer’s instructions (Thermo Scientific). 20µg of total protein was loaded in each pocket of SDS-polyacrylamide gel electrophoresis gel (BioRad) along with Precision Plus Protein Dual Color Standards (BioRad). The gel was then transferred onto a polyvinylidene difluoride membrane. The membrane was blocked in 5% BSA (w/v) dissolved in TBS with 0.1% Tween (v/v) (TBS-T) for 2 hours, washed three times in TBS-T for 5 min, and incubated with primary antibody in 5% BSA in TBS-T at 4 °C overnight on a rotator. On the next day, the membrane was washed three times in TBS-T for 5 min and incubated with a horseradish peroxidase (HRP) antibody in 5% BSA in TBS-T at 4 °C for 2 hours on a rotator. The membrane was then washed three times in PBS-T for 5 min. Chemiluminescence signal was detected using ECL Western Blotting Detection Reagent (Cat. No. 32106; Thermo Scientific). Antibodies used are listed in Supplementary Table 4. Uncropped, unprocessed scan of representative blot is presented in Supplementary Fig. 8.

Quantitative polymerase chain reaction (qPCR)
To quantify gene expression in Supplementary Fig. 6d&h, we performed qPCR. First, µSCLS were lysed using Lysis Buffer RLT from RNeasy Mini Kit (Cat. No. 74104; Qiagen). Three tissues in one device were lysed together. Then, the lysates were homogenized using QIAshredder (Cat. No. 79654; Qiagen). RNA was extracted from the homogenized lysates using RNeasy Mini Kit (Cat. No. 74104; Qiagen) following manufacturer’s instructions. The RNA concentrations of resultant solutions were measured using a spectrophotometer (NanoDropTM One/OneC; Thermo Scientific). Then, the solutions were diluted in nuclease free water to equalize the RNA concentrations across samples within the same batch. Complementary DNA (cDNA) was synthesized from the RNA solutions in a thermal cycler (T100; Bio-Rad Laboratories, Inc.) using iScriptTM cDNA Synthesis Kit (Cat. No. 1708890; Bio-Rad Laboratories, Inc.) following manufacturer’s instructions. Then, qPCR was performed on CFX Connect Real-Time PCR Detection System (Bio-Rad Laboratories, Inc.) by loading mixtures of the cDNA solutions, primers (Supplementary Table 5; Integrated DNA Technologies), and iTaq Universal SYBR Green Supermix (Cat. No. 1725121; Bio-Rad Laboratories, Inc.). The ratio and running conditions were set according to manufacturer’s instructions of iTaq Universal SYBR Green Supermix (Cat. No. 1725121; Bio-Rad Laboratories, Inc.). qPCR measurement was repeated two times for technical repeat. The quantification cycle (Cq) values were used for quantifying the initial mRNA concentrations. First, the Cq values from two technical repeats were averaged. Then, the initial mRNA amount was calculated as 2−𝐶𝑞. Then, the mRNA amount was normalized by the geometric mean of two housekeeping genes, GAPDH and ACTB. The mRNA expression levels normalized by housekeeping genes were further normalized by expression levels in hESCs.

Knock down assay using small interfering RNA (siRNA)
Differentially expressed genes identified in Supplementary Fig. 6d were knocked down using siRNA (Supplementary Fig. 6e-g and Supplementary Table 7). Pre-designed siRNAs were purchased from Integrated DNA Technologies. Three siRNAs were provided per gene. Individual siRNAs were diluted in Nuclease-Free Duplex Buffer (Integrated DNA Technologies) to a concentration of 10 µM. Equal volumes of diluted siRNAs targeting the same gene were mixed. The siRNA mixture was further combined with LipofectamineTM RNAiMAX Transfection Reagent (Cat. No. 13778030; Thermo Fisher Scientific) and Opti-MEMTM (Cat. No. 31985062; Thermo Fisher Scientific) at a volumetric ration of 1:2:50. The mixture was incubated at room temperature for 5 min before being added to culture media, adjusting the final concentration of siRNA in media to be 30 nM. Note that the 30 nM was the total concentration of three different siRNAs targeting the same gene. siRNAs were added to culture media on days 5, 7, 9, and 11.

Dorsally and ventrally located cells
Dorsally located and ventrally located cells in Fig. 4 were determined by the cells’ DV location relative to the NE. Location of the dorsal and ventral poles of three tissues in one device were averaged to produce an averaged dorsal pole and an averaged ventral pole. The DV location of averaged ventral pole was defined as 0 and that of the averaged dorsal pole was defined as 1. Then, the locations of individual cells were projected onto the DV axis defined by averaged dorsal and ventral poles. Cells with DV positions of greater than 1 were defined as dorsally located and cells with DV position of less than 0.5 were defined as ventrally located.

Migration tracking
To track the migration of NCCs, µSCLSs were fabricated with SOX10::T2A-Cre lineage reporter hESC line. On day 12 and onward, the samples were observed daily for emergence of cells expressing ZsGreen with flattened morphology, indicating migrating NCC. When such cell was identified, the sample was imaged every 20 minutes using an inverted epifluorescence microscope (Zeiss Axio Observer Z1; Carl Zeiss MicroImaging) enclosed in an environmental incubator (XL S1 incubator, Carl Zeiss MicroImaging), maintaining the cell culture conditions at 37 °C and 5% CO2. Then, the migration of individual NCCs were tracked using the Fiji plugin TrackMate. A LoG detector was used on the ZsGreen images with estimated object diameter of 25 µm and quality threshold of 0.05 to locate the NCCs. Then, migration was tracked using a simple LAP tracker with linking max distance of 64 pixels, gap-closing max distance of 64 pixels, and gap-closing max frame gap of 2. The polar histogram in Fig. 4e was generated by first calculating the migratory direction. Migratory direction of each cell was represented by the angle between the displacement vector connecting the start and end point of a NCC and the DV axis. Then, the angles of all NCCs were pooled into a polar histogram using a MATLAB function polarhistogram.

Mean square displacement
Mean square displacement at each timepoint t in Fig. 4g was calculated as below:
MSD(𝑡) = 1/𝑁(𝑡) ∑ 𝑁(𝑡) 𝑖=1 |𝐱(𝑖)(𝑡)−𝐱(𝑖)(0) | 2 (1)
Where 𝑁(𝑡) is the number of cells tracked until time t, 𝐱(𝑖)(𝑡) is the position vector of cell i at time t, and 𝐱(𝑖)(0) is the position vector of cell i at time 0.

Random walk simulation
A 2-dimensional random walk simulation was performed using a custom MATLAB script. In the simulation, a particle randomly took a step in either of 4 orthogonal directions in the x-y plane, namely, -x, +x, -y, and +y directions. The number of particles were set to 120, which was the average number of tracks of migrating NCCs in each experimental sample. To determine the step size, the mean square displacement (MSD) of NCCs were fitted to a parabolic curve with the constant term set to 073:
𝑓(𝑡) = 𝑈𝑥2 +4𝐷𝑥. (2)
Then, considering that images were taken every 1/3 hours, step size was determined as below to match the diffusion coefficient observed in migrating NCCs:
𝛿 = √4𝐷/3 (3)
The simulation was performed 5 times.

Quantification of directionality and circular standard deviation
Directionality was quantified by dividing the distance between initial and final position of a cell by the total length of the cell’s migratory path. To calculate the circular standard deviation of migratory directions74, First, the mean resultant length 𝑅̅ was calculated as below:
𝑅̅ = √𝑆̅2 +𝐶̅2. (4)
Here,
𝑆̅ = 1/𝑁 ∑𝑁 𝑖=1 sin𝜃 ,𝐶̅ = 1/𝑁 ∑𝑁 𝑖=1 cos𝜃 . (5)
Then, the circular standard deviation was calculated as:
𝑠𝑑 = √−2ln𝑅̅. (6)

Statistical analysis
All experiments were conducted in at least three biological replicates, except scRNA-seq assays, which were performed in one independent experiment. Sample sizes are indicated in the figure legends. For comparisons between two datasets, P values were calculated using two-sided Student’s t-test. For analysis of significance of SMAD1/5 complex phosphorylation level increase in Fig 3g, a one-sample t-test with a hypothesized mean of 1. For comparison of zero-inflated data that does not follow normal distribution in Fig. 4j&n, a zero-inflated rank test was used75. For comparison of gene expression levels across more than two conditions in Supplementary Fig. 6h, a one-way ANOVA followed by a post-hoc Tukey’s RHSD test were performed. No statistical methods were used to predetermine sample size. Samples were randomly allocated to different experimental groups. Investigators were not blinded to allocation during experiments and outcome assessment. All statistical analyses were performed using R 4.4.0.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, Bok et al. present a microfluidic-based in vitro model of the human spinal cord, termed the microfluidic spinal cord-like structure (μSCLS), to investigate dorsoventral (DV) patterning driven by opposing gradients of morphogens. Building on their previously established human neural tube-like structure platform, the authors apply an antiparallel gradient of BMP4 (dorsalizing) and SAG (ventralizing) to induce DV patterning in μSCLS. Notably, they report the emergence and ventral migration of neural crest cells (NCCs) from the dorsal region in response to a chemoattractant cue, mimicking in vivo behavior. Furthermore, they perform single-cell transcriptomic analysis to demonstrate the presence of all 13 neuroepithelial (NE) subtypes along the DV axis of spinal cord, suggesting that the μSCLS model recapitulates key aspects of human embryonic spinal cord development. An additional and intriguing finding is the dorsalizing effect of retinoic acid (RA) in their system with the potential crosstalk with BMP signaling pathway, which appears to contradict prior in vitro findings but aligns with certain in vivo studies, as discussed in the manuscript. Overall, this study provides a valuable contribution to the field of spinal cord organoid models, with a clear focus on DV patterning. However, to enhance the novelty, rigor, and impact of the study, several critical points should be addressed:\n\nMajor Comment:\nA central claim of the manuscript is that the μSCLS faithfully recapitulates the full DV axis of the human spinal cord, including all 13 NE subtypes. While this is supported by single-cell transcriptomic data, spatial validation of these subtypes remains limited. Although the immunofluorescence results demonstrate spatial patterning, only a subset of markers is presented, and not all 13 NE subtypes along the DV axis are spatially confirmed. As a result, the manuscript’s core claim remains insufficiently substantiated at the spatial level.\nGiven that the μSCLS model exhibits a broader DV width compared to the previously reported human neural tube-like system, it provides a valuable opportunity to demonstrate more refined and extended DV patterning. Capitalizing on this advantage through enhanced spatial validation would further highlight the strength and novelty of the model.\nTo underscore the distinctiveness and impact of the current system, it would be essential to include spatial resolution of cell types across the DV axis using spatial mapping methods/techniques such as spatial transcriptomics, multiplexed or iterative immunofluorescence. The key hallmark of DV patterning in the spinal cord is the precise spatial organization of progenitor domains in response to opposing morphogen gradients. A more detailed spatial analysis would substantially strengthen the manuscript and clearly differentiate it from prior work.\n\nMinor Comments:\n1. In Figure 1h, the authors present regional immunofluorescence images (boxes I–III) showing SOX10, LMX1A, and FOXA2 expression. To aid interpretation and provide context, it would be beneficial to include an overview image of the entire tissue adjacent to these panels, similar to the layout used in Figure 1g. This would help readers appreciate the spatial relationships between regions and better assess the DV organization of the markers shown.\n2. On page 9, the authors note the presence of mesodermal cells, attributed to early WNT activation during neural induction and transient induction of a neuromesodermal progenitor (NMP) state. It would be informative to clarify the spatial distribution of these mesodermal populations although they are very small number of cells within the μSCLS. Specifically, do they localize preferentially to the dorsal or ventral regions, or are they interspersed throughout the tissue?\n3. The manuscript hints that the observed NCCs are of trunk identity but does not directly state the regional identity of the NCCs in AP axis. It would be helpful to specify the regional/axial identity of the NCCs; for example, if they express HOX genes consistent with a trunk-level identity. Inclusion of this information, either through transcriptomic data or immunostaining, would provide a clearer understanding of the axial identity of the NCCs and further validate the model’s relevance to trunk-level spinal cord development.", "reviewer_2": "This study establishes a human neural tube–like model using a microfluidic system to investigate dorsal-ventral (D–V) patterning under retinoic acid (RA) regulation and to examine neural crest cell (NCC) migration. While the work is technically sound and the system well executed, the manuscript would benefit from clarification on its novelty and deeper mechanistic insights to clearly distinguish it from the authors’ prior publications. The findings on RA–BMP crosstalk and NCC migration are intriguing but need further elaboration. Below are major and specific points for improvement.\n\nMajor Concerns：\nThe methodology and some key findings (D–V patterning) closely resemble the authors’ previous publications. The manuscript should explicitly highlight: What new biological insights are provided here (mechanistic understanding of RA–BMP interactions, detailed NCC migration dynamics). Whether the microfluidic system itself is novel, or if this work represents new applications.\nThe authors identify BMP downstream genes with RA binding sites, but it remains unclear:Which of these genes are functionally relevant for D–V patterning?Are these genes direct targets of both RA and BMP, or is their regulation indirect?\nFor the Neural Crest Cell Migration:Do NCCs migrate along the μSCL microfluidic scaffolds, or do they delaminate fully into the Matrigel/gelatin?Are there any fate changes in NCCs during or after migration?\n\nSpecific Concerns to Address:\n1.The abstract mentions “human-specific features,” but these are not clearly defined or highlighted in the results.\n2.More markers for the 11 neural progenitor subtypes should be included via staining to strengthen cell identity claims.\n3.Dextran diffusion is used to mimic morphogen gradients, but it may not fully replicate endogenous conditions—this limitation should be discussed.\n4.RA is known to regulate anterior–posterior (A–P) identity. How do the authors exclude A–P patterning effects in their D–V model?\n5.The “dorsal intermediate” region is defined solely by Pax3/Pax6 overlap. Additional markers would help refine and validate this domain.\n6.NMPs are typically caudal/posterior, yet the model claims thoracic identity. Did the authors observe any T/SOX2 double-positive cells in the model?\n7.The absence of V3 interneurons might reflect either low abundance or a true lack of specification. A specific V3 marker should be used to determine whether these progenitors are simply rare or genuinely missing.\n8.The manuscript does not address whether NCCs undergo fate transitions during migration. Clarifying this would strengthen the interpretation of NCC behavior.", "reviewer_3": "Bok et al. presents a microfluidics-based, human PSC-derived model recapitulating Dorsal-Ventral (DV) patterning of the spinal cord and neural crest cells (NCCs) migration. They then used the platform to study the role of Retinoic Acid (RA) in DV patterning and found a synergistic crosstalk between RA and BMP signaling. They further analyzed cellular dynamics of NCCs and demonstrated a role of CXCL12 in regulating migration.\nIn my view, this system represents the epitome of recent efforts in harnessing the power of stem cells to illuminate development. Not only can this system recapitulate essential features of DV patterning and NCC migration, but also it allows deep mechanistic investigations (live imaging, controllable parameters, etc.). Although many in vitro systems capturing spinal cord development have been published, including a recent one (Luo et al. 2025 Cell Stem Cell) that also focused on DV patterning, the system here offers controllable and quantitative powers that could be very useful for future mechanistic studies. I have a few comments regarding rigor:\n1. Although the authors showed a phenotypic interaction between RA and BMP and demonstrated that RA treatment increases SMAD1/5 phosphorylation, the proposed molecular mechanism regarding RAR binding to regulatory elements of BMP pathway-associated genes is highly speculative. I suggest the authors tone down the language used in the abstract “… clarifies the dual role of RA…”.\n2. Does the RA-BMP crosstalk need input from Shh signaling at all? The authors showed that RA produced a dorsalizing effect under BMP-SAG opposing gradients but had no effect on patterning when BMP gradient was withdrawn. Would RA still dorsalize patterning if SAG gradient is removed? This experiment could further probe RA’s effects on DV patterning.\n3. Could the authors speculate why their system acquires a thoracic identity? Only HOXC9 and HOXC10 were investigated, could the single cell data further support the thoracic identity?\n4. In the SOX10 lineage tracing experiment, the ZsGreen signals emerged after delamination. Could this delay be an artifact of the lineage tracing system? For example, following SOX10 expression, it takes time to express Cre and then fluorescence expression. So the ZsGreen signal might not reflect the onset of SOX10 expression."}, "author_response": "We would like to express our sincere gratitude to the reviewers for their dedicated efforts in evaluating our manuscript. Their thoughtful critiques and constructive suggestions have been instrumental in enhancing the overall quality of our work. We greatly appreciate the time and expertise that each reviewer has contributed, which has collectively had a significant impact on the improvement of our manuscript. In this revision process, we have conducted numerous new experiments to thoroughly address all comments and concerns raised by the reviewers. We sincerely hope that the current version of the manuscript meets the reviewers’ expectations and standards. For ease of review, all edits and revisions in the manuscript are highlighted in YELLOW.\nResponse to Reviewer #1\nGeneral comment: A central claim of the manuscript is that the μSCLS faithfully recapitulates the full DV axis of the human spinal cord, including all 13 NE subtypes. While this is supported by single-cell transcriptomic data, spatial validation of these subtypes remains limited. Although the immunofluorescence results demonstrate spatial patterning, only a subset of markers is presented, and not all 13 NE subtypes along the DV axis are spatially confirmed. As a result, the manuscript’s core claim remains insufficiently substantiated at the spatial level. Given that the μSCLS model exhibits a broader DV width compared to the previously reported human neural tube-like system, it provides a valuable opportunity to demonstrate more refined and extended DV patterning. Capitalizing on this advantage through enhanced spatial validation would further highlight the strength and novelty of the model. To underscore the distinctiveness and impact of the current system, it would be essential to include spatial resolution of cell types across the DV axis using spatial mapping methods/techniques such as spatial transcriptomics, multiplexed or iterative immunofluorescence. The key hallmark of DV patterning in the spinal cord is the precise spatial organization of progenitor domains in response to opposing morphogen gradients. A more detailed spatial analysis would substantially strengthen the manuscript and clearly differentiate it from prior work.\nResponse: We thank the reviewer for their critical and constructive feedback. Spatial transcriptomics would indeed offer a powerful approach to illustrate dorsal-ventral (DV) patterning within the microfluidic spinal cord-like structures (µSCLS). However, the thin geometry of the tissue (~150 µm width) presents a significant technical challenge for sectioning along the DV axis, which is necessary for the application of spatial transcriptomics.\nIn response to the reviewer’s suggestion, we employed additional antibodies and performed RNA fluorescence in situ hybridization (RNA-FISH) to further characterize neural progenitor domains and their spatial distribution in the µSCLS (Extended Data Fig. 2a-e & 2g and Fig. R1). Specifically, to identify dorsal intermediate domains, we conducted co-immunostaining for PAX7 and PAX6, where regions positive for both markers (PAX7+PAX6+) denote the dp4-dp6 domains (Extended Data Fig. 2b and Fig. R1b). Furthermore, co-immunostaining for NKX6-1, PAX6, and OLIG2 revealed a NKX6-1+PAX6+OLIG2– region dorsal to the OLIG2+ region, corresponding to the p2 domain (Extended Data Fig. 2c and Fig. R1c). Spatial mapping of DBX1 and DBX2 expression via RNA-FISH showed a DBX1+ domain within a broader DBX2+ region; co-expression of DBX1 and DBX2 indicates the p1 and dp6 identities, while a dorsal DBX2+DBX1– domain corresponds to the dp5 identity (Extended Data Fig. 2d&e and Fig. R1d&e). Additionally, immunostaining for NKX2-2 together with OLIG2 and FOXA2 identified the p3 domain, positioned between the OLIG2+ pMN domain and the FOXA2+ floor plate (FP) region (Extended Data Fig. 2g and Fig. R1f).\nAlthough simultaneous visualization of all neural progenitor domains within the µSCLS was not achieved in this revision, the observed spatial relationships among neural progenitor markers support the presence of dorsal, ventral, and intermediate progenitor domains in the µSCLS. These domains are correctly organized along the DV axis within the µSCLS, reflecting a spatial arrangement consistent with that seen in vivo in the spinal cord.\n[...]"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed my comments.", "reviewer_2": "All my concerns have been addressed and I have no more questions", "reviewer_3": "The authors have satisfactorily addressed all my concerns. I support publication of the paper in its current form."}, "author_response": ""}]

Decision letter