Loss of TDP-43 may cause abnormal RNA processing in ALS and FTD

This approach allowed us to subdivide our events into three main categories—ALEs, IPAs and 3'Ext (Fig. 1a)—overcoming the limitations of comparable available tools that focus on specific event categories^24–28. APA events were widespread, and we defined cryptic APA events as ones with less than 10% mean usage in controls and more than 10% usage change after TDP-43 knockdown. We identified 227 cryptic APAs to be present in at least one dataset (adjusted P < 0.05; Fig. 1b, Supplementary Fig. 1 and Supplementary Table 2). Cryptic ALEs (n = 92) included previously identified cryptic exons such as STMN2, ARHGAP32 and RSF1 (Fig. 1b and Supplementary Fig. 2). In total, 108 3'UTR cryptics were identified, of which 86 are novel 3'UTR extensions (3'Ext; for example, TLX1; Fig. 1c), and 20 were 3'UTR shortening events at loci with novel 3'Exts (3'shortening). Twenty IPA events were also detected, including CNPY3, which was identified with an independent bioinformatics approach and experimentally validated²⁹. The remaining nine events could not be uniquely assigned to ALEs or IPAs based on annotation and are defined as 'complex'. Multiple non-cryptic APA events were also detected and are reported in Supplementary Table 3.

We experimentally validated strong activation of cryptic APA and confirmed the expression of multiple predicted PASs by performing 3' rapid amplification of cDNA ends (3'RACE) in i3Neurons (Extended Data Fig. 1) and inspecting poly(A)-tail ligation-dependent, oligo-dT primer-free i3Neuron direct RNA nanopore sequencing¹³ (Supplementary Fig. 3a−c). We further evaluated global cryptic polyadenylation site (PAS) precision by pooling across TDP-43 depletion RNA-seq samples poly(A)-tail-containing reads (PATRs; Supplementary Fig. 3d), which allows independent defining of PASs³⁰. Cryptic and expression-matched annotated PASs were similarly identified, further supporting the novel cryptic APA events (Supplementary Fig. 3e). Finally, the commonly used tool DaPars2 (ref. ³¹), when provided with the predicted 3'Ext coordinates, reproduced cryptic 3'Ext activation (Supplementary Fig. 4). These findings collectively support the validity of our cryptic APA discovery pipeline.

Out of 227 cryptic APAs detected by our analysis across datasets, most (138) satisfied cryptic expression criteria (<10% mean usage in controls and >10% usage change after TDP-43 knockdown) when considering the median across datasets. Fifty-one APAs were, instead, consistently below 10% usage threshold in controls but did not sufficiently increase after TDP-43 depletion to meet the cryptic criteria definition across datasets. Twenty-eight APAs showed, instead, a significant increase upon TDP-43 loss across datasets but had more than 10% median usage in controls, therefore placing them outside the cryptic criteria but demonstrating consistent regulation by TDP-43 (Supplementary Fig.). Altogether, these data highlight a widespread presence of cryptic APA upon TDP-43 loss. 5

Next, we investigated TDP-43 binding patterns around cryptic APAs using TDP-43 individual-nucleotide resolution UV crosslinking and immunoprecipitation (iCLIP) data generated in SH-SY5Y cells¹⁰. We focused on ALEs and 3'Ext events as the low number of IPA and 3'shortening events (n = 20 in both cases) did not allow reliable binding profile inferences. TDP-43 binding was enriched around the splice acceptor of cryptic ALEs, as previously described in cryptic splice junctions, and downstream of the cryptic PAS of ALEs (Fig. 1d), supporting TDP-43 acting as a repressor of both splicing and polyadenylation. Intriguingly, TDP-43 binding was also enriched immediately downstream of the annotated proximal PAS of 3'Ext events (Fig. 1e), supporting a role for TDP-43 in enhancing poly(A) usage, consistent with previous reports of TDP-43 binding with respect to regulated PAS⁵.

