A gene expression atlas of a juvenile nervous system

Abstract

Although the fundamental architecture of metazoan nervous systems is typically established in the embryo, substantial numbers of neurons are added during post-natal development while existing neurons expand in size, refine connectivity, and undergo additional differentiation. To reveal the underlying molecular determinants of post-embryonic neurogenesis and maturation, we have produced gene expression profiles of all neuron types and their progenitors in the first larval stage (L1) of C. elegans. Comparisons of the L1 profile to the embryo and to the later L4 larval stage identified thousands of differentially expressed genes across individual neurons throughout the nervous system. Key neuropeptide signaling networks, for example, are remodeled during larval development. Gene regulatory network analysis revealed potential transcription factors driving the temporal changes in gene expression across the nervous system, including a broad role for the heterochronic gene lin-14. We utilized available connectomic data of juvenile animals in combination with our neuron-specific atlas to identify potential molecular determinants of membrane contact and synaptic connectivity. These expression data are available through a user-friendly interface at CeNGEN.org for independent investigations of the maturation, connectivity and function of a developing nervous system.

Figure 1.. Single-cell RNA sequencing profiles of the L1 hermaphrodite nervous system

A) Neurons in the C. elegans L1 larval nervous system. B) UMAP projection of 89,073 L1 neurons. Arrow denotes cell groups selected for sub-UMAP shown in panel D. C) UMAP projection of L1 nervous system colored by the generation time (embryonic or post-embryonic) of the neuron class. The majority of post-embryonic neurons grouped together. D) Sub-UMAP of the cluster denoted in panel B (arrow) shows clear separation of neuronal classes. E) UMAP as in B colored by expression of the post-mitotic neuron marker sbt-1. The two clusters with low sbt-1 expression within the boxed region correspond to HSN and VC. F) Enlargement of the boxed region in E, colored by co-expression of egl-5 and unc-86, which uniquely marks HSN. G) Enlargement of the boxed region in E, colored by expression of the transcription factor ham-2, which is restricted to HSN and VC. H) Confocal micrograph of a ham-2 fosmid reporter in late L1 larvae shows expression in VC neurons in the ventral cord and in the HSN neuron pair. Scale bar = 10 μm.

Figure 2.. Neuronal progenitor cells in the Q and T lineages

A) QL and QR lineages each giving rise to three neurons and two daughter cells that undergo apoptosis (X). B) Sub-UMAP of Q.pa and descendants, SDQ, AVM and PVM colored by neuron class. C) Sub-UMAP of the Q.pa lineages with post-mitotic neurons colored by maturation index. D) The T.p lineage produces two glial cells (phasmid sockets), three neurons (PVW, PHC, and PLN) and a cell death (X) on each side of the worm. E) PHATE plot showing annotations of T.pp and descendants. F) PHATE plot showing expression of lin-32 in the T.pp lineages. G) Confocal micrographs of an endogenous lin-32 reporter and a membrane bound mKate expressed in the T lineage by a seam cell promoter. H) PHATE plot showing expression of vab-7 in PVW (T.ppa) and a subset of progenitor cells. I) Confocal micrographs showing endogenously tagged vab-7 in combination with the membrane bound mKate in the T-lineage. Expression of vab-7 begins in PVW (T.ppa), its sister T.ppp, and the T.ppp descendants, including PHC. All scale bars = 10 μm.

Figure 3.. P-lineage derived progenitors and neurons

(A) Diagram of P-cell migrations to the ventral cord and subsequent cell divisions during the L1 stage. (B) Diagram of the P3-P8 lineages, each of which produce five motor neurons. (C) Quantification of HLH-3::tdTomato expressing cells counted during each collection period (hph, hours post hatch). (D) Confocal micrographs of HLH-3::tdTomato expressing cells in the ventral cord during each collection period. The first two images featuring the bilateral array of P-cells before they migrate to the ventral cord, show two different Z-planes in the same individual. (E) Sub-UMAP of P1-P12 lineage derived progenitors and neurons colored by cell type, (F) colored by the time in which cells were collected, (G) colored by expression of the cell cycle marker cdk-1, (H) colored by expression of the post-mitotic neuron marker sbt-1, (I) colored by expression of the transcription factor ceh-6, and (J) colored by expression of the transcription factor cnd-1. (K) Heatmap showing scaled expression values of 44 transcription factors differentially expressed between cell types in the P.aa sub-lineage.

Figure 4.. Expression of transcription factors, RNA-binding proteins and cell adhesion molecules in the larval nervous system

