This article is mentioned in:

Introduction

The Wnt/β-catenin signaling pathway serves important roles in the tumorigenesis of colorectal cancer (CRC). In colon cancer cell lines, the disruption of the adenomatous polyposis coli protein, caused by loss of heterozygosity or mutation, leads to deregulation of the β-catenin protein (1). Subsequently, it is transported into the nucleus and activates its target genes by recruiting cofactors of the transcription factor/lymphoid enhancer binding factor family (2). It is crucial to identify β-catenin target genes, since they are involved in cellular processes that contribute to proliferation and migration in colorectal carcinoma (3–8). Based on a literature review, the Nusse group (9) revealed 28 β-catenin target genes in human colon cancer, which are presented on their website (https://web.stanford.edu/group/nusselab/cgi-bin/wnt/target_genes). In an experimental approach, RNA interference (RNAi) of β-catenin was implemented by treating DLD1 and SW480 cells with small interfering RNA (siRNA) (10). By comparison with controls, in which samples were treated with mock siRNA, the identified differentially expressed genes are considered to be potential β-catenin target genes (10). The study then incorporated results from another study with a similar design, in which LS174T cells were treated with short hairpin RNA against β-catenin (11), and identified a total of 335 target genes. Recently, a multi-omics approach was used by Ewing et al (12) to decipher the oncogenic β-catenin network in HCT116 cells by comparing transcriptome, expression proteome and interactome data of wild type and β-catenin mutated samples. The results were subsequently integrated into a functional molecular network. However, target genes identified by these studies (10–12) lack functional information. Additionally, expression levels occasionally fail to reflect molecular dysfunction where post-translational modifications are involved (13).

Therefore, the present study aimed to identify β-catenin target genes by leveraging large-scale RNAi and CRISPR-CRISPR associated protein 9 (Cas9) genetic perturbation datasets, since they may provide an opportunity to derive gene-gene functional associations. Initially, in yeast, the gene-gene associations were revealed by epistatic analysis, where the phenotypic readout of a genetic perturbation depends on the status of a second gene (14). Based on this, genetic interaction networks could be constructed to reveal functionally associated genes (15). For human cell lines, the utility of this approach is limited due to the exponential increase in the combinatorial space (16,17). However, genome-wide single gene perturbation screens are more advanced and can be applied in a large collection of cell lines (18,19). For these screens, highly variable genetic dependencies for cellular fitness are observed in cancer cell lines, which may reflect diverse alterations during tumorigenesis (19). Upon perturbation of each gene, the genomic state of each cancer cell exhibits a unique overall fitness response. Notably, Project Acheilles achieved promising results by systematically elucidating genetic vulnerabilities across 501 cancer cell lines using RNAi (18) and in 342 cancer cell lines using CRISPR-Cas9 (19). The data could provide an opportunity to identify potential β-catenin target genes by deriving gene-gene functional associations.

The present study reanalyzed data generated from large-scale RNAi (18) and CRISPR-Cas9 (19) single gene perturbation screen studies and ranked gene candidates by their fitness correlation of β-catenin. The present study focused on genes that were highly correlated with β-catenin, and identified 13 β-catenin target genes. These genes were fitted into statistical models that could successfully predict CRC prognosis. Overall, the present study applied a novel strategy to identify β-catenin target genes, and the findings could establish the utility of this approach to reveal functional associations among Wnt/β-catenin signaling and deepen the understating of CRC pathogenesis.

Materials and methods

Filtering fitness data for analysis

Fitness screening data from RNAi (18) and CRISPR-Cas9 (19) experiments in CRC cell lines (Table SI) were downloaded from the DepMap website (depmap.org/portal/). Subsequently, genes were filtered based on the dependence scores upon genetic depletion across cell lines in both datasets. Therefore, genes that had a minimum dependency score >-0.3, for both the RNAi and CRISPR-Cas9 dataset, and were expressed in the cell lines with transcripts per kilobase million (TPM) >0, were selected for downstream analysis. Gene expression data (20) were retrieved from the Cancer Cell Line Encyclopedia (CCLE) project (portals.broadinstitute.org/ccle).

Computing the β-catenin fitness correlation landscape

Following filtering of the fitness data, the Pearson correlation coefficient (PCC) of fitness profiles of all other genes with β-catenin was computed to generate the fitness profile correlation landscape. The genes were ranked based on the PCC scores (between-1 and 1; Fig. 1A). It was hypothesized that potential targets could be among those genes exhibiting a positive or negative correlation with the fitness profile of β-catenin.

Figure 1. Overview of fitness profile correlations between β-catenin and selected genes in the present study. (A) Scatter plot of ranked genes based on the PCC. Red dots represent β-catenin targets identified in the present study. The ranked genes were significantly enriched in (B) proteasome and (C) L ribosomal proteins. L ribosomal proteins, large ribosomal proteins; PCC, Pearson correlation coefficient; RNAi, RNA interference; Cas9, CRISPR associated protein 9.

Microarray data

All microarray data used in the present study were retrieved from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (ncbi.nlm.nih.gov/geo/) (21) using GEO accession numbers. The siRNA β-catenin and control treatment microarray dataset accession numbers were GSE44097 (10) for the DLD1 and SW480 cell lines, and GSE18560 (12) for Ls174T cells. The CRC dataset accession numbers were GSE68468 (22), GSE14333 (23), GSE17536 (24), GSE17537 (25), GSE24549 (26), GSE24550 (27), GSE31595 (28), GSE37892 (29) and GSE39582 (30). For the GSE24549 and GSE24550 datasets, gene expression profiles were generated using an Affymetrix Human Exon 1.0 ST array, and GSE68468 gene expression profiles were generated using an Affymetrix Human Genome U133A array. The remaining data were generated using an Affymetrix Human Genome U133 Plus 2.0 array. The disease-free survival outcome information was retrieved from the Prediction of Clinical Outcomes for Genomics database (https://precog.stanford.edu/index.php) (31) using the aforementioned NCBI GEO accession numbers.

Raw intensity files (*.CEL) were downloaded from the GEO database and then processed using an in-house bioinformatics pipeline. Briefly, the raw files were loaded into R version 3.5.2 (32) environment using the Bioconductor package oligo (version 1.46.0) (33). The rma algorithms from the oligo package were applied for background correction and normalization. Gene annotation was processed based on the custom chip definition files (version 22.0.0) (34) downloaded from the BrainArray website (http://brainarray.mbni.med.umich.edu/Brainarray/Database/CustomCDF/genomic_curated_CDF.asp).

RNA-sequencing (RNA-seq) data

RNA-seq data were retrieved from the NCBI Sequence Read Archive database using the following accession numbers: SRP029880 (35) for CRC samples and SRP101345 (36) for HCT116 cell lines with either mutant catenin b1 (CTNNB1)/β-catenin allele disrupted or wild-type CTNNB1/β-catenin allele disrupted. The raw reads files (*.fastq) were downloaded and processed using Salmon software (version 0.13.0) (37) to quantify the expression levels of transcripts. Reference transcriptome data (GRCh38 release 94) were downloaded from the Ensembl database (http://www.ensembl.org/)

The cancer genome Atlas (TCGA)-colon adenocarcinoma (COAD) data

RNA-seq expression profiles of colon cancer (n=521) were retrieved from the TCGA-COAD dataset deposited in the NCI Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov) using the Bioconductor package, TCGAbiolinks (version 2.10.5) (38). The raw fragments per kilobase million values were converted to TPM (transcripts per million) values for downstream analysis. The differences between the expression levels between primary solid tumors and normal solid tissues were tested using a Wilcoxon rank sum test.

Differential gene expression analysis

The limma (version 3.38.3) (39) package was implemented to conduct the differential expression analysis. The main strategy of this package is to use linear models to evaluate differential expression levels in the context of a multifactor designed experiment. The P-values were adjusted using Benjamini and Hochberg correction (40).

Cox proportional hazard modeling

The Cox proportional hazards modeling of CRC disease-free survival and log-rank tests was conducted using the R package, survival (version 2.44–1.1; CRAN.R-project.org/package=survival). Using this statistical model, gene expression values were summarized into the unified scores. The median value of these scores was chosen as a threshold, samples with scores above which were defined as high-risk groups and low-risk groups were determined. The Kaplan-Meier plot was generated using the R package, survminer (version 0.4.4; CRAN.R-project.org/package=survminer).

Gene set enrichment analysis (GSEA)

GSEA was conducted using a Bioconductor package, fgsea (version 1.8.0) (41), with gene family annotations from the HUGO database (https://www.genenames.org). In the present study, GSEA determined whether members of the gene families tended to be enriched at the top (positively correlated) or bottom (negatively correlated) of the ranked gene list and measured the probability of such distribution to be generated randomly. The NES (normalized enrichment score) was adjusted from the ES (enrichment score), which reflected to what extent the gene sets from the gene families were overrepresented in the ranked gene list (42).

Results

ß-catenin fitness correlated genes enriched in CRC associated signaling pathways

As demonstrated in a previous study (43), the correlated fitness effects across cancer cell lines can serve as a measure to reveal gene-gene functional associations. In other words, genes that share similar fitness profiles may be functionally associated, which makes them potential regulating targets. Therefore, the present study aimed to identify β-catenin target genes by computing its fitness correlation with other genes perturbed in CRC cell lines (CRISPR-Cas9, n=27; RNAi, n=45). The present study applied the stringent computational pipeline in parallel and retrieved a total of 8,200 genes from the CRISPR-Cas9 dataset and 8,381 genes from the RNAi dataset. From these genes, the present study computed the PCC of fitness profiles with β-catenin. As presented in Fig. 1A, the genes were ranked from the most positive to the most negative correlation with β-catenin and they were significantly enriched in the gene groups of proteasome (Fig. 1B; CRISPR-Cas9, NES=−1.45, P=0.05; RNAi, NES=−2.25, P=0.0024) and large ribosomal proteins (Fig. 1C; CRISPR-Cas9, NES=1.54, P=0.016; RNAi, NES=−1.47, P=0.04). The revealed associations between β-catenin and the proteasome and ribosomal proteins may highlight candidate β-catenin targets.

Identified β-catenin target genes are involved in CRC tumorigenesis

The correlated genes were filtered using stringent criteria (P≤0.03; Pearson correlation test), and 210 genes in the CRISPR-Cas9 dataset and 754 genes in RNAi dataset were identified as potential β-catenin targets, since their fitness profiles were significantly correlated with β-catenin in CRC cell lines. Since the resulting fitness profile correlations exhibited similar trends in terms of gene set enrichment analysis, the overlapping significantly correlated genes from both datasets were selected, and 13 genes were identified as high confidence β-catenin targets (Fig. 2; Table I). Notably, these genes seldom exhibited significant expression alterations in colon cancer cell lines that were treated with siRNA that inhibited β-catenin, or harbored β-catenin mutations (Fig. 3).

Figure 2. Fitness profiles of β-catenin and 13 identified target genes in the CRISPR-Cas9 and RNAi datasets. Rows (genes) were hierarchically clustered. ABHD10, abhydrolase domain containing 10; AFG3L2, AFG3 like matrix AAA peptidase subunit 2; Cas9, CRISPR associated protein 9; COPS8, COP9 signalosome subunit 8; CTNNB1, catenin b1; DHX15, DEAH-box helicase 15; GMPS, guanine monophosphate synthase; IARS2, isoleucyl-tRNA synthetase 2; MED28, mediator complex subunit 28; MYB, MYB proto-oncogene, transcription factor; RNAi, RNA interference; RPL35, ribosomal protein L35; RRAD, RRAD, Ras related glycolysis inhibitor and calcium channel regulator; SRSF6, serine and arginine rich splicing factor 6; SRSF10, serine and arginine rich splicing factor 10; WWP2, WW domain containing E3 ubiquitin protein ligase 2.

Figure 3. Volcano plots revealing differentially expressed genes for siRNA-transfected vs. control samples and β-catenin mutant vs. wild type samples. β-catenin targets identified in the present study were highlighted. Data were obtained from the National Center for Biotechnology Information Gene Expression Omnibus database (accession numbers, GSE18560: top-left, GSE44097: top-right and bottom-left, and SRP101345: bottom-right). ABHD10, abhydrolase domain containing 10; AFG3L2, AFG3 like matrix AAA peptidase subunit 2; COPS8, COP9 signalosome subunit 8; CTNNB1, catenin b1; DHX15, DEAH-box helicase 15; GMPS, guanine monophosphate synthase; IARS2, isoleucyl-tRNA synthetase 2; MED28, mediator complex subunit 28; MYB, MYB proto-oncogene, transcription factor; RRAD, RRAD, Ras related glycolysis inhibitor and calcium channel regulator; siRNA, small interfering RNA; SRSF6, serine and arginine rich splicing factor 6; SRSF10, serine and arginine rich splicing factor 10; WWP2, WW domain containing E3 ubiquitin protein ligase 2.

Table I. β-catenin target genes identified in the present study.

Table I.

β-catenin target genes identified in the present study.

	CRISPR-Cas9		RNAi
Symbol	PCC	P-value	PCC	P-value	Chromosomal location	Description
ABHD10	0.430	0.025	−0.477	0.016	3q13.2	Abhydrolase domain containing 10
AFG3L2	−0.427	0.026	0.371	0.012	18p11.21	AFG3 like matrix AAA peptidase subunit 2
COPS8	−0.465	0.014	0.336	0.024	2q37.3	COP9 signalosome subunit 8
DHX15	−0.587	0.001	−0.440	0.003	4p15.2	DEAH-box helicase 15
GMPS	0.496	0.009	0.349	0.019	3q25.31	Guanine monophosphate synthase
IARS2	−0.421	0.029	−0.341	0.022	1q41	Isoleucyl-tRNA synthetase 2, mitochondrial
MED28	−0.473	0.013	−0.421	0.004	4p15.32	Mediator complex subunit 28
MYBa	0.484	0.011	0.462	0.001	6q23.3	MYB proto-oncogene, transcription factor
RPL35	0.455	0.017	0.559	<0.001	9q33.3	Ribosomal protein L35
RRAD	0.559	0.002	0.487	0.001	16q22.1	RRAD, Ras related glycolysis inhibitor and calcium channel regulator
SRSF10	−0.537	0.004	−0.383	0.010	1p36.11	Serine and arginine rich splicing factor 10
SRSF6	0.420	0.029	−0.366	0.013	20q13.11	Serine and arginine rich splicing factor 6
WWP2	0.423	0.028	0.373	0.013	16q22.1	WW domain containing E3 ubiquitin protein ligase 2

{ label (or @symbol) needed for fn[@id='tfn1-ol-0-0-10724'] } P-values indicate the significance level of the PCC value.

a Previously reported as known β-catenin target (51). Cas9, CRISPR associated protein 9; PCC, Pearson correlation coefficient; RNAi, RNA interference.

Additionally, the present study explored the expression levels of these genes in primary CRC and normal tissue samples from three independent cohorts, including a total of 805 samples. As presented in Fig. 4, all 13 genes were differentially expressed in solid tumor tissues from cohort 1. Among the genes, ten were upregulated and three were downregulated in CRC tumor samples. Although similar trends were observed in the other cohorts, certain genes (ABHD10, IARS2, MYB, RRAD, WWP2 in cohort2; IARS2 and WWP2 in cohor3), were not significantly differentially expressed between tumor and normal samples.

Figure 4. Boxplots demonstrating gene expression levels of β-catenin targets in CRC primary solid tumors vs. normal solid tissues. Analysis was performed using the Wilcoxon rank sum test. *P<0.05, **P<0.01, ***P<0.001, Normal solid tissue compared with respective Primary solid tumor. ABHD10, abhydrolase domain containing 10; AFG3L2, AFG3 like matrix AAA peptidase subunit 2; COPS8, COP9 signalosome subunit 8; CTNNB1, catenin b1; DHX15, DEAH-box helicase 15; GMPS, guanine monophosphate synthase; IARS2, isoleucyl-tRNA synthetase 2; MED28, mediator complex subunit 28; MYB, MYB proto-oncogene, transcription factor; RNA-seq, RNA sequencing; RPL35, ribosomal protein L35; RRAD, RRAD, Ras related glycolysis inhibitor and calcium channel regulator; SRSF6, serine and arginine rich splicing factor 6; SRSF10, serine and arginine rich splicing factor 10; TCGA-COAD, The Cancer Genome Atlas-colon adenocarcinoma; WWP2, WW domain containing E3 ubiquitin protein ligase 2.

Modeling target genes as a CRC prognostic biomarker

There is increasing interest to discover CRC prognostic biomarkers as patients could benefit from earlier diagnosis and a personalized treatment strategy. For this purpose, the present study aimed to establish gene signature prognostic models that could predict disease-free survival outcomes in CRC. The expression profiles of the 13 genes were fitted into a Cox proportional hazards model and a unified score was computed for each sample. For each cohort, the median value of the aforementioned unified score was selected to divide high- and low-risk groups. Log-rank tests indicated that the high- and low-risk groups exhibited significant differences in disease-free survival prognosis except for GSE31595 (Fig. 5).

Figure 5. Kaplan-Meier plots demonstrating disease-free survival analysis of prognostic models based on expression profiles of β-catenin targets. P-values were calculated using the log-rank test.

Discussion

The present study explored large-scale RNAi and CRISPR-Cas9 screening data and identified β-catenin target genes based on fitness profile similarities generated from a large-scale study, CCLE (20) and the DepMap (18,19) project. As illustrated in a previous study, this fitness correlation strategy can recapitulate gene functional modularity by rewiring human protein complexes and establishing a human functional similarity network (43). A total of 79 CRC cell lines were considered in the present study. This set the CRC context, in which β-catenin regulates target genes. For both RNAi and CRISPR-cas9 datasets, the genes where the fitness data were highly correlated with β-catenin were significantly enriched in the proteasome and ribosomal families. Lack of proteasomal degradation of β-catenin facilitates its entry to the nucleus and targeting genes, including cyclin D1 (4) and AKT1 (44). This eventually leads to the proliferation and differentiation of cells. Furthermore, disturbances in ribosomal proteins have been observed in a variety of cancerous tissues, including glioblastoma (45), breast (46), esophagus (47), liver (48) and cervix tissues (49). For CRC, similar trends have been observed. For example, the expression levels of several ribosomal proteins are abnormally regulated in primary (50) and metastatic (51) CRC. It is worth noting that, in the previous strategy, the identified β-catenin targets from different cell lines merely overlapped (10). This limitation was overcome by computing gene fitness correlation and thus yielded reliable results.

Notably, among the 13 identified targets, the MYB proto-oncogene, transcription factor (MYB), has previously been reported as a known β-catenin target (52). Upregulated MYB levels and activated β-catenin may induce robust upregulation of MYC promoter activity in CRC (53). Previous studies have suggested that the dysregulation of MYB is associated with several rare types of cancer, including adenoid cystic cancer (54–60), spiradenocarcinoma (61) and cutaneous cylindroma (62). For the mediator complex subunit 28 gene, the same trend in expression alterations has been observed in epithelial-derived types of cancer (63). Although, to the best of our knowledge, there are no reports indicating the involvement of COP9 signalosome subunit 8 (COPS8) in CRC, a proteomic study revealed that COPS8 is upregulated in prostate cancer (64). The serine and arginine rich splicing factor 6 and serine and arginine rich splicing factor 10 genes are members of the serine-arginine family, which regulates RNA splicing. As a result, their roles in regulating alternative splicing may promote cancer pathogenesis (65–67). WW domain containing E3 ubiquitin protein ligase 2 modulates transforming growth factor b-dependent transcription and epithelial-mesenchymal transition (68,69). It should be noted that some of the identified targets seldom exhibited altered expressions in the siRNA β-catenin approaches. This suggests that the method used in the present study may serve as a complementary alternative for the identification of β-catenin targets, which may be potentially missed by other strategies.

The present study further revealed the associations of the 13 identified targets with CRC disease free survival outcomes in eight independent cohorts encompassing 1,310 individuals. Indeed, the prognostic potential of these genes has been revealed in previous studies (70,71). For instance, loss of the AFG3 like matrix AAA peptidase subunit 2 gene, located in the 18p11.32–21 region, have been associated with a significantly longer progression-free survival in patients with CRC (70). Additionally, MYB is associated with metastasis in pancreatic tumors (71).

However, it worth mentioning that the present study used a stringent criterion to filter the data and reported β-catenin targets with high confidence. While there are two datasets available, CRISPR-Cas9 and RNAi, the present study selected the overlap of both datasets for downstream analysis. This may lead to the misinterpretation of the remaining potential targets. For example, survivin, cyclin-D1 and axin-2 are well established β-catenin targets. The present study could not recapture them, since they exhibited less significant correlations with β-catenin comparing the identified 13 targets in this study. It was anticipated that a more robust strategy, such as machine learning, should be employed to explore the data. Additionally, the prognostic model requires improvements by statistical modeling, so that these findings can be applied to clinical practice using less complicated assays. Finally, experimental studies are required to verify the findings. The identified β-catenin target genes are of high confidence and the pathways are associated with CRC pathogenesis, which provides resources for the research community.

Supplementary Material

Supporting Data

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets generated and/or analyzed during the present study are available in the figshare repository, doi.org/10.6084/m9.figshare.8872769.

Authors' contributions

HZ and BS conceived of and designed the study. HZ and LH performed the gene fitness correlation, gene set enrichment, differential gene experssion and prognostic modelling analysis. DY collected and re-formatted the data. HZ and BS wrote the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References