The data of 24 cancer types from TCGA project, including gene expression profile and clinical information, were downloaded from the UCSC XENA database (https://xenabrowser.net/datapages/) to determine the expression levels of TP53TG1. Prognosis analysis of TCGA data was performed with GEPIA2 (Tang et al., 2017).

Unique molecular identifier (UMI) count matrix of single-cell RNA-seq data from one cervical cancer tissue sample and one normal adjacent tissue sample was downloaded from GSE168652 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE168652). The 10x UMI count matrix was converted into a Seurat object using the R package Seurat (Butler et al., 2018) (version 4.0.4). Cells with UMI numbers <500 or with detected genes <200 or with over 15% mitochondrial-derived UMI counts were considered low-quality, and thus were removed. Genes detected in less than three cells were removed for downstream analysis.

After quality control, the UMI count matrix was log-normalized. Then the top 2,000 variable genes were used to create potential anchors with Find Integration Anchors function of Seurat. Subsequently, Integrate Data function was used to integrate data. To reduce the dimensionality of the scRNA-Seq dataset, principal component analysis (PCA) was performed on an integrated data matrix. With Elbowplot function of Seurat, the top 50 PCs were used to perform the downstream analysis. The main cell clusters were identified by the Find Clusters function offered by Seurat, with resolution set as default (res = 0.6). Finally, cells were clustered into 14 major cell types. Then they were visualized using tSNE or UMAP plots. To identify the cell type for each cluster, gene markers for each cell cluster were detected using the “Find Markers” function in Seurat package (v4.0.4). Then cell types were annotated using the cell markers provided in Gene Expression Omnibus (GEO) dataset (Ho et al., 2019).

pcDNA3.1-TP53TG1 (TP53TG1-overexpressing plasmid) was purchased from Youbio Biotech (Changsha, PRC). Full-length sequence of the plasmid:

ccc​tgt​ctc​cag​tgg​gcg​tct​tgg​gcc​ccg​gct​cta​ttc​tgg​gct​gcg​ggc​ctg​gga​agggct​cgc​cgg​gtg​cca​aat​gag​ctg​tcc​taa​ctc​tgc​ggg​gct​gca​gct​tcc​tgc​atg​atg​ctg​gggagc​ttg​gcg​cct​gac​cca​gga​tct​aga​agg​cac​tct​ggg​cag​gcc​gcg​ctc​cgc​cca​cga​agg​tac​cca​acc​ctc​tgg​gat​aga​tgc​agg​aag​cga​tgg​tta​aga​ccc​att​ttc​acc​caa​ctt​ctc​gcc​gca​ggt​ctg​gct​tac​cac​acg​ctc​ctc​ccc​att​ccc​agt​gag​ccg​ctt​ttt​gca​gca​cca​ggc​gaa​cac​tta​cac​cag​tgc​ttt​gta​aag​gaa​tct​tat​tgt​cca​ccc​cgt​gtc​ttg​gca​aaa​gaa​cag​tga​tca​cac​aga​ttc​cta​ctt​ggg​ctc​ttt​cct​tta​atc​ttc​gga​ggc​tga​gtt​tgc​cca​act​cag​gtttaa​cca​cca​agg​act​ctg​aga​gct​ggc​agg​tct​gag​taa​ccc​tgg​taa​caa​ttc​tct​tca​cct​tat​caa​aac​ctg​agc​taa​aac​caa​tgc​atc​agc​tga​tga​tga​cag​cag​aga​gtg​gca​ggg​ctg​agg​acc​caa​agt​cat​ttc​cca​ggc​tgg​cgg​aga​ata​aac​tgc​cag​gga​gaa​gaa​tga​gaa​gac​agg​agacaaactgtt​tgg​aaa​gct​aaa​tct​tcc​ctc​tta​atg​aat​aaa​ggt​ttt​tgc​ctt​gtc​tta​aaa​aaa​aa.