A) Dot plot showing expression of homeodomain transcription factors across 121 neuron types at the L1 stage. Neuron types are grouped by functional category (bottom). Transcription factors are clustered. Each dot is colored by scaled TPM value. The size of the dot represents the proportion of cells expressing that transcription factor. B) Cumulative distribution of various transcription factor families across 121 neuron types at the L1 stage. The number in parentheses denotes the total number of genes expressed from each category C) Quantification of the total number of genes expressed at threshold 2 in each neuron class at the L1 (left) and L4 (right) stages. Neuron types are divided by birthtime (embryonic vs post-embryonic). D) Quantification of the total number of transcription factors expressed at threshold 2 in each neuron class at the L1 and L4 stages. E) Quantification of the total number of RNA binding proteins expressed at threshold 2 in each neuron class at the L1 and L4 stages. F) Quantification of the total number of cell adhesion molecules expressed at threshold 2 in each neuron class at the L1 and L4 stages. All statistical tests were performed using linear models featuring the birthtime, stage, and number of cells per neuron class as covariates. Between group comparisons were performed on the estimated marginal means using the Tukey p-value adjustment for multiple comparisons. G) Combined box and jitter plots of gene expression stability from L1 to L4 for several gene families with neuronal function. This comparison is restricted to embryonically born neurons. Stability is measured for each gene by the fraction of L1-expressing neurons which still express the gene at L4. H) Combined box and jitterplots of Jaccard similarity between L1 and L4 for neuronal gene families. The Jaccard similarity index accounts for both loss and gain of expression between the L1 and L4 stages. Only a subset of significant comparisons are shown, designated by horizontal lines. For G-H), Kruskal-Wallis tests followed by Dunn’s multiple comparison test. * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001.

Figure 5.. Thousands of genes are differentially expressed across larval development

A) Heatmap showing the average log fold change for genes with significant differential expression between embryonic and L1 stages. Positive log fold change values (orange, red) reflect higher gene expression in the L1 whereas negative log fold change values (blue) reflect higher expression in the embryo. Both genes (columns) and neurons (rows) are clustered. Neuronal functional categories are depicted on the left. B) Heatmap as in A, but for differential expression between the L1 and L4. Positive fold change values (red, orange) reflect higher expression in L1, negative values reflect higher expression in L4 (blue). Y-axis labels include functional categories and developmental birthtime (embryonic vs postembryonic) for neurons. C) Dotplot of scRNA-seq data showing expression of hlh-32 in AVF and VB2 in L1 and L4. D) Confocal micrographs showing expression of endogenously GFP tagged hlh-32 in AVF in L1, but not L4. Note stable expression in the AVF sister cell VB2. Scale bar = 5 μm. E) Quantification of hlh-32 expression in AVF in L1 and L4/adult worms, Fisher’s Exact Test. F) Histogram showing the distribution of neuron classes in which a gene exhibited significant differential expression for the embryo vs L1 comparison. G) Histogram as in F, but for the L1 vs L4 comparison. H) Heatmap showing transcription factors (TFs) with significant TF activity (columns) across neurons (rows) between the embryonic and L1 stages. Warmer colors (red-yellow; positive numbers) denote increased expression of a TF regulon in L1. Cooler colors (blue; negative numbers) denote decreased expression of a TF regulon in L1. Neuronal functional categories are color coded, and TFs with known roles as terminal selectors (in any neuron) are denoted in dark blue. Dotted lines denote the increased activity of hlh-4 in ADL. TFs with significant changes across many neurons are indicated by black arrows. I) Heatmap as in H, showing transcription factors (TFs) with significant TF activity (columns) across neurons (rows) between the L1 and L4 stages. Neurons are color-coded by both functional categories and birthtime, and TFs with known roles as terminal selectors (in any neuron) are denoted in dark blue. TFs with significant changes across many neurons are indicated by arrows. J) Bar graph showing the fraction of instances of differential TF activity which also showed differential TF expression (black) or which lacked significant differential expression (white). Most cases of differential activity of a TF did not show differential expression of the TF at the transcript level.

Figure 6.. Developmental alterations in neuropeptide networks

A) Thresholded neuropeptidergic connectome (mid-range) showing the developmental pattern of connections between sending neurons (y-axis) and receiving neurons (x-axis). Core (conserved) connections and developmentally dynamic connections are color-coded. Only neuropeptide-receptor couples with functionally validated threshold of EC₅₀ < 500 nM binding in vitro^, were included. Neurons are grouped by functional category (Ph = Pharyngeal). B) Percentage distribution of the conserved connections between NPPs (Neuropeptide Precursor genes) and GPCRs (G-Protein Coupled Receptors) forming the core connectome. 13% of those connections do not share a common NPP-GPCR pair across development. The other 87% connections share at least one NPP-GPCR pair (64%) or all NPP-GPCR pairs (23%). C) Weighted neuropeptidergic connectomes, indicating for every cell-cell connection how many pairs of neuropeptide-receptors mediate the connections. D) Total number of degrees (y-axis, mid-range networks) in homologous neurons across development (x-axis) classified by functional categories. Bars are color-filled according to the subsets of core connections and developmentally dynamic neuropeptides (See key in B) of each neuron across development. E) Proportion of developmentally dynamic neuropeptide connections (y-axis) per neuron classified by embryonic development (colors and x-axis; “non-mature” refers to the HSN and VC neurons). Kruskal-Wallis test with Dunn’s correction, ****P < 0.0001, ***P < 0.001, *P < 0.025. F) Variation in mesoscale grouping across development. The same 4 groups (hubs core, sensory core, motor core, periphery) are conserved but some individual neurons change between groups through development. G) Network topologies of the indicated NPP/GPCR pairs. Empty circles represent neuropeptide expression, filled circles represent receptor expression. Local networks display restricted NPP and GPCR expression (≤50 neurons), pervasive networks display broad NPP and GPCR expression (>50 neurons), broadcaster networks display restricted NPP but broad GPCR expression, integrative networks display broad NPP but restricted GPCR expression. NLP-58 and NLP-49 = NPPs; TKR-1 and SEB-3 = GPCRs

