Vinculin influences essential processes in enteric nervous system development and Hirschsprung disease pathogenesis

Vinculin (VCL), a linker between cells and their environment, has rarely been linked to disease. This study examines the role of VCL in the development of the enteric nervous system (ENS) and its relationship to Hirschsprung disease (HSCR). Using whole-genome sequencing and in vitro assays, we identified 4 VCL mutations associated with HSCR, most causing loss of function. Neural crest–specific Vcl knock-out mice (Vcl cKO) displayed ENS defects resembling short-segment HSCR, including partial colonic aganglionosis and abnormal gut musculature. Single-cell transcriptomics revealed dysregulation of genes involved in neuronal differentiation and MAPK signaling. Spatial RNA-seq revealed reduced ENS-mesenchyme interactions in Vcl cKO mice, accompanied by significant disruption of the Pleiotrophin (PTN) pathway; Ptn knock-out mice exhibited phenotypes similar to those of Vcl cKO mice, underscoring the importance of ENS-mesenchyme crosstalk. VCL works as a hub gene crucial for cell connection and signaling pathways essential for ENS formation. VCL deficiency subtly impacts various developmental stages and neighboring cells, cumulatively leading to a phenotype similar to short-segment HSCR. This research highlights the role of VCL in maintaining cellular interactions and signaling pathways, such as MAPK and PTN, which are crucial for ENS development and may inform therapeutic targets for ENS disorders.