iCLIP data, typically generated in control cells, are not sensitive in detecting binding to cryptic 3'Ext regions, as these events can be detected only at very low levels with physiological TDP-43 presence. We, therefore, sought to corroborate our findings by adapting PEKA³² to infer de novo hexamer enrichment relative to cryptic landmarks. Previously defined hexamers enriched around TDP-43 iCLIP binding sites⁶ (Supplementary Fig. 6a) were overrepresented among the most enriched hexamers proximal to all cryptic landmarks, with the strongest signal overall observed at both the 3' splice site (3'ss) and PAS of ALE events (Supplementary Fig. 6b). To assess the concordance with iCLIP binding profiles, we visualized the positional coverage of the hexamer group most strongly associated with TDP-43 binding⁶. For ALEs, we observed a notable peak immediately upstream of splice acceptors and a strong peak downstream of PAS (Fig. 1d), although previous reports of splice-site-dependent STMN2 cryptic ALE repression³³ suggest that the binding at PAS may have secondary effects. Enriched signal was also observed immediately downstream of the distal PAS of 3'Exts (Fig. 1e).

To experimentally validate the direct relationship between TDP-43 binding and cryptic PAS usage, we generated a reporter for the ELK1 3'Ext APA (Fig. 1f). Nanopore sequencing showed a strong upregulation of the distal cryptic PAS upon TDP-43 knockdown in neuronal cells (Fig. 1g), confirming similar behavior to endogenous ELK1. We then focused on 150 base pairs downstream of the constitutive poly(A) site, where iCLIP data show TDP-43 binding to occur, and generated a series of constructs where we removed or increased UG content to disrupt or enhance TDP-43 binding (Fig. 1f). Under normal TDP-43 levels, cryptic PAS usage was enhanced by UG depletion, whereas it was reduced by UG dinucleotide content increase (Extended Data Fig. 2a,b). Increasing levels of TDP-43 knockdown enhanced cryptic PAS usage in constructs with normal, increased or moderately disrupted UGs, whereas constructs with severe UG depletion did not respond to TDP-43 depletion, confirming a direct regulation by TDP-43 (Fig. 1h).

Overall, our data support a direct role for TDP-43 binding in both enhancing and repressing PAS usage, therefore leading to cryptic APA upon TDP-43 loss.

Next, we used the New York Genome Center (NYGC) ALS Consortium RNA-seq dataset to assess cryptic APA in a larger cohort of CNS cases with or without TDP-43 pathology (Supplementary Table 4). Cryptic 3'Exts often demonstrated low basal expression in control samples in our in vitro datasets, confounding the detection in postmortem bulk RNA-seq datasets, in which only a very small proportion of cells is expected to have TDP-43 pathology. IPA detection is further complicated by the fact that normal pre-mRNA reads also map to IPA regions, creating significant noise in bulk RNA-seq. We, therefore, focused on ALEs, where detection of the associated upstream cryptic splice junctions provide direct evidence of expression. As cryptic ALEs are expected to be dependent on nuclear TDP-43 depletion, we defined criteria based on spliced read detection to identify cryptic events with specific expression in tissues and disease subtypes where TDP-43 pathology is present. Of 118 cryptic ALE junctions, 7 fulfilled specificity criteria (Supplementary Table 5), in contrast to 56 out of 313 cryptic splicing events collated from i3Neurons with TDP-43 knockdown¹³ (Fig. 2c and Extended Data Fig. 3). STMN2 was most frequently detected in tissues with expected TDP-43 proteinopathy, and several other ALEs were among the most frequently detected specific cryptic events, including SYNJ2 (third; Fig. 2d) and PHF2 (eighth; Fig. 2e).

Altogether, this suggests that cryptic APAs are detectable in postmortem tissue affected by TDP-43 pathology, highlighting their potential relevance in loss-of-function disease mechanisms and their promising utility as biomarkers.

Cryptic splicing events impact expression, often leading to a reduction in transcript levels^9–11. We, therefore, assessed the effect of cryptic APAs on their own transcripts in i3Neurons¹³ (Supplementary Fig. 8a) and found that the majority of events (86 out of 126) coincide with a significant change in expression, equally split between significant upregulation and downregulation. When subdivided further into cryptic APA categories, no category showed a clear bias for upregulation or downregulation (19 out of 34 3'Ext, 17 out of 37 ALE and 6 out of 10 IPA genes are downregulated). This suggests that cryptic APAs are associated with differential expression but have variable effects on transcript levels.

Regulation of both ALE and 3'Ext usage has been demonstrated to impact protein abundance through distinct mechanisms^37,38, but differential RNA abundance does not necessarily imply a coordinated change in protein levels. To assess whether changes in gene expression were also reflected in translation levels, we performed differential translation analysis of Ribo-seq data generated from i3Neurons with TDP-43 depletion¹³.

Interestingly, the three 3'Ext-containing genes that were most upregulated at both RNA and translation levels (Fig. 3a,b) encode for three transcription factors: ELK1, SIX3 and TLX1. The regulation of these 3'Ext events is reproducible across in vitro datasets (Supplementary Fig. 1). As ELK1 increase was previously associated with neuronal toxicity^43–45 and its levels are consistently higher in mature neurons, compared to SIX3 and TLX1, which are associated with neuronal development^46,47, we decided to focus our investigations on ELK1. We tested whether the increase in Ribo-seq also corresponded to an upregulation of steady-state protein, and western blots confirmed a significant increase in ELK1 protein expression upon TDP-43 knockdown in i3Neurons (Fig. 3c). We next asked whether the activity of ELK1, which functions as a transcription factor in the ternary complex factor (TCF) family⁴⁸, could be altered in the context of TDP-43 loss. We assessed whether ELK1 target genes defined by chromatin immunoprecipitation followed by high-throughput sequencing (ChIP–seq) in HeLa cells were also affected in TDP-43 knockout HeLa cells⁴⁹, in which the cryptic 3'Ext is robustly upregulated (Fig. 3d). Using GSEA, we observed a significant change in ELK1 target gene expression upon TDP-43 knockout (Fig. 3d). This suggests that cryptic 3'Exts can lead to change in function in the context of TDP-43 loss.

We investigated the mechanisms by which cryptic 3'UTRs could mediate increased translation levels of ELK1, SIX3 and TLX1. We revisited differential splicing analysis of i3Neuron RNA-seq datasets^10,13 and confirmed that cryptic 3'Exts are the only differential RNA processing events occurring in these three transcription factor RNAs upon TDP-43 depletion. As alternative 3'UTRs have been linked to differences in RNA stability⁵⁰, we reasoned that increased RNA stability could account for changes in overall RNA abundance and translation levels. To investigate changes in RNA stability in i3Neurons with TDP-43 depletion, we performed SLAM-seq⁵¹, which allows the detection of newly synthesized RNAs through incorporation of a uridine analogue (4SU). Different lengths of 4SU treatment allow the estimation of gene-level RNA half-lives. We observed increased half-lives in cryptic 3'Ext-containing genes ELK1, TLX1 and SIX3 (Fig. 3e). To confirm that the 3'Ext half-life change was due to the cryptic APA event, we performed an isoform-specific analysis for ELK1 in control i3Neurons, where the distal (cryptic) 3'Ext is sufficiently expressed to be analyzed and not prevalent enough to confound the evaluation of the proximal (shared) isoform. We observed elevated ELK1 3'Ext half-life relative to the proximal PAS (Supplementary Fig. 9a). Altogether, this suggests that increased RNA abundance and translation of cryptic 3'Ext genes are mediated by increased RNA stability.

Given that translation depends on extranuclear localization of mRNAs, we tested whether cryptic 3'Ext transcripts localize to the cytoplasm and contribute to the increased translation levels^52–55. Focusing on the ELK1 cryptic 3'Ext, we designed probes to recognize the common proximal sequence and the distal sequence specific to the 3'Ext and performed fluorescence in situ hybridization (FISH) in i3Neurons where we could detect both probes in the nuclei, cytoplasm and neurites (Fig. 3f). Consistent with RNA-seq, we observed a significant increase in total foci for both the total and cryptic-specific probes upon TDP-43 knockdown (Fig. 3g and Supplementary Fig. 9b,c). To specifically discriminate and quantify proximal and distal ELK1 APA subcellular localization, we performed 3'RACE on SH-SY5Y cells after nuclear cytoplasmic fractionation. We found that both isoforms are predominantly localized to the cytoplasm and that, upon TDP-43 knockdown, the proximal canonical APA is reduced, whereas the cryptic 3'Ext is increased (Fig. 3h and Extended Data Fig. 4). Finally, we evaluated ELK1 isoform-specific ribosome recruitment using fractionation and sequencing (Frac-seq) data from neural progenitor cells⁵⁶. We found ELK1 cryptic 3'Ext to be relatively enriched in ribosome-associated fractions, supporting a preferential engagement of the cryptic 3'Ext with the translation machinery (Supplementary Fig. 9d). Overall, these findings show that ELK1 cryptic 3'Ext has increased RNA stability, localizes to the cytoplasm and neurites and is translated, driving the increase in ELK1 protein.

CRISPR interference (CRISPRi) knockdown experiments were performed in the WTC11 iPSC line harboring stable TO-NGN2 and in dCas9−BFP−KRAB cassettes at safe harbor loci⁶². CRISPRi knockdown of TDP-43 in iPSCs was achieved using single guide RNA (sgRNA) targeting the transcription start site of TARDBP (or non-targeting control sgRNA)¹⁰, delivered by lentiviral transduction. sgRNA sequences were as follows: non-targeting control GTCCACCCTTATCTAGGCTA and TARDBP GGGAAGTCAGCCGTGAGACC. iPSCs were differentiated into cortical-like i3Neurons as described previously^10,63 and fixed 9 days after re-plating for RNA-FISH.

For RNA-seq experiments ('Humphrey i3 cortical'), i3Neurons were induced as previously described⁶³ with the addition of SMAD and WNT inhibitors⁶⁴ (SB431542 10 µM; LDN-193189 100 nM; XAV939 2 µM, all from Cambridge Bioscience). After induction, cells were cultured in BrainPhys Media (STEMCELL Technologies) with 20 ng ml⁻¹ BDNF (PeproTech), 20 ng ml⁻¹ GDNF (PeproTech), 1× N2 supplement (Thermo Fisher Scientific), 1× B27 supplement (Thermo Fisher Scientific), 200 nM ascorbic acid (Sigma-Aldrich), 1 mM dibutyryl cyclic-AMP (Sigma-Aldrich) and 1 µg ml⁻¹ laminin (Thermo Fisher Scientific), as previously described⁶⁵, and harvested 30 days after differentiation. The 'Zanovello i3 Cortical' samples were generated as previously described for the dual TDP-43/UPF1 knockdown experiments¹⁰. Only the TDP-43/Control and Control/Control transfection conditions were used for RNA-seq. See the 'RNA-seq' section for library preparation details.

An iPSC line with an N-terminal HaloTag on both endogenous copies of TDP-43 (Halo-TDP-43 i3Neurons) was generated by CRISPR−Cas12 gene editing⁶¹. The parental cell line used was the WTC11 cell line with integrated dCas9−Krab and NGN2 cassettes as mentioned previously⁶². The homology-directed repair (HDR) template used was Addgene plasmid 178131. Editing was done with Cas12 CRISPR RNA (crRNA) (Integrated DNA Technologies) with GGAAAAGTAAAAGATGTCTGAAT as the targeting sequence. Recombinant Cas12 (Cpf1 ultra; Integrated DNA Technologies) was electroporated with HDR template and Cas12 crRNA using the P3 Primary Cell 4-D Nucleofector Kit (Amaxa, V4XP-3024). iPSCs were then single-cell plated, and positive colonies were selected with HaloTag TMR dye (Promega) and verified by polymerase chain reaction (PCR) of genomic DNA.

For proteolysis-targeting chimera (PROTAC)-mediated knockdown of Halo-TDP-43, i3Neurons were treated with HaloPROTAC-E⁶⁶ (30 nM) on days in vitro 14 (DIV14) and harvested on DIV28. This protocol allows to avoid incurring in maturation alterations caused by loss of TDP-43, as this occurs at a later step; we, therefore, used this approach to validate ELK1 protein increase as transcription factor levels can be sensitive to maturation stages.

Cortical-like i3Neurons were cultured on 13-mm glass coverslips and fixed in 4% paraformaldehyde (PFA)/sucrose on day 9. RNA-FISH was performed using the QuantiGene ViewRNA ISH Cell Assay Kit (Invitrogen, QVC0001), according to the manufacturer's instructions. Protease was used at 1:1,000 dilution. Two probe sets were used to detect the canonical ELK1 transcript (TYPE 4 probe, 488-nm) or specifically the distal 3'UTR cryptic extension (TYPE 1 probe, 550-nm). Confocal images were acquired with an LSM 980 laser scanning confocal microscope with Airyscan 2 (Zeiss), using a ×40 oil immersion objective.

For each biological replicate, 10 images were acquired for the control and TDP-43 knockdown conditions. For each image, foci for both probes were counted within the 106.07-µm × 106.07-µm field of view on FIJI/ImageJ using the maximum intensity z-projection function to flatten the 2-µm-thick z stack. The 'Find Maxima' function using the same prominence setting between conditions was performed to quantify total numbers of RNA foci. To separately count nuclear and cytoplasmic foci, the Cell Counter plugin was used. For each probe and field of view, the total number of foci was divided by the number of DAPI-stained nuclei to give the average number of foci per cell. To calculate the nuclear:extranuclear ratio for the 'Total ELK1' probe, the number of nuclear foci was divided by the number of extranuclear foci in each field of view. For each probe and condition, the mean number of foci per cell and the nuclear:extranuclear ratio were calculated from the 10 images and normalized, for each biological replicate, to the respective control condition. Statistical significance was evaluated using a one sample t-test with a log transformation and the Benjamini− Hochberg false discovery rate procedure, testing the null hypothesis that mean = log(1).

Halo-TDP-43 i3Neurons were homogenized in lysis buffer (25 mM Tris-HCl, 150 mM NaCl, 1% NP-40, 1% glycerol, 2 mM EDTA, 0.1% SDS, protease inhibitor (cOmplete EDTA-free protease inhibitor cocktail; Roche) and phosphatase inhibitor (PhoSTOP; Roche)). Samples were loaded on a NuPAGE 4−12% Bis-Tris protein gel (Invitrogen), which was run in NuPAGE MOPS buffer. Proteins were transferred onto PVDF blotting membrane (Amersham) through wet transfer for 1 hour and 30 minutes at 200 mA in transfer buffer (25 mM Tris, 192 mM glycine and 20% methanol). The membrane was blocked in 5% milk in TBST (20 mM Tris, 150 mM NaCl and 0.1% Tween 20) and incubated overnight with primary antibodies diluted in 5% milk in TBST (anti-ELK1 (Abcam, ab32106) 1:500, anti-TDP-43 (Abcam, ab104223) 1:2,000 and anti-tubulin (Sigma-Aldrich, MAB1637) 1:5,000). After 1-hour incubation with horseradish peroxidase (HRP)-conjugated secondary antibodies diluted in 5% milk in TBST (anti-mouse HRP (Bio-Rad, 1706516) 1:10,000 and anti-rabbit HRP (Bio-Rad, 1706515) 1:10,000), the membrane was developed using Immobilon Classico HRP substrate (Sigma-Aldrich) and the Bio-Rad ChemiDoc system.

For the fractionation experiments, SH-SY5Y cells were treated for 10 days with 25 ng ml⁻¹ doxycycline hyclate (Sigma-Aldrich) to induce the short hairpin RNA (shRNA) against TDP-43. After 10 days, cells were trypsinized, pelleted and resuspended in 1× PBS. Before re-pelletting them, a fraction for each sample was saved for protein analysis to assess TDP-43 depletion. The other fraction was used for the subcellular fractionation with the Ambion PARIS Kit (Life Technologies), according to the manufacturer's instructions. RNA from the nuclear and cytosolic fractions was extracted with the Direct-zol kit (Zymo Research) with on-column DNase I treatment. For each experimental condition, 2 μg of cytoplasmic RNA and an equal volume of nuclear RNA fraction were reverse transcribed with the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific) according to the manufacturer's instructions and analyzed by RT−qPCR with PowerUp SYBR Green Master Mix (Thermo Fisher Scientific). DNA amplification was monitored on a QuantStudio 5 Real-Time PCR system (Applied Biosystems). GAPDH and pre-GAPDH transcripts were used as cytosolic and nuclear controls, respectively. The oligonucleotides used for the analyses are reported in Supplementary Table 9.

For each condition, equal amounts of total RNA were reverse transcribed in a 20-µl reaction with the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific), according to the manufacturer's instructions, using 1 µl of 50 µM oligo dT-anchor RT primer. cDNAs were diluted to 1 ng µl⁻¹, and the expression of each target was evaluated through RT−qPCR with PowerUp SYBR Green Master Mix (Thermo Fisher Scientific) using a gene-specific forward primer and the PCR universal reverse primer. DNA amplification was monitored on the QuantStudio 5 Real-Time PCR system (Applied Biosystems). Unless otherwise specified in the figure legend, relative RNA quantity was calculated as the fold change (2^−ΔΔCt) with respect to the experimental control sample set as 1 and normalized over GAPDH, used as an endogenous control. The oligonucleotides used for the analyses are reported in Supplementary Table 8.

For the initial test, we cloned the ELK1 proximal 3'UTR and the first 800 bp of ELK1 cryptic 3'Ext into the region downstream of the mGreenLantern coding sequence in a dual-fluorescent (mScarlet and mGreenLantern), dual-promoter reporter plasmid. We then transfected two groups of SK-N-BE(2) cells with our construct: one treated with 1,000 ng ml⁻¹ doxycycline and one untreated, and each group had triplicates. This cell line contains the SMARTvector, which enables TDP-43 knockdown upon doxycycline treatment. One day after transfection, we combined triplicates together for RNA extraction and performed 3'RACE to generate DNA samples. Subsequently, we submitted these samples for nanopore sequencing and analyzed the sequencing data to assess APA site usage. First, we used minimap2 version 2.28 (ref. ⁶⁷) to perform alignment. Subsequently, we determined the polyadenylation site for each read by locating the sequence of 10 consecutive adenosines and their corresponding position in the alignment reference.

For the subsequent UG replacement experiment, we constructed a plasmid library. This cloning included three steps. (1) We inserted a restriction site between the mGreenLatern coding sequence and ELK1 3'UTR within the construct. (2) We digested the construct with AfeI (New England Biolabs, R0652) and AccI (New England Biolabs, R0161), whose cutting sites are located at proximal 3'UTR and cryptic 3'Ext, respectively. Next, using Gibson assembly⁶⁸, we assembled the digested plasmid backbone with the inserts (described below) to produce the library. Plasmids with different inserts were referred to as variants. (3) We used the restriction site inserted in the first step to incorporate a 15-mer random barcode into each variant. After this, each variant in the library corresponded to one or more unique barcodes, which could be used to identify inserts during sequencing data analysis.

Each insert consisted of three distinct fragments: the first fragment comprised the last 192 bp of ELK1 proximal 3'UTR, whereas the second and third fragments comprised the first 350 bp of ELK1 cryptic 3'Ext. Moreover, the first and last 28 bp of each fragment were conserved to enable Gibson assembly with adjacent fragments and the plasmid backbone. To emphasise the importance of the first 150 bp of cryptic 3'Ext within the second fragment, we focused on it in our results.

Before transfection, we conducted nanopore sequencing to identify each variant's corresponding unique barcodes. We followed the protocol described above to transfect the plasmid library into SK-N-BE(2) cells with SMARTvector in four treatment groups (0, 30 ng ml⁻¹, 60 ng ml⁻¹ and 1,000 ng ml⁻¹ doxycycline). The protocol was performed in triplicate for each variant, and replicates were not combined before RNA extraction. After obtaining the nanopore sequencing results, we used a custom script to extract the barcode sequence from each read to identify which insert the read should be aligned to. Reads were aligned, and the APA site usage was determined by using the method described above.

Variant design and analysis code are available at https://github.com/MaxChien1996/replace_UG_in_first_800_bp_of_ELK1_extended_3_prime_UTR↗.

SH-SY5Y and SK-N-BE(2) cells were transduced with a SMARTvector lentivirus (V3IHSHEG_6494503) containing a doxycycline-inducible shRNA cassette for TDP-43. Transduced cells were selected with puromycin (1 μg ml⁻¹) for 1 week, before being plated as single cells and expanded to obtain a clonal population. Cells were grown in DMEM/F12 + GlutaMAX (Thermo Fisher Scientific) supplemented with 10% FBS (Thermo Fisher Scientific) and 1% penicillin−streptomycin (Thermo Fisher Scientific). For induction of shRNA against TDP-43, cells were treated with the following amounts of doxycycline hyclate (Sigma-Aldrich) and collected after 10 days:For experiments in SH-SY5Y cells (curves), 75 ng ml⁻¹For experiments in SH-SY5Y cells (cycloheximide), 25 ng ml⁻¹For experiments in SK-N-BE(2) cells, 1,000 ng ml⁻¹

Strand-specific, poly(A)-enriched sequencing libraries for the 'Humphrey i3 cortical' dataset were prepared using the KAPA mRNA Hyper Prep Kit. One hundred total nanograms of RNA was used as input material for poly(A)-positive mRNA capture. Fragmentation was performed for 6 minutes at 85 °C to obtain a target fragment size of 300−400 bp, and 13 cycles of PCR amplification were performed. The resulting libraries were sequenced 2 × 150 bp on an Illumina NextSeq 2000 machine.

RNA was extracted from i3Neurons ('Zanovello i3 Cortical') and SH-SY5Y and SK-N-BE(2) cells using the RNeasy Mini Kit (Qiagen) following the manufacturer's protocol including the on-column DNA digestion step. RNA concentrations were measured by NanoDrop, and 1,000 ng of RNA was used for reverse transcription. Samples undergoing RNA-seq were furthermore assessed for RNA quality on a TapeStation 4200 (Agilent), resulting in an RNA integrity number (RIN) higher than 9.4 for all samples. Sequencing libraries were prepared with poly(A) enrichment using the TruSeq Stranded mRNA Prep Kit (Illumina) and sequenced on an Illumina HiSeq 2500 or NovaSeq 6000 machine at UCL Genomics with the following specifics: SH-SY5Y cells: 2 × 100 bp, depth >40 million per sample SK-N-BE(2) and 'Zanovello i3 Cortical' cells: 2 × 150 bp, depth >40 million per sample

'Humphrey i3 Cortical' samples were processed as previously described⁶⁹ using the RAPiD-nf Nextflow pipeline. In brief, adapters were trimmed from raw reads using Trimmomatic⁷⁰ version 0.36, and reads were aligned to the GRCh38 genome build using gene models from GENCODE version 30 (ref. ⁷¹) with STAR⁷² version 2.7.2a. The RAPiD-nf pipeline is available at https://github.com/CommonMindConsortium/RAPiD-nf/↗.

The 'Brown' SH-SY-5Y, SK-N-BE(2) and i3Neuron datasets were processed as previously described¹⁰. Unless otherwise stated, all short-read RNA-seq datasets were processed using the following pipeline. Raw reads in FASTQ format were quality trimmed for a minimum Phred score of 10 and otherwise default parameters using fastp⁷³ (version 0.20.1). Quality trimmed reads were aligned to the GRCh38 genome build using gene models from GENCODE version 40 (ref. ⁷¹) with STAR⁷² (version 2.7.8a). Quality trimmed reads are used as input for any tools that require FASTQ files as input (for example, PAPA and Salmon). Our alignment pipeline is implemented in Snakemake⁷⁴ and is available at https://github.com/frattalab/rna_seq_snakemake↗.

SLAM-seq was performed on cortical-like i3Neurons following protocols adapted from Herzog et al.⁵¹. Samples were treated with 100 µM 4SU on day 7 for 0, 1, 4, 8, 12 and 24 hours before immediate wash with PBS. Each timepoint had two replicates for both control and TDP-43 knockdown, excluding 4 hours where one of the control replicates did not pass RNA quality controls and so was not submitted for sequencing.

RNA was extracted using the Qiagen RNA isolation and purification kit. RNA concentration was estimated using a NanoDrop Microvolume Spectrophotometer (Thermo Fisher Scientific). After ensuring an adequate amount of RNA in each sample, iodoacetamide (IAA) treatment was applied to each, facilitating the thiol modification of incorporated 4SU.

Sequencing libraries were prepared with the KAPA RiboErase RNA Hyper Kit and sequenced (2 × 250 bp) on an Illumina NovaSeq SP. Using the 'rna_seq_snakemake' alignment pipeline (https://github.com/frattalab/rna_seq_snakemake↗), raw FASTQ files were quality trimmed using fastp⁷³ with the parameter 'qualified_quality_phred: 10' and aligned without soft clipping to the GRCh38 genome build using STAR⁷² (version 2.7.0f) with gene models from GENCODE version 34 (ref. ⁷¹). GRAND-SLAM (version 2.0.7b) was run on the aligned data using gene models from GENCODE version 34 (ref. ⁷¹) using the '-trim5p 10 -trim3p 10' parameter to ignore mismatches at the ends of reads. The output files containing the estimated new-to-total RNA ratios (NTRs) of each gene were used to estimate the half-life of each gene using the recommended workflow in grandR⁷⁵.

For analyses on specific isoform stability, the reads were aligned to a custom general transcription factor (GTF) containing all 3'UTR isoforms quantified by PAPA (see the 'Identification of cryptic last exons with PAPA' section) using the fastq2EZbakR pipeline (https://github.com/isaacvock/fastq2EZbakR↗, version 0.2.0). Half-lives for the bins aligning to the ELK1 long and short UTR were calculated using the 'EstimateFractions' function from EZbakR⁷⁶ version 0.0.0.9000 to retrieve the fraction of old RNA. Decay constants and 95% confidence intervals for each bin were calculated using a custom script ('isoform_specific_analysis.Rmd' in the 'tdp43-apa' repository) using weighted nonlinear regression. In brief, for each bin and condition, fraction old RNA estimates were inversely weighted proportional to the squared s.e. estimate, and nonlinear least-squares regression was performed to model the fraction remaining as an exponential decay function. We note that this method is used here to detect relative changes in RNA half-lives between conditions and not to provide the exact half-life estimates.

Although there are many tools for de novo alternative polyadenylation detection within 3'UTRs from RNA-seq data, all suffer from poor performance with respect to matched 3' end sequencing approaches^77,78. These tools also cannot detect upstream poly(A) sites or define complete last exon structure. Aptardi is a deep-learning-based approach to refine predicted 3' ends of reference or assembled transcriptomes⁷⁹ but was excluded from a recent benchmarking study due to compute times and resource requirements⁷⁸. TECtool (version 0.4) trains a machine learning model on annotated last exons to classify novel intronic last exons defined upstream of poly(A) sites from the PolyASite atlas²⁴ but can only define ALEs and only supports single-end RNA-seq data, substantially impacting sensitivity. Inspired by findings that general purpose transcript assemblers can sufficiently define individual exons⁸⁰ and a previous workflow combining matched short-read and 3' enriched sequencing¹⁸, our approach extracts last exons from StringTie¹⁹ assembled transcripts and filters based on proximity to 3' end sequencing-derived poly(A) sites. Additionally, we rescue events with poly(A) signal hexamers near the 3' end, an important feature in discriminating 3'UTRs from other transcriptomic regions²⁸ that can also mitigate incomplete coverage of cellular contexts and experimental conditions by 3' sequencing databases.

Transcript assemblies for individual samples were generated using StringTie 2.1.7 (annotation-guided mode). Grouping by experimental condition, a redundant assembly was generated using GffCompare⁸¹ 0.11.2. Next, condition-wise, transcript-level mean transcripts per million (TPMs) were calculated, assigning 0 TPM if absent in a sample. Transcripts were filtered for >1 mean TPM to improve global assembly accuracy⁸². Next, we extracted last exons from sample-wise assembled transcripts and identified novel events that satisfy the following criteria:Predicted PAS does not overlap annotated exons.ALEs—last intron is contained within annotated introns with exactly matching 5'ssIPAs—last exon overlaps annotated exon with a matching 5' end (exact for internal exons, within 100 nucleotides (nt) for first exons due to known imprecision of assembled transcript start sites)3'Ext—overlaps annotated last exon with exactly matching 5' ends and extends the longest exon at the locusIPA and 3'Ext—extends annotated exon by minimum distance (default 100 nt)

Filtered novel last exons were then merged by condition into single GTFs to select a condition-wise representative prediction based on 3' end precision. Last exon 3' ends within 100 nt of PolyASite 2.0 database²⁰ PASs were retained and updated to database coordinates. Alternatively, last exons containing any of the 18 poly(A) signal hexamers²¹ in the final 100 nt were retained, selecting the exons with hexamers closest to the expected 21-nt upstream position.

We then combined the filtered novel and annotated last exons into a combined transcriptome reference. We then defined 'last exon identifiers' based on overlapping regions. Overlapping last exons of each gene were assigned a common identifier, with 3'Exts receiving a unique identifier to the exons the annotated last exons they extend. Regions overlapping annotated first or internal exons were removed to retain only unique last exon sequences. Last exons with 3' ends overlapping annotated first/internal exons were excluded.

Transcript sequences were extracted using GffRead⁸¹ 0.12.1 and used to construct a decoy-aware transcriptome index using Salmon²² 1.5.2 (GRCh38 genome build as decoys). Samples were subsequently quantified using Salmon²² 1.5.2 ('–gcBias' and '–seqBias' flags enabled). TPM values were summed by the last exon identifier, and estimated counts were generated with tximport⁸³ 1.26.0 ('countsFromAbundance=lengthScaledTPM') for differential isoform usage testing with DEXSeq²³ 1.44.0. PAS usage was calculated by dividing last exon isoform expression (TPM) by total gene isoform expression.

PAPA 0.2.0, available at https://github.com/frattalab/PAPA↗, is implemented as a Snakemake⁷⁴ pipeline using PyRanges⁸⁴ 0.0.115 for interval operations and pyfaidx⁸⁵ 0.6.2 and BioPython⁸⁶ 1.79 for genomic sequence operations. Conda environments are used for dependency management.

We ran PAPA in 'identification' mode to predict novel last exons in the i3Neuron, 'Zanovello' SH-SY5Y and SK-N-BE(2) datasets. We provided GENCODE version 40 (ref. ⁷¹) annotations filtered for protein-coding and lncRNA gene transcripts with a 'transcript support level' value ≤ 3 and without the 'mRNA_end_NF' tag⁸⁷.

Predicted last exon GTF files were combined into a single GTF using PAPA's 'combine_novel_last_exons.py' script. All datasets were then quantified and assessed for differential usage using a unified transcriptome reference combining novel and annotated last exons from the filtered GTF. Differential usage was performed using the standard DEXSeq workflow, with the differentiation date added as a covariate for the 'Klim i3 motor' dataset⁴⁰. We defined cryptic APAs as DEXSeq adjusted P < 0.05, mean control usage < 10% and change in mean usage > 10% (TDP-43 knockdown, control). We further manually curated cryptic IPAs, as manual inspection suggested frequent artifacts at regions of reduced coverage in intron retention loci.

TDP-43 knockdown samples from all in vitro datasets were used. Soft-clipped alignments were extracted and 3' ends inferred based on the strandedness of the RNA-seq protocol (reads with soft clips at both ends were excluded from unstranded protocols). PATRs were defined as soft-clipped regions ≥6-nt with ≥80% tail nucleotide content³⁰ (A for rightmost/plus strand; T for leftmost/minus strand) or 3−5-nt overhangs with 100% tail content, with the 3' most-aligned coordinate defining the putative PAS.

PATRs were pooled across datasets and clustered using an iterative approach approximating PolyASite's algorithm²⁰. PASs were extended ±12 nt and overlapped, selecting the position with highest read support as representative. Reads within 12 nt of the representative coordinate were collapsed into a cluster, with the process repeated until all PATRs were assigned.

For cryptic PAS validation, we generated 1,000 covariate-matched annotated PAS samples by stratified sampling without replacement using 'matchRanges' from nullranges⁸⁸ version 1.8.0. We matched for expression (log₂(median TPM + 1)) and the number of unique PASs (separated by ≥12 nt), assessing covariate balance using the 'bal.tab' method from cobalt⁸⁹ version 4.5.5.

We then computed distances between annotated/cryptic PASs and nearest PATR clusters, assigning 0 for overlaps. We reported overlap if one or more PASs passed the distance threshold (10, 25, 50, 100 and 200 nt). At each threshold, we computed a group-wise fraction of overlapping events (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{{p}_{i}}$$\end{document}pi^) and computed two-sided empirical P values to assess whether cryptic and annotated PASs arose from the same distribution as follows:\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=\frac{1}{N}\mathop{\sum }\limits_{i=1}^{N}I(|{\hat{p}}_{i}-\hat{\mu }|\ge |\widehat{{p}_{\rm{obs}}}-\hat{\mu }|)$$\end{document}p=1N∑i=1NI(∣p^i−μ^∣≥∣pobs^−μ^∣)where N = 1,000 (total annotated samples), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\mu }$$\end{document}μˆ = annotated distribution mean \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{{p}_{i}}$$\end{document}pi^, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{{p}_{\rm{obs}}}$$\end{document}pobs^ = cryptic PAS \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widehat{{p}_{i}}$$\end{document}pi^ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I(.)$$\end{document}I(.) is an indicator function.

The PATR extraction pipeline, available at https://github.com/SamBryce-Smith/bulk_polyatail_reads↗ (version 0.1.0), is implemented using Snakemake⁷⁶ version 7.32.4, Python 3.10.13, PyRanges⁸⁶ version 0.0.129, pysam version 0.22.0, pandas version 2.1.4, NumPy version 1.26.3, pyarrow version 15.0.0 and fastparquet version 2024.2.0. Cryptic PAS validation scripts are available under the 'preprocessing' directory at https://github.com/frattalab/tdp43-apa/↗.

Transcript models for ELK1, SIX3 and TLX1 were extracted from National Center for Biotechnology Information RefSeq version 110 annotation. 3'UTR and last exons were overlapped with 3'Ext intervals. If any overlap was detected, the 3' end coordinate of the annotated interval was updated to the 3'Ext 3' end. Upstream transcript intervals were otherwise unmodified. We then analyzed the 'Seddighi i3 Cortical' dataset with DaPars2 (ref. ³¹) using a Snakemake pipeline developed for the APAeval project⁷⁸ (available at https://github.com/iRNA-COSI/APAeval↗). Two separate runs with the original or updated transcript models were performed. BED files of predicted PASs and their relative usages parsed from the DaPars2 output file were used for downstream analysis, extracting the distal events to represent cryptic 3'Ext predictions.

The SH-SY5Y TDP-43 iCLIP data (ArrayExpress: E-MTAB-11243↗) were generated and processed as previously described¹⁰. iCLIP peaks from the two independent replicates were merged into non-redundant intervals for all subsequent analysis.

Cryptic events were defined as last exon isoforms passing cryptic thresholds in any in vitro dataset. The probability of detecting TDP-43 binding events via iCLIP is influenced by the abundance of target RNAs, but, by pooling cryptic events across datasets, we cannot control for the confounding influence of RNA expression between groups. We, therefore, defined background events as isoforms that were assessed for differential usage in all SH-SY5Y datasets and had an adjusted P > 0.05 across all datasets, which biases against observing enriched binding in the cryptic group.

For 3'Ext events, the most distal annotated poly(A) site is selected to represent the proximal site, and background events represent loci with a predicted novel 3'UTR extension. For other event categories, background events include annotated and novel events. Our approach to define a common last exon reference across datasets can result in non-redundant intervals being predicted for the same last exon isoform. We, therefore, implemented a collapsing strategy to define a single representative interval for each event.

First, we filtered for novel predictions matching a PolyASite reference PAS. If distinct reference PASs are reported for the same isoform, the site predicted in the most independent datasets is selected as representative. If distinct sites are detected in the same number of independent datasets, the most proximal site is arbitrarily selected. PolyASite PAS intervals represent clusters. If distinct 3' end predictions overlap with the same PAS cluster, the prediction closest to the PolyASite representative coordinate is selected (most distal prediction is arbitrarily selected in case of ties).

If no isoforms matched a PolyASite PAS, we selected a representative prediction whose poly(A) signal motif minimizes the deviance from the characteristic position 21 nt upstream of the PAS. In case of ties, the most proximal prediction was arbitrarily selected. As distinct intervals still remained for background ALEs and IPAs after 3' end collapsing, we arbitrarily selected the most distal 3' end for nine background IPAs and the most proximal 5' end for four background ALEs.

We constructed TDP-43 binding metaprofiles by extending genomic landmarks by 500 nt in both directions and computing per-position coverage by iCLIP peaks using BEDTools⁹⁰ version 2.31.0. We then calculated mean coverage (fraction of events with an overlapping peak) and s.e. for each position relative to the landmark. We plotted LOESS-smoothed ('span' = 0.1) coverage and confidence intervals (±1 s.e.).

To perform de novo motif enrichment, we adapted PEKA³², which identifies kmers with positional enrichment at iCLIP peaks relative to background crosslink sites while normalizing to the general occurrence in the surrounding genomic context. Therefore, we can substitute iCLIP peaks and global crosslink sites for cryptic and background landmarks, respectively, to identify positionally enriched kmers with respect to cryptic landmarks. For all comparisons, we ran PEKA to search for enriched 6-mers in the proximal window of interest set to 250 nt (the broad window in which iCLIP peaks were observed), and the distal window was set to 500 nt (to maintain consistency with the overall search space for iCLIP peaks). The 'percentile' flag was set to 0 to switch off thresholding of background regions based on read count, and the 'relpos' flag was set to 0 to consider all positions in the proximal window when calculating the enrichment score.

Preferred TDP-43 binding 6-mers were extracted from Halleger et al.⁶. In brief, the 6-mers were defined using PEKA as the top 20 most enriched kmers around intronic iCLIP crosslinks across all wild-type, A326P, G294A, G335A, M337P and Q331K and a 316del346 GFP−TDP-43 in HEK293 cells. The 20 were subsequently separated into the following three groups based on a gradient of enrichment in wild-type and G335A TDP-43 with respect to A326 and 316del346 variants and their consensus sequence:YG-containing [UG]n 6-mers: UGUGUG, GUGUGU, UGUGCG, UGCGUG, CGUGUG, GUGUGCYA-containing [UG]n 6-mers: AUGUGU, GUAUGU, GUGUAU, UGUGUA, UGUAUG, UGCAUGAA-containing [UG]n 6-mers: GUGUGA, AAUGAA, GAAUGA, UGAAUG, AUGAAU, GUGAAU, GAAUGU, UUGAAUwhere 'Y' corresponds to a pyrimidine nucleotide. To assess their overrepresentation among enriched 6-mers relative to cryptic landmarks, we performed a one-sided GSEA using fgsea⁴² version 1.24.0 with default settings for each cryptic landmark. The three 6-mer groups and the union of all three groups were provided as input pathways, and kmers were ranked by their PEKA score. After independent runs for each landmark, Benjamini−Hochberg adjusted P values were calculated with respect to all tested landmarks and 6-mer sets and used to evaluate statistical significance.

To generate maps of coverage of specific kmers, we used cv_coverage⁹¹ version 1.1.0 (https://github.com/ulelab/cv_coverage↗) to scan for occurrences of the YG-containing [UG]n 6-mers in a 500-nt window around cryptic and background landmarks, disabling weighting the occurrence by cDNA count. For coverage plots, the percentage occurrences of each 6-mer were summed separately for the cryptic and background regions. The percentage occurrences were converted to mean coverages and visualized as described for iCLIP maps.

The adapted PEKA code is available at the 'output_mods' branch of the following forked copy of the PEKA repository: https://github.com/SamBryce-Smith/peka↗. A Snakemake pipeline to run PEKA and cv_coverage is available in the 'motifs/peka_snakemake' directory of the 'tdp43-apa' repository.

Sequenced reads from FACS-sorted frontal cortex neuronal nuclei³⁴ were processed as described in Brown et al.¹⁰. The data are available in the Gene Expression Omnibus (GEO) at GSE126543↗.

Nuclear RNA-seq libraries contain both nascent and processed RNA. We, therefore, constructed decoy transcript models that reflect alternative processing decisions at ALE and IPA loci (for example, intron retention) to limit the confounding effect of nascent RNAs on transcript quantification²².

First, we extracted cryptic ALE and IPA coordinates from the unified transcript reference used to quantify cell culture datasets. We then generated decoy transcript models separately for each event type. For IPA events, the unique cryptic IPA region was extended to incorporate the adjacent upstream annotated internal exon. Then, a 'spliced' decoy transcript that traverses the annotated internal exon to the downstream annotated internal exon was generated, alongside an 'intron retention' decoy transcript that contains the same pairs of internal exons merged with the intervening intron. For ALEs, a 'retained intron' decoy transcript was generated that corresponds to the complete intronic region in which the ALE is contained. No decoy transcript models were generated for 3'Ext, 3'shortening and 'complex' events or for ALEs that are the most distal annotated isoform of their gene. Decoy transcript and gene identifiers were appended with suffixes to differentiate from cryptic APAs and annotated transcripts. Finally, the decoy transcripts and cryptic APAs were returned to the unified transcript reference to generate a decoy-augmented last exon reference for quantification.

The decoy-augmented reference was quantified with Salmon version 1.8.0 (ref. ²²) using the 'salmon' sub-pipeline available at https://github.com/frattalab/rna_seq_single_steps↗. As with PAPA, samples are quantified against a decoy-aware transcriptome index with full genome sequence (GRCh38 build) used as decoys⁹² and the '–gcBias' and '–seqBias' flags enabled.

Calculation of percent poly(A) usage (PPAU) was performed using a copy of the 'tx_to_polyA_quant.R' script from the PAPA repository. Sample-wise differences in PPAU were calculated by subtracting PPAU in the TDP-43-positive population from the TDP-43-negative population (that is, a positive difference indicates enrichment in the TDP-43-depleted population). Cryptic APAs with a median sample-wise enrichment of more than 5% were considered as enriched. Scripts to construct decoy transcripts and analyze quantifications are available under the 'postmortem' subdirectory at https://github.com/frattalab/tdp43-apa↗.

The sequencing libraries were generated^35,93 and processed¹³ as previously described. Samples were classified into disease subtypes as previously described¹³. In brief, FTD subtypes were classified by pathology according to the presence of TDP-43 inclusions (FTLD-TDP), FUS or Tau aggregates. Patients with ALS were subcategorized based on presence (ALS-non-TDP) or absence (ALS-TDP) of reported SOD1 or FUS mutations. The following samples were considered as regions where TDP-43 pathology (and specific cryptic junction expression) is expected: motor (ALS-TDP), frontal and temporal cortex samples (FTLD-TDP and ALS-TDP) and cervical, lumbar and thoracic spinal cord samples (ALS-TDP).

We opted to quantify ALE events using junction reads, which provide direct quantification of the occurrence of a splicing event. As of version 0.2, PAPA does not directly report splice junctions associated with ALE events. However, as the filtering criteria applied by PAPA require putative ALE events to have a terminal splice junction with a direct match to an annotated 5'ss, it is possible to infer splice junctions from reference annotation using just the reported last exon coordinates. For ALEs fully contained within annotated introns, the splice junction is defined from the intron start to the start of the ALE. If last exons are distal to the annotated gene, then the closest upstream annotated intron is found. The splice junction is subsequently defined as the region from the intron start to the start of the ALE. Finally, for annotated ALEs, all annotated introns that terminate at the ALE are reported as splice junctions for the event. The above steps are implemented in a custom script, 'last_exons_to_sj.py', available at the 'tdp43-apa' GitHub repository.

Splice junctions for cryptic ALEs and cryptic splice junctions identified in cortical-like i3Neurons¹³ were quantified across the NYGC RNA-seq cohort by extracting counts for provided junctions from the '.SJ.out.tab' files produced by STAR⁷². The code is implemented in the 'bedops_parse_star_junctions' version 0.1.0 Snakemake pipeline and is available at https://github.com/SamBryce-Smith/bedops_parse_star_junctions↗.

We defined detection criteria to prioritize cryptic splice junctions that are specifically in tissue types and samples with expected TDP-43 pathology. Junctions are considered expressed if at least two spliced reads are detected in a sample. Junctions are considered selectively expressed if expressed in at most 0.5% of all samples where TDP-43 pathology is not expected and in at least 1% of samples where TDP-43 pathology is expected. We note that such criteria will exclude events with enriched expression in tissues with expected TDP-43 proteinopathy but that have basal expression in unknown cell types not represented in our in vitro compendium. Such events may still have relevance in mechanisms of disease in specific cell types but are less suitable for discriminating samples with TDP-43 proteinopathy.

i3Neuron Ribo-seq data were generated and processed as previously described¹³. Uniquely mapped reads were assigned to genes based on the union of annotated 'CDS' entries in the GENCODE version 34 standard annotation released using featureCounts⁹⁴ version 2.0.1. Differential expression between TDP-43 knockdown and control was performed using DESeq2 (ref. ⁹⁵) version 1.38.3, and differentially translated genes were defined based on a Benjamini−Hochberg adjusted P value threshold of 0.05. Any last exon passing our cryptic criteria in at least one of the i3 Neuron datasets (Brown i3 cortical, Seddighi i3 cortical, Humphrey i3 cortical) was considered for intersection with differentially translated genes.

GSEA was performed using fgsea⁴² version 1.24.0 with default settings. Cryptic 3'Ext, IPA and ALE containing genes were provided as input pathways, and moderated fold changes were calculated with the 'lfcShrink' function from the DESeq2 package using the default apeglm⁹⁶ method as the shrinkage estimator to rank genes. A threshold of 0.05 Benjamini−Hochberg adjusted P value was used to determine statistical significance.

Read counting was performed using the 'feature_counts' sub-pipeline available at https://github.com/frattalab/rna_seq_single_steps↗. Custom scripts used to perform differential expression and pathway analysis are available at https://github.com/frattalab/tdp43-apa↗.

For cross-referencing with differential RNA expression, we used differential expression analysis from cortical-like i3Neurons performed as previously described¹³. Cryptic last exon-containing genes were highlighted if they passed the statistical significance threshold in the Ribo-seq differential expression analysis.

ELK1 target genes in HeLa cells were accessed from the ChIP-Atlas⁹⁷ on 15 November 2023. We used the 'Target genes' module to obtain a list of target genes that have a ChIP–seq peak within ±1 kb of transcription start sites. The resulting list contained two HeLa datasets (GSM608163↗ and GSM935326↗) and was filtered to target genes identified in both datasets. Given a reported redundancy of function between ELK1 and other members of the TCF family⁹⁸ (ELK3 and, particularly, ELK4), we also attempted to define a unique set of ELK1 target genes. ELK4 target genes in HeLA cells were accessed from ChIP-Atlas on 29 November 2023 using the same parameters. The resulting list contained three HeLa datasets (GSM608161↗, GSM608162↗ and GSM935351↗), and we again filtered for target genes identified in all datasets. ELK3 HeLa ChIP–seq data were not available through ChIP-Atlas at the time of publication and were not considered for further redundancy. ELK3 RNA levels are 10× lower than ELK3 and ELK4 in HeLa TDP-43 knockout cells⁴⁹, so we anticipate that this is unlikely to affect our conclusions. ELK1 and ELK4 target gene lists were intersected to define common and unique target genes for each transcription factor. Final target gene lists used are reported in Supplementary Table 5.

RNA-seq data from HeLa cells with TDP-43 knockout⁴⁹ were accessed from GSE136366↗. The data were processed and differential expression was performed as previously described¹⁰. Genes were ranked by DESeq2's test statistic (log₂ transformed fold change divided by the s.e. of the fold change) after removing genes with differential splicing upon TDP-43 knockout, where we can expect to attribute any changes in gene expression to TDP-43 loss of function. Differentially spliced genes were defined using MAJIQ⁹⁹, considering any genes with a probability greater than 0.95 as differentially spliced. The target gene sets described above were used as input pathways to fgsea⁴² version 1.24.0 using default settings.

The neural progenitor cell short-read Frac-seq data⁵⁶ were downloaded from the GEO at accession number GSE244655↗. RNA-seq quality control and processing was performed as previously described (see 'RNA-seq data processing' section). The PAPA index was used to quantify ELK1 isoform expression with Salmon version 1.8.0, using the 'salmon' sub-pipeline available at https://github.com/frattalab/rna_seq_single_steps↗. TPM values for the ELK1 3'Ext were pooled across ribosome-associated fractions (monosome, light polysome and heavy polysome), and PPAU was recalculated for each fraction and replicate. All ELK1 3'Ext PPAU values were then normalized to the cytosol PPAU within each replicate for subsequent visualization. Statistical significance was evaluated using a two-sided one-sample t-test after log transforming the PPAU ratios, testing the null hypothesis that the mean is equal to log(1).

Our study design involved multiple stages. First, we used transcriptome-wide hypothesis testing of high-throughput RNA-seq datasets to identify a panel of TDP-43-sensitive cryptic polyadenylation events. We performed this screen in neuronal cell models, where we could reliably deplete TDP-43 levels to mimic nuclear loss in disease. We then screened this panel in specialized and bulk postmortem tissue datasets to highlight events whose expression patterns were consistent with disease and TDP-43 pathology status. Finally, we performed targeted experimental assays to validate observations from high-throughput sequencing and to investigate the molecular consequences of specific cryptic polyadenylation events.

Sample sizes for postmortem tissue analysis (Fig. 2b) were determined by the availability of samples at the time of analysis. Sample size for the NYGC ALS Consortium was determined by the number of available samples at the time of analysis (corresponding to a subset of the 21 February 2023 data freeze) as data collection is still ongoing. Sample sizes for novel omics datasets and experimental validation were determined based on previous studies succeeding with similar aims to identify novel isoforms, perform targeted validation and assess their downstream effects on RNA and protein expression¹⁰.

All statistical tests were performed two-sided. One-sample t-tests were performed using log-transformed ratios of within-replicate, control-normalized values (mean count for FISH experiments and percent PAS usage for Frac-seq). Log transformation is a standard transformation to bring a distribution closer to a normal distribution, but the assumption of normally distributed transformed data was not formally tested. For Student's unpaired t-test (3'RACE experiments), equal variances were assumed, and the data distribution was assumed to be normal, but this was not formally tested. Unless otherwise stated, the Benjamini−Hochberg multiple-testing correction method was used to compute 'adjusted' P values.

Randomization was not used in this study, as most of the analyses (experimental and omics-based) were carried out in cell lines that are inherently homogenous. Randomization was not applicable in postmortem analyses as the variable of interest (disease status and expected TDP-43 pathology) is an observed variable, and no intervention was performed. No data were excluded from analysis. FISH images were analyzed blinded to TDP-43 depletion status. For all other experiments, the investigators were not blinded to experimental condition or disease status during experimentation and analysis.

Further information on research design is available in thelinked to this article. Nature Portfolio Reporting Summary