GAPDH (glyceraldehyde-3-phosphate dehydrogenase) was used as a control gene for assessing the effects of TP53TG1 overexpression. cDNA was synthesized according to standard procedures and RT-qPCR was performed on the Bio-Rad S1000 with Hieff qPCR SYBR^® Green Master Mix (Low Rox Plus; YEASEN, China). The information of primers is presented in Additional File 1. The concentration of each transcript was then normalized to GAPDH mRNA level using 2^−ΔΔCT method (Livak and Schmittgen 2001).

The cell proliferation assay was conducted using a Cell Counting kit-8 (CCK-8; Dojindo Molecular Technologies, Inc., Shanghai, China). Briefly, TP53TG1-overexpressing Hela cells and control Hela cells were seeded with 6,000 cells/well in 96-well culture plates. Cells treated with an equal volume of phosphate-buffered saline (PBS) served as controls and vials without cells were used as blank controls. After incubated for 24, 48, and 72 h, 20 μl CCK-8 solution was added to the culture medium and incubated for an additional 3 h. Subsequently, the optical density (OD) of the cells was measured with a PerkinElmer/envision at an absorbance of 450 nm. The cell proliferation rate was calculated using the formula: proliferation rate = (experimental OD value − blank OD value)/(control OD value − blank OD value) × 100%.

Total RNA was extracted using TRIzol reagent and was purified twice in phenol-chloroform. To remove DNA, the purified RNA was then treated with RNase-free RQ1 DNase (Promega Corp.) and its quality and quantity were determined by measuring the absorbance at 260 nm/280 nm (A260/A280) using a Smartspec Plus (Bio-Rad Laboratories, Inc.). The integrity of RNA was then verified by 1.5% agarose gel electrophoresis.

A total of 10 μg total RNA from each sample was used to prepare a directional RNA-seq library. First, the polyadenylated mRNAs were concentrated with oligo (dT)-conjugated magnetic beads (Invitrogen; Thermo Fisher Scientific, Inc.). The concentrated mRNAs were then iron-fragmented at 95°C, end-repaired and ligated to a 5′ adaptor. Reverse transcription (RT) was performed with RT primer harboring a 3′ adaptor sequence and randomized hexamer. The purified cDNAs were amplified and stored at −80°C for subsequent sequencing. Following the manufacturer’s instructions, the libraries were prepared for high-throughput sequencing. The Illumina Novaseq 6000 (Illumina, Inc.) was used to collect data from 150-bp paired-end sequencing (Illumina, Inc.).

