Long read sequencing reveals novel isoforms and insights into splicing regulation during cell state changes

TL;DR: In this article, a custom annotation of a well-studied human neuroblastoma cell line was used to identify isoform expression and usage across differentiation, including a novel transcript of the voltage-gated calcium channel subunit gene.

Abstract: Alternative splicing (AS) is a key mechanism underlying cellular differentiation and a driver of complexity in mammalian neuronal tissues. However, understanding of which isoforms are differentially used or expressed and how this affects cellular differentiation remains unclear. Long read sequencing allows full-length transcript recovery and quantification, enabling transcript-level analysis of AS processes and how these change with cell state. Here, we utilise Oxford Nanopore Technologies sequencing to produce a custom annotation of a well-studied human neuroblastoma cell line and to characterise isoform expression and usage across differentiation. We identify many previously unannotated features, including a novel transcript of the voltage-gated calcium channel subunit gene, CACNA2D2. We show differential expression and usage of transcripts during differentiation, and identify a putative molecular regulator underlying this state change. Our work highlights the potential of long read sequencing to uncover previously unknown transcript diversity and mechanisms influencing alternative splicing.

Summary

INTRODUCTION

• The complex suite of processes that occur during transcription gives rise to a staggering diversity of protein structures, molecular interactions and cell fates.
• The authors identified novel transcriptomic features and performed differential expression and usage analyses to identify transcripts that show variation during differentiation, as well as identifying a novel putative molecular regulator underlying this state change.

RESULTS AND DISCUSSION

• ONT reads accurately detect differential isoform expression Using the Oxford Nanopore GridION platform, the authors generated on average 10,691,538 QC-passed reads per sample (± 1,751,518.6 SD).
• It is therefore important to assess the performance of long read vs short read sequencing in both transcript quantification and its application to differential expression studies (Sessegolo et al. 2019).
• This suggests that RBM5 may play a role in splicing regulation during differentiation of SH-SY5Y cells.

Sampling and Sequencing

• Cell culture and neuronal differentiation A total of 10 technical replicates of human neuroblastoma SH-SY5Y cells were cultured in neurobasal media (Gibco 21103-049) supplemented with B-27 Plus .
• Retinoic acid was added to five replicates to a final concentration of 10mM, to induce cell differentiation to a neuronlike state; whilst five replicates were cultured to confluence in standard media.
• Cells were washed with phosphate buffered saline and harvested in QIAzol to preserve RNA, before being stored at - 80°C until RNA extraction.

RNA extraction and spike-in control

• Total RNA was purified from the 10 replicate cell cultures using a Direct-zol RNA Miniprep Plus kit (Zymo Research), according to the manufacturer’s instructions.
• The whole second-strand reaction was then mixed and incubated at 42°C for 90 minutes.
• The cDNA was quantified using High Sensitivity Qubit assays (ThermoFisher, Q32854) and sized using the 2100 Bioanalyzer instrument (Agilent Technologies, cat. no. G2939BA) High Sensitivity DNA assay (Agilent, 5067-4626).
• The TALON custom gtf contains only features detected with reads present in the dataset, so a complete custom transcriptome annotation was compiled by merging the reference and TALON gtfs.

Differential expression analyses

• Sequin spike-in detection & ONT DE sensitivity Sensitivity in detecting isoform DE using ONT was assessed by a) finding the threshold of detection for each Sequin mix, b) comparing observed vs expected logFC and c) comparing with short read data.
• This was also performed for both the full short read data and a version downsampled to equivalent ONT average nucleotide coverage using bedtools.
• The authors then utilised a standard differential expression pipeline (detailed below).
• The differential expression regression model was specified by splitting the data into Sequin MixA and MixB accordingly.
• Transcript-level counts were then obtained by importing Salmon results with the EdgeR function catchSalmon, using the bootstrap replicates to calculate and apply an overdispersion correction for each count.

Differential usage analyses

