Abstract

   Isoform switching events with predicted functional consequences are
   common in many cancers, but characterization of switching events in
   esophageal adenocarcinoma (EAC) is lacking. Next-generation sequencing
   was used to detect levels of RNA transcripts and identify specific
   isoforms in treatment-naïve esophageal tissues ranging from
   premalignant Barrett’s esophagus (BE), BE with low- or high-grade
   dysplasia (BE.LGD, BE.HGD), and EAC. Samples were stratified by
   histopathology and TP53 mutation status, identifying significant
   isoform switching events with predicted functional consequences.
   Comparing BE.LGD with BE.HGD, a histopathology linked to cancer
   progression, isoform switching events were identified in 75 genes
   including KRAS, RNF128, and WRAP53. Stratification based on TP53 status
   increased the number of significant isoform switches to 135, suggesting
   switching events affect cellular functions based on TP53 mutation and
   tissue histopathology. Analysis of isoforms agnostic, exclusive, and
   shared with mutant TP53 revealed unique signatures including
   demethylation, lipid and retinoic acid metabolism, and glucuronidation,
   respectively. Nearly half of isoform switching events were identified
   without significant gene-level expression changes. Importantly, two
   TP53-interacting isoforms, RNF128 and WRAP53, were significantly linked
   to patient survival. Thus, analysis of isoform switching events may
   provide new insight for the identification of prognostic markers and
   inform new potential therapeutic targets for EAC.

   Keywords: MT: Bioinformatics, esophageal adenocarcinoma, Barrett’s
   esophagus, isoform switching, transcriptomics, TP53, Mortality-linked
   isoforms

Graphical abstract

   graphic file with name fx1.jpg
   [43]Open in a new tab
     __________________________________________________________________

   Kresty and colleagues performed the first comprehensive
   characterization of isoform switching events in an esophageal
   adenocarcinoma progression series. Results were stratified by
   histopathology and TP53 mutation status revealing unique signatures,
   including demethylation, lipid and retinoic acid metabolism, and
   glucuronidation. Importantly, TP53-interacting isoforms, RNF128 and
   WRAP53, were linked with patient survival.

Introduction

   Esophageal cancer is the seventh most common cancer worldwide and one
   of the most lethal cancers.[44]^1 The two main subtypes of esophageal
   cancer include adenocarcinoma (EAC) and squamous cell carcinoma (ESCC).
   Increasing incidence of EAC has been observed throughout Westernized
   countries, including the United States, United Kingdom, and Australia
   in recent years, whereas ESCC has been declining with reduced tobacco
   use.[45]^2 EAC is the predominant subtype of esophageal cancer in the
   United States and the seventh leading cause of cancer-related deaths
   among US men.[46]^3 Moreover, the 5-year survival rate following a
   diagnosis of EAC remains poor at 19.9%.[47]^4 EAC is strongly
   associated with gastroesophageal reflux disease (GERD), obesity, and
   male gender, and frequently develops subsequent to long-standing
   Barrett’s esophagus (BE).[48]^5 BE is a metaplastic replacement of the
   normal squamous esophageal cells by specialized columnar epithelium,
   usually resulting from chronic GERD.[49]^5^,[50]^6 BE represents the
   only known precursor lesion to EAC.[51]^7 Importantly, cancer risk and
   mutational burden increase with progression beyond BE metaplasia to BE
   with low-grade dysplasia (LGD) and high-grade dysplasia
   (HGD).[52]^8^,[53]^9 For BE patients with LGD, the incidence rate of
   EAC is as low as 1.70% per patient-year, whereas the incidence rate of
   EAC dramatically increases to 6.58% per patient-year among BE patients
   with HGD.[54]^10 In addition, many patients diagnosed with BE and HGD
   have occult EAC, further supporting the elevated risk associated with
   BE accompanied by advanced pathologic changes.[55]^11 For these
   reasons, increasing our understanding of the transition from low- to
   high-risk BE pathology or BE with LGD to HGD and EAC may identify
   patterns of progression and in turn facilitate early detection and
   potentially a window for intervention.

   EAC is a cancer with one of the highest mutational
   burdens,[56]^12^,[57]^13 but unlike other similar cancers, the
   mutational burden is spread across a large number of genes with
   relatively low mutational frequencies, except for TP53.[58]13, [59]14,
   [60]15, [61]16 This coupled with the high spatial and temporal
   heterogeneity of EAC has hampered clinical applications targeting
   mutational drivers.[62]^17 Despite improved understanding of the
   genomic landscape contributing to progression from BE to EAC in recent
   years, we still have a poor understanding of how to leverage this
   information for improved risk stratification, or to inform targeted
   strategies for esophageal cancer prevention and treatment. Thus, in the
   current study we sought to move beyond single gene transcript analysis
   to characterize differential alternatively spliced isoforms of genes in
   esophageal tissues based on histological progression and TP53 mutation
   status. This approach permits an investigation of the individual splice
   variants or isoforms that contribute to gene expression and in some
   cases reveals differentially expressed isoforms masked by unchanged
   gene expression.

   Historically, the lack of functional annotation of gene isoforms and
   the dearth of sophisticated isoform analysis tools resulted in
   underutilization of RNA-sequencing (RNA-seq) data.[63]^18 More
   recently, decreased sequencing costs coupled with increased
   availability of next-generation sequencing platforms and advanced
   RNA-seq analysis tools have made transcriptomics quantification with
   isoform resolution possible.[64]^19 Several recent investigations have
   performed isoform switch analysis using TCGA data revealing that
   isoform switches are highly prominent in many cancers and can have
   profound biological impacts.[65]^18^,[66]^20^,[67]^21 To date,
   comprehensive analysis of isoform switching has not been conducted in
   EAC or BE precursor lesions. Herein, we investigated isoform switching
   events using RNA-seq data collected from 57 esophageal tissue samples
   derived from 46 treatment-naïve patients undergoing esophagectomy
   following a diagnosis of HGD or EAC. Significant isoform switching
   events with predicted functional consequences were identified comparing
   patients stratified by histopathology and TP53 mutation status.
   Multiple studies link aberrant TP53 expression with malignant
   progression of BE, particularly BE with HGD.[68]22, [69]23, [70]24,
   [71]25, [72]26 Overall, TP53 mutation remains the strongest known
   driver of EAC progression,[73]^6 with a mutation frequency over 70% in
   EAC patients.[74]^13 Incorporation of TP53 mutation status in the
   analyses markedly increased detection of isoform-switched genes and
   expression of specific isoforms linked to TP53 were significantly
   associated with patient survival. Gene Ontology (GO) enrichment
   analysis of genes associated with differential isoforms revealed both
   shared and unique cancer relevant alterations when comparing BE with
   LGD to BE with HGD histology. In turn, targeting identified processes
   may offer new leads for inhibiting BE progression. As
   proof-of-principle, we selected two novel agents informed by the study
   results and evaluated efficacy as inhibitors of EAC viability. Results
   indicate isoform switching events may provide new molecular insights
   and potentially identify prognostic markers or new therapeutic targets
   for EAC.

Results