Raw sequencing reads containing more than 2-N bases were first discarded. Subsequently, adaptors and low-quality bases were trimmed off the raw reads using a FASTX-Toolkit (v.0.0.13; http://hannonlab.cshl.edu/fastx toolkit/). Short reads of less than 16 nt were dropped. Then clean reads were subsequently aligned to the GRch38 genome by HISAT2 (Kim et al., 2015). Uniquely mapped reads were ultimately used to calculate read number and paired-end fragments per kilobase of exon per million fragments mapped (FPKM) for each gene.

The gene expression levels were measured using FPKM. The DEGs were screened using the software DEseq2 (Love et al., 2014), which analyzes DEGs. Using the fold change (FC ≥ 2 or ≤0.5) and false discovery rate (FDR<0.05) as the cutoff, the results were analyzed to determine whether a gene was differentially expressed.

The alternative splicing events (TP53TG1-AS) and regulated alternative splicing events (TP53TG1-RAS) between TP53TG1 overexpression (TP53TG1-OE) samples and control samples were defined and quantified using the SUVA pipeline described previously (Cheng et al., 2021a). The frequency and proportion of reads of SUVA AS events (pSAR) were calculated.

TP53TG1-mediated ceRNA network. Miranda (score ≥ 150) and Rnahybrid (p-value ≤ 0.05) were used to predict the target relationship between miRNA and TP53TG1. Finally, the results of the two methods were intersected. miRNA-mRNA target pairs from miRDB (http://mirdb.org) and targetscan (http://www.targetscan.org/vert_80/) databases were used to predict the target relationship between miRNA and TP53TG1-regulated DEGs. The regulatory network was established using Cytoscape (https://cytoscape.org/).

Using the KOBAS 2.0 server (Xie et al., 2011), GO analyses and enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were adopted to predict the functions of genes and calculate the distribution frequency in each functional category. The enrichment of each pathway (corrected p < 0.05) was defined using hypergeometric tests and the Benjamini-Hochberg FDR controlling procedure.

One-way analysis of variance and GEPIA2 were used for comparison of expression levels of TP53TG1 among 24 TCGA cancer types, with disease state (Tumor or Normal) as the variable for calculating differential expression. Student’s t-test was used for all other comparisons between the TP53TG1-OE and control groups. For each assay, the results are presented as the mean ± standard error of the mean values of three experiments. The data were analyzed using R software (v3.5.3, https://www.r-project.org/). p < 0.01 or FDR<0.05 was set as the threshold to indicate a statistically significant difference.

Inspired by previous discoveries that TP53TG1 may play oncogenic or anti-oncogenic roles in different cancers (Zhang et al., 2019b; Islam et al., 2021), TP53TG1 expression was compared between tumor and normal tissues in 24 cancer types from TCGA. TP53TG1 expression in tumors was significantly up-regulated (p-value < 0.05) in 15 of the 24 cancer types and was significantly down-regulated only in two cancer types, including colorectal and gastric cancer (Figure 1A). This is consistent with that p53 transcriptionally activates TP53TG1 when the DNA-damaging agent is used except that tumor-specific promoter CpG island hypermethylation-associated silencing of the lncRNA TP53TG1 occurs in colorectal and gastric cancer cells (Diaz-Lagares et al., 2016). Furthermore, GEPIA2 was used to comprehensively analyze the association of TP53TG1 expression with survival rates in 2,376 patients. Pan-cancer analysis shows that higher TP53TG1 expression was positively correlated with a higher overall survival rate (Figure 1B), which confirmed that TP53TG1 exhibits tumor suppressor-like feature in cancer patients (Diaz-Lagares et al., 2016). Furthermore, the association of TP53TG1 expression with survival rates in various types of cancer was also depicted using the survival map (Figure 1C). Surprisingly, we found that TP53TG1 showed positive or negative correlations with prognosis in different cancers, and that this correlation was significant in only five of them. Taken together, these results suggest that TP53TG1 may have a dual role in cancer development, with a predominant oncogenic function. Its function is associated with regulatory networks in cancerous tissue.

For cervical cancer (CESC), the expression of TP53TG1 in tumor tissue was higher than that in adjacent normal tissue (Figure 1A; Supplementary Figure S1A). Prognostic analysis showed that high expression of TP53TG1 was associated with a good prognosis (Figure 1B; Supplementary Figure S1B). Since bulk RNA-seq obtains the average gene expression of various cells, to know whether TP53TG1 is highly expressed in malignant cells or in other cell types, we downloaded and analyzed the single cell data of cervical cancer tissue and adjacent tissue, respectively. (GSE168652) (Li et al., 2021) (Figure 1D). We found that TP53TG1 was highly expressed in cell types C0 (cancer cells), C2 (cancer cells), C6 (endothelial cells), C9 (cancer cells) and C13 (cancer cells) (Figure 1E). These results reveal that TP53TG1 is activated in many malignant tumors and has a potential inhibitory effect on cancer development, resulting in a better prognosis.

To further investigate the function of TP53TG1 in cervical cancer, a TP53TG1-overexpression (TP53TG1-OE) cell model was constructed by transfecting a TP53TG1-overexpression plasmid into HeLa cells, with empty vector as negative control (NC). The overexpression of TP53TG1 was verified by RT-qPCR. The result showed that TP53TG1 gene expression significantly increased in TP53TG1-OE group (p < 0.001; Figure 2A). To characterize the role of TP53TG1 in regulating cell apoptosis and proliferation in HeLa cells, flow cytometric analyses and CCK8 assay were performed, respectively. Cell proliferation significantly decreased in the TP53TG1-OE group after 24 h (p < 0.01; Figure 2B; Supplementary Figure S2A). Meanwhile, TP53TG1-OE significantly promoted cell apoptosis (p < 0.001; Figure 2C). Together, these results indicated that TP53TG1 plays a tumor suppressor role in the progression of Hela cells.

LncRNAs participate in the regulation of gene expression and play vital roles in various biological and pathological processes (Liao et al., 2016). To investigate TP53TG1-mediated transcriptional or post-transcriptional regulation in HeLa cells, cDNA libraries for TP53TG1-OE cells and control cells were constructed for RNA-seq using the Illumina HiSeq Novaseq 6000 platform. Three biological replicates for each group were used and a total of 113–135 Million (M) 150-nucleotide paired-end raw reads per sample were obtained. After removing adaptors and low-quality reads, 103.5–125.4 M clean reads were aligned to the human GRCH38 genome using HISAT2, of which 86.27%–90.24% were uniquely aligned for further analysis (Supplementary Table S1). The level of gene expression was calculated in the units of FPKM (Supplementary Table S2). A total of 26,951 expressed genes were assessed by RNA-seq. Effective overexpression of TP53TG1 was further confirmed by the use of RNA-seq analysis (Figure 3A). The differentially expressed genes (DEGs) between the TP53TG1-OE and control cells were determined using an absolute FC ≥ 2 and a 5% FDR as criteria coupled with the DESeq2 package (Love et al., 2014). A total of 704 up-regulated and 470 downregulated DEGs were identified (Supplementary Table S3). The DEGs associated with TP53TG1-OE are displayed in a volcano plot (Figure 3B).

To further explore the potential biological roles of these DEGs, GO (Supplementary Table S4) and KEGG (Supplementary Table S5) enrichment analyses were performed. The top 10 GO terms in the category biological process were shown in Figure 3C and Figure 3D. And the top 10 KEGG terms in the category biological process were shown in Supplementary Figures S3A,B. The upregulated genes in the TP53TG1-OE cells are mainly enriched in type I interferon signaling pathways, response to virus, and positive regulation of cell migration terms (Figure 3C). Previously, type I interferon was reported to inhibit the growth of tumor stem cells and promote the apoptosis of cancer cells (Bekisz et al., 2010; Aricò et al., 2019). The downregulated genes were mostly associated with DNA damage responses, such as cellular response to DNA damage stimulus and DNA repair. Genes involved in RNA splicing and regulation of cell migration were also significantly enriched in downregulated genes (Figure 3D). DNA damage is closely related to apoptosis in cancer (Zhang et al., 2007; Figueroa-González and Pérez-Plasencia, 2017). These results indicated that TP53TG1 could active the type I interferon signaling pathways and suppress DNA damage responses in HeLa cells, which ultimately promotes the apoptosis of cancer cells.

To confirm the regulatory function of TP53TG1 in gene expression in the HeLa cell line, several DEGs were involved in type I interferon signaling pathways, cellular response to DNA damage stimulus and DNA repair, which were validated by RT-qPCR (Figure 3E). And these genes might be related to prognosis (Supplementary Figure S3C).

LncRNAs were shown to interact with transcriptional factors and thus affect transcription. They are also able to regulate the alternative splicing (AS) by interacting with splicing factors (SFs) or other RNA-binding proteins (RBPs) (Zhang et al., 2019a). Inspired by these previous discoveries, we compared the DEGs with reported RBP gene lists (Castello et al., 2012; Gerstberger et al., 2014; Castello et al., 2016; Hentze et al., 2018). We found that 72 RBP genes were upregulated and 40 RBP genes were downregulated in TP53TG1-OE cells (Figure 4A). Herein we propose that TP53TG1 could regulate alternative splicing by modulating the expression of RBPs. We used the recently published AS analysis software SUVA to identify AS events significantly altered between TP53TG1-OE and NC group (TP53TG1-RAS).

As shown in Figure 4B, by applying a stringent cut-off of p ≤ 0.05 and an AS ratio ≥0.2, 470 high-confidence RASEs were identified (Supplementary Table S8) (Cheng et al., 2021a). There are five RAS types: 249 alt3p, 156 alt5p, 13 contain, 14IR and 38 Olp. The splicing events were corresponding to classical splicing events, in which A3SS events accounted for a large proportion (Figure 4C), which may be one of the characteristics of TP53TG1-RAS. As shown in Supplementary Figure S4B, Alt3p and alt5P events, A3SS and A5SS events dominate the alternative splicing events. Further, about 60% of TP53TG1-RAS events were complex splicing events (Supplementary Figures S4C,D), indicating the complexity of AS regulation by TP53TG1. A splicing event involves two transcripts, and these two transcripts may only account for a very small part of the expression of the whole gene. We hope to find a more dominant transcript undergo AS which was quantified as “pSAR” value by SUVA. Here principal component analysis (PCA) was conducted by using the splicing ratio of different variable splicing events, and it could be seen that the principal components of the two groups could be well separated. And the samples’ correlation is strong (Supplementary Figure S4E). The number of splicing events accounting for different pSAR was shown (Figure 4C). Finally, 470 main splicing events with pSAR ≥ 50% were selected for further analysis (Figure 4D). As shown in Supplementary Figures S4C,D, the TP53TG1-OE group and the NC group can be clearly distinguished by the splicing ratio of TP53TG1-RAS (pSAR ≥ 50%).

Overall, genes linking with TP53TG1-RAS were primarily enriched in the proliferation, apoptosis (Supplementary Table S7) and signal transduction by P53 mediator pathways (Figure 4E). As shown in Figure 4F, the expression of differentially expressed RBPs and the splicing ratio of TP53TG1-RAS events associated with proliferation and apoptosis were used for correlation analysis (Pearson’s correlation coefficient ≥0.8, p ≤ 0.01). These RBPs might potentially regulate TP53TG1-RAS (Supplementary Figure S4F). The size of nodes in the figure represents the number of gene/splicing events associated with them, including two intron retention (IR), one A5SS, two ES and 1 CE which were located on BCLAF (Supplementary Figure S4G).

Here represents the number of gene/splicing events associated with them, including two intron retention (IR), one A5SS, two ES and 1 CE which were located on MVD, CNBP (Figures 5A,B).

CeRNAs have emerged as an important mechanism for lncRNA and miRNA regulatory network (Su et al., 2021). To investigate TP53TG1-mediated ceRNA network, miranda (Enright et al., 2003) and rnahybrid (Rehmsmeier et al., 2004) were used to predict the target relationship between miRNA and TP53TG1, respectively. Then, seven overlapping miRNAs, containing has-miR-6799-5p, has-miR-1273h-5p, has-miR-6779-5p, has-miR-6807-3p, has-miR-6510-3p, has-miR-6732-5p, and has-miR-1972, were screened. The seven miRNAs were mapped into miRDB (http://mirdb.org) and targetscan (http://www.targetscan.org) databases to explore their target genes, accompanied with TP53TG1-regulated DEGs. Finally, integrating the TP53TG1/miRNA interactions with the miRNA/DEGs interactions, a ceRNA network was constructed (Figure 6). The result reveals that lncRNA-TP53TG1 could adsorb miRNAs as a sponge RNA. The predicted ceRNA network comprised 426 nodes (13 miRNAs and 413 mRNAs) and 743 interactions. In this network (Figure 6), several nodes and interactions involving the feature genes play essential roles in CESC, such as hsa-miR-6779-5p (target genes: OAS3, HDAC2 and IFI6), hsa-miR-1972 (target genes: PINX2, OAS2, ZC3H12C), and hsa-miR-6807-3p (target genes: PDCD6IP, KIN, and EZH2), and all these target genes were upregulated. MiR-6779-5p was also identified as a prognosis-related miRNA (Yang et al., 2022).