• Differential transcript usage (DTU) was assessed using the R package IsoformSwitchAnalyzeR v.1.11.3 (Vitting-Seerup & Sandelin 2019) on the same transcript quantification input used for DTE and DGE.
• TPM abundances were imported using the scaledTPM function in tximport and imported into IsoformSwitchAnalyzeR.
• The DTU analysis was run in two parts; first non-expressed isoforms were removed, and switches calculated for each gene using DEXseq (Anders et al. 2012) and nucleotide and peptide outputs for each gene were created for protein assessment.
• Transcripts were assessed for coding potential with CPAT (Wang et al. 2013), protein domain assignment with PFam (Punta et al. 2012), signal peptide prediction with SingalP v.5.0 (Armenteros et al. 2019) and intrinsically disordered regions and binding regions with IUPred2A (Mészáros et al. 2018), using default parameters according to the IsoformSwitchAnalyzeR workflow.
• The second part of the IsoformSwitchAnalyzeR DTU analysis then leveraged these data to identify isoforms switches with potential functional consequences and provide visualisation using default functions.

Hypergeometric enrichment tests

• The set of putative functionally consequential DTUs was checked for RBPs by intersection with a set of known RBPs assayed as part of the ENCODE project (Van Nostrand et al. 2020), revealing the presence of RBM5.
• Corresponding narrow-peak eCLIP bed data for RBM5 were accessed using the ENCODE portal (Davis et al. 2018) for both HepG2 isogenic replicates in the ENCODE repository (RBM5 accessions ENCFF176RGG and ENCFF998ACW downloaded 29/09/2020) and intersected using bedtools to find the most supported subset of binding targets.
• The intersection for both significant DTUs (N=104) and for total genes assessed (N=32325) with the eCLIP binding targets were then obtained by intersection with this subsetted list of ENCODE targets.
• A hypergeometric test for enrichment was performed using the phyper functionality in the R core package ‘stats’ v.4.0.2 (R Core Team 2016).

Ontology and functional association

• To interpret the differentially expressed or used gene sets, the authors assessed gene ontology and known associations with neurologically relevant biology.
• The authors used the GENE2FUNC function in FUMA (Watanabe et al. 2017) to annotate the gene sets within a biological context.
• For transcripts, the corresponding Ensembl gene ID was used.
• In each case, the default thresholds of significance and ontology enrichment were applied.
• Analyses focused on tissue specificity analyses in GTEx v.8 30 tissue types and Gene Ontogeny (GO) Biological Processes.

Figures

Fig 1 Comparison of Sequin spike-in sensitivity of detection between ONT (long read) and short read sequencing (Illumina paired-end, SRS) sequencing for Sequin synthetic spike-ins MixA (plot A) and MixB (plot B). Labelled sequin (red point) in each plot is at ONT limit of quantification (LOQ) in each mix; mix A 0.059 attomol/μl and mix B 0.27 attomol/μl. Plot C shows correlation of ONT differential isoform observed vs observed expression (log2 fold change +1) of synthetic Sequin spike-ins of known concentration. Pearson correlation coefficient is displayed along with a linear regression trend line with standard error in pale grey.

Fig 2 Breakdown of novel transcripts identified using ONT long reads and TALON custom transcriptome annotation. Cassette exons are previously unannotated positions. *CPAT assessment of CDS coding probability (CP ≥ 0.364). **Novel junctions are previously unannotated junctions identified between existing exonic parts. All novel assessments are relative to Gencode v29 human transcriptome annotation.

Fig 4 Panel of A) known RBM5 isoforms with relative size and exon composition and B) differential transcript usage (DTU) of known RBM5 isoforms across differentiation displaying a switch from a truncated, non-coding RBM5 isoform ENST00000474470 to full-length coding isoform ENST00000347869. Figures modified from IsoformSwitchAnalyzeR (Vitting-Seerup, 2017) output.

Table 1 Showing number of genes and transcripts differentially expressed between undifferentiated and differentiated cells, after multiple testing correction using ONT long read counts (FDR < 0.05). Bracketed numbers refer to the portion of total that are TALON-identified novel transcripts. U = undifferentiated and D = differentiated cells, with arrows displaying expression directionality.

Fig 3 Panel of A) gene-level differential expression (DE) smear plot, solid red lines highlighting ± 1.5 logFC threshold, B) gene-level DE of CACNA2D2 during differentiation, C) isoform-level DE smear plot with ± 1.5 logFC threshold, D) CACNA2D2 isoform expression, showing novel TALON isoform with highest read count. Red points on smear plots indicate significant differential expression (FDR < 0.05). Boxplots display median and IQR.