Isoform switching events and enrichment analysis based on histologic
progression and TP53

   Isoform switching analysis was performed comparing patient tissues
   categorized as BE with metaplasia and low-grade dysplasia (BE.LGD) to
   BE high-grade dysplasia (BE.HGD), alone or with stratification based on
   TP53 mutation status as summarized in [75]Figure 1. Patient demographic
   data and frequency of TP53 mutation for each group are summarized in
   [76]Table S1. When comparing BE.LGD with BE.HGD, 75 genes show
   significant isoform switching events with predicted functional
   consequences ([77]Figure 1A). Moreover, with the inclusion of TP53
   mutation status, significant isoform switching events increased to 135
   genes when comparing BE.LGD TP53 wild-type (WT) with BE.HGD TP53 mutant
   (MUT). Next, GO pathway enrichment analysis was performed
   ([78]Figures 1B, [79]Tables S2, and [80]S3) revealing that most of the
   top enriched GO processes are shared between the two comparisons, with
   differences noted based on the significance of flavonoid metabolic
   process observed in BE.LGD versus BE.HGD and retinoic acid metabolic
   process observed in BE.LGD TP53 WT versus BE.HGD TP53 MUT. Dominant
   shared pathways included processes related to glucuronidation and
   metabolism. For BE.LGD versus BE.HGD, 14.7% (11 of 75) of top
   significant isoform-switched genes have an absolute differential
   isoform fraction (dIF) value greater than 0.25 ([81]Figure 1C) and
   nearly all of the top 20 significant genes involved in isoform
   switching are protein-coding genes with the exception of MIR4458HG
   ([82]Figure 1D). Similarly, when comparing BE.LGD TP53 WT and BE.HGD
   TP53 MUT, 17.0% (23 of 135) of top significant isoform-switched genes
   have an absolute dIF value greater than 0.25 ([83]Figure 1E) and again
   protein-coding genes dominate, except for MIR4458HG and CENPBD1P1
   ([84]Figure 1F). Among all identified isoform-switched genes, only
   isoform switching of RNF128 was previously reported in HGD and
   EAC.[85]^27 Despite well-characterized gene-level changes in KRAS with
   EAC progression, isoform switching events were not previously
   reported.[86]28, [87]29, [88]30, [89]31 Additional isoform switch
   events reported previously in other cancer targets were also identified
   in our results, including TPM4 and MINDY1.[90]^27^,[91]32, [92]33,
   [93]34, [94]35, [95]36 Isoform switching of TP53 is not observed in
   BE.LGD TP53 WT versus BE.HGD TP53 MUT, potentially because 84.6%
   (22/26) of TP53 mutations are missense mutations failing to result in
   further transcriptional level changes in TP53 isoforms ([96]Figure S1).
   The high frequencies of missense TP53 mutations and mutation hotspots
   observed in our cohort align with previously published studies in
   EAC.[97]37, [98]38, [99]39 Significant functional consequences
   associated with dominant isoforms identified in BE.HGD and BE.HGD TP53
   MUT include complete open reading frame (ORF) loss, signal peptide
   loss, and the loss of protein-coding isoforms ([100]Figure 1G). Global
   splicing enrichment analysis also reveals that isoform switching events
   are significantly associated with changes in alternative transcription
   start sites (ATSS) and alternative transcription termination sites
   (ATTS) in BE.LGD versus BE.HGD, and changes in alternative 5′ splice
   site (A5) and ATTS in BE.LGD TP53 WT versus BE.HGD TP53 MUT
   ([101]Figure 1H).

Figure 1.

   [102]Figure 1
   [103]Open in a new tab

   Identification of significant isoform switching events with predicted
   functional consequences based on histopathology and TP53 status

   (A) Number of isoform-switched genes detected in Barrett’s esophagus
   with metaplasia and low-grade dysplasia (BE.LGD) compared with BE with
   high-grade dysplasia (BE.HGD) and with inclusion of TP53 mutation
   status, wild-type (WT) or mutant (MUT). (B) Top 10 GO processes
   identified for each comparison (p < 0.001 and FDR <0.001), ranked by
   FDR value. (C) Volcano plot showing the differential isoform fraction
   (dIF) values comparing BE.LGD with BE.HGD. Each dot represents a
   transcript involved in isoform switching. In red are isoforms with
   |dIF| ≥ 0.1 and FDR ≤0.05. (D) Dot plot showing top 20 genes with
   significant isoform switching events for BE.LGD versus BE.HGD. The
   shape of each dot represents gene category and the color of each dot
   indicates significance level (|dIF| ≥ 0.1 and FDR ≤ 0.05). (E) Volcano
   plot showing dIF values in BE.LGD TP53 WT versus BE.HGD TP53 MUT. Each
   dot represents a transcript involved in isoform switching. In red are
   isoforms with |dIF| ≥ 0.1 and FDR ≤0.05. (F) Dot plot showing top 20
   genes with significant isoform switching events BE.LGD TP53 WT versus
   BE.HGD TP53 MUT (|dIF| ≥ 0.1 and FDR ≤ 0.05). (G) Dot plot showing
   functional consequences analysis associated with isoform switching. (H)
   Dot plot showing global alternative splicing event analysis. The size
   of the dot indicates the number of isoform-switched genes in each
   category and the line indicates 95% confidence interval. Significant
   changes in functional consequences and global alternative splicing
   events are colored in red (p ≤ 0.05).

   In addition, gene set enrichment analysis (GSEA) on RNA- and
   ribonucleoprotein (RNP)-related processes were performed on the
   gene-level data, as well as differential gene expression analysis on
   RNA-binding proteins (RBPs), because of the role that RNPs and RBPs
   play on alternative splicing.[104]^40^,[105]^41 GSEA analysis indicates
   significant enrichment of alternative-splicing-associated pathways,
   such as ribonucleoprotein complex biogenesis, spliceosomal small
   nuclear RNP (snRNP) assembly, mRNA splice site selection, and
   regulatory RNA binding ([106]Figure 2A, [107]Table S4, and [108]S5).
   Sixty-one significant pathways were identified in BE.LGD versus BE.HGD,
   whereas 37 significant pathways were identified in BE.LGD TP53 WT
   versus BE.HGD TP53 MUT; thus, enriched pathways were more specific in
   BE.LGD TP53 WT versus BE.HGD TP53 MUT ([109]Table S4 and [110]S5). A
   strong upregulation of leading-edge genes in spliceosomal snRNP
   assembly, mRNA splice site selection, spliceosomal complex assembly,
   and regulatory RNA-binding pathways, as shown in [111]Figures 2B is
   observed in BE.HGD compared with BE.LGD patient samples. Moreover, 54
   and 73 RBPs are significantly differentially expressed in BE.LGD versus
   BE.HGD and BE.LGD TP53 WT versus BE.HGD TP53 MUT, respectively
   ([112]Figures 2C and 2D). Among all identified RBPs, 18.5% (10 of 54)
   and 28.8% (21 of 73) of RBPs had predicted protein interactions with
   isoform-switched genes in BE.LGD versus BE.HGD and BE.LGD TP53 WT
   versus BE.HGD TP53 MUT, respectively. Taken together, both GSEA and
   differential gene expression analysis revealed that expression changes
   of RNA regulatory proteins and RBPs may contribute to observed isoform
   switching events.

Figure 2.

   [113]Figure 2
   [114]Open in a new tab

   Changes associated with RNA regulatory proteins and RNA-binding
   proteins on the gene-level data

   (A) GO pathway enrichments associated with RNA regulatory proteins
   identified in BE.LGD versus BE.HGD and BE.LGD TP53 WT versus BE.HGD
   TP53 MUT (p-value ≤ 0.05 and FDR ≤0.25). Pathways are colored and
   ranked by FDR. (B) Heatmap of the leading-edge genes identified in
   select pathways in the GO pathway enrichment analysis. Dot plot on the
   right shows the pathways that the gene belongs. (C) Volcano plot of the
   expression of RNA-binding proteins (RBPs) in BE.LGD versus BE.HGD. (D)
   Volcano plot of the expression of RBPs in BE.LGD TP53 WT versus BE.HGD
   TP53 MUT. Red dot represents significantly differentially expressed
   RBPs (|log[2]FC| > 0.585 and FDR ≤0.05). Gene names of the genes that
   interact with isoform-switched genes predicted by the STRING database
   are shown.

   Last, qRT-PCR was used to validate isoform switching of RNF128, KRAS,
   and TRIM29 using a subset of patients in the cohort ([115]Figure S2).
   Isoform fractions using qRT-PCR results were comparable to isoform
   fractions using RNA-seq data by the analysis tool.

Significant isoform switching events without differential gene expression

   Interestingly, 42.7% (32 of 75) of significant isoform-switched genes
   in BE.LGD versus BE.HGD ([116]Table 1) and 45.9% (62 of 135) of genes
   in BE.LGD TP53 WT versus BE.HGD TP53 MUT ([117]Table 2) are not
   differentially expressed on the gene level between the two conditions.
   GO enrichment analysis was performed for isoform-switched genes that
   occur in the absence of gene-level expression changes ([118]Tables 1
   and [119]2), with results closely paralleling those shown in
   [120]Figure 1B (data not shown). In brief, isoform switch analysis
   independent of gene-level changes in BE.LGD versus BE.HGD supported top
   enriched biologic processes focused on hormone-related processes,
   cellular glucuronidation, and histone and protein demethylation and
   dealkylation processes. Similarly, isoform-switched genes identified in
   the comparison of BE.LGD TP53 WT versus BE.HGD TP53 MUT without
   gene-level alterations show GO biologic processes overlapping with
   those in [121]Figure 1B, but with less statistical significance likely
   due to reduced power with decreased sample numbers (data not shown).

