Alternative splicing across the C. elegans nervous system

Abstract

Alternative splicing is a key mechanism that shapes transcriptomes, helping to define neuronal identity and modulate function. Here, we present an atlas of alternative splicing across the nervous system of Caenorhabditis elegans. Our analysis identifies novel alternative splicing in key neuronal genes such as unc-40/DCC and sax-3/ROBO. Globally, we delineate patterns of differential alternative splicing in almost 2000 genes, and estimate that a quarter of neuronal genes undergo differential splicing. We introduce a web interface for examination of splicing patterns across neuron types. We explore the relationship between neuron type and splicing, and between splicing and differential gene expression. We identify RNA features that correlate with differential alternative splicing and describe the enrichment of microexons. Finally, we compute a splicing regulatory network that can be used to generate hypotheses on the regulation and targets of alternative splicing in neurons.

Fig. 1. Overview of data collection and splicing analysis.

A Schematic of the experimental procedure. B–E Four methods to analyze alternative splicing, applied to the gene ric-4/SNAP25 in the neurons NSM and PVM. B Raw data visualization. Top track: Gene model of ric-4 (along with non-coding RNAs 21ur-13262 and Y22F5A.10). Bottom tracks: Read coverage and junction counts. Numbers denote junction-spanning reads for the splice junctions of interest. Blue and orange boxes indicate alternative first exons. C Local Splicing Variation (LSV) visualization. Top: Gene model of ric-4, highlighting the alternative first exons corresponding to ric-4a (blue) and ric-4b (orange). Bottom: Posterior mPSI (MAJIQ-defined Percent Selected Index) estimates displayed as violin plots. Numbers indicate the posterior expected mPSI (summing to one in each neuron, as the orange and blue junctions are mutually exclusive). D Transcript-level quantification. Top: Average transcript TPM across all sequenced neurons. Middle: Relative transcript usage (in proportion to the total TPM for the gene, in each neuron), in NSM and PVM. Bottom: Transcript quantification, indicating the mean +/- standard deviation of TPM across samples for each transcript in each neuron. E Neuron-wise comparison between PVM and NSM. Top: The application returns a list of splice junctions with differential usage between the two neurons, including two junctions belonging to the same LSV in the gene ric-4. Bottom: Visualization of junction usage for this LSV in neurons PVM and NSM.

Fig. 2. Comparison to previous splicing data.

A Exon 5 skipping in elp-1 in AVM neuron. Left: Raw data visualization in the neuron AVM. Right: Local quantification of the two LSVs corresponding to the cassette exon. The violin plots represent the posterior distribution computed by MAJIQ. B–D Cassette exon 11.5 in daf-2 in individual neuron types. B daf-2 gene structure. The cassette exon 11.5 (red box) is unique to the transcript daf-2c. C Local quantification of the upstream (top) and downstream (middle) events flanking the cassette exon 11.5. Left: Schematic representation of splice junctions constituting the local event. Right: Splice junction quantifications in 16 neuron types. Bottom: Inclusion pattern of exon 11.5 from Tomioka et al., and agreement with our data. D Raw data visualization of exon 11.5 alternative splicing in 5 neuron types. Top: Gene model. Bottom: Exonic and junction-spanning counts. Right: Exon 11.5 inclusion pattern from Tomioka et al. in the same neurons. E Transcript quantification of daf-2 in the same neurons as (C). The transcripts daf-2a, d, e, and f, (other) which do not include exon 11.5, were grouped together. Right: Exon 11.5 inclusion pattern from Tomioka et al. .

Fig. 3. Detection of novel cassette exon in unc-40.

A Alternative expression of a novel cassette exon in unc-40. Top: Schematic of the gene structure, showing the novel alternative exons 8.5 and 14.5 (green) and the annotated constitutive exons (gray). Bottom: Genome browser representation around exon 14.5, displaying the maximum track, the AVM track, and samples collected from whole animals. B Local quantification of exon 14.5 inclusion. The violin plots represent the posterior distribution computed by MAJIQ. C Validation by RT-PCR of the novel exon 14.5 in cDNA from whole animals (N2 wild type). Top: Location of the primers relative to the alternative exon 14.5. Left, the primer pair only results in amplification from a transcript with exon 14.5 included; right, the primer pair flanking exon 14.5 is able to amplify from both transcripts. Arrowheads indicate the expected product sizes, green arrowhead: transcript with exon 14.5 included (left: 754 bp; right: 837 bp); blue arrowhead: transcript with exon 14.5 skipped (582 bp). For each experiment, the amplification was performed on three biological replicates along with a control without reverse transcriptase (RT). Molecular weights reported in base-pairs. Source data are provided as a Source Data file. D Protein structure of UNC-40, indicating the cassette exons positions vs known protein domains. E Structure of the unc-40 orthologous genes in Caenorhabditis briggsae (Cbr-unc-40) and C. brenneri (Cbn-unc-40). Top: Annotated gene models (unannotated features in green). Bottom: RNA-Seq data. Red arrowheads denote orthologs of exon 8.5 and 14.5. F Alignment of the protein sequences of C. elegans UNC-40, C. briggsae UNC-40, and C. brenneri UNC-40 around exon 8.5 orthologs (top) and around exon 14.5 orthologs (bottom).

Fig. 4. Detection of novel alternative first exon in sax-3.