Multiple loss-of-function mutations in VCL were found in patients with both syndromic and isolated HSCR. To elucidate the role of VCL in the pathogenesis of Hirschsprung disease (HSCR), we reanalyzed a whole-genome sequencing dataset comprising 94 patients with S-HSCR with comprehensive clinical records, which included 13 patients presenting with both HSCR and VSD, representing approximately 13% of the cohort (19). In total, we identified 7 de novo heterozygous mutations distributed across various exons of the VCL gene in 9 patients with HSCR (Figure 1A and Supplemental Data 1; supplemental material available online with this article; https://doi.org/10.1172/JCI198531DS1). Among these, A977P mutation is located within the functional domain of the muscle-specific isoform (isoform 1), which is not expressed in NCCs (Supplemental Figure 1), it is plausible that this mutation does not contribute to NCC-associated defects.

Figure 1

HSCR-associated mutations in VCL interrupt focal adhesion assembly. (A) The schematic shows that mutations affect different functional domains of the VCL protein. (B) Immunofluorescent stains show the focal adhesions (FAs) as marked by colocalization of Paxillin and Vinculin in HeLa cells overexpressing various VCL mutants. Scale bars: 10 μm. (C) Bar graph shows the average FA size found in cells overexpressing various VCL mutants (mean ± SEM, > 50 cells were analyzed; ***: P < 0.001, 1-way ANOVA). co-IP assay indicates that M209L substitution in VCL (D) enhances head and tail interaction, (E) leading to reduced binding affinity to phosphorylated paxillin.

We then evaluated the potential impacts of the remaining 6 VCL mutations on focal adhesion (FA) assembly by measuring FA sizes, as described in previous studies (17). In brief, expression constructs carrying either WT or mutant VCL were overexpressed in the human cervical cancer cell line HeLa, and FAs were detected using immunocytochemistry with the FA marker Paxillin. The sizes of the FAs were measured and compared. Notably, 4 out of 6 mutations were found to be deleterious, leading to a reduction in FA size, primarily located within the head or tail domains of VCL (Figure 1, B and C, and Table 1). The M209L substitution was identified in 3 nonconsanguineous syndromic patients with HSCR presenting with VSD and exhibited a prominent effect on FA assembly. Therefore, we further investigated the mechanism by which the M209L substitution disrupts FA assembly. Specifically, this mutation is situated in head domain 1 (Vh1) of VCL, and it may promote autoinhibitory head-to-tail interactions. To test this hypothesis, we generated expression constructs incorporating either the head (amino acids 1–258) or tail (amino acids 879–1,066) regions of VCL, tagged with FLAG or c-MYC, respectively. Both WT and M209L mutant head constructs, along with the c-MYC–tagged tail construct, were coexpressed in HeLa cells, followed by coimmunoprecipitation assays. Our findings revealed that the M209L substitution enhances the head-to-tail interaction, likely resulting in VCL inactivation (Figure 1D). Moreover, VCL mutants exhibited a diminished binding affinity to phosphor-Paxillin (Figure 1E), which may account for the observed FAs defects. Importantly, impaired FA assembly may reduce the ability of cells to respond to external stimuli and interact with their environment, thereby disrupting various cellular processes.

Table 1

A summary of VCL mutations identified in patients with HSCR and their impacts on FA size.

Neural crest–specific ablation of VCL caused partial colonic aganglionosis. To directly investigate the roles of VCL in the development of ENCCs, both copies of the Vcl gene were specifically deleted from NCCs using the Wnt1-Cre mouse line, resulting in Vcl conditional knock-out (cKO) mutants. These mutants exhibited both cardiac and enteric nervous system (ENS) defects and did not survive beyond a few hours after birth. In terms of ENS phenotypes, nearly all E18.5 mutants analyzed (n > 20) showed varying degrees of colonic aganglionosis, as determined by acetylcholinesterase (AChE) staining. As demonstrated in Figure 2A, the ENS network did not fully cover the distal colon in the Vcl cKO mutants, whereas the ENS in the proximal and middle colon appeared unaffected and closely resembled that of the control group. We subsequently performed IHC analyses on whole-mount colon preparations, consistently revealing aganglionic regions in the distal colon of the Vcl cKO mutants (Figure 2B). Interestingly, the neurons and glial cells adjacent to these aganglionic areas were not properly organized nor distributed as ganglia, in contrast to the arrangement observed in the control (Figure 2C). We then quantified the number of HuD⁺ neurons and SOX10⁺ glial cells in the distal colon. In the control colon, both neurons and glial cells were evenly distributed with comparable density across 4–8 random views (600 mm² each) from 3 embryonic guts. The counts of neurons and glia were normalized against the total number of cells per view, as indicated by DAPI staining. In summary, the control group exhibited 41.84% ± 1.39% neurons and 22.65% ± 1.07% glia, resulting in a neuron-to-glia ratio of 1.89 ± 0.09. In contrast, in the comparable regions, the mutants displayed a lower percentage of neurons (34.56% ± 3.80%) and a higher percentage of glia (29.29% ± 3.07%), along with greater variability. The neuron-to-glia ratio in the mutants was significantly reduced to 1.20 ± 0.12 (Figure 2C).

Figure 2

Malformation of enteric ganglia in Vcl cKO at E18.5. (A) Three regions (proximal: C1; middle: C2 and distal: C3) of colon were collected from control and Vcl mutants (Vcl KO^Wnt1–Cre) for whole-mount acetylcholinesterase (AchE) staining. Scale bars: 250 μm. (B) Whole-mount immunofluorescence of TUJ1, HuD and SOX10 shows ENS network in the distal colon of E18.5 Vcl cKO was aberrantly organized. Scale bars: 100 μm.(C) Magnified images of the regions highlighted in the red dotted boxes in B. Scale bars: 25 μm. The percentages of neurons (HuD⁺) and glia (SOX10⁺) normalized to the total number of cells (DAPI⁺) and the relative neuron-to-glia ratio in each region are shown in the bar charts. (mean ± SEM, n: number of embryonic guts analyzed, ***: P < 0.001, Student t test, 2-tailed). 6–8 randomly selected regions were analyzed.

We then generated Vcl cKO with a YFP background (Wnt1-Cre; Rosa26^YFP;Vcl^f/f), allowing for the labeling of ENCCs with YFP. Additional IHC analyses were performed on sections of the distal colon from E18.5 control and mutant specimens. Consistent with the data obtained from whole-mount staining, we observed a significant reduction in the percentage of committed neurons (PHOX2B⁺SOX10^–YFP⁺), accompanied by an increase in the number of uncommitted ENCCs (PHOX2B⁺SOX10⁺YFP⁺) (Figure 3A). In sum, Vcl mutants displayed partial aganglionosis and exhibited an immature ENS.

Figure 3

Immature ENS in Vcl cKO at E18.5. Immunofluorescence of (A) PHOX2B & SOX10; and (B) SMA and YFP counterstained with DAPI on cross sections of E18.5 distal colon, and the quantitative data are shown in bar graphs. (mean ± SEM; n, number of embryonic guts analyzed; P < 0.05 was considered significantly different; ***P < 0.001, Student t test, 2-tailed). 6–8 randomly selected regions were analyzed. Mes, mesenchyme; e, endoderm; lu, lumen; ci, circumferential muscle; ol, outer longitudinal layer; il, inner longitudinal layer; m.p., myenteric plexus. Scale bars: 10 μm (A); 50 μm (B).

Upon analyzing the crosssections of the mutant guts at E18.5, we observed that the smooth muscle in the Vcl cKO colon was less densely packed and less organized, with a notable reduction in the thickness of the mesenchymal layer (Figure 3B). Since Vcl was specifically deleted in NCCs, this is likely a secondary consequence of the ENS defect, suggesting that the ENS plays a crucial role in the development and patterning of the mesenchyme and smooth muscle.

Aberrant migration of ENCCs in Vcl mutants. We reasoned that the incomplete colonization of the distal colon in mutants would be the result of migration defects of ENCCs. Thus, we collected embryonic guts from E11.5 to E13.5, which is the key developmental window for gut colonization, and analyzed the migration of ENCCs in the ex vivo guts during these stages. A slight but significant delayed migration of ENCCs was observed in E11.5 Vcl cKO embryonic gut (Figure 4A).

Figure 4

Failure to form the migration chain in Vcl mutants. (A) The migration of ENCCs in hindgut and colonization in cecum were examined at E11.5. White dashed lines indicate the migration distance. Bar graph shows significantly delayed migration of ENCCs in the mutant mice. P < 0.05 was considered significantly different and marked by *, Student t test, 2-tailed. Scale bar: 100 μm. (B) Time-lapse images from live imaging of ENCC migration in E12.5 embryonic guts, where YFP-labelled ENCCs located at the migratory wavefront (asterisks) were tracked and their migratory paths were indicated by arrows. Bar graph shows the net speed of migration. Scale bar: 200 μm. (C) Whole mount immunofluorescence of YFP-labeled ENCCs in E13.5 hindguts. Failure in formation of migratory chains and presence of solitary ENCCs (arrowheads) were observed in Vcl cKO. P < 0.5 was considered significantly different and marked by “*”, Student t test, 2-tailed. Scale bar: 100 μm.

To investigate their migration patterns, live imaging was conducted over a period of 12 to 15 hours using E12.5 hindguts to observe the behavior of YFP-labeled ENCCs in both control and mutant conditions. The control ENCCs exhibited a distinctive migration pattern, moving in chains that often converged and diverged as they elongated distally to invade the uncolonized regions of the gut. The nodes formed at the intersections of these ENCC chains, along with the interconnecting chains themselves, are proposed to play a role in guiding the arrangement of ganglia into a lattice-like neuronal network (20). ENCCs in the control gut remained interconnected as continuous strands throughout their migration (Figure 4B). In contrast, the migration of ENCCs was disrupted in Vcl cKO mutants. Although these mutant cells could migrate toward the distal end of the hindgut, they displayed an inconsistent migration trajectory. More significantly, the mutant cells failed to establish a migratory chain at the migratory front, resulting in a number of solitary ENCCs (Figure 4B and Supplemental Video 1). We then further monitored the behavior of ENCCs at the migratory front. While the migration direction of ENCCs in both Vcl cKO and control groups was similar, with a tendency to move toward the distal end of the hindgut, the mutant ENCCs consistently deviated from the primary migratory path (Supplemental Figure 2A). Additional analysis of their migration tracks revealed that the average speed of cell migration was somewhat reduced in the mutants (Supplemental Figure 2B), whereas their persistence (defined as the ratio of net distance traveled to total distance traveled) remained comparable with that of control cells (Supplemental Figure 2C). Alongside the irregular migratory patterns observed in VCL-deficient ENCCs, the overall net migration speed was also diminished (Figure 4B). By E13.5, a distinct neuronal meshwork had formed in the control gut; however, the mutant ENCCs were unable to maintain connections with neighboring cells, leading to a significant presence of solitary ENCCs at the distal end of the hindgut (Figure 4C). This observation indicates that the connections between ENCCs were disrupted during migration, which likely interferes with subsequent ganglionogenesis.

Reconstruction of the differentiation paths of ENCCs reveals that loss of VCL delays their differentiation along the Branch A. In addition to the migration defect, the mutant exhibited a reduced number of mature neurons and an increased presence of uncommitted ENCCs, as illustrated in Figure 2C and Figure 3A. This suggests a potential differentiation defect in the mutant ENCCs. Notably, delayed differentiation was observed as early as E13.5, with the differentiation defect becoming more pronounced by E15.5 (Supplemental Figure 3). Therefore, we analyzed the transcriptomes of ENCCs that were isolated and enriched through fluorescence-activated cell sorting (FACS) from 7 guts of control (Wnt1-Cre; Rosa26^YFP) and Vcl mutant (Wnt1-Cre; Rosa26^YFP;Vcl^f/f) embryos at E13.5 and E15.5 (Supplemental Figure 4), each from 2–3 litters. This was performed using single-cell RNA-seq with 10X Genomics to explore further how the loss of Vcl influences the molecular dynamics of ENCCs along their differentiation trajectories. Additionally, we carried out 10X Visium spatial transcriptomics (ST) analysis on sagittal sections of E13.5 embryos to investigate how alterations in signaling among various ENS cells and their neighboring cells contribute to the observed defects in the mutants. By cross referencing our scRNA-seq and ST-RNA-seq datasets, we aimed to clarify the molecular mechanisms by which VCL governs ENCC development and their interactions with adjacent cells. Human induced pluripotent stem cell (hiPSC) ENS-based and mouse models were then used for functional validation. A schematic summarizing our analysis pipeline is presented in Figure 5A.

Figure 5

scRNA-seq analysis revealed a delayed progression of Vcl mutant cells along the neuronal lineage differentiation trajectory. (A) The workflow summarizing scRNA-seq data analysis and corresponding functional validations. (B) The integrated UMAP (Uniform Manifold Approximation and Projection) projection of 30,157 cells ENCCs at E13.5 and E15.5, colored by cell types. (BP, biopotent progenitor; GP, glial progenitor; Branch A & Branch B neurons, annotated by integrating the data from Morarach et al (46); ENMFB, enteric mesothelial fibroblast.) (C) The expression of marker genes and proliferative markers across cell types. (D) The proportion of cell types across different samples. Empirical Bayes moderated t test is used for the cell proportion comparison, with Benjamini-Hochberg–adjusted P value < 0.05 considered significant. The proportion of E13.5 BP, E13.5 Branch A, and E15.5 Branch A neurons is affected. (E) Comparison of RNA velocity in control and Vcl mutant. (F) The density plot shows the distributions of pseudotime of cell states inferred by Slingshot, revealing a delayed differentiation of neurons. (G) Expression of Branch A markers in control and Vcl mutant. 2-sided Wilcoxon rank-sum test, ****P < 0.0001.

After performing quality control on the scRNA-seq data, we identified a total of 30,157 cells, with each cell exhibiting an average of 3,677 genes and 13,860 unique molecular identifiers (UMIs) (refer to Supplemental Figure 5 and Supplemental Data 2). To refine the differentiation trajectories, we integrated our dataset with a previously published dataset of E18.5 ENS cells 21 and reannotated them based on the expression of canonical marker genes and lineage-specific transcriptional factors identified in the original study (Supplemental Figure 6). The Uniform Manifold Approximation and Projection (UMAP) plots revealed 5 cell clusters comprising 3 distinct differentiation branches: neuronal, which includes inhibitory (Branch A) and excitatory (Branch B) neurons, and glial lineages (Figure 5B). All paths originated from the highly proliferative bipotent progenitors (BP) characterized by high expression of Mki67, Ube2c, Cdc20, and Ccnb1. The cells on the neuronal differentiation trajectory progressed through a neuronal intermediate stage, Neuroblast, which is marked by early neuronal markers (Tubb3^high/Elavl4^high/Cartpt⁻/Prph⁻), before diversifying into 2 branches distinguished by the complementary expression of Etv1 and Bnc2, representing Branch A and Branch B neurons, respectively. By E15.5, Branch A constituted the predominant neuronal population, coexpressing markers of inhibitory neurons like Nos1 and Vip, while a smaller subset expressed Bnc2 with a reduction in Etv1 expression, corresponding to Branch B excitatory neurons. In terms of the glial lineage differentiation trajectory, glial progenitors (GPs) displaying high levels of Sox10 and Fabp7 began to emerge at this stage. Additionally, a distinct population of progenitors that expressed unique marker sets referring to the enteric mesothelial fibroblasts (ENMFBs) was identified (Figure 5, B and C, and Supplemental Data 3).

We subsequently examined the impact of Vcl loss on cell composition. Both the control and mutant groups contained all 5 clusters of neural cells and the ENMFB cluster. Interestingly, the Vcl mutant showed an increased proportion of BP cells at E13.5, coupled with a decrease in Branch A neurons at both E13.5 and E15.5 (Figure 5D). Further RNA velocity analyses revealed that the differentiation trajectories in both control and Vcl mutant populations were similar (Figure 5E), indicating that the absence of Vcl does not disrupt the differentiation of ENCCs toward a specific lineage, nor does it introduce a differentiation bias that would lead to a reduction in Branch A.

To investigate the heterogeneity of cells within the BP-to-Branch A lineage, we aimed to order the cells in pseudotime and infer the trajectories along this differentiation path. We began by performing a Principal Component Analysis (PCA) to reorganize all the cells in the BP-to-Branch A lineage. The distribution of cells projected to principal components 1 (PC1) and 3 (PC3) illustrated a continuous pseudotemporal trajectory from the BP to Branch A. Consequently, we inferred the pseudotime based on these 2 dimensions (Supplemental Figure 7, A–C). As shown in Figure 5F, both control and mutant cells from E13.5 and E15.5 exhibited a well-ordered arrangement along the pseudotime of differentiation. In the Vcl mutant, the BP cells at E13.5, Neuroblasts at E15.5, and Branch A neurons demonstrated significantly lower pseudotime values compared with the control cells, indicating a delayed differentiation along the Branch A neuronal lineage and immaturity of the Branch A neurons (Figure 5F and Supplemental Figure 7D). Consistently, a significant reduction in the expression of Branch A marker genes, coupled with elevated expression of the proliferative marker Mki67, was observed, suggesting fewer mature Branch A neurons in Vcl mutant (Figure 5G).

In accordance with the RNA-seq data, a significantly greater number of ENCCs remained proliferative (Ki67⁺) in the mutants at E18.5, while only postmitotic ENCC derivatives were observed in the control group at this stage (Figure 6A). This finding further suggests the immature nature of the ENS in the mutants. Consistent with this observation, there were significantly fewer mature neurons (HuD⁺) and a larger population of immature neurons (HuD⁺SOX10⁺) in the mutants compared with the controls (Figure 6B). Notably, the inhibitory lineage (nNOS⁺) appeared to be the most severely impacted (Figure 6C), which corresponds to neurons in Branch A.

Figure 6

Increased progenitor cells and delayed differentiation in Vcl mutants. (A) Mitotic ENS neurons (Ki67⁺PHOX2B⁺YFP⁺); (B) mature (HuD⁺SOX10^–) and immature (HuD⁺SOX10⁺) neurons; and (C) inhibitory (nNOS⁺) neurons. Bar graphs show the quantitative data (mean ± SEM; n, number of embryonic guts analyzed; *P < 0.05; ***P < 0.001, Student t test, 2-tailed; scale bars: 10 μm (A); 25 μm (B, C)).

The overarching roles of VCL in the neuronal lineage differentiation of ENCCs. We aimed to elucidate the mechanisms responsible for the delayed neuronal differentiation and maturation. To that end, we conducted further analyses on the 4 most impacted cell states: BP at E13.5, Neuroblast at E15.5, and Branch A neurons at both E13.5 and E15.5 (Figure 7A). We identified a total of 3,016 differentially expressed genes (DEGs) (FDR < 0.01), which included 35 and 14 consistently downregulated and upregulated genes, respectively, along the BP-to-Branch A trajectory, while 2,967 genes showed a dynamic expression profile. To identify the most biologically relevant DEGs, we employed a 2-tier approach that considered both changes in expression levels and the biological roles of the genes. First, we conducted Gene Ontology (GO) enrichment analysis on the DEGs across each cell state, categorizing them according to biological processes. We then focused on the top 50 GO pathways in each cell state and calculated the corresponding pathway scores for both control and mutant samples (Supplemental Data 4). The 10 most significantly disrupted biological processes within each cell state, comprising 640 DEGs, were selected for further examination based on changes in pathway scores. Subsequently, we examined how these 640 DEGs could influence cellular progression by reclustering them alongside their potential driving genes, specifically transcription factors (TFs), utilizing a gene regulatory network (GRN) inference strategy (Figure 7B and Supplemental Data 5). Among the 640 DEGs, 447 were identified as potential target genes of E2f1, Egr1, and Klf7 based on expression correlation and motif binding analyses. These target genes were utilized to construct a GRN, which was organized into 4 distinct modules based on their regulatory relationships. We also recorrelated the genes within these modules with their respective cell states by analyzing the dynamic expression profiles of 3 hub transcription factors: E2f1, Egr1, and Klf7, in control and Vcl mutant cells along the differentiation trajectory (Supplemental Figure 8). Notably, E2f1 exhibited significant differential expression primarily in the BP, while Klf7 showed disruption specifically in the Branch A neurons. In contrast, expression of Egr1 decreased consistently across the entire trajectory from BP to Branch A in mutant cells, with the greatest difference from the control cells in Neuroblast (Supplemental Figure 8, A and B). Similar dynamic regulons were also observed consistently (Supplemental Figure 8C). Consequently, we designated Module 1 to BP, Module 2 to Neuroblast, and Modules 3 and 4 to the Branch A neurons (Figure 7C). To assess the impact of Vcl loss on these module genes, we then calculated the module scores for both control and mutant cells. The results indicated that all modules exhibited significantly diminished signals in their corresponding cell states in Vcl mutants (Figure 7C), suggesting that the absence of Vcl disrupts the cellular processes governed by these module genes.

Figure 7

Disrupted gene regulatory network and pathways in BP-to-Branch A lineage. (A) The density plot shows the pseudotime distributions of 4 disrupted cell states. E13.5 BP, E15.5 Neuroblast, E13.5 Branch A neurons, and E15.5 Branch A neurons. (B) Gene regulatory network of the DEGs along the BP-to-Branch A lineage, colored by regulatory modules defined by regulatory relationship (Module1, genes regulated by E2f1 and Egr1; Module2, genes regulated by Egr1 only; Module3, genes regulated by Egr1 and Klf7; Module4, genes regulated by Klf7 only). Lines connecting Vcl represent the potential regulation and known protein-protein interactions from the STRING database. (C) The boxplot shows the module score in the corresponding disrupted cell states. (D) The dot plot shows the Spearman correlation on the expression of Vcl and genes in modules. The top 10 genes ranked by correlation are labeled and shown in B. The dot plots show (E) the biological processes enriched by the genes in the modules and (F) KEGG pathways enriched by the DEGs identified in Neuroblast and BranchA neurons.

In line with these observations, genes within modules showed stronger expression correlation with Vcl compared with other genes (Supplemental Figure 8D). Moreover, VCL demonstrated protein-protein interactions and coexpression relationships with numerous genes across all 4 modules (Figure 7, B and D). These results suggested that VCL mediates the genes in these modules directly or indirectly to govern various cellular processes along BP-to-Branch A differentiation. For instance, integrin β1 (Itgb1) within Module 2 is essential for ENCC migration and subsequent ganglionogenesis (20–22). Another gene, transgelin 2 (Tagln2), encoding an actin-binding protein that is linked to the migration and proliferation of tumor cells (23), also interacts directly with VCL at the protein level, which may suggest that the cell migration process is affected in Vcl mutants. In Module 3, VCL interacts with catenin α-2 (CTTNA), which links cadherin adhesion receptors and the cytoskeleton to regulate cell-cell adhesion and differentiation in the nervous system (24). In addition, Alpha actinin-1 (ACTN1) is a cross-linking protein that interacts with F-actin, playing a crucial role in anchoring actin to various intracellular structures. Morphological changes are essential for initiating neuronal differentiation; thus, VCL likely interacts with these proteins, governing the morphological changes necessary to support subsequent neuronal differentiation.

Additionally, the GO enrichment analysis of the genes within the modules highlighted several disrupted biological processes. Notable findings include neurogenesis within BP (Module 1), cell junction assembly in both BP and Neuroblast modules (Modules 1 and 2), morphology-related regulation and actin filament organization affecting both neuroblasts and Branch A neurons (Modules 2 and 3), and synaptic vesicle function specifically within Branch A neurons (Module 4) (Figure 7E). According to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, the upregulated genes were associated with neuronal disorders while the significantly downregulated genes were involved in MAPK and Rap1 signaling pathways, which are critical for neuronal differentiation (25) (Figure 7F and Supplemental Data 6). In summary, the results indicate that VCL deficiency disrupts multiple cellular processes involved in neuronal lineage differentiation. This disruption begins at the early stages, affecting cell morphogenesis and cell-cell adhesion, which, in turn, impacts the BP and Neuroblast stages. Additionally, the maturation of Branch A neurons is notably hindered, with the MAPK and RAP1 signaling pathways being the most significantly affected. Ultimately, these changes contribute to a delay in the formation of mature Branch A neurons.

To evaluate the direct involvement of VCL in the neuronal lineage differentiation of ENCCs, we established an in vitro model using hiPSCs. First, an inducible CRISPR/Cas9–hiPSC (iCas9-hiPSC) line was used to specifically knock down VCL in ENCCs or committed neuronal progenitors (NPs), where the expression of Cas9 protein was induced by the addition of doxycycline (DOX) (Figure 8A). A 2-step differentiation protocol was used to model ENCC-to-neuron differentiation. iCas9-hiPSCs were first directed to the NCC lineage by dual-SMAD inhibitors (LDN193189 and SB431542) with a WNT agonist (CHIR99021) and then caudalized with retinoic acid (RA) to obtain posterior/vagal NCCs (hereafter referred to as hENCCs), which were further enriched using FACS and characterized based on the expression of the NC-specific surface markers (HNK1, p75^NTR, CD49 and SOX10) (Supplemental Figure 9), as described previously (17, 26). The neuronal lineage differentiation of hENCCs was then induced by culturing the hENCCs in neuronal differentiation medium supplemented with various neurotrophic factors (see Supplemental Methods). The VCL gene was knocked down in hENCCs or committed NPs by transfecting cells with VCL-targeting gRNAs and in the presence of doxycycline on day 12 or day 15, respectively (Figure 8A). By day 30 of differentiation, the neuronal differentiation capacity of the cells was monitored based on the expression of neuronal markers (neurofilament, NF and Protein gene product 9.5, PGP9.5). When VCL was knocked down (KD) at NCC stage on day 12, the mutant cells showed a weaker ability to aggregate together, and this severely abolished the subsequent neuronal lineage differentiation. As a result, significantly fewer cells were obtained at day 30 of differentiation in VCL KD group, and, among them, there was a lower percentage of cells expressing neuronal markers (NF⁺ or PGP9.5⁺, a marker for clinical diagnosis) compared with the control group (Figure 8B and Supplemental Figure 10A). Similarly, even if we bypassed the neuronal initiation step and KD VCL on day 3 after culturing ENCCs in neuronal differentiation medium, the percentage of neurons was significantly lower in the KD group, while the total number of cells was more comparable (Figure 8C and Supplemental Figure 10B). Our scRNA-seq analysis revealed the MAPK pathway as the most significantly disrupted signaling pathway in VCL-deficient cells, particularly in Neuroblast and Branch A neurons. This disruption of the pathway was consistently observed in hENCC-derived neurons when VCL was knocked down (Figure 8D). These findings suggest that VCL is essential for the activation of the MAPK pathway, which is crucial for initiating neuronal lineage differentiation of ENCCs and their subsequent maturation.

Figure 8

VCL is implicated in multiple stages of neuronal lineage differentiation of ENS cells. (A) Schematic shows the human iPSC-based model of ENS development. An inducible-iPSC line was used, where the addition of doxycycline induced the expression of Cas9. Cotransfection of VCL-specific gRNA drove the deletion of VCL gene in FACS-enriched ENCCs (Day 12) or ENCCs committed to the neuronal lineage (Day 15). Immunocytochemistry with antibodies against pan-neuronal markers (NF and PGP9.5) on cells with VCL knockdown on (B) day 12 and (C) day 15. Bar charts show the percentage of neurons marked by NF or PGP9.5 over the total number of DAPI⁺ cells. Data shown in the bar graphs are mean ± SEM from 3 independent experiments. Student t test, 2-tailed. (D) Western blot analyses of ITGB1, MAPK, phospho-MAPK and VCL. ACTIN was used as the loading control. The relative expression levels were shown in the bar graphs (mean ± SEM from 3 independent experiments). Student t test, 2-tailed. Scale bars: 200 μm.

Spatial transcriptomic analysis reveals disruption of cell-cell interactions among ENCCs and with the gut mesenchyme in Vcl cKO. VCL is essential for cell-cell interactions and our scRNA-seq data showed perturbation of integrins and cell junction assembly in Vcl mutant ENCCs. Thus, we also examined whether loss of Vcl in ENCCs perturbs the communications between ENS cells and their neighborhoods, which may interrupt the ENS development. To this end, we conducted spatial-RNA-seq (ST-seq) analysis on sagittal sections of E13.5 embryos, utilizing replicate tissue sections spaced approximately 16 microns apart. The sequencing of these samples was performed to a median depth of 171,916,081 reads (with an interquartile range of 153,220,340 to 200,680,799), resulting in an average of 6,845 genes and 26,110 unique molecular identifiers (UMIs) per spot (Supplemental Figure 11A). The gut regions were manually delineated based on histological images in each section. We then integrated and jointly analyzed the replicate sections from the control and mutant groups to cover more independent regions of the gastrointestinal tract. Subsequently, deconvolution was employed to estimate the cell type composition for each ST data spot, using a published scRNA-seq dataset of embryonic guts at the same developmental stage as a reference. Through this analysis, ST spots were categorized into specific gut regions according to the relative expression of genes characteristic of the large intestine (Fxyd4, Muc2, Ntm, Fabp1, Cdx2, Satb2), small intestine (Tff3, Tdo2, Lum, Gpr50, Agr2, Sulf1), and stomach (Barx1, Sox2, Gata4, Igf1, Nkx2-5) (Supplemental Figure 11, B and C). Using a similar approach, the spots were further annotated as epithelial (EPI), mesenchymal (MES), and neural crest (NC) cells. Only cells that expressed YFP (YFP⁺) within these regions were classified as NC (Supplemental Figure 11, D and E).

Among the 4 sections, we focused on Control 1 and Mutant 2, which had more comparable numbers of various cell types and covered a larger region of the large intestine compared with the other 2 sections. In particular, the ENS phenotypes of Vcl cKO were only observed in the distal colon, so the subsequent analyses were restricted to those regions identified as the large intestine in these 2 sections at similar spatial locations (Figure 9A). Within the selected regions, comparable interactions among EPI and MES were found in the control and mutant (Supplemental Figure 12, A and B), while no common communication between NC and MES was detected. It is likely attributed to the specific deletion of Vcl in NCCs, abolishing the communications between NC and their neighboring cells in the mutant. We therefore focused on the interactions within NC and other cell types (Supplemental Data 7). The overall interactions related to NCCs in the Vcl mutant exhibited a reduction in both quantity and interaction strength (Figure 9B). Among the 20 interrupted pathways detected, NCAM and PTN pathways exhibited the highest probability of contributing to the reduced interactions found in the mutant (Figure 9C and Supplemental Figure 12C). To further analyze the disrupted cell-cell interactions, we categorized communications by cell-type pairs and identified the top 3 interactions ranked by probability in each cell-type pair (Figure 9D). Among the signal flows from NC to MES, 2 interactions related to PTN (Ptn-Sdc2/3) were significantly affected, which aligns with the observed reduction in PTN expression in ENCCs at E18.5 (Figure 10A). The interaction between Fn1 and the Itga5+Itgb1 pair was found to be the most impacted signal flow from MES to NC, exhibiting the highest probability. Additionally, a decrease in ITGB1 expression was detected in ENS cells at E18.5 (Figure 10B). Limited NC-NC interactions were detected from the ST data, likely due to the sparse number of NC spots available. Therefore, we investigated NC-NC interactions using scRNA-seq data from E13.5 (Supplemental Figure 13). Notably, we consistently observed reduced Cdh2-Cdh2 interactions in Vcl cKO mutants across both the scRNA-seq and ST-seq datasets, alongside diminished N-cadherin (CADH2 encoded by Cdh2) expression levels in NC cells (Figure 10, C and D). In summary, Vcl deficiency disrupts a crucial signaling pathway PTN (27, 28) and affects adhesion molecules (NCAM and CADH2) (20–22) that are vital for cell-cell interactions, cell migration, and subsequent gangliogenesis. The compromise of these signaling and adhesion molecules impairs cell-cell communication and interactions between NC and neighboring cells, ultimately hindering the development of the ENS.

Figure 9

Disrupted cell-cell interactions in Vcl cKO revealed by spatial transcriptomics of embryonic guts at E13.5. (A) The spatial annotation of the selected regions of the large intestine at E13.5. (EPI, epithelial; MES, mesenchymal cells; NC, neural crest cells.) (B) Bar graphs illustrate the total number and strength of NC-related cell-cell interactions. (C) The number and strength of NC-related interactions across different pathways are presented. (D) The top 3 disrupted cell-cell interactions in each group, ranked by probability. All of these interactions are entirely abolished in the mutant within the corresponding cell-type pairs.

Figure 10

Downregulation of PTN, ITGB1 and CADH2 in Vcl-deficient ENS at E18.5. IHC using antibodies against (A) PTN; (B) ITGB1; and (C) N-CADHERIN in the E18.5 distal colon of control and Vcl cKO mutants. ENS cells were YFP-labeled. (D) A Western blot was conducted with whole gut lysate from control and mutants. The relative CADH2 expression levels are shown compared to the control. GAPDH served as the loading control. Scale bars: 25 μm (A and B, left); 10 μm (A and B, right); 15 μm (C).

Loss of Ptn perturbed the formation of circumferential smooth muscle cell layer in Vcl cKO and Ptn KO. Our spatial RNA-seq data revealed multiple disruptions in cell-cell communication between ENCCs and MES cells in the Vcl cKO colon. While the functions of NCAM and CADH2 in the ENS are well established, the role of PTN remains less defined. Our findings suggest that PTN is the most significantly affected pathway mediating communication between ENCCs and MES cells. Further expression analysis demonstrated that PTN is expressed in both ENCCs and gut MES cells; notably, its expression in ENCCs increases with developmental stages, whereas in MES cells, the highest expression is found at E15.5 (Figure 11A and Supplemental Figure 14). Notably, both ENCCs and MES cells exhibited downregulation of PTN in Vcl cKO mutants across all developmental stages examined. We next sought to understand the biological functions of PTN signal in ENS and gut development. The Ptn^–/– mutants consistently exhibited ENS defects resembling the Vcl cKO mutants, but with less severity. They included hypoganglionic colon (Figure 11B), a disorganized ENS network of reduced neuron-to-glia ratio (Figure 11, C and D) and immature ENS with reduced percentage of postmitotic neurons (Figure 11E). Intriguingly, PTN is essential for the smooth muscle lineage differentiation of MES cells, which contribute to the formation of the 2 circumferential smooth muscle layers of the colon: the lamina muscularis interna and lamina muscularis externa. Ptn^–/– mutants at E15.5 exhibited a thinner lamina muscularis interna and a delayed spatial separation between the 2 smooth muscle layers. A similar phenotype was observed in Vcl cKO mutants, where only a few Calponin-positive (CNN1-positive) smooth muscle cells were present in the lamina muscularis externa at this stage (Figure 11F). This phenotype likely results from impaired smooth muscle differentiation.

Figure 11

Dynamic expression of PTN in ENS and mesenchyme is essential for the formation of the muscle layer of the gut. (A) IHC of PTN in E13.5, E15.5 & E18.5 of control (Wnt1-Cre; Rosa26^YFP) and Vcl cKO (Vcl KO^Wnt1–Cre) guts. (B) Whole-mount acetylcholinesterase (AchE) staining of control and Ptn KO colon. (C) Whole-mount immunofluorescence of TUJ1, HuD, and SOX10 shows ENS network aberrantly organized in the distal colon of E18.5 Ptn mutant. The average neuron-to-glia ratio in the 4–5 selected regions of distal colon is shown in the bar graphs. Immunofluorescence of (D) PHOX2B and SOX10 and (E) Ki67 counterstained with DAPI on cross sections of E18.5 distal colon. (F) Immunofluorescence staining for Calponin (red) and counterstained nuclei with DAPI (blue) in E15.5 Control (Ptn^+/+) and Ptn^−/− colons. m.p., myenteric plexus; Lu, lumen; lamina muscularis interna (yellow arrowhead) and the lamina muscularis externa (green arrowhead). Bar graphs show the quantitative data (mean ± SEM; n, number of embryonic guts analyzed, *P < 0.05; ***P < 0.001, Student t test, 2-tailed). Scale bars: 25 um (A); 250 um (B); 20 um (C, E, and F); 10 um (D).

LL and ML performed bioinformatics analyses under the supervision of ESWN and AJ, respectively. XW analyzed the mouse mutants. PHYC and KKYW provided the clinical data. MMGB conducted the genetic analyses. PL and CYY performed in vitro studies. ESWN supervised the whole project. ESWN and LL prepared the manuscript.

• Seed funding from the University of Hong Kong.
• GRF grants from the Research Grants Council of Hong Kong Special Administrative Region, China Hong Kong (HKU17111824 &17107322) to ESWN.

View Supplemental data

View Supplemental data set 1

View Supplemental data set 2

View Supplemental data set 3

View Supplemental data set 4

View Supplemental data set 5

View Supplemental data set 6

View Supplemental data set 7

View Unedited blot and gel images

View Supplemental video 1

View Supporting data values

We thank Frank Lai, Pang Li, Kathy Lui, Reeson Wang, Kasey Lau, and Ophelia Liu for establishing the iCas-iPSC line, conducting some of the in vitro and live-imaging studies. Confocal imaging and RNA sequencing were performed using equipment maintained by the Li Ka Shing Faculty of Medicine Faculty Core and the Centre for PanorOmic Sciences, the University of Hong Kong, respectively.

Address correspondence to: Elly SW Ngan, Department of Surgery, School of Clinical Medicine, The University of Hong Kong, Faculty of Medicine Building, 21 Sassoon Road, Pokfulam, Hong Kong, China. Phone: 852.2831.5410; Email: engan@hku.hk.

Conflict of interest: The authors have declared that no conflict of interest exists.

Copyright: © 2025, Liu et al. This is an open access article published under the terms of the Creative Commons Attribution 4.0 International License.

Reference information: J Clin Invest. 2026;136(3):e198531.https://doi.org/10.1172/JCI198531.

1. Montalva L, et al. Hirschsprung disease. Nat Rev Dis Primers. 2023;9(1):54.
   View this article via: CrossRef PubMed Google Scholar
2. Heanue TA, et al. Enteric nervous system development and Hirschsprung’s disease: advances in genetic and stem cell studies. Nat Rev Neurosci. 2007;8(6):466–479.
   View this article via: CrossRef PubMed Google Scholar
3. Nagashimada M, et al. Autonomic neurocristopathy-associated mutations in PHOX2B dysregulate Sox10 expression. J Clin Invest. 2012;122(9):3145–3158.
   View this article via: JCI CrossRef PubMed Google Scholar
4. Anderson RB, et al. The cell adhesion molecule l1 is required for chain migration of neural crest cells in the developing mouse gut. Gastroenterology. 2006;130(4):1221–1232.
   View this article via: CrossRef PubMed Google Scholar
5. Faure C, et al. Gangliogenesis in the enteric nervous system: roles of the polysialylation of the neural cell adhesion molecule and its regulation by bone morphogenetic protein-4. Dev Dyn. 2007;236(1):44–59.
   View this article via: CrossRef PubMed Google Scholar
6. Fu M, et al. BMP signaling regulates murine enteric nervous system precursor migration, neurite fasciculation, and patterning via altered Ncam1 polysialic acid addition. Dev Biol. 2006;299(1):137–150.
   View this article via: CrossRef PubMed Google Scholar
7. Bays JL, et al. Vinculin in cell-cell and cell-matrix adhesions. Cell Mol Life Sci. 2017;74(16):2999–3009.
   View this article via: CrossRef PubMed Google Scholar
8. Thievessen I, et al. Vinculin-actin interaction couples actin retrograde flow to focal adhesions, but is dispensable for focal adhesion growth. J Cell Biol. 2013;202(1):163–177.
   View this article via: CrossRef PubMed Google Scholar
9. Plotnikov SV, et al. Force fluctuations within focal adhesions mediate ECM-rigidity sensing to guide directed cell migration. Cell. 2012;151(7):1513–1527.
   View this article via: CrossRef PubMed Google Scholar
10. Case LB, et al. Molecular mechanism of vinculin activation and nanoscale spatial organization in focal adhesions. Nat Cell Biol. 2015;17(7):880–892.
    View this article via: CrossRef PubMed Google Scholar
11. Johnson RP, et al. An intramolecular association between the head and tail domains of vinculin modulates talin binding. J Biol Chem. 1994;269(17):12611–12619.
    View this article via: CrossRef PubMed Google Scholar
12. Johnson RP, et al. F-actin binding site masked by the intramolecular association of vinculin head and tail domains. Nature. 1995;373(6511):261–264.
    View this article via: CrossRef PubMed Google Scholar
13. Chen H, et al. Coincidence of actin filaments and talin is required to activate vinculin. J Biol Chem. 2006;281(52):40389–40398.
    View this article via: CrossRef PubMed Google Scholar
14. Biswas R, et al. Mechanical instability of adherens junctions overrides intrinsic quiescence of hair follicle stem cells. Dev Cell. 2021;56(6):761–780.
    View this article via: CrossRef PubMed Google Scholar
15. Zemljic-Harpf AE, et al. Heterozygous inactivation of the vinculin gene predisposes to stress-induced cardiomyopathy. Am J Pathol. 2004;165(3):1033–1044.
    View this article via: CrossRef PubMed Google Scholar
16. Marg S, et al. The vinculin-DeltaIn20/21 mouse: characteristics of a constitutive, actin-binding deficient splice variant of vinculin. PLoS One. 2010;5(7):e11530.
    View this article via: CrossRef PubMed Google Scholar
17. Lai FP, et al. Correction of hirschsprung-associated mutations in human induced pluripotent stem cells via clustered regularly interspaced short palindromic repeats/Cas9, restores neural crest cell function. Gastroenterology. 2017;153(1):139–153.
    View this article via: CrossRef PubMed Google Scholar
18. Gui H, et al. Whole exome sequencing coupled with unbiased functional analysis reveals new Hirschsprung disease genes. Genome Biol. 2017;18(1):48.
    View this article via: CrossRef PubMed Google Scholar
19. Tang CS, et al. Identification of genes associated with hirschsprung disease, based on whole-genome sequence analysis, and potential effects on enteric nervous system development. Gastroenterology. 2018;155(6):1908–1922.
    View this article via: CrossRef PubMed Google Scholar
20. Druckenbrod NR, et al. Behavior of enteric neural crest-derived cells varies with respect to the migratory wavefront. Dev Dyn. 2007;236(1):84–92.
    View this article via: CrossRef PubMed Google Scholar
21. Chevalier NR, et al. A neural crest cell isotropic-to-nematic phase transition in the developing mammalian gut. Commun Biol. 2021;4(1):770.
    View this article via: CrossRef PubMed Google Scholar
22. Lake JI, et al. Enteric nervous system development: migration, differentiation, and disease. Am J Physiol Gastrointest Liver Physiol. 2013;305(1):G1–24.
    View this article via: CrossRef PubMed Google Scholar
23. Shi J, et al. Transgelin-2 contributes to proliferation and progression of hepatocellular carcinoma via regulating Annexin A2. Biochem Biophys Res Commun. 2020;523(3):632–638.
    View this article via: CrossRef PubMed Google Scholar
24. Schaffer AE, et al. Biallelic loss of human CTNNA2, encoding αN-catenin, leads to ARP2/3 complex overactivity and disordered cortical neuronal migration. Nat Genet. 2018;50(8):1093–1101.
    View this article via: CrossRef PubMed Google Scholar
25. RD Y, et al. Rap1 mediates sustained MAP kinase activation induced by nerve growth factor. Nature. 1998;392(6676):622–626.
    View this article via: CrossRef PubMed Google Scholar
26. Li Z, et al. Transcriptomics of Hirschsprung disease patient-derived enteric neural crest cells reveals a role for oxidative phosphorylation. Nat Commun. 2023;14(1):2157.
    View this article via: CrossRef PubMed Google Scholar
27. Lu KV, et al. Differential induction of glioblastoma migration and growth by two forms of pleiotrophin. J Biol Chem. 2005;280(29):26953–26964.
    View this article via: CrossRef PubMed Google Scholar
28. Polykratis A, et al. Characterization of heparin affin regulatory peptide signaling in human endothelial cells. J Biol Chem. 2005;280(23):22454–22461.
    View this article via: CrossRef PubMed Google Scholar
29. Olson TM, et al. Metavinculin mutations alter actin interaction in dilated cardiomyopathy. Circulation. 2002;105(4):431–437.
    View this article via: CrossRef PubMed Google Scholar
30. Vasile VC, et al. A missense mutation in a ubiquitously expressed protein, vinculin, confers susceptibility to hypertrophic cardiomyopathy. Biochem Biophys Res Commun. 2006;345(3):998–1003.
    View this article via: CrossRef PubMed Google Scholar
31. Graneli C, et al. Histopathological dimensions differ between aganglionic and ganglionic bowel wall in children with Hirschsprung’s disease. BMC Pediatr. 2022;22(1):723.
    View this article via: CrossRef PubMed Google Scholar
32. Broders-Bondon F, et al. N-cadherin and β1-integrins cooperate during the development of the enteric nervous system. Dev Biol. 2012;364(2):178–191.
    View this article via: CrossRef PubMed Google Scholar
33. Wang RX, et al. Vinculin mediates a dynamic signal flow underlying cardiac outflow tract development as revealed by single-cell and spatial transcriptomic analysis. [published online October 28, 2025]. J Adv Res. https://doi.org/10.1016/j.jare.2025.10.062.
    View this article via: CrossRef Google Scholar
34. Papadimitriou E, et al. On the role of pleiotrophin and its receptors in development and angiogenesis. Int J Dev Biol. 2022;66(1-2-3):115–124.
    View this article via: CrossRef PubMed Google Scholar
35. Polykratis A, et al. Characterization of heparin affin regulatory peptide signaling in human endothelial cells. J Biol Chem. 2005;280(23):22454–22461.
    View this article via: CrossRef PubMed Google Scholar
36. Zhou BY, et al. Identification of signaling pathways that specify a subset of migrating enteric neural crest cells at the wavefront in mouse embryos. Dev Cell. 2024;59(13):1689–1706.
    View this article via: CrossRef PubMed Google Scholar
37. Li G, et al. PTEN deletion leads to up-regulation of a secreted growth factor pleiotrophin. J Biol Chem. 2006;281(16):10663–10668.
    View this article via: CrossRef PubMed Google Scholar
38. Li YS, et al. Pleiotrophin gene expression is highly restricted and is regulated by platelet-derived growth factor. Biochem Biophys Res Commun. 1992;184(1):427–432.
    View this article via: CrossRef PubMed Google Scholar
39. Tang CS, et al. Uncovering the genetic lesions underlying the most severe form of Hirschsprung disease by whole-genome sequencing. Eur J Hum Genet. 2018;26(6):818–826.
    View this article via: CrossRef PubMed Google Scholar
40. Li MX, et al. A comprehensive framework for prioritizing variants in exome sequencing studies of Mendelian diseases. Nucleic Acids Res. 2012;40(7):e53.
    View this article via: CrossRef PubMed Google Scholar
41. Gonzalez F, et al. An iCRISPR platform for rapid, multiplexable, and inducible genome editing in human pluripotent stem cells. Cell Stem Cell. 2014;15(2):215–226.
    View this article via: CrossRef PubMed Google Scholar
42. DeKelver RC, et al. Functional genomics, proteomics, and regulatory DNA analysis in isogenic settings using zinc finger nuclease-driven transgenesis into a safe harbor locus in the human genome. Genome Res. 2010;20(8):1133–1142.
    View this article via: CrossRef PubMed Google Scholar
43. Lau ST, et al. Activation of hedgehog signaling promotes development of mouse and human enteric neural crest cells, based on single-cell transcriptome analyses. Gastroenterology. 2019;157(6):1556–1571.
    View this article via: CrossRef PubMed Google Scholar
44. Li P, et al. αE-catenin inhibits a Src-YAP1 oncogenic module that couples tyrosine kinases and the effector of Hippo signaling pathway. Genes Dev. 2016;30(7):798–811.
    View this article via: CrossRef PubMed Google Scholar
45. Liu JA, et al. Identification of GLI mutations in patients with hirschsprung disease that disrupt enteric nervous system development in mice. Gastroenterology. 2015;149(7):1837–1848.
    View this article via: CrossRef PubMed Google Scholar
46. Morarach K, et al. Diversification of molecularly defined myenteric neuron classes revealed by single-cell RNA sequencing. Nat Neurosci. 2021;24(1):34–46.
    View this article via: CrossRef PubMed Google Scholar