Table 1.

   Significant isoform switching events without altered gene expression in
   BE.LGD versus BE.HGD
   Gene name Isoform ID dIF FDR Gene name Isoform ID dIF FDR
   MIR4458HG ENST00000502001 0.326 <0.001 KDM6B ENST00000254846 −0.123
   0.020
   ENST00000652260 −0.355 <0.001 ENST00000575521 0.146 0.039
   RPL22L1 ENST00000463836 −0.149 <0.001 [122]AC103810.1 ENST00000585187
   −0.201 0.021
   ENST00000475836 0.124 <0.001 PCSK5 ENST00000376752 0.134 0.022
   [123]AC234782.4 ENST00000674488 −0.259 <0.001 ENST00000674117 −0.147
   0.023
   ENST00000674271 0.248 0.024 [124]AC073896.1 ENST00000548360 0.125 0.023
   SLC25A35 ENST00000581320 0.152 0.009 ENST00000549318 −0.124 0.024
   ENST00000577745 −0.131 0.025 HSD11B1 ENST00000615289 0.195 0.023
   ZNF350 ENST00000243644 −0.138 0.009 ENST00000367027 −0.199 0.029
   ENST00000599258 0.125 0.013 RBP2 ENST00000511956 0.267 0.024
   SNHG18 ENST00000508179 0.157 0.010 [125]AC009509.2 ENST00000536922
   0.229 0.026
   ENST00000655411 −0.246 0.015 UGT2B7 ENST00000305231 −0.217 0.026
   TMEM104 ENST00000335464 −0.182 0.011 GHR ENST00000537449 −0.147 0.026
   ENST00000584246 0.119 0.016 SLC2A5 ENST00000484798 −0.219 0.031
   COX19 ENST00000466146 0.117 0.012 [126]AP000866.2 ENST00000504932 0.163
   0.033
   WRAP53 ENST00000316024 0.126 0.014 CADM3-AS1 ENST00000415675 −0.128
   0.033
   ENST00000396463 −0.154 0.040 ENST00000609696 0.128 0.035
   DISP2 ENST00000559721 0.230 0.014 RSPH1 ENST00000291536 0.101 0.038
   LINC02542 ENST00000660534 −0.195 0.019 CFAP46 ENST00000486104 −0.132
   0.041
   ENST00000668266 0.120 0.042 CA4 ENST00000585705 0.147 0.041
   NEU4 ENST00000407683 −0.239 0.019 CNTN1 ENST00000347616 −0.167 0.047
   [127]AC093866.1 ENST00000513572 0.153 0.019 RTN2 ENST00000245923 −0.105
   0.048
   ENST00000612706 −0.102 0.028 NUDT18 ENST00000613958 0.108 0.049
   GDF15 ENST00000252809 −0.145 0.020 ENST00000611621 −0.108 0.049
   [128]Open in a new tab

Table 2.

   Significant isoform switching without altered gene expression in BE.LGD
   TP53 WT versus BE.HGD TP53 MUT
   Gene name Isoform ID dIF FDR Gene name Isoform ID dIF FDR
   MIR4458HG[129]^a ENST00000502001 0.326 <0.001 NEU4[130]^a
   ENST00000407683 −0.239 0.019
   ENST00000652260 −0.355 <0.001 FSIP2-AS1 ENST00000429929 −0.152 0.019
   MINDY1 ENST00000361936 −0.135 <0.001 [131]AC015802.4 ENST00000592622
   −0.200 0.019
   ENST00000470877 0.100 <0.001 ALDH1A2 ENST00000558239 0.110 0.019
   EVI2B ENST00000330927 0.108 <0.001 PTRH1 ENST00000543175 0.202 0.019
   ENST00000577894 −0.108 <0.001 ZNF324 ENST00000196482 0.136 0.020
   CENPBD1P1 ENST00000493504 −0.103 <0.001 ENST00000593925 −0.126 0.034
   ENST00000651608 0.117 <0.001 GDF15[132]^a ENST00000252809 −0.145 0.020
   [133]AC234782.4[134]^a ENST00000674488 −0.259 <0.001
   [135]AC103810.1[136]^a ENST00000585187 −0.201 0.021
   ENST00000674271 0.248 0.024 PCSK5[137]^a ENST00000376752 0.134 0.022
   LINC02615 ENST00000505133 0.110 <0.001 ENST00000674117 −0.147 0.023
   ENST00000513851 −0.157 0.024 DCDC2 ENST00000378454 0.174 0.022
   TSPAN6 ENST00000494424 −0.114 <0.001 ENST00000378450 −0.177 0.025
   ENST00000614008 0.108 0.021 LTB4R2 ENST00000527924 0.108 0.023
   SMIM6 ENST00000556126 0.116 0.001 CRYBA2 ENST00000392096 −0.202 0.024
   ENST00000579469 −0.116 0.001 GHR[138]^a ENST00000537449 −0.147 0.026
   SESN1 ENST00000436639 −0.133 0.001 NEURL2 ENST00000545238 −0.240 0.028
   TMED1 ENST00000588289 0.103 0.002 ENST00000372518 0.240 0.029
   ENST00000214869 −0.113 0.019 MSH5-SAPCD1 ENST00000476085 −0.137 0.028
   [139]AC026254.2 ENST00000666005 0.220 0.003 ENST00000493662 0.106 0.031
   LINC00964 ENST00000657752 −0.128 0.005 TNNI1 ENST00000622580 −0.180
   0.029
   ENST00000663940 0.107 0.009 ENST00000555340 0.164 0.043
   TMEM241 ENST00000577448 0.100 0.005 SFTA2 ENST00000634371 −0.113 0.030
   ENST00000475185 0.128 0.029 ENST00000359086 0.106 0.037
   PLA2G15 ENST00000562966 0.103 0.006 CCDC170 ENST00000537358 −0.269
   0.030
   ENST00000219345 −0.115 0.044 TLCD1 ENST00000292090 0.153 0.031
   ANKRD54 ENST00000609706 0.125 0.006 SOCS4 ENST00000395472 0.121 0.031
   ZNF350[140]^a ENST00000243644 −0.138 0.009 [141]AP000866.2[142]^a
   ENST00000504932 0.163 0.033
   ENST00000599258 0.125 0.013 C1orf53 ENST00000436652 0.147 0.033
   GAS2 ENST00000278187 0.155 0.009 TRIM29 ENST00000532195 −0.122 0.036
   ATP9B ENST00000586722 0.115 0.009 EFHC2 ENST00000420999 0.221 0.037
   SPATA17 ENST00000470448 0.118 0.011 ENST00000343571 −0.221 0.038
   ENST00000491809 −0.186 0.048 RSPH1[143]^a ENST00000291536 0.101 0.038
   TMEM104[144]^a ENST00000335464 −0.182 0.011 CARNS1 ENST00000307823
   −0.190 0.040
   ENST00000584246 0.119 0.016 CFAP46[145]^a ENST00000486104 −0.132 0.041
   COX19[146]^a ENST00000466146 0.117 0.012 CA4[147]^a ENST00000585705
   0.147 0.041
   WRAP53[148]^a ENST00000316024 0.126 0.014 ONECUT1 ENST00000305901
   −0.219 0.043
   ENST00000396463 −0.154 0.040 LINC00278 ENST00000652068 −0.135 0.043
   DISP2[149]^a ENST00000559721 0.230 0.014 TMEM107 ENST00000437139 0.106
   0.045
   BCL7A ENST00000432926 −0.197 0.015 ZBTB18 ENST00000622512 −0.102 0.046
   ENST00000261822 0.160 0.039 ENST00000358704 0.102 0.047
   [150]AC117453.1 ENST00000651749 −0.182 0.016 RTN2[151]^a
   ENST00000245923 −0.105 0.048
   ENST00000656199 0.151 0.021 LINC01186 ENST00000670979 0.115 0.049
   SYT15 ENST00000449358 −0.138 0.017 CCDC190 ENST00000524710 −0.139 0.050
   ENST00000374321 0.131 0.026
   S100A5 ENST00000368717 −0.247 0.017
   ENST00000368718 0.247 0.020
   [152]Open in a new tab
   ^a

   Genes shared between BE.LGD versus BE.HGD and BE.LGD TP53 WT versus
   BE.HGD TP53 MUT. BE.LGD: Barrett’s esophagus with metaplasia or
   low-grade dysplasia; BE.HGD: Barrett’s esophagus with high-grade
   dysplasia.

   Moreover, several identified isoform-switched genes that do not show
   overall differential gene expression have documented interactions with
   TP53, including WRAP53, KDM6B, and TRIM29 ([153]Figures 1D and 1F,
   [154]Tables 1 and [155]2, and [156]Figure S3). These data highlight the
   potential impact of stratifying the isoform data by TP53 mutational
   status, as multiple isoform changes in TP53-related genes were
   identified that otherwise would go undetected.

