CstF64-Induced Shortening of the BID 3′UTR Promotes Esophageal Squamous Cell Carcinoma Progression by Disrupting ceRNA Cross-talk with ZFP36L2

Alternative polyadenylation (APA) is an RNA-processing event allowing for RNA polymerase II termination to occur at distinct locations, thus generating diverse 3′-ends (1). Nearly 70% of annotated human genes have more than one polyadenylation site (PAS) suggesting that APA is a major engine to drive transcriptomic diversity (2). Although a subset of APA occurs to alter the coding potential within C-terminal regions of proteins, the majority of APA events result in altered length and content of the 3′UTR. Although the coding sequence (CDS) is unchanged, through modulation of 3′UTR length, the impact of APA is profound as it can influence mRNA stability, mRNA translation, mRNA localization, and protein localization (1, 3). Moreover, APA can also re-route miRNA, resulting in altered expression of many genes that are not even subject to APA (4).

The importance of APA is only beginning to be appreciated and has been found to be involved in multiple biological processes and disease states. As a general observation, 3′UTRs show progressive lengthening during development and cellular differentiation (5–8), whereas present as broadly shortened in proliferating cells (9, 10). Dysregulated APA contributes to human diseases through mutations or germline variations that either create/destroy specific PASs (11–13) or through altered expression of APA regulators, resulting in more broad changes in PAS usage (1, 3). Among these disease states, cancer has been found to be highly sensitive to 3′UTR status where shortening is observed in tumors relative to normal tissues (14, 15). The preponderant view is that 3′UTR shortening within tumors is a mechanism to activate oncogenes through lost inhibition of miRNA (16). In addition, we and others have shown that broad 3′UTR shortening can disrupt competing endogenous RNA (ceRNA) networks leading to the repression of tumor suppressor genes through miRNA release and enhanced targeting (4, 17). Taken together, these studies suggest that 3′UTR shortening may promote tumorigenesis in cis, in trans or in combination through multiple mechanisms, and underscore the complexity of APA in regulating gene expression.

Mortality due to esophageal cancer ranks sixth all over the world (18) and is particularly problematic in China where it is the fourth leading cause of cancer-related death (19). The two primary histologic subtypes are esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC). In high-risk areas like China, ESCC comprises over 90% of all esophageal cancer cases (18). Through analysis of single-nucleotide variations (SNV) and copy-number variations (CNV) along with gene expression, similar profiles and patterns have been found among ESCC, head and neck squamous cell carcinoma (HNSC) and lung squamous cell carcinoma (LUSC), but not between ESCC and EAC (20, 21). This further supports the characterization of ESCC as a squamous-like subtype (22), suggesting that ESCC, HNSC, and LUSC may share some tumorigenic mechanisms and molecular signatures.

Despite the emerging evidence showing that APA plays a critical role in tumorigenesis, the importance of APA in ESCC has not been well described. In our previous work, we analyzed paired tumor and nontumor tissues from 94 patients with ESCC using whole-genome sequencing and high-depth RNA-seq and reported an integrated analysis of genomic alterations, including SNVs, CNVs, and structural variations, in these Chinese individuals with ESCC (20). BRCA1-associated protein and TP53-induced glycolysis and apoptosis regulator were shown to be overexpressed in ESCC due to copy-number gain and have been demonstrated to be involved in ESCC progression or chemoresistance (23, 24). This dataset (herein referred as “Chang ESCC”) represents a robust and distinct resource to further investigate molecular signatures in ESCC.

To explore the mechanism of aberrant regulation of gene expression at posttranscriptional level, we leveraged this unique resource in ESCC, and used the computational algorithm termed “Dynamic analyses of APA from RNA-Seq (DaPar),” to identify and quantify PAS usage (14). Through APA analysis of paired tumor/nontumor esophageal tissues of 94 individuals in the Chang ESCC RNA-seq dataset, we find a recurrent shortening of the BID 3′UTR mRNA, which leads to overexpression of BID in tumor tissues relative to nontumor controls. This shift of PAS selection is most likely due to a broad upregulation of CstF64 in tumors. Consequently, shortening of the BID 3′UTR inhibits the expression of tumor suppressor ZFP36L2 through disturbing the ceRNA cross-talk between the BID and ZFP36L2 3′UTR. Finally, we demonstrate that 3′UTR shortening of BID and low expression of ZFP36L2 are correlated with worse survival in ESCC thereby providing two new prognostic indicators of this aggressive disease. Overall, our results reveal a critical and previously undetected ceRNA network that consists of BID and ZFP36L2 3′UTR that is disrupted in ESCC and leads to increased tumor progression.