Full text

Long read sequencing reveals novel isoforms and insights into splicing regulation during cell state
changes
David J Wright
1
, Nicola Hall
2,3
, Naomi Irish
1
, Angela L Man
1
, Will Glynn
1
, Arne Mould
2,3
, Alejandro
De Los Angeles
2,3
, Emily Angiolini
1
, David Swarbreck
1
, Karim Gharbi
1
, Elizabeth M Tunbridge
2,3
,
Wilfried Haerty
1
*
1
Earlham Institute, Norwich Research Park, Norfolk NR4 7UZ, UK
2
Department of Psychiatry, Medical Sciences Division, University of Oxford, Oxfordshire OX3 3JX,
UK
3
Oxford Health NHS Foundation Trust, Oxford, Oxfordshire OX3 7JX, UK
* corresponding author: Wilfried.Haerty@earlham.ac.uk
ABSTRACT
Alternative splicing (AS) is a key mechanism underlying cellular differentiation and a driver of
complexity in mammalian neuronal tissues. However, understanding of which isoforms are
differentially used or expressed and how this affects cellular differentiation remains unclear. Long read
sequencing allows full-length transcript recovery and quantification, enabling transcript-level analysis
of AS processes and how these change with cell state. Here, we utilise Oxford Nanopore Technologies
sequencing to produce a custom annotation of a well-studied human neuroblastoma cell line and to
characterise isoform expression and usage across differentiation. We identify many previously
unannotated features, including a novel transcript of the voltage-gated calcium channel subunit gene,
CACNA2D2. We show differential expression and usage of transcripts during differentiation, and
identify a putative molecular regulator underlying this state change. Our work highlights the potential
of long read sequencing to uncover previously unknown transcript diversity and mechanisms
influencing alternative splicing.
INTRODUCTION
The complex suite of processes that occur during transcription gives rise to a staggering diversity of
protein structures, molecular interactions and cell fates. Alternative splicing (AS) allows different
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted April 28, 2021. ; https://doi.org/10.1101/2021.04.27.441628doi: bioRxiv preprint

transcripts to be generated from a single gene. Differential transcript expression (the overall abundance
of a given transcript) or transcript usage (the abundance of a given transcript relative to that of others
produced from the same gene) are key mechanisms for regulating cell lineage commitment and function
(Breschi et al. 2020; Chepelev & Chen 2013). In vertebrates, AS is particularly prominent in the brain,
and regulates multiple aspects of neurodevelopment including neurogenesis, synaptogenesis, cellular
migration and axon guidance (Grabowski 2011; Ule et al. 2005; Raj & Blencowe 2015) in a temporally
precise manner (Weyn-Vanhentenryck et al. 2018; Liu et al. 2018; Burke et al. 2020). These
neurodevelopmental processes are defined by ordered switches in exon usage and expression across a
spectrum of genes, controlled by a suite of highly specific RNA-binding proteins (RBPs) such as
NOVA2 (Saito et al. 2019), PTBP1 and PTBP2 (Boutz et al. 2007; Linares et al. 2015; Keppetipola et
al. 2012). A number of more ubiquitous RBPs may also help regulate these neuronal AS events, though
which ones and what specific roles they play remain poorly understood (Jackson et al. 2020; Gallego-
Paez et al. 2017).
Since AS can give rise to mRNAs that encode protein isoforms that exhibit distinct, or even opposing
effects, it is essential to understand an individual gene’s products at transcript-level resolution (Clark et
al. 2007; Yi et al. 2018; Liu et al. 2018; Yuste et al. 2020). However, the diversity of full-length
transcripts remains poorly understood, as exemplified by the recent study of the L-type voltage gated
calcium channel (VGCC) gene, CACNA1C (Clark et al. 2020). Furthermore, many unknowns remain
as to the nature and regulation of changes in transcript expression during differentiation and
development. For example, are there pronounced switches in primary transcript expression in a few key
genes, or more nuanced expression differences across the transcriptome? Furthermore, although some
of the molecular mechanisms that drive the observed ‘switches’ in transcriptional profiles occurring
during lineage commitment have been identified, many of these processes remain to be determined.
As well as being of importance for understanding normal developmental processes, AS is also of clinical
relevance, since aberrant transcriptional processes are implicated in many diseases (Scotti & Swanson
2016). Disease-associated mutations can directly affect AS by disrupting existing splice sites and/or
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted April 28, 2021. ; https://doi.org/10.1101/2021.04.27.441628doi: bioRxiv preprint