Isoform switching events unique or shared based on stratification by
pathology and TP53 status

   Although a large overlap was observed in GO enrichment analysis in
   BE.LGD versus BE.HGD (n = 75) compared with BE.LGD TP53 WT versus
   BE.HGD TP53 MUT (n = 135), the number of isoform switching events
   nearly doubled when TP53 stratification was incorporated. Thus, a Venn
   diagram was generated to identify shared and unique isoform-switched
   genes in the analysis as shown in [157]Figure 3A. Among all significant
   genes in the analysis, 60 isoform-switched genes are shared between the
   two comparisons ([158]Figures 3A and [159]Table S6), with 15 unique
   genes identified in BE.LGD versus BE.HGD ([160]Figures 3A and
   [161]Table S7), and 75 genes uniquely identified in BE.LGD TP53 WT
   versus BE.HGD TP53 MUT ([162]Figures 3A and [163]Table S8). As shown in
   [164]Figure 3A, 80% (60 of 75) of the isoform-switched genes in BE.LGD
   versus BE.HGD are shared, whereas only 44.4% (60 of 135) of the
   isoform-switched genes in BE.LGD TP53 WT versus BE.HGD TP53 MUT are
   shared, supporting that greater numbers of unique alterations are
   associated with mutant TP53. Top significant genes in each section of
   the Venn diagram include many protein-coding genes, but also several
   non-protein-coding genes including SNHG18, [165]AC009509, CADM3-AS1,
   MIR4458HG, CENPBD1P1, and LINC02615 ([166]Figures 3B–3D). Previously
   reported genes involved in isoform switching (RNF128, TPM4, and KRAS)
   were shared between the two comparisons, with the exception of MINDY1,
   which was uniquely identified in BE.LGD TP53 WT versus BE.HGD TP53
   MUT.[167]^27^,[168]32, [169]33, [170]34, [171]35, [172]36 Next, GO
   enrichment analysis of the non-overlapping isoforms in BE.LGD versus
   BE.HGD and BE.LGD TP53 WT versus BE.HGD TP53 MUT comparisons were
   conducted. For genes uniquely identified in BE.LGD versus BE.HGD, top
   enriched processes are related to histone and protein demethylation,
   dealkylation, and GDP catabolic process ([173]Figure 3E and
   [174]Table S9). Top enriched processes based on genes unique in BE.LGD
   TP53 WT versus BE.HGD TP53 MUT are more strongly associated with lipid
   and alcohol metabolism, retinoic acid metabolism and biosynthesis, and
   midgut development ([175]Figure 3E). For genes shared between the two
   comparisons in the overlapping section of the Venn diagram, GO enriched
   processes are similar to the previous analysis ([176]Figure 1B) with
   top processes related to glucuronidation and metabolism
   ([177]Figures 3E and [178]Table S10). Although enriched processes had a
   significant p-value (≤0.05), none of the enriched processes unique to
   BE.LGD TP53 WT versus BE.HGD TP53 MUT reached a significant false
   discovery rate (FDR) ≤0.05, but instead ranged from 0.08 to 0.16,
   suggesting greater heterogeneity of isoform-enriched genes in that
   grouping.

Figure 3.

   [179]Figure 3
   [180]Open in a new tab

   GO processes associated with significant isoform switching based on
   histopathology and TP53 mutation status

   (A) Venn diagram showing 15 unique genes involved in isoform switching
   in BE.LGD versus BE.HGD, 75 unique genes involved in isoform switching
   in BE.LGD TP53 WT versus BE.HGD TP53 MUT, and 60 genes shared. (B–D)
   Top 10 most significant isoform-switched genes in each section of the
   Venn diagram (|dIF| ≥ 0.1 and FDR ≤ 0.05). (E) GO enrichment analyses
   showing the top 10 enriched biological processes in each part of the
   Venn diagram (p < 0.01).

TP53-linked isoform changes correspond to patient outcomes

   TP53 is the strongest known genetic driver of EAC with a mutation
   frequency of over 70% in EAC patients.[181]^8^,[182]^13 In turn,
   isoform-switched genes with direct interaction with TP53 were
   identified using the STRING database ([183]Figure S3). Patient survival
   was also determined to investigate whether there is an association
   between expression of specific isoforms of TP53-linked genes and
   survival among EAC patients.

   Isoform switching of RNF128 was recently reported to show decreased
   expression of RNF128-202 and increased expression of RNF128-201
   contributing to the stabilization of mutant TP53 and potentially EAC
   progression.[184]^27 Our analysis confirms the RNF128 isoform-switching
   event with similar trends observed ([185]Figure 4A). The overall
   expression of RNF128 is significantly lower in patients with BE.HGD
   compared with patients with BE.LGD and isoform usage of each isoform is
   significantly altered, supporting the isoform switching event.
   Functional domain prediction using Pfam did not reveal any major
   changes in protein domains. However, a significant difference in
   patient survival was observed in the analysis ([186]Figure 4B).
   Patients with higher expression levels of RNF128-202 had a
   significantly higher survival probability compared with patients with
   lower expression of the same isoform (Wilcoxon p-value = 0.031).
   Patient samples with high expression of RNF128-202 were largely
   composed of BE.LGD samples with only one BE.HGD and one EAC patient
   sample, whereas tissues with low expression levels of RNF128-202 were
   contributed by either BE.HGD or EAC samples.

Figure 4.

   [187]Figure 4
   [188]Open in a new tab

   TP53-linked isoform switches predict patient survival

   (A) Isoform switching of RNF128. Transcripts involved in isoform
   switching of RNF128 and their predicted functional domain Error bars
   indicate 95% confidence intervals and significance levels of isoform
   usage were determined by the Mann–Whitney U test (∗ FDR ≤ 0.05). (B)
   Kaplan-Meier survival curve of patients stratified by expression level
   of RNF128 transcript ENST00000324342 and the histology proportion of
   patients in each group. (C) Isoform switching of WRAP53. Transcripts of
   WRAP53 and their corresponding isoform usage. Error bars indicate 95%
   confidence intervals and significance levels of isoform usage were
   determined by the Mann–Whitney U test, ns indicates no significant
   change, (∗ FDR ≤ 0.05). (D) Kaplan-Meier survival plot of patients
   stratified by expression of WRAP53 transcript ENST00000316024 and
   histology distribution of patients in each group. (E) Survival
   probability of BE.LGD or BE.HGD patient samples (n = 32). (F) Patient
   survival probability of BE.LGD or BE.HGD + EAC (n = 41). Only the most
   severe pathological sample for each patient was kept in the analysis if
   a patient contributed multiple biopsy samples. Significance of survival
   differences was calculated by the Gehan-Breslow-Wilcoxon test.

   WRAP53 is another isoform-switched gene with direct TP53 interactions
   ([189]Figure S3). WRAP53 transcripts with exon regions overlapping with
   the TP53 exon are believed to be the antisense of TP53, which can
   result in TP53 induction.[190]^42 On the gene level, significant
   expression differences were not observed in WRAP53 between BE.LGD and
   BE.HGD patient samples. However, on the isoform level, a significant
   increase of WRAP53-201 and a decrease of WRAP53-202 are observed
   ([191]Figure 4C). The 5′-UTR of WRAP53-201 contains a shared region
   with the first exon of TP53, whereas WRAP53-202 does not have any
   overlapping regions. Moreover, patients with high expression levels of
   WRAP53-201 have a significantly lower survival probability compared
   with patients with low expression levels of WRAP53-201 ([192]Figure 4D;
   Wilcoxon p-value = 0.050). All high-expressing patient samples were
   either BE.HGD or EAC, whereas low-expressing samples were contributed
   by BE.LGD tissues. Still, two BE.HGD patient samples and one EAC sample
   had low WRAP53-201 expressions, suggesting isoform expression level is
   not solely determined by histopathological progression.

   In addition, survival differences based on histopathological
   categorization were calculated to evaluate the possibility that
   significant differences in patient survival were solely driven by the
   histopathological distribution. The survival probability of patients
   with BE.HGD is marginally lower than the survival probability of BE.LGD
   patients ([193]Figure 4E; Wilcoxon p-value = 0.113) and when EAC
   samples are included in the analysis, the survival probability is also
   marginally lower and borderline statistically significant
   ([194]Figure 4F; Wilcoxon p-value = 0.058). Ultimately, these results
   support that isoform switching events may reveal specific gene isoforms
   that impact patient survival to a greater extent than histopathological
   changes alone and specific isoforms may be targetable as novel targets
   for prevention or treatment.

