We obtained The Cancer Genome Atlas (TCGA) pan-cancer somatic mutation data (MC3 MAF v0.2.8 file) from Genomic Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). The previously reported multiple filtering steps were used to eliminate artifacts and reduce false-positive rates (Thorsson et al., 2018). To include one of wga, native_wga_mix, and PASS, only the non-silent mutations were retained (the values in Variant_Classification column should be one of Frame_Shift_Del, Frame_Shift_Ins, In_Frame_Del, In_Frame_Ins, Missense_Mutation, Nonsense_Mutation, Nonstop_Mutation, Splice_Site, and Translation_Start_Site). Mutations calls were required to be called by two or more mutations callers (the values in NCALLERS column >1). The mutation hotpots were visualized using the MutationMapper tool in cBioPortal (https://www.cbioportal.org/mutation_mapper).

We obtained TCGA thresholded somatic copy-number alteration scores (SCNAs) (ISAR-corrected GISTIC2.0 all_thresholded.by_genes file) for 9,125 patient samples from Genome Data Commons. Genes assigned a positive value of +2 were considered as high-level amplification events, and genes with a negative value of −2 were considered as deep deletion events, as previously described (Stewart et al., 2018; Lim and Cech, 2021). OncoPrint plots of mutation and SCNA were generated by using the ComplexHeatmap R package.

The mRNA expression data from 10,304 samples of tumor and 719 samples of normal tissues across 33 cancer types were obtained from the TCGA database using the University of California Santa Cruz cancer genome browser (http://xena.ucsc.edu/). We collected transcriptome and clinical data and survival information. In this study, the log2-transformed RSEM was used for all mRNA analyses. We obtained 20 cancer types with a sufficient number of tumor-normal matched pairs (n ≥ 3). To assess the gene expression pattern of the CST complex, we used the Wilcox test to contrast the differences in CST gene expression levels between tumor and matched normal tissues and considered p value <0.05 as statistically significant.

To measure overall CST activity and CS activity in cancer, we derived a CST score and a CS score based on CST expression. First, the expression levels of CST genes were normalized within each cancer type according to previous publications (Korkut et al., 2018; Wang et al., 2018). Then, the CST score and CS score were obtained by calculating the first Z, normalizing the values for the CST genes. The mean across the genes was calculated for each sample to yield the CST score and CS score per sample.

Telomerase activity was estimated from the expression of a 13-gene signature, named EXTEND score (Expression-based Telomerase ENzymatic activity Detection score) (Noureen et al., 2021). Cell stemness, an indicator of cell proliferation, was represented as the mRNAsi score calculated by using an innovative one-class logistic regression machine-learning algorithm (OCLR) (Malta et al., 2018). A previous pan-cancer study provided a list of genomic instability scores, including the mutation burden score, the homologous recombination deficiency (HRD) score, the SCNA burden score, the aneuploidy score, and the loss of heterozygosity (LOH) score (Knijnenburg et al., 2018). We assessed the level of the genetic instability score among the three expression clusters using the Kruskal test.

The hallmark gene sets associated with the CST complex were identified by using the GSEA R package. We calculated the Spearman correlation coefficient between mRNA expression of all protein-coding genes and the CST complex. All the genes were ranked according to their correlation coefficients, and the pre-ranked gene lists were run against the MSigDB’s Hallmark Gene Set (HGS). All gene sets were annotated by their functional categories.

The patients’ clinical data were retrieved from Genome Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). Survival analysis included 10,223 patients across 33 cancer types with complete transcriptome data and survival information. We assessed the correlation between the gene expression of the CST complex and the patient’s overall survival times. Kaplan–Meier survival curves of overall survival and the univariate Cox proportional hazards model were generated using the Survival and Survminer R packages. We used 1) The log-rank test to compare patient survival curves among the three clusters; 2) Within each cancer type, we determined the optimal cutpoint for the CTC1-STN1 score (CS score) and TEN1 expression using the maximally selected rank statistics from the maxstat R package. Then, the statistical significance of survival differences between the low and high groups was determined using the log-rank test; and 3) the univariate Cox proportional hazards model to assess the correlation between CST expression and the patient’s overall survival time.

Gene expression data and overall survival data were obtained from four independent cancer datasets, including metastatic melanoma (GSE78220, GSE91061, and GSE100797) and metastatic urothelial cancer (IMvigor). The CS score of each sample was calculated in each dataset. Kaplan–Meier survival curves of overall survival were generated using the Survival and Survminer R packages. In metastatic melanoma, the correlation between CS score and PD-L1 expression in each dataset was calculated by linear regression in R. In the IMvigor dataset, we analyzed the correlation between CS score and TGFBR2 expression.