Figure 7.. Differential expression of secreted and cell surface molecules based on patterns of membrane contact and synaptic connectivity

A) Heatmap showing enrichment of hmr-1 and syg-1 in neurons contacting AVK (left of vertical red line) compared to neurons not contacting AVK (right of vertical red line). T13C2.6 and tsp-21 (below horizontal red line) are enriched in cells not contacting AVK. Red arrows indicate genes with known binding partners expressed in AVK. B) Volcano plot of secreted and cell surface enrichment in neurons based on membrane contact. Each dot represents a gene tested for a given neuron. Positive log2 Fold Changes indicate enrichment in cells with membrane contact. Light gray dots represent instances with Benjamini-Hochberg adjusted p-values > 0.05. Dark grey dots represent instances with Benjamin-Hochberg adjusted p-values < 0.05. Black dots represent significant instances in which at least one known binding partner is expressed in the tested neuron (i.e., syg-1 is enriched in neurons that contact AVK, and the syg-1 binding partner sax-3 is detected in AVK). C) Scatterplot showing Benjamini-Hochberg corrected p-values for each gene/neuron combination testing expression enrichment in membrane contacted vs non-contacted cells in L1 (x-axis) and in L4 (y-axis). Genes only detected at one age are not shown. Pink: significant enrichment only in the L1, blue: significant enrichment only in the L4, green: significant enrichment in both stages. D) Smoothed density plot of the fraction of membrane-contacted neurons with synaptic connections in the L1 nerve ring. The blue curve shows the fraction of membrane-contacted neurons with any synaptic connections (median of 0.2). The green curve represents the fraction of neurons that receive postsynaptic output (median: 0.12), and the pink curve represents the fraction of membrane-contacted neurons that provide presynaptic input (median: 0.09). E) Top: Graphical representation of cell adhesion molecule (CAM) expression in cells with membrane contact only (left) and cells with synaptic connections (right). Bottom: CAM binding partners confirmed by biochemical assays. F) Cartoon representations of the SMB neurons that innervate the nerve ring, from WormAtlas. G) 3D reconstruction from NeuroScan of late L1 nerve ring featuring three neurons with membrane contact sites and synapses. SAADL (orange) has both membrane contact (red patches) and synaptic contacts (red triangles) with SMBVL (gray), whereas SIAVL (light blue) only has membrane contact (dark blue) with SMBVL. H) Heatmap showing normalized expression (TPM) of 7 CAMs enriched in SMB-non-synaptic adjacent cells (below horizontal line) compared to in SMB postsynaptic outputs (above horizontal red line). Red arrows denote genes encoding proteins with known binding partners expressed in SMB. I) Heatmap showing normalized expression of genes enriched in SAA-non-synaptic, adjacent neurons (14 genes, to the left of vertical red line) and genes enriched in SAA-synaptic partners (2 genes, to the right of the vertical red line). Red arrows denote genes encoding proteins with known binding partners expressed in SAA. J) Volcano plot showing enrichment of CAMs based on synaptic connectivity. Each dot represents one gene/neuron comparison. Negative log fold change values represent enrichment in non-synaptic adjacent cells, whereas positive log fold change denote enrichment in synaptically-connected cells. Light gray dots correspond to Benjamini-Hochberg adjusted p-values > 0.05. Dark grey dots represent significant enrichment (Benjamini-Hochberg adjusted p-value < 0.05), and black dots denote significant enrichment with known binding partners expressed in the tested neuron. K) Scatterplot showing -log10 transformed p-values for each gene-neuron combination in L1 (x-axis) and in L4 (y-axis). Genes detected at only one age are not shown. Significant cases (Benjamini-Hochberg adjusted p-value <0.05) of enrichment only in L1 (pink), significantly enriched only in L4 (blue), and cases enriched at both L1 and L4 (green). The p-value adjustments were made for multiple testing within each neuron.