Isoform switching informs targetable processes for EAC inhibition

   Our results revealed that mutant TP53 increases isoform switching
   events and modulates numerous lipid and retinoic acid-linked processes
   including unique changes in 9-cis-retinoic acid, leading us to evaluate
   the inhibitory potential of two Food and Drug Administration
   (FDA)-approved and mechanistically relevant pharmacological agents for
   targeting EAC. TP53 MUT EAC cells (OE33 and JHAD1) were treated with
   adapalene,[195]^43 an agonist of the retinoic acid receptor (RAR) β and
   γ, and the mutant TP53 targeting agent APR-246.[196]^44 These two drugs
   were selected as proof-of-concept agents given identification of
   9-cis-retinoic acid metabolic processes and retinoic acid biosynthetic
   processes in the comparison of BE.LGD TP53 WT versus BE.HGD TP53 MUT
   ([197]Figure 3E). Gene-level GSEA reveals significant downregulation of
   retinol metabolism in BE.HGD TP53 MUT patients compared with BE.LGD
   TP53 WT ([198]Figure S4), suggesting that activation of retinoic acid
   signaling may elicit anti-cancer activity targeting EAC. Therefore,
   adapalene was used to test whether stimulating retinoic acid signaling
   in TP53 MUT EAC cell lines would induce EAC cell death. Both OE33 and
   JHAD1 cells show dose-dependent inhibition of viability following
   adapalene treatment with the lethal dose 50 (LD[50]) in the 1.5- to
   2.0-μM range ([199]Figures 5A–5D). At 72 h post-treatment, 1.5 μM
   adapalene significantly decreases OE33 cell viability to 42.59%
   ([200]Figures 5A and 5B). Similarly, at 72 h, 2.0 μM adapalene
   treatment significantly decreases JHAD1 cell viability to 41.50%
   ([201]Figures 5C and 5D). Importantly, OE33 cells treated with
   adapalene do not recover following treatment removal ([202]Figures S5A
   and S5B). However, JHAD1 cells do show recovery upon adapalene removal,
   suggesting adapalene is only inhibitory in JHAD1 cells when present,
   with little durable effect once removed ([203]Figures S5C and S5D).
   Thus, a combination of agents may prove more efficacious in JHAD1
   cells.

Figure 5.

   [204]Figure 5
   [205]Open in a new tab

   Enrichment analysis of isoform switching events supports retinoic acid
   and TP53 signaling pathway targeting to inhibit EAC growth

   (A) Viability of adapalene-treated OE33 cells measured using Calcein-AM
   staining 48 and 72 h post-treatment with (B) representative fluorescent
   images of adapalene-treated Calcein-AM stained OE33 cells. (C)
   Viability of adapalene-treated JHAD1 cells assessed using Calcein-AM 48
   and 72 h post-treatment with (D) representative fluorescent images of
   adapalene-treated stained JHAD1 cells. (E) Viability of OE33 cells
   treated with APR-246 at 48 and 72 h post-treatment with
   (F) representative fluorescent images of stained OE33 cells treated
   with APR-246. (G) Viability of JHAD1 cells treated with APR-246 at 48
   and 72 h post-treatment with (H) representative fluorescent images of
   stained JHAD1 cells treated with APR-246. Significant differences of
   viability were assessed by ANOVA with Tukey’s post hoc test for
   multiple comparisons between treatments. Within each time point,
   treatments were significantly different from a = vehicle (VEH), b =
   1.0 μM adapalene or 10 μM APR-246, c = 1.5 μM adapalene or 20 μM
   APR-246, and d = 2.0 μM adapalene or 40 μM APR-246. Data are shown as
   mean ± standard error. Note: VEH-48h and VEH-72h images are the same in
   5B and 5F because adapalene and APR-246 treatment of OE33 cells
   occurred in the same plate to permit direct comparisons. Similarly,
   VEH-48h and VEH-72h images are the same in 5D and 5H for JHAD1 cells.
   Scale bars, 500 μm.

   Next, APR-246, a small molecule agent that restores TP53 WT functions
   in TP53 MUT cells by covalently binding to cysteine residues in TP53
   MUT cells,[206]^44 was evaluated. TP53 mutation is a well-known driver
   for EAC progression and several isoform-switched genes have direct
   interactions with TP53 ([207]Figure S3). OE33 and JHAD1 cells were
   treated with various concentrations of APR-246. APR-246 treated OE33
   cells show a dose-dependent response with an LD[50] above 20 μM, but
   clearly less than 40 μM, which permanently induces cell death
   ([208]Figures 5E, 5F, [209]S5E, and S5F). At 72 h, 20 μM APR-246
   treatment significantly decreased OE33 cell viability to 62.84% and
   40 μM APR-246 treatment significantly decreased OE33 cells viability to
   8.57%. OE33 cells treated with 40 μM APR-246 do not recover upon
   treatment withdrawal ([210]Figures S5E and S5F). For JHAD1 cells, a
   sharp decrease in cell viability was observed between 20 μM and 40 μM
   treatments. At 72 h, 20 μM APR-246 treatment significantly decreases
   cell viability to 64.61%, whereas 40 μM treatment completely eradicates
   JHAD1 cells ([211]Figures 5G and 5H). JHAD1 cells did recover following
   drug removal when lower concentrations of APR-246 were utilized
   (<40 μM; [212]Figures S5G and S5H). Representative fluorescent images
   from the viability assays performed in OE33 and JHAD1 cells show
   decreases in viability from adapalene ([213]Figures 5B and 5D) and
   APR-246 ([214]Figures 5F and 5H). Collectively, these results confirm
   that targeting the retinoic acid and TP53 pathways using
   pharmacological agents is successful in inducing EAC cancer cell death,
   and that evaluations of mechanistically driven combinations are
   warranted to screen for enhanced durable effects. These data further
   validate the identification of pathways based on isoform-switching
   analysis using RNA-seq datasets.