The robust estimation data of immune infiltration levels for TCGA using six algorithms (CIBERSORT, EPIC, MCPCOUNTER, QUANTISEQ, TIMER, and XCELL) were downloaded from the Tumor Immune Estimation Resource database (TIMER 2.0, http://timer.cistrome.org/). We analyzed the correlation between CS score and immune cell infiltration levels within each cancer type. The correlations of CS score with immune infiltration were evaluated by tumor purity-adjusted partial Spearman’s correlation and statistical significance.

TIDE (tumor immune dysfunction and exclusion) is a computational method to model two primary mechanisms of tumor immune evasion: T cell dysfunction and T cell exclusion. The TIDE score can predict immune checkpoint blockade (ICB) clinical response based on pre-treatment tumor profiles. The TIDE score of each sample in the TCGA cohort was analyzed by using the web application (http://tide.dfci.harvard.edu). We further calculated Spearman correlation coefficients between CST score, CS score, Ten1 expression, and TIDE score in each cancer type.

TCGA DNA methylation data files were downloaded from Genomic Data Commons. We mapped the Illumina methylation array probes to individual genes using the Illumina Human Methylation 450 k R annotation data package. We retained those probed mapped to the promoter region. For genes with multiple probes, mean beta values were used. To estimate the overall methylation level, we then calculated the median beta value for CST genes in each sample. To examine the regulation of CST expression by DNA methylation, we calculated the Spearman correlation between the DNA methylation beta value and mRNA expression for each gene.

A recent pan-cancer study reconstructed lncRNA regulatory networks (lncNETs) using the molecular tumors of TCGA (Feng et al., 2017). We retrieved a list of lncRNA that is associated with CST expression. Next, we retained those lncRNAs which had a Spearman correlation coefficient < −0.15 and p < 0.05 with CTC1 and STN1. The gene sets associated with that lncRNA were identified by using the GSEA R package.

The batched-corrected, normalized miRNA expression data (pancanMiRs_EBadjOnProtocolPlatformWithoutRepsWithUnCorrectMiRs_08_04_16. xena) were downloaded from the TCGA database using the University of California Santa Cruz cancer genome browser (http://xena.ucsc.edu/). We calculated Spearman correlation coefficients between expression for miRNA mature strands and mRNA expression levels of the three CST genes in individual cancer types. We filtered the result by requiring correlations to have a coefficient < −0.2 and p > 0.05. Then, we filtered by retaining miRNA with experimentally validated miRNA–mRNA interactions (from TargetScan 7.2 database and miRDB database).

Connectivity Map (CMap, https://clue.io/) is a chemical genomics database that collects gene-expression profiles from cultured human cells treated with small molecules (Huang et al., 2012). In each cancer type, there is a list of top 150 positively or negatively correlated genes associated with CST score, CS score, and TEN1 expression level. Compounds with enrichment score >90 or < −90 in at least 17 cancer types (>1/2 of all cancer types) are significantly positive or negative compounds, respectively.

The CST database is a resource tool open to all users. It is constructed in the standard MVC (Model-View-Controller) pattern, which consists of a server side and a client side. The outputs of the database are plots and tables. All plots are created by R (version R ×64 4.0.3).

Heatmaps for gene expression data and GSEA results were depicted using the pheatmap R package. The OncoPrint plots of mutation and SCNA were generated using the ComplexHeatmap R package. The differences among the three clusters were plotted as a boxplot by using the ggplot2 R package. Cytoscape software is a free open-source platform that enables biological network analysis and two-dimensional visualization. The connective network among CTC1, STN1, lncRNA, and pathways was visualized using Cytoscape (version 3.7.2). The Kaplan–Meier survival curves of overall survival were plotted using the survminer and survival R packages. Cox regression was performed using the survival R package and visualized using the ggplot2 R package. The miRNA-CST regulation and lncRNA-CST regulation were plotted as heatmaps by using the pheatmap R package.

CTC1 and TEN1 conditional knockout HCT116 cells were cultured in McCoy’s 5A Medium (modified) supplemented with 10% FBS, with antibiotics as previously described (Feng et al., 2017; Feng et al., 2018). Tamoxifen (Sigma, H7904) was added to 10 nM for 7 days to induce Cre expression to disrupt the CTC1 or TEN1 gene.

The total RNA was extracted from cells using the Trizol Reagent (Invitrogen) according to the manufacturer’s protocol. To generate cDNA, total RNA by reverse transcription using the PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa) was determined. cDNA was analyzed by qPCR with QuantStudio (Thermo Fisher) using PowerUp™ SYBR™ Green Master Mix (Thermo Fisher). The relative expression levels of the target genes were determined by the 2^−ΔΔCt method.

CTC1 forward 5′−CTA​GAC​TCC​TAT​CCT​ATT​CGG​GA−3′ and reverse 5′−ACA​CAC​TCC​AGG​TCG​CAT​C−3′, TEN1 forward 5′−GGC​CAA​GTT​CCT​GAT​GGG−3′ and reverse 5′− CAG​TGT​TAC​TCT​GGA​CTG​AAT​CAT−3′, C-MYC forward 5′− CGT​CTC​CAC​ACA​TCA​GCA​CAA−3′ and reverse 5′− CAC​TGT​CCA​ACT​TGA​CCC​TCT​TG−3′.

To investigate the genomic characteristics of the CST complex in cancers, we analyzed the mutation frequency and SCNA of patients across 33 cancer types in the pan-cancer cohort. Genetic alterations include gene mutations (truncating and missense) and SCNA (amplification and deep deletion). The proportion of CST alterations was relatively low; about 4.79% of tumors had at least one type of the alterations. The OncoPrint shows a panoramic view of CST genomic aberrations (Figure 1B). The overall genetic alteration of CST genes ranged from 1 to 2%. Within the CST complex, the predominant alteration pattern varied for each gene. The alterations of CTC1 and STN1 were dominated by mutations and deletions, while TEN1 was more common by amplifications (Figures 1B,C). The lollipop plots show the mutation patterns and some potential recurrent hotspot mutations in CST (Supplementary Figure S1A). We found two significant mutation sites, which were missense mutations of STN1 (H317Y) and TEN1 (R119P/Q/W) (Supplementary Figure S1A). Examining the expression level of the corresponding genes among different alteration groups, we found that the amplification group displayed the highest expression, while the deletion group showed the lowest expression in all the three genes (Supplementary Figure S1B).

When examining the genetic alteration status of the CST complex across 33 cancer types, we found that various cancer types have different patterns of genetic alterations. A high-amplification frequency of TEN1 was observed in breast invasive cancer (BRCA), liver hepatocellular carcinoma (LIHC), mesothelioma (MESO), and uterine carcinosarcoma (UCS), while a relatively high deletion rate of CTC1 and STN1 was found in colon adenocarcinoma (COAD), lymphoid neoplasm diffuse large B-cell lymphoma (DLBC), and LIHC and prostate adenocarcinoma (PRAD) (Figure 1C). CTC1 mutation predominately occurred in uterine corpus endometrial carcinoma (UCEC) (6.59%), skin cutaneous melanoma (SKCM) (4.91%), and COAD (3.69%) (Figure 1C). Some cancer types were found with both alterations of TEN1 amplification/CTC1 deletion (LIHC) or both of CTC1/STN1 deletions (PRAD) or all of the TEN1 amplifications and CTC1/STN1 deletions (BRCA and UCEC). To investigate whether the two/three alternations simultaneously existed in these cancers, we did OncoPrint for each cancer type. However, we found few overlaps of different alternation types in one patient (Figure 1D). Particularly, in BRCA, TEN1 amplification was the most frequent type of genetic alteration. In LIHC, the proportion of CTC1 deletion or mutation and TEN1 amplification was 4%. The alteration pattern of the CST complex in PRAD is mainly the deletion of CTC1 and STN1, but CTC1 mutations are more common in UCEC (Figure 1D). These data illustrate a diverse alteration pattern for the CST complex, with mutation and deletion for CTC1 and STN1, and amplification for TEN1 as major events. The major events are not simultaneous in any cancer type, indicating that the individual component of the CST complex has different expression profiles and functions in cancer.

When comparing the expression of CST between tumor and paired adjacent normal samples, we found that CTC1 and STN1 were largely downregulated, whereas TEN1 was upregulated in tumors (Figure 2A). CTC1 and STN1 presented consistently lower expression in tumor samples for the majority of cancers, including bladder urothelial carcinoma (BLCA), BRCA, COAD, kidney chromophobe (KICH), kidney renal clear cell carcinoma (KIRC), lung adenocarcinoma (LUAD), and lung squamous cell carcinoma (LUSC) (p < 0.05). In contrast, TEN1 was upregulated in LUSC, KIRC, PRAD, thyroid carcinoma (THCA), and UCEC (p < 0.05) (Figure 2A).

To gain a comprehensive view of CST expression, unsupervised consensus clustering was performed on patient samples across 33 cancer types based on the expression of the CST complex using the K-means clustering method. Three distinct clusters emerged from cluster analysis, with cluster A displaying the highest expression of CTC1 and STN1 and low expression of TEN1, while cluster B had the most abundant expression of TEN1 (Figure 2B). Analyzing the distribution of samples, we found that acute myeloid leukemia (LAML) (99%), THCA (93%), and thymoma (THYM) (94%) were mainly located in cluster A, whereas glioblastoma multiforme (GBM), brain lower grade glioma (LGG), LIHC, and UCS were nearly depleted in this cluster (the proportion of samples was less than 5%) (Figure 2C). Adrenocortical carcinoma (ACC) (80%), BLCA (75%), and uveal melanoma (UVM) (83%) were largely located in cluster B. GBM (85%) and LGG (96%) are common central nervous system tumors, predominately resided in cluster C (Figure 2C). The samples of some cancer types were almost equally distributed among the three clusters, including DLBC, esophageal carcinoma (ESCA), stomach adenocarcinoma (STAD), and UCEC (Figure 2C). These data indicate that there is a dysregulated expression profile of the CST complex in a cancer-type-dependent manner.

The CST complex has been demonstrated to regulate telomerase, cell growth, and genome integrity. To assess the functionality of the expression pattern, we examined a group of scores measuring telomerase activity (represented as EXTEND score), cell stemness (an indicator for cancer progression), and genome instability (scores for aneuploidy score, SCNA burden, HDR, LOH, and mutation burden). We first found that three clusters showed different survival rates (Figure 2D). Cluster A, which had the highest expression of CTC1 and STN1, showed the best survival (Figure 2D). Additionally, cluster A had the lowest levels of scores for telomerase activity, stemness, and genome instability (Figures 2E,F). Collectively, these data not only confirm that the CST complex is a terminator of telomerase and that it maintains genome integrity but also implicate that CTC1 and STN1 can still execute these functions when TEN1 is low.

Having established an initial link between the CST complex and survival, we further investigated this link in greater detail. To measure the activity of the CST complex, we derived the CST score (based on mRNA expression of three CST genes) for each cancer type. We also derived the CS score (based on mRNA expression of CTC1 and STN1) as CTC1-STN1 had distinct expression profiles when compared to TEN1.

To obtain the predictive value of the CS score and TEN1 expression in overall survival, we performed Kaplan–Meier survival curve analysis by grouping the patients into the high and low groups based on optimal cutoff values of the CS score and TEN1, respectively. The survival analyses reveal two cancer groups in which both CS score and TEN1 are significantly associated with survival (Figures 3A,B). The first cancer group contained nine cancer types showing a consistent survival pattern for CS score and TEN1 (Figure 3A). In contrast, the second cancer group had six cancer types exhibiting an opposed survival pattern for CS score and TEN1 (Figure 3B). Among 15 cancer types, the CS score predicted better survival in 13 cancer types including ACC, BLCA, cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), ESCA, head and neck squamous cell carcinoma (HNSC), kidney renal papillary cell carcinoma (KIRP), LUAD, pancreatic adenocarcinoma (PAAD), KICH, KIRC, LGG, SKCM, and THYM. However, TEN1 was associated with better survival in nine cancer types and with worse survival in six cancer types (Figures 3A,B). To confirm these observations, we thus performed Cox regression analysis. Cox regression revealed that CTC1, STN1, the CS score, and the CST score were associated with better survival in at least four cancer types, while TEN1 was associated with worse survival only in LGG (Figure 3C). Taking all these data together with data in Figure 2A, we conclude that CTC1 and STN1, but not TEN1, may prevent tumor progression and protective factors for survival in many cancer types.

To understand the underlying mechanisms by which the CST complex regulates survival, cell proliferation, and genome integrity, we performed gene set enrichment analysis using the CST score, CS score, and TEN1 expression. This analysis identified several significant gene sets (Supplementary Figures S2A,B). Among them, the CST and CS score were negative, whereas TEN1 was positively associated with DNA repair, upregulated genes in response to UV (UV_RESPONSE_UP), and MYC targets (Figures 4A,B). Experimentally, we confirmed that knockout of CTC1, but not TEN1, increased the c-MYC mRNA level (Figure 4C). Altogether, these data indicated that CTC1 may be a negative regulator of c-MYC.

We noticed that CST had both positive and negative correlations with several immune-related gene sets dependent on the cancer type (Supplementary Figures S2). To further study how CST connects with cancer immunity, we calculated the correlation between CST and tumor immune dysfunction and exclusion (TIDE) score. TIDE score is a measurement of tumor immune escape, reflecting T cell dysfunction and exclusion and predicting cancer immunotherapy response (Jiang et al., 2018). We observed that CTC1 and STN1 (i.e., the CS score) were negatively associated with the TIDE score, suggesting that the CS score predicted a better response to immunotherapy (Figure 5A). Using six immune cell estimation methods, we found that the CS score was positively correlated with CD8 T cells and B cells in several TCGA cancer types (Figure 5B). CD8 T cells and B cells are important biomarkers for immune checkpoint blockade (Riaz et al., 2017; Helmink et al., 2020). The abundance of CD8 T cells and B cells was positively correlated with a favorable response to immune checkpoint blockade. Thus, all these data indicate that CTC1-STN1 might serve as a potential predictive biomarker for immune checkpoint blockade.

To confirm the predictive value of CTC1-STN1, we curated three datasets of expression data for melanoma patients treated with anti-PD-1, anti-PD-L1, or anti-CTLA4 and one dataset of expression data for metastatic urothelial cancer patients treated with anti-PD-L1 (Hugo et al., 2016; Riaz et al., 2017; Lauss et al., 2017; Mariathasan et al., 2018). For all these datasets, the CS score predicted a better response. Patients with a high CS score showed better survival (Figure 5C,D). Supporting this observation, we also found that CS score was positively correlated with PD-L1 in melanoma and with TGF-beta receptor TGFBR2 in metastatic urothelial cancer (Figures 5E,F). PD-L1 and TGFBR2 had been shown to be higher in the good responder group (Patel and Kurzrock, 2015; Mariathasan et al., 2018). Taken all these data, we conclude that CTC1 and STN1 predict the better response to immune therapy.

To investigate mechanisms that were responsible for dysregulated expression of the CST complex, we focused on DNA methylation, miRNAs, and lncRNAs as they are key regulators of gene expression. Surprisingly, Spearman correlation analyses between CST expression and the corresponding methylation level showed that the DNA methylation level of TEN1 was negatively correlated with its mRNA expression, while those of CTC1 and STN1 were positively correlated (Figure 6A). This positive correlation data suggest that DNA methylation was not responsible for the downregulation of CTC1 and STN1 expression in tumor samples.

Therefore, we tested whether non-coding RNAs could be an important regulator. Interestingly, by mining the miRNA-target interaction database and analyzing the correlation between miRNA and CST expression, we found that CTC1 and STN1, but not TEN1, could be repressed by multiple miRNAs in various cancer types (Figure 6B). In addition, we observed that some lncRNAs negatively associated with CTC1 and STN1 expression. Remarkably, those lncRNAs positively associated with those pathways had negative connections with CTC1 and STN1 (Figure 6C). For example, the negatively correlated pathway of CTC1/STN1 and MYC_TARGETS_V1 (Figure 4B) was positively correlated with CTC1/STN1-associated lncRNAs including CSAG2, CSAG3, H19, DSCR9, and PART1 (Figure 6C), suggesting that lncRNAs might enhance MYC expression by downregulating CTC1 and STN1. Collectively, these findings suggest that non-coding RNAs might be responsible for the suppression of CTC1 and STN1 in cancer.

CST has been reported to promote PARP inhibitor sensitivity in BRCA1-deficient cancer cells, highlighting that it can be a potential therapeutic target. Thus, we searched for chemical compounds that may modulate CST expression by connectivity map analysis. CS score was used since CTC1 and STN1 were co-expressed in most cancers. By setting a stringent cutoff as having a correlation with CS or TEN1 in at least 10 from the 33 cancer types, we discovered both positively (in red) and negatively (in blue) correlated compounds (Figures 7A,B). Several cell growth inhibitors, such as inhibitors for AKT, EGFR, mTOR, and PI3K, were positively correlated with CS score, implying that inhibition of cell growth might induce CTC1 and STN1 expression (Figure 7A). In contrast, chemical compounds affecting the metabolism were potential compounds regulating CTC1 and STN1 negatively. These compounds include the nitric oxide synthase inhibitor, bile acid, and pyruvate dehydrogenase kinase inhibitor (Figure 7A).

To make the data available and benefit the research community, we have built a user-friendly, open-access database named CST database (http://bioinfo-sysu.com/CSTDatabase/) (Figure 8A). Users can search and download the multi-omic data (Figures 8B–D). These data repositories of the secondary analysis are allowed to upload and release publicly.