Dissection of molecular and histological subtypes of papillary thyroid cancer using alternative splicing profiles

Abstract

Despite growing evidence of the relevance of alternative splicing (AS) to cancer development and progression, the biological implications of AS for tumor behaviors, including papillary thyroid cancer (PTC), remain elusive. With the aim of further understanding the molecular and histological subtypes of PTC, we in this study explored whether AS events might act as new molecular determinants. For this purpose, AS profiles were analyzed in RNA-sequencing data from The Cancer Genome Atlas (TCGA) and from a Korean patient dataset. A total of 23 distinct exon-skipping (ES) events that correlated significantly with PTC oncogenic activity and differentiation scores were identified. The two top-ranked ES events, NUMA1_17515 in exon 18 of NUMA1 and TUBB3_38175 in exon 6 of TUBB3, showed high correlations with oncogenic activities and discriminated histological and molecular subtypes of PTC. Furthermore, two novel intron-retention (IR) events for TUBB3 were uncovered. All ES and IR events for the TUBB3 gene were predicted to induce nonsense-mediated mRNA decay. The relative abundances of intron reads in the PTC dataset from TCGA showed IR levels to differ significantly among PTC subtypes, possibly reflecting their different tumor behaviors. This study provides a landscape of AS changes among PTC subtypes and identified two significant AS events, NUMA1_17515 and TUBB3_38175, as potential AS biomarkers for PTC subclassification and characterization. The AS events identified in this study may be involved in the development of phenotypic differences underlying the functional characteristics and histological differentiation of PTCs.

Fig. 1. Overall features of alternative splicing (AS) events in molecular subtypes of papillary thyroid cancer (PTC).

a Numbers of AS events that differed significantly among molecular subtypes of PTC (by ANOVA). AP alternate promoter, AT alternate terminator, ES exon skip, RI retained intron, AD alternate donor site, AA alternate acceptor site, ME mutually exclusive exons. b Top five molecular subtype pairs that exhibited differences in ES events (by Tukey’s test). c Motif enrichment profiles of BRUNOL4 and BRUNOL5 for significant ES events in THCA.1 vs. THCA.4. The X-axis indicates the intron position up/downstream of the target exon (red box). Solid red and blue lines indicate motif scores for 37 included and 78 excluded ES events between THCA.1 and THCA.4, respectively. d Comparison of CELF4 mRNA expression levels between THCA.1 and THCA.4 cells. The p value was determined using the t test.

Fig. 2. Association of the NUMA1_17515 exon skipping event with biological features of papillary thyroid cancer (PTC).

a Correlations of the NUMA1_17515 event with thyroid differentiation scores (TDSs) and ERK and RAS activity scores. Each dot represents an individual patient sample. Pearson’s correlation coefficient (r) and the linear regression line are shown for each plot. b Percent spliced-in (PSI) values for molecular (left) and histological (right) subtypes. The p values were determined by ANOVA. c Correlation between NUMA1_17515 PSI values and NUMA1 expression levels. Histological subtypes are indicated by different colors. d Association of NUMA1_17515 PSI values with progression of extrathyroidal extension (left) and regional lymph node metastasis (right). Extrathyroidal extension: 0 = none; 1 = minimum (T3); 2 = advanced (T4a and T4b). Regional lymph node metastasis: 0 = negative; 1 = positive. e Correlations of the NUMA1_17515 event with TDSs, histological subtype, and NUMA1 expression level in the Seoul National University Hospital (SNUH) thyroid cancer dataset.

Fig. 3. Validation of the NUMA1_17515 exon-skipping event in thyroid cancer cell lines.

a Isoforms of the NUMA1 gene in the exon expression panel of the GTEx portal. The location of NUMA1_17515 is marked by a red arrowhead (chr11:72012401–72012442, hg38). Only one transcript (ENST00000393695, red asterisk) had an NUMA1_17515-regulated exon. Exon numbers below the NUMA1 gene structure are in accordance with TCGA SpliceSeq. b Primers were designed to amplify sequences in neighboring exons of exon 18 (blue arrows). The exon targeted by NUMA1_17515 is colored in red. A primer (right) across exons 19 and 20 was designed. Numbers below each exon represent sizes produces using the primer set. The expected lengths of the RT–PCR products are 169 and 211 bp, excluding or including exon 18 (42 bp), respectively. c Gel electrophoresis of RT–PCR products from three thyroid cancer cell lines (TT2609-C02, FTC-133, and ML-1). d Confirmation of the inclusion of exon 18 in the 211-bp amplicon and its exclusion from the 169-bp amplicon based on Sanger sequencing.

Fig. 4. TUBB3_38175 exon skipping (ES) event as a classifier for molecular subtypes of papillary thyroid cancer (PTC).

a Random forest variable importance plot showing significant ES events for classification using the Gini importance score as an indicator of the subtype-discriminating ability. b TUBB3_38175 percent spliced-in (PSI) value profiles for molecular (left) and histological (right) subtypes. The p values were determined by ANOVA. c Isoforms of the TUBB3 gene from the exon expression panel of the GTEx portal. The location of TUBB3_38175 is marked by the red box (chr11:72012401–72012442, hg38) (ENST00000555609). Sequence-based annotation of exon 6 in the ASpedia database indicates the sites of nonsense-mediated decay (NMD) (blue dots). The transcript with a shorter form of exon 6 (exon 6b) is indicated by a solid arrow (ENST00000557490). d Correlation between TUBB3_38175 PSI values and TUBB3 expression levels. Histological subtypes are indicated by different colors. e Association of TUBB3_38175 PSI values with progression of extrathyroidal extension (left) and regional lymph node metastasis (right).

Fig. 5. Validation of exon-skipping events and identification of novel intron retention (IR) events in TUBB3.

a Validation of exon 6 skipping. Primers were designed to amplify sequences in neighboring exons of exon 6 (blue arrows). The exon targeted by TUBB3_38175 is colored in red. Numbers below each exon represent sizes from the primer. Agarose gel electrophoresis of RT–PCR products from three thyroid cancer cell lines showed a transcript of the expected size (160 bp), which was confirmed by Sanger sequencing. b Identification of new IR events between exons 5 and 7 of the TUBB3 gene (details available in Supplementary Additional file 1). Sequence alignment of three RT–PCR amplicons (m1–3) to the human genome (hg38) confirmed two types of novel IR events between exons 5 and 7, namely, retention of a full-length intron and an intron that was 231 bp shorter. For GENCODE transcripts (left), blue lines indicate coding transcripts, and red lines display problem transcripts. In the right plots, red colors indicate IR regions. c IR ratios among papillary thyroid cancer (PTC) histological subtypes. The p values were calculated using the t test.

Fig. 6. Model of dynamic alternative splicing regulation of the TUBB3 gene.

Five kinds of transcripts (T0–T4) were identified between exon 5 and exon 7. Three (T0, T1, and T2) were generated by exon-skipping events, and two (T3 and T4) were produced by intron-retention events. Red boxes represent unique sequences compared with that of the T0 transcript. The positions of premature stop codons are shown as red stop signs. The green Pac-Man denotes isoforms subject to nonsense-mediated decay (NMD). Only the T0 transcript is predicted to avoid NMD, whereas the others are likely to be subjected to transcriptional degradation due to premature stop codons.