Discussion

   In this study, we performed isoform switch analyses using an RNA-seq
   dataset composed of 57 esophageal tissues across a continuum of
   pathologies ranging from BE with metaplasia, to BE with dysplasia (both
   low and high grade), and EAC. Isoform switch analyses were conducted
   based on histopathological progression and TP53 mutation status given
   its key role in EAC progression. Our goal was to characterize isoform
   switching events in EAC and understand whether factors associated with
   increased risk for progression to EAC, namely the presence of HGD or
   mutant TP53, show unique isoform switching events that may inform
   future risk stratification, prevention, or therapeutic efforts.

   Seventy-five genes involved in isoform switching were identified
   comparing low-risk with high-risk histopathologies or BE.LGD with
   BE.HGD groups. Inclusion of TP53 mutation status coupled with
   histopathological stratification markedly increased the number of
   isoform-switched genes identified, resulting in 135 genes. Previously
   published studies suggest that alternative splicing is performed by
   spliceosomes and could also be controlled by RBPs, both of which are
   crucial components of various ribonucleoprotein (RNP)
   complexes.[215]^40^,[216]^41 In our data, gene-level GSEA analysis
   revealed many significantly enriched processes related to alternative
   splicing in BE.HGD compared with BE.LGD and in BE.HGD TP53 MUT compared
   with BE.LGD TP53 WT patient samples, such as ribonucleoprotein complex
   biogenesis, spliceosomal snRNP assembly, mRNA splice site selection,
   and spliceosomal complex assembly as evidenced by strong upregulation
   of genes in these pathways in BE.HGD compared with BE.LGD patient
   samples ([217]Figure 2). For example, small nuclear ribonucleoprotein
   D1/2 polypeptide (SNRPD1 and SNRPD2) were two of the leading-edge genes
   identified in both spliceosomal snRNP assembly and spliceosomal complex
   assembly pathways. Previously published studies suggest that
   overexpression of SNRPD1 and SNRPD2 are common in many cancers, leading
   to abnormal alternative splicing and regulating cell cycle and
   autophagy.[218]^45^,[219]^46 In EAC, TP53 mutation is the most commonly
   observed mutation, with a mutation frequency over 70%.[220]^13 TP53
   mutation is not only a driver mutation contributing to the EAC
   progression, but TP53 mutations also contribute to alternative splicing
   activity in cancer.[221]^47 In pancreatic ductal adenocarcinoma (PDAC),
   TP53 mutations lead to significant exon usage changes and upregulate
   expression of RNA processing factors important for RNA
   splicing.[222]^47 Therefore, aberrant RNA splicing machinery activity
   resulting from TP53 mutation may explain why more isoform-switched
   genes are observed in the comparison of BE.LGD versus BE.HGD with the
   incorporation of TP53 mutation status. Aside from significant pathway
   enrichment of alternative-splicing-related processes, multiple RBPs are
   differentially expressed on the gene level. The number of significantly
   altered RBPs identified are markedly increased when TP53 mutation
   status is incorporated in the analysis and several RBPs have protein
   interactions with identified isoform-switched genes ([223]Figures 2C
   and 2D). Our data support the merit of delineating specific
   interactions between RBPs and isoform-switched genes in future research
   samples using cross-linking-immunoprecipitation-sequencing. Global
   functional consequence analysis revealed three shared significant
   changes between BE.LGD versus BE.HGD and BE.LGD TP53 WT versus BE.HGD
   TP53 MUT, which are the loss of ORFs, signal peptide, and coding
   transcript, suggesting that dominant isoforms observed in BE.HGD,
   regardless of TP53 mutation status, will have profound impact on
   related cellular functions, which might also be contributing factors
   for EAC progression.

   Comprehensive analysis of isoform switching events has not previously
   been performed in BE or EAC. A previous pan-cancer analysis utilizing
   the Pan-Cancer Analysis of Whole Genomes included seven EAC cases
   reporting isoform-related changes in only one EAC patient when compared
   with non-patient-matched normal esophageal tissues from Genotype-Tissue
   Expression (GTEx), presenting findings as “most dominate transcripts,”
   not isoform switching events making generalizability uncertain.[224]^48
   In addition, isoform switching of RNF128 was previously reported in HGD
   and EAC.[225]^27 RNF128 encodes an E3 ubiquitin ligase responsible for
   the degradation of TP53. Similarly, we show that RNF128 undergoes
   isoform switching in patients with BE.LGD compared with BE.HGD, with
   increased expression of RNF128-201 and decreased expression of
   RNF128-202. RNF128-201 has limited ubiquitin ligase activity and fails
   to degrade mutant TP53, whereas RNF128-202 is the main isoform
   responsible for TP53 degradation.[226]^27 Importantly, we observed a
   significant patient survival advantage among RNF128-202 high expressors
   compared with low expressors. In addition, patients with high
   expression levels of RNF128-201 had a non-significantly lower survival
   probability compared with patients with low expression levels (data not
   shown).

   With the exception of RNF128, isoform switching events have been
   largely unexplored in BE and EAC. Herein, we identify large numbers of
   events linked to advanced pathology, as well as TP53 mutation (as
   detailed in [227]Tables S6, [228]S7, and [229]S8). Two reports
   characterizing isoform switching events in ESCC identified numerous
   unique isoforms with only MINDY1 as a gene in common with our isoform
   switch analysis, supporting divergent molecular changes based on
   esophageal cancer subtype.[230]^36^,[231]^49 MINDY-1, also known as
   FAM63A, is a newly identified deubiquitinating protein that
   preferentially removes K48-linked ubiquitin molecules.[232]^50 In BE
   and EAC, MINDY1 was identified in a set of 90 genes significantly
   predicting disease progression by distinguishing EAC progressors from
   patients with non-dysplastic BE.[233]^51 Moreover, MINDY-1 was found to
   promote bladder cancer progression by stabilizing YAP,[234]^52 a key
   tumorigenesis pathway member that can also be induced by conjugated
   bile acids in EAC.[235]^53 In our analysis, increased usage of the
   MINDY1 transcript that lacks deubiquitinase (DUB) domains and a
   decreased usage of the MINDY1 transcript that contains five DUB domains
   are observed, potentially suggesting impaired MINDY-1 function during
   EAC progression. Despite the fact that the MINDY1 transcripts involved
   in isoform switching are different in our analysis than the data
   published in ESCC,[236]^36 a similar switching trend is observed in
   that the isoform transcript that does not contain DUB domains is
   increased while the usage is decreased for a DUB domain-containing
   transcript. These data suggest that MINDY1 isoforms lacking
   deubiquitinase domains may be important for both EAC and ESCC
   progression.

   In ESCC, reported isoform switching events were dominated by cell
   regulation of cell motility, cell-to-cell junction organization,
   regulation of cell migration, and adhesion-linked GO processes.[237]^49
   Select isoform-switched genes in our BE progression dataset hold some
   similar functions. For example, isoform switching of TPM4 was
   originally observed in breast cancer with loss of TPM4.1 associated
   with increased migration, disruption of cell-cell adhesions, and cancer
   invasiveness.[238]^32 Moreover, TPM4.1 inhibited cellular invasion by
   regulating the Rac1-myosin IIB signaling.[239]^32 Significant decreases
   in TPM4.1 are observed when comparing BE.LGD with BE.HGD patients, and
   thus loss of TPM4.1 during EAC progression may contribute to invasive
   cell behavior.[240]^32

   Amplification and mutation of KRAS have been previously reported in
   EAC.[241]^30^,[242]^54^,[243]^55 Despite the mutation frequency of KRAS
   in EAC being less than 20%, it is considered an EAC driver and proposed
   targeting of KRAS mutations may sensitize a subgroup of EAC patients to
   targeted treatment.[244]^30 Alternative splicing of KRAS is well
   documented in other cancers,[245]^33^,[246]^35^,[247]56, [248]57,
   [249]58, [250]59, [251]60 resulting in two isoforms, KRAS4A and
   KRAS4B,[252]^33 both of which can carry KRAS mutations.[253]^33
   Although expression of both KRAS4A and KRAS4B has been observed in
   cancer,[254]^58 KRAS4B is thought to be more tumorigenic than
   KRAS4A.[255]^33 KRAS4A expression inhibits apoptosis, whereas KRAS4B
   does not.[256]^57 KRAS4B also regulates the matrix metalloproteinase
   MMP-2, which is linked to esophageal cancer development, via the
   PI3K-AKT signaling pathway.[257]^56^,[258]^61 Moreover, in colorectal
   carcinoma, expression of KRAS4B was found to be significantly
   associated with larger tumor size while expression of KRAS4A was
   significantly associated with better patient survival.[259]^35 In our
   isoform switching analysis and qRT-PCR validation, we observed
   increased expression of KRAS4B and decreased expression of KRAS4A,
   suggesting that KRAS4B may exert a similar oncogenic role in EAC
   progression and warrant investigation as a therapeutic target
   ([260]Figures S2C and S2D).

   Interestingly, isoform switch analysis in this study also identified
   significant isoform-switched genes without significant differences in
   gene-level expression ([261]Tables 1 and [262]2). Thirty-two isoform
   switching events (32 of 75, or 42.7%) were identified in comparing
   BE.LGD with BE.HGD, whereas BE.LGD TP53 WT compared with BE.HGD TP53
   MUT resulted in identification of 62 isoform switching events (62 of
   135, or 45.9%) without overall gene expression changes. GO enrichment
   analysis revealed similar enriched biological processes for those genes
   compared with enrichment analysis of all isoform-switched genes (data
   not shown). Isoform-switched genes that occur in the absence of altered
   gene-level expression include KDM6B, UGT2B7, WRAP53, and TRIM29. KDM6B,
   also known as lysine-specific demethylase 6B, is a histone demethylase
   that can act as either a tumor suppressor or oncogene in
   cancer.[263]^62 KDM6B is reportedly overexpressed in pancreatic
   premalignancy with its expression decreasing with progression to PDAC,
   supporting a role for KDM6B in pancreatic carcinogenesis.[264]^63
   Although overall expression differences were not significant between
   patients with BE.LGD compared with BE.HGD, expression of the KDM6B
   coding transcript decreased while the expression of the non-coding
   transcript increased, suggesting KDM6B might also be important during
   early EAC progression. Conversely, KDM6B has been shown to promote ESCC
   progression and KDM6B levels significantly increased in patients with
   lymph node metastasis.[265]^64

   As highlighted in [266]Figure 3, the Venn diagram depicts all
   significant isoform-switched genes and reveals common and unique
   isoform switching events based on histopathology alone or both
   histopathology and TP53 mutation status. Subsequent enrichment analysis
   revealed common and unique changes in GO biologic processes based on
   stratification by histopathology and TP53 mutation status. Multiple
   biologic processes related to glucuronidation were among the top
   significantly enriched GO processes shared between two comparisons,
   primarily caused by significant isoform switching events observed in
   UGT1A1 and UGT2B7. In recent years, multiple studies have shown that
   UDP-glucuronosyltransferase enzymes (UGTs) are linked to increased
   cancer risk, drug resistance, and cancer progression.[267]^65^,[268]^66
   Multiple case-control studies have been published identifying
   associations between UGT polymorphisms and cancer
   risk.[269]^65^,[270]^66 It is hypothesized that UGT polymorphisms can
   decrease glucuronidation of carcinogens and molecules that promote
   cancer, such as estrogens, dietary carcinogens, and tobacco
   carcinogens, resulting in cancer progression.[271]^65^,[272]67,
   [273]68, [274]69, [275]70, [276]71, [277]72 However, studies have not
   found a significant association between UGT polymorphisms and EAC
   risk.[278]^73 Based on our analysis, we identified a significant
   increase in the usage of a UGT1A1 isoform that lacks a UDP-glucuronosyl
   and UDP-glucosyl transferase (UDPGT) domain, suggesting there is
   decreased glucuronidation activity of UGT1A1. Moreover, isoform
   switching results of UGT2B7 showed a significant decrease in expression
   of the isoform that encodes a UDPGT domain, which also suggests
   impairment of UGT2B7 glucuronidation activity in EAC progression.

   In addition, enrichment of processes related to demethylation from
   isoform-switched genes uniquely observed in the comparison of BE.LGD
   versus BE.HGD were detected. Demethylation processes in EAC have been
   extensively studied, and previous studies identified a key role for
   demethylation during disease progression, including promoter
   demethylation for chemokine genes and upregulation of histone
   demethylation proteins (KDM4C and KDM6A).[279]^74^,[280]^75 Isoform
   switch analysis and GO biological process enrichment analysis
   identified KDM5A and KDM5B as the main genes contributing to the
   observed enrichment of demethylation process, both of which have
   oncogenic roles in ESCC development but have not been assessed in
   EAC.[281]^64^,[282]^76 For GO processes uniquely enriched in BE.LGD
   TP53 WT versus BE.HGD TP53 MUT, processes related to retinoic acid
   biosynthesis and 9-cis-retinoic acid biosynthesis were identified
   ([283]Figure 3E). Gene-level analysis revealed downregulation of
   associated genes with GSEA showing downregulation of retinol metabolism
   in patients with BE.HGD TP53 MUT ([284]Figure S4). Isoform switching
   analysis of retinol dehydrogenase 5 (RDH5) and aldehyde dehydrogenase 1
   family member A2 (ALDH1A2), key proteins in the 9-cis-retinoic acid
   biosynthesis, revealed downregulation of both RDH5 and ALDH1A2 coding
   transcripts, while expression levels of non-coding transcripts were
   increased, further suggesting decreased activity of such processes
   ([285]Table S8). Downregulation of RDH5 has been reported in multiple
   types of cancer, including hepatocellular carcinoma, colon adenomas and
   carcinomas, and thyroid carcinoma.[286]77, [287]78, [288]79 In
   hepatocellular carcinoma, RDH5 was found to suppress proliferation and
   metastasis by reversing epithelial-mesenchymal transition.[289]^79
   Downregulation of ALDH1A2 was also reported in several cancer types and
   suggest to act as a tumor suppressor.[290]^80^,[291]^81 Studies have
   also pointed out that cancer cells treated with 9-cis-retinoic acid
   induce apoptosis and inhibit cancer cell growth both in vitro and
   in vivo.[292]^82^,[293]^83 Moreover, retinoic acid receptor RARβ, was
   shown to act as a tumor suppressor in lung cancer.[294]^84 In BE
   patients with LGD or HGD, expression of both RARβ and RARγ are
   significantly lower compared with normal esophageal tissues, and the
   retinoic acid-regulated nuclear protein regulation pathway is the most
   significantly enriched pathway in EAC progressors compared with
   patients with non-dysplastic BE, suggesting the potential role of
   retinoic acid-related processes in EAC progression.[295]^51^,[296]^85

   Significant isoform-switched genes with direct interaction with TP53
   were also identified in this study, including RNF128, KDM6B, and WRAP53
   ([297]Figure S3). WRAP53 has three separate alternative start exons
   (1α, 1β, and 1γ), which generate three different transcripts of WRAP53
   (α, β, and γ).[298]^42 The precise function of WRAP53 is largely
   unknown; however, it was previously reported that the 5′UTR region of
   WRAP53α, located in exon 1α, overlaps with the first exon of TP53 and
   is a natural antisense of TP53, inducing expression of both WT and MUT
   TP53.[299]^42 The other two forms of WRAP53 (β and γ) do not have exon
   regions that overlap with exons of TP53.[300]^42 In the isoform
   switching analysis, we observed increased expression of a WRAP53
   transcript with an overlapping exon region with TP53 and decreased
   expression of a WRAP53 transcript that does not contain an overlapping
   region. Moreover, protein functional domains predicted by Pfam are the
   same between these two isoforms. Based on this, we speculate that
   increased expression of a transcript with an overlapping region of TP53
   may induce expression of mutant TP53 contributing to disease
   progression; however, such switching events are not likely to have a
   profound impact on WRAP53 function since the functional domains remain
   unchanged.

   The Venn diagram and subsequent GO biological process analysis in
   [301]Figure 3 show differences in pathway enrichment based on pathology
   and TP53 mutation status, suggesting that patients stratified into
   different subgroups may potentially benefit from different treatment
   regimens. As a proof-of-principle, we used two FDA-approved drugs that
   target the 9-cis-retinoic acid pathway and TP53 mutation to test this
   hypothesis. Two EAC cell lines, OE33 and JHAD1, were treated with the
   RARβ and RARγ agonist, adapalene, and results indicated that adapalene
   was a potent inducer of cell death in both cell lines. Considering RARβ
   and RARγ are two of the receptors of 9-cis-retinoic acid, these data
   suggest activation of 9-cis-retinoic acid pathway may hold potential
   for targeting EAC.[302]^86 Moreover, to test the potential benefit of
   blocking mutant TP53 expression in EAC, OE33 and JHAD1 cells (both TP53
   mutant) were treated with APR-246, a small molecule that restores TP53
   WT functionality.[303]^44 Results showed that post-treatment with
   APR-246, both EAC cell lines experienced significant loss of viability,
   which is in alignment with previously published results in EAC and
   ESCC, although a discrepancy in effective treatment concentration of
   OE33 cell was noticed.[304]^87^,[305]^88 Therefore, treatment with
   APR-246 suggests that inhibition of mutant TP53 expression and isoform
   expression linked to mutant TP53 may be a viable preventive or
   treatment approach for inhibiting EAC. Moreover, enrichment of the
   lipid metabolism process is also observed in our analysis for BE.HGD
   patients that carry the TP53 mutation. Metformin, a lipid metabolism
   modulating drug, was also shown to inhibit proliferation of EAC cell
   lines,[306]^89 further supporting the plausibility of targeting
   specific pathways identified in the analysis and the importance of
   patient stratification based on TP53 mutational status. Finally,
   epidemiological and limited preclinical data point to statins as
   esophageal cancer inhibitors, in alignment with agents impacting lipid
   metabolism holding promise in a subset of EACs.[307]90, [308]91,
   [309]92, [310]93, [311]94

   Treatment of EAC cells with adapalene and APR-246 shows promising
   results inhibiting EAC cell viability in vitro, suggesting that
   targeting specific pathways or isoforms might offer novel therapeutic
   insights. Genetic manipulation of specific isoforms is required to
   further explore the possibility of inhibiting EAC cell growth by
   modulating isoform expression levels. Development of synthetic
   splice-switching oligonucleotides may also offer a future therapeutic
   direction as specific isoforms are better characterized and
   isoform-specific inhibition showed potent results as targeted cancer
   therapies.[312]^95^,[313]^96 Finally, survival analysis of patients
   stratified by expression levels of RNF128 and WRAP53 isoforms showed
   that expression of specific isoforms could be used as prognostic
   markers to better predict patient survival probability. Taken together,
   isoform switching analysis may provide novel insights for the
   identification of prognostic markers and inform new potential
   therapeutic targets for EAC.