forming novel or cryptic sites, as observed in the VGCC CACNA1A gene in Episodic Ataxia Type 2
(Jaudon et al. 2020). Alternatively, AS can alter disease presentation, as is seen in the case of Timothy
Syndrome where the localisation of the disease-causing mutation in one of two mutually exclusive
exons of CACNA1C determines syndrome severity (Splawski et al. 2004). Global changes in differential
isoform expression are also associated with psychiatric conditions (Gandal et al. 2018).
Transcriptome profiling and annotation are essential first steps in investigating gene, isoform and exon
expression or usage differences during cell differentiation. Until recently profiling was hampered by
technological constraints, relying on short read sequencing technology (Stark et al. 2019; Wang et al.
2009). Whilst short read technologies provide cheap, accurate and high-coverage reads, with good
differential expression analysis power (Wang et al. 2009), their ability to resolve and quantify full-
length transcripts is inherently limited (Byrne et al. 2017). In this context, the advent of long read
technologies has rapidly improved our ability to characterise the transcriptome (Byrne et al. 2017, 2019)
revealing, for example, the complexity of the transcriptional landscape of the mammalian brain (Wang
et al. 2019; Sessegolo et al. 2019).
Here, we use long read sequencing to identify and quantify isoforms during a cellular state change;
specifically, during the differentiation of the well-validated SH-SY5Y neuroblastoma line into neuron-
like cells. SH-SY5Y cells exhibit a stable genomic structure and have been widely used to investigate
AS mechanisms and cellular differentiation from a neuroblast-like state (Kovalevich & Langford 2013;
Shipley et al. 2016; Agholme et al. 2010; Truckenmiller et al. 2001), into a neuronal-like state (Forster
et al. 2016; Mendsaikhan et al. 2018). We generated a custom high-coverage long read transcriptome
annotation (using Oxford Nanopore Technology [ONT] cDNA sequencing), validated with orthogonal
short read sequencing (Illumina paired-end short read) data. We identified novel transcriptomic features
and performed differential expression and usage analyses to identify transcripts that show variation
during differentiation, as well as identifying a novel putative molecular regulator underlying this state
change.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted April 28, 2021. ; https://doi.org/10.1101/2021.04.27.441628doi: bioRxiv preprint

RESULTS AND DISCUSSION
ONT reads accurately detect differential isoform expression
Using the Oxford Nanopore GridION platform, we generated on average 10,691,538 QC-passed reads
per sample (± 1,751,518.6 SD). We also generated an average of 105,349,119 lllumina read pairs per
sample (± 17,312,599.92 SD, Table S1). Whilst the utility of long read sequencing for recovering full
length transcripts is widely accepted, there remains uncertainty as to the sensitivity of this technology
for differential expression analysis. It is therefore important to assess the performance of long read vs
short read sequencing in both transcript quantification and its application to differential expression
studies (Sessegolo et al. 2019). We investigated the ability of the ONT data to detect Sequin spike-ins
(Hardwick et al. 2016) of known concentration in a set of two different concentration mixes. We found
that ONT limit of quantification (minimum transcript concentration) was 0.059 attomol/μl for mixA
and 0.27 attomol/μl for mixB (Figs. 1A, 1B). By downsampling the short read data to the ONT average
nucleotide coverage, we show there is similar power to detect transcripts by ONT and downsampled
short read, although ONT scores lower than the complete short read data (Table S2). Next, we directly
assessed the ability of the ONT reads to detect differential isoform expression. We calculated expected
log
2
fold-change (logFC) from the differences in concentrations between Sequins in mixA and mixB
and compared this with observed logFC from our differential expression pipeline (see Methods). There
was a strong correlation between expected and observed logFC of R
2
= 0.973, p-value = 2.2e
-16
(Fig.
1C), demonstrating that the ONT data can be used to detect differential expression over a broad range
of logFC values. Collectively, our results indicate that our ONT data is sufficiently sensitive and
powered to detect all but the lowest concentration transcripts and, further, that the ONT reads are
suitable for differential isoform expression analysis. These results add to the growing literature
highlighting the importance of assessing suitability of long read sequencing for differential expression.
TALON custom long read annotation reveals novel features of the human transcriptome
The TALON custom annotation provided a total of 3,274 novel transcripts prior to validation using
short read sequencing. We found short read support for 2,567 of the 3,274 (78.41%) novel transcripts
recovered from the ONT read data (Fig. 2) by stringent removal of transcripts that contained a novel
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted April 28, 2021. ; https://doi.org/10.1101/2021.04.27.441628doi: bioRxiv preprint