A Structure of the gene sax-3/ROBO indicating the novel exon 5.5 and alternative 5’ splice site at the end of exon 9 (green), and the annotated alternative 3’ splice site at the start of exon 11 (red). Bottom: Enlarged 3’ end of the locus. B Local quantification of LSVs in the representative neurons AVL and AVM for the novel alternative first exon 5.5 (left), the novel alternative splice site (middle), and the annotated alternative splice site (right). Violin plots represent the posterior distribution computed by MAJIQ. C, D Validation by RT-PCR of the novel exons in cDNA from whole animals (N2 wild type). Alternative first exon 5.5 (C), top: Location of primers relative to the exon. Left, the primer pair only amplifies the annotated transcript including exons 1-5; right, the primer pair only amplifies the novel transcript using exon 5.5 as first exon. Arrowheads: the expected product sizes (307 bp, 286 bp). Alternative 5’ splice site in exon 9 (D), top: Location of the primers relative to the event. Left, the primer pair only amplifies the novel transcript with shorter exon 9; right, the primer pair only amplifies the annotated transcript with longer exon 9. Arrowheads: expected product sizes (242, 248 bp). Three biological replicates along with a control without reverse transcriptase. Molecular weights reported in base-pairs. Source data are provided as a Source Data file. E Protein structure of SAX-3, indicating the impact of the novel splice variants. Top: Overall structure of annotated isoform SAX-3A. Bottom: Overall structure of the novel short isoform starting at novel exon 5.5 (SAX-3S). The position of the novel splice site in exon 9 is indicated. F Structure of the orthologs of sax-3 in C. briggsae (Cbr-sax-3) and C. brenneri (Cbn-sax-3). Top: Gene models. Bottom: RNA-Seq data. Red arrowheads denote the orthologs of exon 5.5. G Protein sequence alignment of the short isoforms of SAX-3 orthologs from C. elegans, C. briggsae, and C. brenneri.

Fig. 5. Differential AS between neuron types.

A Euler plot of the number of genes presenting differential alternative splicing (DAS) in this analysis, and from a literature review (Supp. Data 3). Source data are provided as a Source Data file. B Proportion of genes DAS within families of neuronally significant genes. The number of genes in the family is indicated inside each bar. Hypergeometric test: * Significantly enriched, # significantly depleted. One-sided hypergeometric test and Benjamini-Hochberg correction for multiple comparisons. Exact p-values reported in Source Data file. C Proportion of DAS genes per neuron pair, among genes co-expressed in the neurons of this pair. D Number of DAS genes plotted against number of differentially expressed genes for each neuron pair. The line corresponds to a linear fit, the shaded area to the 95% confidence interval. E Estimate of the number of genes in which DAS can be detected relative to the number of neuron types sequenced. Each subsampling was repeated 10 times; the box plots are centered on the median and between the first and third quartiles, the whiskers extend to the furthest value no further than 1.5 times the interquartile range. Source data are provided as a Source Data file.

Fig. 6. Landscape of alternative splicing event types.

A Number of events by types. Bar length indicates the number of canonical events of each type predicted from genome annotation. Dark gray: events that did not display DAS between neuron types in our dataset; dark red: events DAS in at least one pair of neurons; dark blue: events DAS between tissues, based on published tissue-specific TRAP-Seq data, but not between individual neuron types in our data; light gray: events that could not be measured in either dataset. B Distribution (density plot) of length or conservation score of particular sequences within splicing events, for DAS events (red) vs non-DAS events (gray). Only features with a statistically significant difference between DAS and non-DAS events are represented here. Vertical dashed lines represent the median of each group. Lengths in basepairs (bp). C Histogram of cassette exon lengths, separating DAS exons (red) and non-DAS exons (gray). The vertical dashed bar at 27 bp delimits microexons. D For each exon (each dot), proportion of neuron pairs where the exon is differentially AS between the two neuron types (proportion calculated among the neuron pairs in which the exon-containing gene is expressed in both neurons of the pair). p = 0.0035 (two-sided Wilcoxon test). E Potential impact of nucleotides addition on open reading frame. For alternative splice sites and cassette exons overlapping a CDS, the color indicates whether inclusion of the alternative event results in a nucleotide addition that is a multiple of 3 (lilac) or is not a multiple of 3 (green), separating events DAS between neuron types. Alt. 3’ ss, p = 0.023; alt. 5’ ss, p = 0.053; cassette exon, p = 0.41 (Chi-square test, comparison of the proportion of PFS events between DAS and non-DAS events). F Location of events within the gene body. For alternative splice sites and cassette exons, the color indicates non-PFS (lilac) and PFS events (green). The density curve indicates the distance from the event to the closest extremity (start or end) of the gene. ***: p < 0.001 (two-sided Wilcoxon test, comparison of distance from extremity between PFS and non-PFS events, test statistics reported in Source Data). Source data are provided as a Source Data file.

Fig. 7. Splicing Regulatory Network.

A Schematic representation of the approach: The splicing and expression measurements in each samples are grouped in a single data matrix, to compute a covariance matrix. The precision matrix is estimated, whose side quadrants represent the weights of a splicing regulatory network linking Splice Factors (SF) to cassette exons. B Comparison of metrics for selected parameters. The four metrics are represented across all methods. The Frobenius loss was inverted such that for all metrics, higher is better. The arrow represents the selected “best” method. C The four metrics plotted against sparsity for a range of penalties. The shaded area highlights the selected penalty of 0.25, corresponding to a sparsity of 97%. The red triangles correspond to values significantly different from the permuted data. D, E Subnetworks centered on exon 5 of C07A12.7 (D) and exon 11.5 of daf-2 (E). Blue diamond: cassette exon, red ellipse: putative splice factor displaying a non-zero weight in the network, green ellipse: splice factor identified from the literature to regulate this exon. Source data are provided as a Source Data file.