Materials and methods

RNA-seq and analysis of esophageal patient tissues

   RNA was isolated from 57 esophageal tissues derived from 46 patients
   undergoing esophagectomy following a diagnosis of EAC or HGD at the
   University of Michigan as previously described.[314]^27^,[315]^97
   Signed informed consent was obtained from each patient and all
   procedures were consistent with the protocol submitted and approved by
   the institutional review board. Tissue samples were collected from the
   tumor or within a 6 cm of the surrounding area. Histopathology of
   patient samples was characterized and confirmed by two independent
   pathologists. Patient tissues were categorized based on histopathology
   into the following groups: BE with metaplasia (n = 6) or LGD (n = 19)
   combined and referred to as BE.LGD (n = 25), BE with HGD based on >40%
   HGD and referred to as BE.HGD (n = 21), and EAC (n = 11) prior to RNA
   isolation. Total RNA was extracted using standard TRIzol methodology
   and all isolates had RNA integrity number greater than 7.0.
   Strand-specific RNA-seq library preparation and sequencing were
   performed at the University of Michigan Advanced Genomics Core, using
   the Illumina HiSeq 100-base pair paired-end sequencing platform.
   Quality control and adapter trimming were performed using Qiagen CLC
   Genomics Workbench 21.0 ([316]https://digitalinsights.qiagen.com/).
   After processing, the average read count per sample was 72 million
   (range: 40–149 million, median: 63 million). RNA quantification and
   differential gene expression analysis were performed using Kallisto
   (version 0.46.2, 100 bootstraps with pseudoalignment saved) and Sleuth
   (version 0.30.0).[317]^98^,[318]^99 On average, 83% to 90% of read
   fragments were pseudoaligned per sample. Gencode Human Release 34
   (GRCh38.p13) was used as the gene and transcript annotation source in
   the analysis, which contains 59,667 genes and 228,116 gene transcripts.
   Likelihood-ratio tests were used to determine significantly differently
   expressed genes by comparison groups and the
   Benjamini-Hochberg-adjusted FDR was applied to the analysis. Genes with
   p-value and FDR less than 0.05 were considered significant. Sequence
   files are deposited to the NCBI Gene Expression Omnibus
   (GEO:[319]GSE193946) using the protocol approved by the institutional
   review board. TP53 mutation analysis on patient samples was performed
   using the PlexSeq process by PlexSeq Diagnostics (Cleveland, OH). The
   lollipop plot of TP53 mutation was generated using MutationMapper in
   cBioPortal.[320]^100^,[321]^101 Differentially expressed RBPs were
   identified using the result from Sleuth and the list of known human
   RBPs was acquired from a previously published study.[322]^102

Identification of isoform switching events

   Identification of isoform switching events was performed in R (version
   4.1.1) using the IsoformSwitchAnalyzeR package (version
   1.14.1).[323]^18^,[324]^103 Kallisto was used to quantify transcript
   abundance, which was imported to the package, followed by ORF
   prediction and isoform switch testing using an embedded function using
   the DEXSeq package.[325]104, [326]105, [327]106, [328]107 Isoform
   switching testing was performed by calculating the isoform fraction
   (IF) defined as the fraction of gene expression originating from each
   associated isoform, and the differential isoform fraction (dIF) value
   was subsequently determined by calculating the difference of IF values
   between two conditions of interest. Several external analyses were
   performed to predict functional consequences associated with isoform
   switching, including prediction of coding potential (CPAT),[329]^108
   coding domains (Pfam),[330]^109 signal peptides (SignalP),[331]^110
   intrinsically disordered regions (IUPred2A),[332]^111 and sensitivity
   to nonsense-mediated decay.[333]^18^,[334]^105^,[335]^112^,[336]^113
   Results from external analyses were combined with isoform switching
   testing to identify isoform switch events with predicted functional
   consequences.[337]^18 Results were filtered using an absolute dIF
   cutoff at 0.1 and FDR cutoff at 0.05. Alternative splicing analysis and
   genome-wide enrichment analysis were also
   performed.[338]^18^,[339]^103^,[340]^105 Volcano plots were generated
   using an embedded function in the package and the dot plots were
   generated using ggplot2 (version 3.3.5) in R.[341]^18^,[342]^114 The
   Venn diagram was generated using InteractiVenn.[343]^115 Isoform
   switching plots were generated directly by the IsoformSwitchAnalyzeR
   package. ENSEMBL transcript IDs on the isoform switching plots were
   matched to the corresponding transcript names. The overlapping region
   between WRAP53 isoforms and TP53 was identified using NCBI
   BLAST.[344]^116

Gene and isoform functional analysis

   Comprehensive expression analysis was performed using the list
   of genes/transcripts that showed significant isoform switching
   events. The gene names of statistically significant (|dIF| ≥ 0.1 and
   FDR ≤0.05) isoform switches for comparison groups of interest were
   uploaded and analyzed in Metacore and Cortellis Solution software
   ([345]https://clarivate.com/products/metacore/, Clarivate Analytics,
   London, UK). The Enrichment Analysis, Metabolic Network, GO Molecular
   Function, and GO Localization tools were used to identify enriched
   pathways, networks, and processes that were altered in each comparison
   group. Significance was calculated by using the hypergeometric test and
   the default Metacore database with Homo sapiens as the background for
   each analysis. All p-value and FDR significant Pathway Maps, Process
   Networks, Diseases, GO Processes, Metabolic Networks, GO Molecular
   Functions, and GO Localizations for each comparison were exported for
   further analysis. Protein interactions between identified genes and
   TP53 were predicted using the STRING database (version 11.5) with
   minimum required interaction score set to medium confidence (0.400) and
   max number of interactors to show set to query proteins only.[346]^117
   TP53 was added to the list of proteins identified from isoform-switched
   genes and used as input for STRING protein interaction prediction.

Gene-level GSEA

   Gene-level GSEA analysis was performed as previously
   described.[347]^118 In short, ranking scores of each gene were
   calculated using the results of differential gene expression analysis
   and saved as a ranked list file. The file was imported to GSEA software
   (version 4.2.2; Broad Institute, Cambridge, MA) and analysis was
   performed using GO gene sets in the Molecular Signatures Database with
   1,000 permutations, max and min size set to 500 and 15,
   respectively.[348]119, [349]120, [350]121, [351]122, [352]123, [353]124
   Leading-edge genes of each identified pathway, which were genes that
   contributed the most to the enrichment, were extracted from GSEA
   results. The heatmap and bar plots were generated using ComplexHeatmap
   (version 2.12.0) and ggplot2, respectively in R.[354]^114^,[355]^125
   The heatmap was plotted using z-score-transformed transcripts per
   million (TPM) values of leading-edge genes.

Isoform expression determination via quantitative reverse-transcriptase PCR

   cDNA for the qRT-PCR experiment was generated using the High-Capacity
   cDNA Reverse Transcription Kit (ThermoFisher Scientific, Waltham, MA),
   following the manufacturer’s instructions, from RNA extracted from
   isolated patient samples. qRT-PCR was performed on 50 ng cDNA from
   BE.LGD (n = 9) and BE.HGD (n = 11) patient samples using SsoAdvanced
   Universal SYBR Green Supermix (Bio-Rad Laboratories, Hercules, CA) with
   CFX Connect Thermo Cycler (Bio-Rad Laboratories), following the
   manufacturer’s instructions. Primer sequences for KRAS isoforms and
   GAPDH were acquired from previously published studies and primers for
   RNF128 isoforms, TRIM29 isoforms, and HPRT were designed using NCBI
   Primer-BLAST ([356]Table S11).[357]126, [358]127, [359]128 Primers were
   ordered from Eurofins Scientific (Louisville, KY). Relative isoform
   expression was quantified against housekeeping genes (GAPDH and HPRT)
   and isoform fraction was calculated by dividing relative isoform
   expression by the total relative expression of both isoforms in a
   condition.

Cell culture

   OE33 and JH-EsoAd1 (JHAD1) cells were cultured in RPMI 1640 medium
   containing 2.0 mM L-glutamine, 10^4 units/mL penicillin, 10^4 μg/mL
   streptomycin, 1 mM sodium pyruvate, and 10% fetal bovine serum (FBS).
   Cell culture reagents were acquired from either ThermoFisher Scientific
   or Sigma-Aldrich (St. Louis, MO). All cells were incubated at 37°C with
   5% CO[2] and maintained as monolayers. The OE33 cell line was obtained
   from the European Collection of Authenticated Cell Cultures (ECACC,
   Wiltshire, UK) and the JHAD1 cell line was kindly shared by Dr. James
   R. Eshleman (Johns Hopkins University, Baltimore, MD).

Cellular viability assay

   OE33 and JHAD1 cells were seeded in black-walled clear bottom 96-well
   plates at a concentration of 6,000 and 8,000 cells/well, respectively.
   Following overnight incubation, cells were treated with adapalene
   (1.0–2.0 μM; Tocris Bioscience, Bristol, UK) or APR-246 (10–40 μM;
   Cayman Chemical Company, Ann Arbor, MI). Treatments were prepared by
   diluting the stock solution to the desired experimental concentrations
   in RPMI 1640 medium containing 2.0 mM L-glutamine, 10^4 units/mL
   penicillin, 10^4 μg/mL streptomycin, 1 mM sodium pyruvate, and 5% FBS.
   Vehicle control was 0.08% DMSO diluted in the same medium. To determine
   cell viability, OE33 and JHAD1 cells were stained using Calcein-AM
   (ThermoFisher Scientific) at 48 and 72 h post-drug treatment as
   previously reported.[360]^129 Briefly, fluorescent images and readings
   were acquired using the SpectraMax MiniMax Imaging Cytometer (Molecular
   Devices, San Jose, CA). Following the initial fluorescence reading for
   determining viability, cells were replenished with fresh media and
   incubated for another 48 h to determine whether cells recovered upon
   drug removal. Calcein-AM staining was similarly used to quantify
   recovery of viability across drug concentrations.

Statistical analysis

   Survival analysis and calculation of statistical significance for
   viability data were performed using GraphPad Prism v9.1.0 (GraphPad
   Software, San Diego, CA). p-values ≤ 0.05 were considered significant
   unless otherwise noted. Survival analysis was performed on patients
   stratified by expression levels of specific transcripts. Patients
   were divided into tertiles based on TPM expression values for each
   transcript with upper compared with lower tertiles used to investigate
   survival differences using the Gehen-Breslow-Wilcoxon test.
   One-tail two-sample t-tests were applied to determine significant
   changes in isoform fraction between BE.LGD and BE.HGD patient samples
   in qRT-PCR experiments and statistically significant differences
   in cellular viability assay were determined by using a one-way ANOVA
   with Tukey’s post hoc multiple comparisons test.

Data availability statement

   RNA-sequencing files generated and analyzed during this study are
   available in the NCBI Gene Expression Omnibus (GEO: [361]GSE193946).

Acknowledgments