We used algorithm DaPars (https://github.com/ZhengXia/DaPars; ref. 14) to characterize APA events using standard RNA-seq data. In brief, DaPars uses a linear regression model to identify the de novo proximal poly(A) site (pPAS) and estimates the abundance of long and short form of 3′UTRs, then calculates the percentage of distal PAS usage index (PDUI) as relative PAS usage. The difference in PAS usage between tumor and nontumor samples is then quantified as a change in PDUI (ΔPDUI), whose negative or positive value represents shortening or lengthening of 3′UTRs, respectively. We used the Wilcoxon signed rank test to analyze paired tumor-normal sample sets and Wilcoxon rank sum test to analyze unpaired sample sets. Genes with FDR < 0.05 and ΔPDUI < −0.15 or > 0.15 in tumors comparing with nontumor tissues were considered significantly 3′UTR shortened or lengthened genes, respectively.

Samples were collected from patients histopathologically defined ESCC from Cancer Hospital, Chinese Academy of Medical Sciences (Beijing) and Cancer Hospital of Zhejiang Province (Hangzhou) between 2010 and 2014. A total of 94 patients underwent esophagectomy without chemotherapy or radiotherapy before surgery were involved in this study. Surgically removed tumors and paired adjacent nontumor esophageal tissues (≥5 cm from tumor site) were collected from each patient. Clinical information was obtained from medical records and survival time of patients with ESCC was measured from the date of diagnosis to the date of last follow-up or death (Supplementary Table S1). Written informed consents were obtained from all patients. This study was performed in accordance with the Helsinki Declaration and was approved by the Institutional Review Board of Cancer Hospital, Chinese Academy of Medical Sciences and Zhejiang Cancer Hospital.

To compare APA usage of ESCC with other types of cancers, PDUI data of HNSC (426 tumor and 42 nontumor), LUSC (220 tumor and 17 nontumor), cervical squamous cell carcinoma (CESC, 160 tumor and 3 nontumor), and stomach adenocarcinoma (STAD, 285 tumor and 33 nontumor) samples were obtained from a published article using algorithm DaPars to analyze APA usage (25). To compare the RNA level expression, RNA data of HNSC (521 tumor and 44 nontumor), LUSC (501 tumor and 51 nontumor), CESC (252 tumor and 3 nontumor), and STAD (415 tumor and 35 nontumor) samples were obtained from The Cancer Genome Atlas (TCGA) project (http://gdac.broadinstitute.org, RRID: SCR_003193).

Human ESCC cell lines KYSE30 (RRID: CVCL_1351), KYSE150 (RRID: CVCL_1348), KYSE180 (RRID: CVCL_1349) and KYSE450 (RRID: CVCL_1353) were kind gifts from Dr. Y. Shimada of Hyogo College of Medicine (Japan). Cell lines were authenticated by DNA finger printing analysis and confirmed to be free of Mycoplasma infection by qRT-PCR. All cell lines were cultured in RPMI-1640 medium with 10% FBS (Hyclone) at 37°C in a humidified incubator with 5% CO₂. All experiments were performed within 2–8 passages after thawing. We measured protein expression of BID and CstF64 in eight ESCC cell lines to exclude cell lines with too high/low level of BID and CstF64 expression. In addition, considering that shorter doubling time of cells is faster to get enough cells to perform experiments, we finally chose KYSE30 and KYSE450 to perform experiments.

We randomly selected 4 paraffin-embedded tissue sections for immunohistochemical analysis of BID expression in ESCC. The sections were incubated with primary antibody against BID (1:100, Abcam Cat# ab32060, RRID: AB_725689) at 4°C overnight and then detected with the ABC Kit (Pierce).

We used 3′-Rapid-amplification of cDNA ends (3′-RACE) to determine the poly(A) sites (PAS) of BID with SMARTer RACE 5′/3′ kit (TaKaRa). The sequences of gene-specific primers used for 3′-RACE are listed in Supplementary Table S2. The 3′-RACE products were then cloned and Sanger sequenced to accurately identify the BID PASs.

The CDS together with 3′UTR of BID was amplified and cloned into pcDNA3.1 vector at EcoRI and XbaI sites, yielding pcDNA3.1-BID. The resulting plasmid was used to test qRT-PCR primers for accuracy in measuring common and extended 3′UTRs of BID.

Reporter vectors in the psiCHECK-2 backbone (Promega) were generated containing the 218- short or 1,440-bp long BID 3′UTR using the restriction enzymes XhoI and NotI. To obtain expression of only the long 3′UTR of BID, the proximal poly(A) signal sequence was mutated from AAUAUA to CCCCCC using the Muta-direct Kit (SBS Genetech). The reporter vectors containing short or mutated long 3′UTR were designated as psiCHECK2–BID-short or psiCHECK2–BID-long, respectively. Reporter vector containing the 1917-bp ZFP36L2 3′UTR with the PmeI and NotI sites was also constructed and designated as psiCHECK2-ZFP36L2. Plasmids with mutations at the putative binding sites of miR-124–3p and miR-506–3p were also constructed and designated as psiCHECK2-BID-long-mut or psiCHECK2-ZFP36L2-mut.

eGFP was PCR-amplified and used to replace Renilla luciferase in psiCHECK2, psiCHECK2-BID-short, and psiCHECK2-BID-long at NheI and XhoI sites. The resulting plasmids were named as psiCHECK2-eGFP, psiCHECK2-eGFP-BID-short, and psiCHECK2-eGFP-BID-long and used to transfect KYSE450 cells together with mCherry construct as cotransfected control. The fluorescence of eGFP or mCherry were detected by fluorescence microscopy.

The CDS of CSTF2 (CstF64) was amplified and cloned into pcDNA3.1 vector at KpnI and NotI sites, yielding transient overexpressing plasmid pcDNA3.1-CstF64.

The CDS of ZFP36L2 was synthesized by Tsingke Biological Technology and cloned into pcDNA3.1 backbone, yielding pcDNA3.1-ZFP36L2. To construct plasmid containing both CDS and 3′UTR of ZFP36L2, linearized pcDNA3.1-ZFP36L2 vector was generated by PCR using CloneAmp HiFi PCR Premix (Takara). The 3′UTR of ZFP36L2 was amplified and cloned downstream of ZFP36L2 CDS into pcDNA3.1-ZFP36L2 vector using In-Fusion HD Cloning Kit (Takara), yielding pcDNA3.1-ZFP36L2–3UTR.

The authenticity of all the constructs was verified by Sanger sequencing. All primers used in plasmid construction are shown in Supplementary Table S2.

Total RNA extraction and reverse transcription were performed using TRizol reagent (Invitrogen) and PrimeScript RT master mix (TaKaRa), respectively. qRT-PCR was accomplished in triplicate using SYBR green reagent (Life Technologies). The accurate quantity of total and long isoforms of BID was calculated by using an absolute qRT-PCR method. We formulated standard curves with serial dilution approaches using the plasmid pcDNA3.1-BID, containing both the common and extended amplicons, as standard template. Each C_t value in the standard curve represents an exact quantity of common or extended 3′UTRs. The primer sequences used for qRT-PCR are shown in Supplementary Table S2.

Luciferase reporter assays was performed according to the manufacturer's instructions (Promega). Briefly, ESCC cells were seeded at 5∼10 × 10⁴ cells per well in 48-well plates, and 200 ng of reporter plasmid was transfected into cells or cotransfected into cells with 10 pmol of siRNA (Sigma-Aldrich) or 0.2, 1 or 5 pmol of miRNA mimics or 10 pmol of miRNA inhibitors (JTS scientific) using Lipofectamine 2000 (Invitrogen). Cells were collected 24 hours after transfection. Firefly luciferase activity was detected and used to normalize Renilla luciferase activity.

ESCC cells were treated with actinomycin D (Selleck) and RNA samples were collected at 0, 4, 8, and 12 hours after actinomycin D treatment. The expression levels of total or long BID mRNA were detected by qRT-PCR with primers located in the common or extended 3′UTR, respectively. The expression of short BID mRNA was calculated as the difference between total and long expression. The relative RNA levels were normalized against GAPDH RNA.

siRNAs targeting CstF64, CstF77, CPSF100, BID or ZFP36L2 were purchased from Sigma-Aldrich. Transfections of each siRNA were performed with Lipofectamine 2000 (Invitrogen). The siRNA sequences used for RNA-interfering were given in Supplementary Table S2.

Protein extracts from cells or tissue samples were subjected to SDS-PAGE and transferred to polyvinylidene difluoride membrane (Millipore). Antibody against BID (Cat# 2002, RRID: AB_10692485) was from Cell Signaling Technology. Antibody against ZFP36L2 (Cat# sc-365908, RRID: AB_10847357) was from Santa Cruz Biotechnology. Antibody against CstF64 (Cat# A301–092A, RRID: AB_873014), CstF77 (Cat# A301–096A, RRID: AB_873019), and CPSF100 (Cat# A301–582A, RRID: AB_1078863) was from Bethyl Laboratories. Membranes were incubated overnight at 4°C with primary antibody and visualized with a phototope–horseradish peroxidase Western Blot Detection kit (Thermo Fisher Scientific).

The miRNAs that target BID or ZFP36L2 3′UTRs were analyzed in silico using miRmap (https://mirmap.ezlab.org, RRID: SCR_016508; ref. 26). MiRmap combines thermodynamic, probabilistic, evolutionary, and sequence-based features to predict the miRNA repression strength as the “miRmap score.” We chose the overlapping miRNAs potentially targeting both BID and ZFP36L2 and those with “miRmap score” >85 for further experiments.

The CRISPR/Cas9 system was used to generate genomic deletion of extended region of BID 3′UTR in KYSE450 cells. Single guide RNAs (sgRNA) were designed and cloned into the sgRNA scaffold of plasmid pU6-(BbsI)_Myc-Cas9-T2A-mCherry. The resulting plasmids were cotransfected with a plasmid resistant to blasticidin into KYSE450 cells at 10:1. 24 hours after transfection, cells were selected with blasticidin (Thermo Fisher Scientific) for 48 hours and seeded into a 96-well plate for one cell per well. The cells were then genotyped and Sanger sequenced to select clones that successfully deleted extended region of BID 3′UTR. Sequences of sgRNAs and primers for genotyping were shown in Supplementary Table S2.

Δext and control KYSE450 cells were seeded in 6-well plates and each well contains 5 × 10⁵ cells. For induction of apoptosis, cells were treated with 20 ng/mL TNFα (Sigma-Aldrich) for 24 hours. Cells were then stained by the Annexin V-FITC/PI Apoptosis Detection Kit (KeyGEN BioTECH) and subsequently analyzed by flow cytometer.

For cell viability assays, Δext and control KYSE450 cells were plated in 96-well plates and each well contains 4,000 cells in 100 μL cell suspension. Cell viabilities were measured with the CCK-8 kit (Dojindo) at 24, 48, 72, and 96 hours after cultivation. For colony formation assays, 1,000 cells were seeded in 6-well plates and allowed to grow until visible colonies formed. Cell colonies were then fixed with methanol, stained with crystal violet, and counted.

Male NOD-scidIL2Rg^null (NSG) mice ages 4–6 weeks (Beijing IDMO Co., RRID: IMSR_JAX:005557) were subcutaneously injected with Δext and control KYSE450 cells (2 × 10⁶) and raised for 30 days. Tumor volume was measured every 5 days and calculated by 0.5 × length × width². Animal experiments were carried out in compliance with approved protocols and guidelines from the Institutional Animal Care and Use Committee of the Chinese Academy of Medical Sciences.

Spearman's correlation was used to calculate the correlations between the PDUI of BID and expression of BID or cleavage and polyadenylation (CPA) members. Pearson correlation was used to measure the correlations between expressions of two genes. Correlations were considered significant when P < 0.05 and |r| > 0.30. Kaplan–Meier analysis and log-rank test were used to compare survival by different levels of PDUI or expression, which were divided by quartile. Hazard ratios (HR) and 95% confidence intervals (CI) were calculated using Cox proportional hazards models with age, sex, tumor stage, smoking status, and drinking status as covariates. The Student t test was used to compare the difference of two groups when the data show a normal distribution; otherwise, the Wilcoxon rank sum test was used. A P value of <0.05 was considered significant for all statistical analyses.

The data generated in this study are available within the article and its Supplementary Data files or available upon reasonable request from the corresponding author.

APA has been commonly observed in cancers (14, 15) but the importance in ESCC is poorly understood. To systematically describe the APA landscape in ESCC, DaPars (14) was used to analyze RNA-seq data from 94 paired tumor and nontumor tissues of the Chang ESCC dataset (20). APA usage was quantified as the PDUI. The difference between tumor and nontumor samples was then described as ΔPDUI. Initially, we analyzed the RNA-seq data in the Chang ESCC dataset and found 90 genes had their 3′UTRs shortened in tumors and 24 genes exhibited 3′UTR lengthening (mean |ΔPDUI| > 0.15, FDR < 0.05). Among the 3′UTR shortened genes, BH3 interacting domain death agonist (BID) displayed the greatest difference in 3′UTR length between tumor and nontumor tissues with the most significant FDR, whereas 3′UTR of SPRR2G was the most lengthened (Fig. 1A). The 3′UTR has been shown to be the primary target for miRNA inhibition (27), thus genes that have shorter 3′UTRs are posited to exhibit higher expression. Among the 90 genes exhibiting 3′UTR shortening within our paired-patient sample cohort, we found that 20 genes also displayed higher expression in the tumor samples (Fig. 1B; Supplementary Table S3). Among these 20 candidate genes with increased expression and 3′UTR shortening, we further observed that 98.9% (93/94) of patients with ESCC had at least one 3′UTR shortening event (ΔPDUI < −0.15) within this group. Importantly, the most commonly shortened gene was BID where it could be detected in 82 out of 94 (87.2%) patient samples (Fig. 1C).

To further verify our findings using an independent collection of RNA-seq datasets, we performed a similar analysis in HNSC and LUSC datasets from TCGA. In all three datasets, BID was the only gene that exhibited both 3′UTR shortening and higher expression in tumors relative to normal controls (Fig. 1D; Supplementary Table S4). To investigate whether this phenomenon is specific in squamous cell carcinoma (SCC), TCGA RNA-seq datasets for both CESC and STAD were also analyzed using DaPars. Surprisingly, we found the BID 3′UTR to be significantly shortened in tumor versus nontumor tissues in all four SCC sample sets, but not in STAD (Fig. 1E; Supplementary Fig. S1A). We further found that BID mRNA level was also expressed at higher levels in all SCCs (Fig. 1F; Supplementary Fig. S1B). Moreover, we observed that 3′UTR length was inversely correlated with BID mRNA levels in the Chang ESCC dataset, as well as the LUSC and CESC datasets from TCGA (Fig. 1G; Supplementary Fig. S1C), suggesting that overexpression of BID may result from 3′UTR shortening in tumors. Finally, we conducted Western blot analysis in 12 pairs of randomly selected ESCC and their matched nontumor controls to test the protein level of BID in ESCC. The results showed that 11 of 12 (91.7%) tumor tissues had higher BID expression than the matched, nontumor counterparts (Fig. 1H). By IHC in 4 sections that containing ESCC and nontumor tissue in a single section, we observed a remarkably high BID staining signal in tumor tissues compared with the nearby nontumor tissues (Fig. 1I; Supplementary Fig. S1D). Overall, these results suggest that 3′UTR shortening of BID and higher expression of BID may represent a pan-SCC–specific phenomenon.

Given the recurrent shortening of the BID 3′UTR in SCCs, we decided to focus on characterizing the molecular mechanism of this phenomenon. To verify whether upregulation of BID is caused by its 3′UTR shortening, we first performed BID 3′-RACE using mRNA from the ESCC cell line KYSE450 and found that two distinct 3′UTR isoforms of BID were present (Fig. 2A and B). The two 3′-RACE products were cloned and Sanger sequenced to accurately identify and confirm the annotated BID pPAS and distal PAS (dPAS) are being used (Fig. 2C). On the basis of this information, we designed PCR amplicons to measure the common and extended 3′UTRs of BID (Fig. 2A; Supplementary Fig. S2A) and quantified the relative usage of the two BID 3′UTRs in four different ESCC cell lines (Fig. 2D; Supplementary Fig. S2B). In both KYSE30 and KYSE450 cell lines, we observed less expression of the lengthened BID 3′UTR relative to the total mRNA, indicating that BID 3′UTR shortening is occurring in tumor cells (Fig. 2D). To verify shortening of the BID 3′UTR, we performed qRT-PCR in 11 pairs of randomly selected ESCC and the matched nontumor tissues and found 10 out of 11 (90.0%) tumor tissues had shorter 3′UTRs compared with nontumor tissues (Fig. 2E). To determine whether the BID 3′UTR shortening can impact BID expression, we first constructed a heterologous dual luciferase reporter system where the Renilla luciferase ORF is placed upstream of either the short or long 3′UTR of BID (Supplementary Fig. S2C). Notably, the construct containing the long BID 3′UTR also had a mutated proximal poly(A) signal sequence to ensure that only the dPAS is used. Within the same plasmid, a control firefly luciferase ORF using a vector-derived PAS was constitutively expressed and was used to normalize transfection efficiency. Transfection of these reporter plasmids revealed that the short 3′UTR showed significantly higher luciferase activity relative to the long 3′UTR (Fig. 2F; Supplementary Fig. S2D). To further confirm the relationship between the expression and 3′UTR length using an orthogonal reporter system, we used GFP to fuse to the short or long BID 3′UTR and monitored GFP expression (Supplementary Fig. S2E). Cells transfected with the short form of the BID 3′UTR showed stronger GFP signal than cells expressing GFP containing the long BID 3′UTR (Fig. 2G). Through RNA decay assays, we found that BID mRNAs with short 3′UTRs were significantly more stable than those with long 3′UTRs (Fig. 2H). Overall, these results indicate that 3′UTR shortening is sufficient to induce expression of BID by increasing the stability of mRNA.

To find a candidate upstream regulator of BID poly(A) site selection, we analyzed the correlation between mRNA level expression of the 15 core members of the CPA machinery and PDUI of BID in the Chang ESCC dataset. Within this comparison, we found that CPSF100, CstF64, and CstF77 showed significant negative correlation with the PDUI of BID (Fig. 3A; Supplementary Table S5). We expanded our analysis to include not only measuring the expression of the CPA members within the Chang ESCC dataset but also the three SCC datasets from TCGA and found that CPSF100, CstF64, and CstF77 were overexpressed in at least three out of the four datasets. Notably, among these three CPA factors, CstF64 expression was most elevated in LUSC and CESC (Fig. 3B; Supplementary Fig. S3A). To determine which of these three CPA proteins is capable of functionally regulating BID PAS selection, we transfected KYSE30 and KY450 cells with siRNAs targeting CstF64, CstF77, or CPSF100 as well as a control siRNA. The knockdown efficiency of each siRNA was confirmed by Western blot (Fig. 3C; Supplementary Fig. S3B). We quantified BID poly(A) site usage using qRT-PCR analysis and found that knockdown of CstF64 led to the most significant lengthening of the BID 3′UTR (Fig. 3D; Supplementary Fig. S3C), and additional siRNA was used to confirm the on-target effect (Supplementary Fig. S3D and S3E). As a consequence of BID 3′UTR lengthening in CstF64 knockdown cells, expression of BID was significantly reduced (Supplementary Fig. S3D). Conversely, overexpression of CstF64 both shortened the BID 3′UTR and increased BID expression (Fig. 3E and F). Overall, these data suggest that CstF64 expression is a critical determinant of BID 3′UTR length.

BID has not been experimentally demonstrated to function as an oncogene and is not present on recent cataloguing of oncogenes. In light of this, we explored the hypothesis that BID 3′UTR shortening may promote tumor cell growth through repression of tumor-suppresser genes in trans by disruption of their ceRNA cross-talk (4). To find a potential ceRNA partner of BID in the Chang ESCC dataset, we designed an analysis pipeline (Fig. 4A). We posited that candidate genes would demonstrate coexpression with BID in nontumor tissues but would lose this coexpression in tumor tissues as BID undergoes 3′UTR shortening. Moreover, we predicted that the genes demonstrating disruption of coexpression in tumors will also be downregulated in tumor samples. These initial filtering criteria resulted in the identification of 93 genes as potential ceRNA for BID. To further refine the list, we added two more filters where a gene must have been predicted to be a tumor suppressor by Tumor Suppressor and Oncogene Explorer (28) and should have common miRNA-binding sites with BID. These additional criteria yielded the identification of ZFP36L2 as the single candidate ceRNA partner of BID (Fig. 4A). mRNA decay activator protein ZFP36L2 is an RNA-binding protein and binds to the AU-rich elements in the 3′UTRs of target genes to promote deadenylation and degradation (29). Numerous studies have showed that ZFP36L2 plays a tumor suppressor function, such as ZFP36L2 has been suggested as a candidate target gene of p53 and participate in p53-dependent apoptosis (30). Deletion of ZFP36L2 perturbs thymic development and leads to T lymphoblastic leukemia (31). Loss-of-function mutations in ZFP36L2 are found to be enriched in metastatic colorectal cancer and associated with poor overall survival. Moreover, loss of ZFP36L2 has been shown to promote metastasis of colorectal cancer cells (32). Altogether, these previous studies indicate the ZFP36L2 could be a tumorigenesis-promoting target for ceRNA disruption though BID 3′UTR shortening.

To characterize the relationship between ZFP36L2 and BID, we found a positive correlation in their expression in nontumor tissues (r = 0.40, P < 0.0001; Fig. 4B, left), whereas ZFP36L2 was uncorrelated with BID in tumor tissues, in agreement with 3′UTR shortening of BID (r = –0.04, P = 0.69; Fig. 4B, right). Moreover, in the Chang ESCC dataset, we observed that ZFP36L2 mRNA levels were significantly lower in tumor than in their nontumor counterparts (Fig. 4C, left). We further investigated the expression of ZFP36L2 in HNSC, LUSC, and STAD from TCGA. We observed that ZFP36L2 was downregulated in HNSC, LUSC, and CESC, but not in STAD (Fig. 4C; Supplementary Fig. S4), consistent with the results of BID 3′UTR length in SCCs and STAD, which is similarly suggestive of a pan-SCC–specific phenomenon.

To directly test whether BID and ZFP36L2 exist as part of a ceRNA network, we first examined the regulatory relationship between the two genes. We designed four siRNAs to knockdown BID, the first two target the CDS of BID and the other two target the extended region of BID 3′UTR (Fig. 4D). Notably, siRNA#1 and siRNA#2 can target both BID isoforms whereas siRNAs #3 and siRNA#4 only target the BID long 3′UTR. We transfected KYSE30 and KYSE450 cells with each of the four siRNAs and found that, as expected, siRNAs targeting CDS led to more significant reduction of BID expression than siRNAs targeting extended 3′UTR. Despite the observed difference in BID knockdown efficiency, all four siRNAs could decrease the expression of ZFP36L2 comparably, suggesting that the BID long 3′UTR isoform is the more potent ceRNA to ZFP36L2 (Fig. 4E). Moreover, transfection of siRNAs targeting ZFP36L2 results in the downregulation of ZFP36L2 as well as BID (Fig. 4F). Furthermore, we cotransfected KYSE450 cells with a reporter plasmid fusing the 3′UTR of ZFP36L2 and along with two siRNAs targeting BID. Knockdown of BID significantly repressed the luciferase activity of ZFP36L2 reporter (Fig. 4G). We also conducted the reciprocal experiment and found that that luciferase activity of a reporter with the long BID 3′UTR was inhibited after knockdown of ZFP36L2 in KYSE450 cells (Fig. 4H). These results indicate that endogenous protein expression of BID and ZFP36L2 is interdependent, and their respective 3′UTRs are sufficient to mediate this coexpression.

To identify miRNAs that target both BID and ZFP36L2, we conducted an in silico analysis with miRmap (26). Using this approach, we identified miR-124–3p, miR-143–5p, miR-506–3p, miR-3944–5p, and miR-4796–3p as potential shared miRNAs (Supplementary Fig. S5A; Supplementary Table S6). Putative binding sites of each miRNA are shown in distinct colors (Fig. 5A). We cotransfected BID short, long or ZFP36L2 3′UTR luciferase reporters with miRNA mimics of each of these candidate miRNAs and found that miR-124–3p and miR-506–3p could potently inhibit luciferase activity from both the BID long and ZFP36L2 reporters but not from the BID short reporter (Fig. 5B). Notably, the repression of the BID long and ZFP36L2 reporters is dose dependent (Fig. 5C). These results suggest that shortening of the BID 3′UTR abrogates the binding sites of miR-124–3p and miR-506–3p, inducing BID expression by escaping the inhibition of miRNAs. In addition, transfection of inhibitors of miR-124–3p and of miR-506–3p together could increase the luciferase activity of both BID long and ZFP36L2 reporters (Supplementary Fig. S5B). To further validate the specific inhibitory effect of miR-124–3p and miR-506–3p, the predicted binding sites of miR-124–3p and miR-506–3p were mutated in the reporter plasmids of the BID long and ZFP36L2 3′UTR (Supplementary Fig. S5C). Mutation of the putative binding sites diminished the effect of miR-124–3p and miR-506–3p mimics, resulting in reduced repression (Fig. 5D). Together, these results suggest that BID and ZFP36L2 act as ceRNA partners by competitive binding of miR-124–3p and miR-506–3p.

Our results, thus far, indicate that BID 3′UTR shortening can enhance tumorigenesis through disruption of a ceRNA network with ZFP36L2. To directly determine the biological consequences of BID 3′UTR loss, we performed CRISPR/Cas9 genome-editing to construct a KYSE450 cell line that lacks the extended region of the BID 3′UTR. Our rationale was that although KYSE450 cells preferentially express the short BID 3′UTR, there are still significant levels of the long 3′UTR BID isoform, thus deletion of this extended region could result in a greater release of miRNA causing further downregulation of ZFP36L2. Using two sgRNAs targeting regions flanking the extended region of BID 3′UTR, we isolated and validated two independent BID 3′UTR deletion clonal lines, which we named “Δext#1 and Δext#2” (Fig. 6A; Supplementary Fig. S6A). Comparing with control (parental) cell lines, luciferase activity of reporter with ZFP36L2 3′UTR was reduced in Δext cells (Fig. 6B). Furthermore, protein levels of BID were upregulated, whereas ZFP36L2 were downregulated in Δext cells relative to control cells (Fig. 6C). These results further indicate that BID and ZFP36L2 expressions are coregulated and that this is mediated through the BID 3′UTR.

We then examined the effects of 3′UTR shortening on tumor cell growth and observed that the proliferation of Δext cells was significantly enhanced relative to control cells (Fig. 6D). The increased proliferation was unlikely to be caused by alterations in cell death as we observed no difference at the level of apoptosis between Δext and control cells (Supplementary Fig. S6C). Moreover, we also observed that colony formation was markedly increased in Δext cells compared with control cells (Fig. 6E). Finally, we transplanted Δext and control cells subcutaneously into NOD-scidIL2Rg^null (NSG) mice and found that xenograft tumors derived from Δext cells had significantly higher tumor volumes (Fig. 6F and G) and tumor weights (Fig. 6H) than control cells.

To explore whether the phenotype of Δext cells is induced by downregulation of ZFP36L2, we transfected Δext cells with a plasmid bearing the CDS of ZFP36L2 and found that overexpression of ZFP36L2 could inhibit proliferation of Δext cells back to the level of control cells (Fig. 6I; Supplementary Fig. S6D). To further investigate whether the ceRNA disruption contributes to the phenotypic changes in Δext cells, we constructed plasmids bearing the ZFP36L2 CDS together with its 3′UTR (ZFP36L2–3UTR) or the extended BID 3′UTR (BID-ext). We observed that transfection of ZFP36L2–3UTR raised the level of ZFP36L2 protein and could also reduce cell proliferation (Fig. 6J). Importantly, we also observed that cotransfection of the inhibitors of miR-124–3p and miR-506–3p could further increase ZFP36L2 expression and cause an even greater reduction in cell growth (Fig. 6J). We also found that cotransfection of the BID-ext plasmid could increase the expression of ZFP36L2 protein and further reduce cell proliferation (Fig. 6K). Collectively, these results show that 3′UTR shortening of BID is sufficient to drive ESCC tumorigenesis, which is mediated, at least in part, by trans inhibition of ZFP36L2.

To investigate the clinical relevance of this ceRNA regulatory network, we examined the relationships between survival of all the 94 patients with ESCC and 3′UTR length of BID, expression of ZFP36L2, or expression of BID. The results showed that only the BID 3′UTR length and ZFP36L2 mRNA level expression were significantly correlated with the 94 patients with ESCC survival time, whereas the expression of BID was not correlated (Fig. 7A). We found that the 3′UTR shortening of BID, quantified by lower PDUI value, had an excess HR for death compared with BID long 3′UTRs in Chang ESCC dataset (HR, 2.07; 95% CI, 1.10–3.88; P = 0.024; Fig. 7A, left). Furthermore, low expression of ZFP36L2 also had an excess HR (HR, 1.91; 95% CI, 1.14–3.19; P = 0.014; Fig. 7A, middle). These data suggest that the 3′UTR shortening of BID (as opposed to BID overexpression) may play a role in tumorigenesis and this is likely mediated through misregulation of ZFP36L2 in trans.

We also conducted a stratified analysis of the Chang ESCC patients' dataset and found that correlations between short survival time and 3′UTR shortening of BID or low repression of ZFP36L2 were not detected in early patients with ESCC (stage II; short BID: HR, 1.53; 95% CI, 0.41–5.78; P = 0.53; Fig. 7B, left; low ZFP36L2: HR, 1.09; 95% CI, 0.35–3.38; P = 0.89; Fig. 7B, middle), but only occurred in progressed patients with ESCC (stage III and IV; short BID: HR, 3.23; 95% CI, 1.39–7.55; P = 0.0066; Fig. 7C, left; low ZFP36L2: HR, 2.60; 95% CI, 1.43–4.74; P = 0.0018; Fig. 7C, middle). Overall, these results indicate that 3′UTR shortening of BID and low expression of ZFP36L2 can predict worse survival, acting as new clinical predictors. Furthermore, the prediction power becomes more significant in patients with ESCC of stage III and IV, suggesting that 3′UTR shortening of BID may be a more important feature in advanced tumors.

Our model, on the basis of the ceRNA theory, shows that BID may compete with ZFP36L2 for miR-124–3p and miR-506–3p inhibition (Fig. 8). In tumor cells, upregulated CstF64 may contribute to the shortening of the BID 3′UTR. This shortening, in turn, allows BID mRNA to escape the inhibition of miRNAs binding on its extended 3′UTR, releasing miR-124–3p and miR-506–3p, thereby breaking the balance between BID and ZFP36L2. The lost ceRNA network ultimately results in downregulation of ZFP36L2 expression. Decreased ZFP36L2 promotes the proliferation of tumor cells and also predicts poor survival of patients with ESCC.

RNA-seq has become a common method for gene expression analysis and has been involved in numerous large-scale genomic studies. RNA-seq data available from TCGA have become especially valuable to retrospective cancer investigations. However, one potential limitation to the TCGA repository is a relative dearth of comparative nontumor tissue RNA-seq relative to what is available for tumors. Our Chang ESCC dataset is particularly useful to analyze in parallel with TCGA as it includes 94 tumor tissues and 94 paired nontumor tissues, which enabled our study here to effectively discover 3′UTR shortening genes in esophageal tumors compared with nontumor tissue. Our analysis in the Chang ESCC and TCGA datasets revealed a pattern of BID 3′UTR length and ZFP36L2 expression similar in ESCC, HNSC, LUSC, and CESC, but differing with STAD (Figs. 1E and 4C; Supplementary Fig. S1A and S4), consistent with our previous finding that gene expression pattern in ESCC is similar to that in HNSC and LUSC, but differs from that in STAD (20), suggesting that SCCs may share some molecular mechanisms of pathogenesis.

The polyadenylation machinery is mainly composed of four protein subcomplexes, CPA specificity factor (CPSF), cleavage stimulation factor (CstF), cleavage factor I (CFI) and cleavage factor II (CFII). Modulation of the expression levels of core polyadenylation machinery components has been revealed as an important mechanism of APA regulation (33). Among these factors, the CstF complex includes subunits of 50, 64, and 77 kDa with the CstF64 subunit directly contacting the G/U-rich region typically located downstream of the poly(A) signal (3). Modulation of CstF64 expression has been found to alter PAS selection and this was first demonstrated where induction of CstF64 is sufficient to increase the use of upstream intronic PAS of the IgM heavy chain transcript during B-cell differentiation (34). Consistently, others have found that depletion of CstF64 leads to an increase of relative use of distal PAS (35). In our studies, we find that CstF64 is overexpressed in Chang ESCC dataset and LUSC and CESC from TCGA, and the CstF64 expression levels are particularly elevated in LUSC and CESC (Fig. 3B). Our results mirror those from others where CstF64 is reported to be upregulated in tumor tissues relative to nontumor in five out of seven measured cancer types (14). Moreover, CstF64 was shown to have the greatest correlation between gene expression fold change and the number of 3′UTR shortening events per sample in tumors (14). Furthermore, mRNA level expression of CstF64 has been reported to be negatively correlated with 3′UTR lengths (r = −0.72). This observation contrasts with what was observed for the paralog of CstF64, CstF64tau, which was uncorrelated (r = −0.01; ref. 15). These studies are consistent with our results where CstF64 but not CstF64tau is negatively correlated with 3′UTR length of BID in ESCC (Supplementary Table S5). Collectively, our data suggest that upregulation of CstF64 alone may be powerful enough to shorten 3′UTRs without any additional changes in CstF64tau expression in tumors.

Initially, it was not expected that BID 3′UTR shortening would play such a prevalent role in ESCC tumor progression. BID is well-described as a pro-apoptotic member of BCL2 family and is cleaved by caspase 8 into its active form, truncated BID (tBID), linking the extrinsic, and intrinsic apoptosis pathway (36, 37). Although resistance to apoptosis is considered as a hallmark of cancer (38), increasing studies have found that proapoptotic proteins can also contribute to the development of tumor, notably both CD95 and PUMA (39, 40). This is also the case for BID as its depletion in mice significantly delayed tumor development, suggesting a tumor-promoting role of BID (41–44). Furthermore, in a diethylnitrosamine-induced hepatocellular carcinoma model, hepatocyte-specific BID depletion is thought to protect mice from tumorigenesis by suppressing inflammation-related compensatory proliferation (42). This phenotype is likely related to a nonapoptotic function of BID in full-length form, which produces inflammatory cytokines in response to nucleotide-binding and oligomerization domain (NOD) activation (45). Overall, these studies may explain our finding that BID is overexpressed in ESCC (Fig. 1H and I) but does not promote apoptosis (Supplementary Fig. S6C), supporting an oncogenic potential in ESCC through its non-apoptotic function.

Despite being assigned a Gene Ontology term and pathway analysis indicating ZFP36L2 as involved in cell cycle regulation, its role in tumorigenesis remains less well characterized (28). In a study characterizing ESCC and EAC, ZFP36L2 has been recognized as an important tumor suppressor gene specific to ESCC (21). Both the downregulation and the antiproliferation phenotypes are shown to be specific to ESCC but not EAC. Furthermore, the downregulation of ZFP36L2 in SCCs has been attributed to SCC-specific hypermethylation in its super-enhancer (21). Here, we provide a distinct mechanism of ZFP36L2 repression in SCCs. We find that downregulation of ZFP36L2 (Fig. 4C; Supplementary Fig. S4) and 3′UTR shortening of BID (Fig. 1E; Supplementary Fig. S1A) coordinately occur in a variety of SCCs, but not in STAD, suggesting a SCC-specific ceRNA cross-talk between BID and ZFP36L2. Moreover, both 3′UTR shortening of BID and low expression of ZFP36L2 are correlated with worse survival of ESCC, whereas expression of BID is not correlated (Fig. 7A), supporting the hypothesis that shortening of the BID 3′UTR promotes tumorigenesis primarily through inhibiting ZFP36L2 rather than upregulating its own expression.

Revealing the proproliferation effect and prognostic power of BID 3′UTR shortening, our study provides the potential of a novel therapeutic strategy through restoring 3′UTR length of BID. Approaches promoting distal poly(A) signal usage of specific genes by blocking access of the CPA machinery to proximal poly(A) signal are emerging. In particular, antisense oligonucleotides (ASO) could target proximal poly(A) signal by complementary base pairing, leading redirection of CPA complex to distal poly(A) signal (46). Indeed, others have shown that ASOs can treat prostate cancer by blocking an intronic poly(A) signal to restore full-length androgen receptor expression and inhibit androgen-independent proliferation (47). Thus, further functional experiments could contribute to the development of a new therapeutic strategy targeting poly(A) signal of BID.

In summary, our data reveal a recurrent shortening of the BID 3′UTR in SCCs and a previously unidentified ceRNA network that connects the BID and ZFP36L2 3′UTRs. As ESCC undergoes tumor progression, we propose that shortening of the BID 3′UTR leads to both the upregulation of BID and downregulation of ZFP36L2. Deletion of 3′UTR of BID promotes the proliferation activity of ESCC cells through disrupting the ceRNA cross-talk of BID and ZFP36L2. In addition, our study identified a significant association between poor survival and 3′UTR shortening of BID or downregulation of ZFP36L2 in patients with ESCC. Although overexpression of CstF64 is partially responsible for the 3′UTR shortening of BID, the complete set of factors that determine BID PAS selection remains undetermined and could represent novel therapeutic targets. Overall, our findings provide evidence for the first time that 3′UTR shortening of BID and downregulation of ZFP36L2 are widespread phenomenon specifically in SCCs and may serve as prognostic markers and therapeutic targets in ESCC.

D. Lin reports grants from Natural Science Foundation of China during the conduct of the study. E.J. Wagner reports a patent for Poly-Click-seq issued and licensed to Eric Wagner. No disclosures were reported by the other authors.

A. Lin: Conceptualization, data curation, formal analysis, validation, investigation, visualization, writing–original draft, writing–review and editing. P. Ji: Formal analysis, funding acquisition, project administration, writing–review and editing. X. Niu: Formal analysis, validation, investigation. X. Zhao: Data curation, software, formal analysis. Y. Chen: Data curation, software, formal analysis. W. Liu: Investigation. Y. Liu: Investigation. W. Fan: Formal analysis. Y. Sun: Formal analysis. C. Miao: Investigation. S. Zhang: Project administration. W. Tan: Funding acquisition, project administration. D. Lin: Supervision, funding acquisition. E.J. Wagner: Conceptualization, resources, supervision, funding acquisition, writing–review and editing. C. Wu: Conceptualization, resources, supervision, funding acquisition, writing–review and editing.

This research was supported by National Science Fund for Distinguished Young Scholars (81725015 to C. Wu), Medical and Health Technology Innovation Project of Chinese Academy of Medical Sciences (2016-I2M-3–019 to D. Lin; 2016-I2M-4–002 to C. Wu; 2019-I2M-2–001 to D. Lin and C. Wu; 2016-I2M-1–001 and 2019–12M-1–003 to W. Tan), Beijing Outstanding Young Scientist Program (BJJWZYJH01201910023027 to C. Wu) and National Natural Science Foundation of China (81988101 to D. Lin and C. Wu). E.J. Wagner and P. Ji acknowledge startup funds from the Department of Biochemistry and Molecular Biology at UTMB.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.