exon lacking at least 15 reads depth across 75% of its length (see methods). The supported novel
transcripts collectively include a total of 49 novel cassette exons (18 frame-conserving) along with 928
and 1046 novel 5’ and 3’ splice sites respectively, with 464 instances of exons exhibiting both a novel
5’ and 3’ splice site. Additionally, we identified 92 novel junctions between previously annotated splice
sites. In total 929 (36.19%) of the validated novel transcripts were putatively coding; either frame-
conserving, or assumed to be coding via CPAT (Wang et al. 2013) assessment, whilst 1638 (63.81%)
were assumed noncoding due to either induction of a frameshift, a noncoding parent gene or via
noncoding classification from CPAT (fig. 2 for full breakdown).
Our stringent filtering criteria and validation likely results in an underestimate of the true quantity of
novel features present. Despite this, our long read sequencing approach still identifies >3000 novel
transcripts, nearly a thousand of which are putatively coding. Collectively, our data highlight the extent
of previously undescribed transcriptome diversity, even within a highly specialised (and well-studied)
cell model. Our work concurs with the growing body of other studies using long reads for transcriptome
assessment (Gleeson et al. 2020; Sessegolo et al. 2019, Soneson et al. 2019); relying on short reads
substantially underestimates transcriptome diversity.
ONT differential gene expression supports neuron-like characteristics of differentiated SH-SY5Y
cells
Differential gene expression analysis revealed 4,239 genes differentially expressed (FDR q < 0.05)
between differentiation states, with 2,041 and 2,198 genes overexpressed in undifferentiated and
differentiated cells, respectively (Table 1, Fig 3A). We performed FUMA analyses to explore the
functional significance of these genes. The upregulated genes in differentiated cells showed greatest
overlap with those upregulated in brain compared with other tissue types (p
adj
= 3.4 x 10
-41
), whilst for
those more highly expressed in undifferentiated cells, there was overlap with those downregulated in
brain (p
adj
= 2.8 x 10
-45
, second only to pancreas: p
adj
= 1.8 x 10
-45
). Gene Ontogeny biological pathway
terms showing differential expression across differentiation included neurogenesis, neuron
development, cell differentiation and regulation of nervous system development (Table S5).
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted April 28, 2021. ; https://doi.org/10.1101/2021.04.27.441628doi: bioRxiv preprint

Frequently asked questions

Q1. What contributions have the authors mentioned in the paper "Long read sequencing reveals novel isoforms and insights into splicing regulation during cell state changes" ?

Here, the authors utilise Oxford Nanopore Technologies sequencing to produce a custom annotation of a well-studied human neuroblastoma cell line and to characterise isoform expression and usage across differentiation. The authors show differential expression and usage of transcripts during differentiation, and identify a putative molecular regulator underlying this state change. Alternative splicing ( AS ) allows different ( which was not certified by peer review ) is the author/funder. Their work highlights the potential of long read sequencing to uncover previously unknown transcript diversity and mechanisms influencing alternative splicing. 

Q2. What are the future works mentioned in the paper "Long read sequencing reveals novel isoforms and insights into splicing regulation during cell state changes" ?

Future work needs to further investigate how changes in RBM5 regulation impacts isoform expression through experimental confirmation of RBM5 binding targets in SH-SY5Y cells, and mutagenesis of RBM5 and its binding sites. Finally, their findings indicate that changes in RBM5 expression profiles may act as a molecular mechanism for the coordination of these changes, paving the way for future functional studies.