CSNK1E/CTNNB1 are synthetic lethal to TP53 in colorectal cancer and are markers for prognosis.

Two genes are called synthetic lethal (SL) if their simultaneous mutations lead to cell death, but each individual mutation does not. Targeting SL partners of mutated cancer genes can kill cancer cells specifically, but leave normal cells intact. We present an integrated approach to uncovering SL pairs in colorectal cancer (CRC). Screening verified SL pairs using microarray gene expression data of cancerous and normal tissues, we first identified potential functionally relevant (simultaneously differentially expressed) gene pairs. From the top-ranked pairs, ~20 genes were chosen for immunohistochemistry (IHC) staining in 171 CRC patients. To find novel SL pairs, all 169 combined pairs from the individual IHC were synergistically correlated to five clinicopathological features, e.g. overall survival. Of the 11 predicted SL pairs, MSH2-POLB and CSNK1E-MYC were consistent with literature, and we validated the top two pairs, CSNK1E-TP53 and CTNNB1-TP53 using RNAi knockdown and small molecule inhibitors of CSNK1E in isogenic HCT-116 and RKO cells. Furthermore, synthetic lethality of CSNK1E and TP53 was verified in mouse model. Importantly, multivariate analysis revealed that CSNK1E-P53, CTNNB1-P53, MSH2-RB1, and BRCA1-WNT5A were independent prognosis markers from stage, with CSNK1E-P53 applicable to early-stage and the remaining three throughout all stages. Our findings suggest that CSNK1E is a promising target for TP53-mutant CRC patients which constitute ~40% to 50% of patients, while to date safety regarding inhibition of TP53 is controversial. Thus the integrated approach is useful in finding novel SL pairs for cancer therapeutics, and it is readily accessible and applicable to other cancers.

Introduction

Colorectal cancer (CRC) is the third most common cancer worldwide with nearly 1 million annual incidences and 400,000 yearly deaths [1]. Multiple genes are mutated during tumorigenesis, with those contributing to tumor formation and growth being called cancer genes [2], which can be categorized into oncogenes, tumor suppressor genes and stability genes. Mutant oncogene and tumor suppressor gene frequently implicated in colorectal cancer are CTNNB1 and TP53, respectively, and are involved in cancer cell proliferation, while stability genes such as MSH2 are responsible for the safeguarding of genomic integrity. Although to date some therapeutics directed against oncogenes have led to increases in patient survival, many fail due to intrinsic or adaptive resistance of cancer cells to the therapeutics. For example, when CRC patients are treated with cetuximab (EGFR inhibitor), resistance arise due to crosstalk between Src family kinases and EGFR [3].

Two genes are called synthetic lethal (a type of genetic interaction [4]) when their simultaneous mutations lead to cell death, while each single mutation does not. Genome-scale mappings of SLs in S. cerevisiae were obtained through high-throughput synthetic genetic array analyses [5], [6], [7]. Wong and colleagues successfully predicted synthetic sick or lethal (SSL) interactions in S. cerevisiae by integrating multiple types of data, e.g., gene expression, protein-protein interaction and properties of network topology of gene triples [8]. About 80% of SSL interactions therein were discovered successfully via validating < 20% of the predicted interactions. In C. elegans, Zhong and Steinberg [9] computationally integrated interactome, gene expression and phenotype data to predict genome-wide genetic interactions; they further experimentally verified the predictions for two human disease-associated genes.

The concept of synthetic lethality can be applied to exploit cancer-cell specific mutations for cancer therapeutics [10]. Targeting synthetic lethal (SL) partners of mutated cancer genes will specifically kill cancer cells bearing the mutations but spare normal cells. Therefore, synthetic lethality strategy offers an elegant alternative for non-druggable mutant tumor suppressor genes and difficult-to-target-directly oncogenes, such as TP53 and KRAS, respectively, by targeting their SL partners. Lately, pioneering studies of SL partners in BRCA1 and BRCA2-deficient cancer cells identified PARP1 as a promising drug target [11], [12]. The clinical relevance of synthetic lethality has been rapidly recognized. Multiple phase II and III clinical trials of PARP inhibitors were conducted for breast and ovarian cancer patients with BRCA1 or BRCA2 mutation [13], [14].Over the last few years, genes having SL interactions (SLs) in cancers have been actively studied using individual RNAi experiments, or by large scale RNAi screenings that uncover multiple SL gene pairs [15], [16], [17], which in general were centered on one gene. Astsaturov and colleagues identified SLs in human by combining both computational and experimental approaches [18].They first combined pathway maps, protein-protein interactions, gene expression data and human orthologs of Drosophila Egfr genetic interaction partners, to predict 2689 SL candidates of EGFR. They then selected 683 candidates for RNAi screening, by their appearance in at least two of these information sources or by prior biological knowledge, which resulted in 61 SLs. Wang and Simon [19] developed a promising computational approach to prefilering p53 SL genes before RNAi screening. They analyzed five gene expression datasets from NCI-60 panel, The Cancer Genome Atlas and Cancer Cell Line Encyclopedia. Kinase-encoding genes highly expressed in functional p53-mutant samples versus samples without functional p53 mutation were first identified, and were further prioritized by inferred significant networks, functional annotation and drug sensitivity analysis.

To augment these RNAi-based experimental methods, we have developed a cost-effective approach, which integrates computational and experimental methods to uncover SL gene pairs for colorectal cancer (CRC). Our approach incorporates SL pairs in literature, microarray gene expression data, protein expression (immunohistochemistry [IHC] of CRC tissues) and phenotypic (clinicopathological factors) data. Let A-B denote a gene (protein) pair of A and B henceforth. Screening across 600 + previously verified SL pairs using microarray gene expression data of cancerous and normal tissues (a pilot study), we found that some SL gene pairs, such as BRCA1-PARP1 and TP53-SGK2, were simultaneously differentially expressed in high percentages in several cancerous tissues including CRC. Studies in S. cerevisiae and C. elegans have also shown that co-expression of gene pairs may be a relevant feature for predicting genetic interactions genome-wide [6], [8], [9]. Therefore, these verified SL pairs with differential gene expression patterns were classified as initial candidate tumor-dependent (TD) gene pairs (Table 1), from which we selected 17 genes to conduct IHC staining at multiple cellular locations. Next, to find novel SL pairs, we combined all possible pairs from the ~ 20 individual IHCs giving a total of 169 protein pairs. Because protein pairs related to tumor cell viability and malignancy may be indicated by their correlation with clinicopathological features obtained from CRC patients, we then tested each of these protein pairs for synergistic correlation with clinical features of the same set of patients to identify TD pairs, which are our predicted SL pairs. The cancer phenotypes observed (e.g., overall survival) suggested that the tumor cells may depend on each protein pair for viability and, therefore, simultaneously mutating each TD pair may kill tumor cells. RNAi, small-molecule inhibitor and in vivo experiments were then performed to validate the top two predicted SL partners of TP53, which is involved in multiple pathways, leading to challenges in therapeutic development. A graphical illustration of our approach is presented in Supplementary Figure W1. Finally, we identified protein pairs whose abnormal IHCs were correlated with poor patient survival, and analyzed the IHCs of the 169 protein pairs and clinicopathological features to find independent prognostic markers of CRC patients.

Table 1

Initial panel of candidate TD pairs for CRC, showing the percentages of the four gene expression patterns.

TD pairs		Percentages of gene pairs from 70 Asian CRC versus 12 non-cancerous tissues that were expressed 2-fold or more
Gene 1	Gene 2	(up, up) pattern	(up, down) pattern	(down, up) pattern	(down, down) pattern	Permutation P value*	FDR†
FEN1	RAD54B	0.67‡	0.00	0.00	0.00	0.0001	1×10- 5
BRCA2	PARP1	0.34	0.00	0.00	0.00	0.0001	1×10- 5
MSH2	POLB	0.33	0.00	0.00	0.00	0.0001	1×10- 5
BRCA1	PARP1	0.26	0.00	0.00	0.00	0.0001	1×10- 5
MYC	AURKB	0.17	0.00	0.00	0.00	0.0001	1×10- 5
RB1	SKP2	0.06	0.00	0.00	0.00	0.0002	2×10- 5
CSNK1E	CTNNB1	0.03	0.00	0.00	0.00	0.0001	1×10- 5
BCR	WNT5A	0.03	0.01	0.00	0.00	0.0001	1×10- 5
KRAS	SLKRAS§	0.01	0.00	0.00	0.00	0.0001 ~ 0.7200	1×10- 5 ~ 0.0303
TP53	SGK2	0.00	0.33	0.00	0.20	0.0004	3×10- 5
KRAS	SLKRAS¶	0.00	0.01	0.00	0.00	0.0001 ~ 0.8540	1×10- 5 ~ 0.0351
PTEN	PARP1	0.00	0.00	0.01	0.00	0.0001	1×10- 5
ABL1	WNT5A	0.00	0.00	0.01	0.00	0.0043	0.0003

This P value was for the highest percentage of the four patterns computed by permutation test with 10,000 repeats.

FDR was estimated by q-value.

The four fractions are computed from gene pairs that were 2-fold differentially expressed, thus they might not sum up to 100%.

57 verified KRAS SL pairs were identified in the (up, up) pattern.

19 verified KRAS SL pairs were identified in the (up, down) pattern.

Materials and Methods

Computation of Gene Expression Profiles of Colorectal Cancer Versus Non-Cancerous Tissues

We filtered gene expression data sets according to the following parameters: both cancer and non-cancerous tissues, no treatments, no metastasis and affymetrix chips (up to Nov. 2010). The CRC gene profiles satisfying the above criteria were downloaded from the Gene Expression Omnibus database [20]. Because mutated genes involved in the oncogenesis of a given cancer are known to vary between patients of different ethnic backgrounds [21], we collected gene expression data from patients of Han Chinese origin (GSE 9348 in [22]), the same ethnic background as the IHC and clinicopathological data used in later sections. Gene expression profiles from all 70 colorectal cancer and 12 non-cancerous tissues in the dataset were first quantile-normalized (expresso in R), then for a given gene the log ratio of its expression in each cancer tissue versus that of the averaged non-cancerous tissues was computed. The dataset used is shown in Supplementary Dataset S1.

Inferring the Initial Panel of TD Gene Pairs for Colorectal Cancer Using Microarray Gene Expression Data

We first collected 663 gene pairs whose SL interactions were either validated using various human cancer cell lines [11], [23] or genome-wide RNAi knockdown [15], [16] (Supplementary Table W1). For each SL pair collected, the percentages of tissues in which both genes were (up, up), (up, down), (down, up), or (down, down) were computed, respectively. Next, the top-ranking 11 gene pairs and several KRAS pairs were selected to make up the panel of initial candidate TD pairs for CRC.

Permutation Test and False-Positive Rates of the Percentages of Paired Gene Expression in Table 1

To assess the statistical significance (P value) of the percentages of the (up, up) or (up, down) patterns of each gene pair in Table 1, for each percentage we performed a permutation test to generate its nonparametric distribution. The total rearrangements of the labels of (70) cancer and (12) non-cancerous tissues was equal to , from which we randomly chose 10,000 rearrangements. For each rearrangement we computed the percentage of a pair’s pattern to form its distribution, from which we assessed the significance of an observed percentage. Next, we applied q-value (qvalue in R; [24]) to estimate false discovery rate (FDR) of the significance of the gene pairs and the predicted SL pairs in Table 1, Table 3, respectively.

Table 3

The 11 predicted SL pairs which had synergistic effects when correlated with clinicopathological features.

Protein 1	Protein 2	P value			FDR*	Log-rank (P value)
		Protein 1	Protein 2	Protein pairs
A. Overall survival
CSNK1E	P53	0.065	0.120	0.015	0.467	0.005
CTNNB1	P53	0.346	0.120	0.025	0.467	0.018
CSNK1E	MYC	0.065	0.369	0.027	0.467	0.064
ABL1	CSNK1E	0.169	0.065	0.033	0.467	0.066
AURKB	CSNK1E	0.524	0.065	0.050	0.467	0.066
B. Tumor size
MSH2(N)	POLB	0.069	0.123	0.001	0.168	N/A†
BRCA1	POLB	0.123	0.123	0.022	0.865	N/A
CSNK1E	POLB	0.115	0.123	0.045	0.865	N/A
C. Metastasis
MSH2(C)	RB1	0.255	0.698	0.001	0.179	N/A
BRCA2	MSH2(C)	0.400	0.255	0.048	0.989	N/A
D. Lymph node metastasis
BRCA1	CSNK1E	0.127	0.268	0.049	1.000	N/A

FDR was estimated by q-value.

N/A denotes not applicable.

Preparation of Tissue Microarrays

One hundred seventy-one representative cancer specimens (Supplementary Table W2) were randomly chosen from hematoxylin and eosin–stained sections and confirmed by pathologists. Four tissue cores (2 mm in diameter) were obtained from each paraffin block from which three cores of cancerous and one core of noncancerous tissues were cut longitudinally. The tissue cores were set into new paraffin blocks using a fine steel needle to produce tissue microarrays. This study was approved by the institutional review board and ethics committee of Changhua Christian Hospital and the institutional review board of Academia Sinica.

Immunohistochemistry

The sections (4 μm) of tissue microarray were deparaffinized in xylene and hydrated in serial dilutions of alcohol. Endogenous peroxidase activity was quenched by 3% H2O2. The sections were performed by treatment with boiling citrate buffer (10 mmol/L) for 20 minutes for antigen retrieval. The tissues were then incubated with 17 primary antibodies (Supplementary Table W3). Incubation with the primary antibody was carried out for 30 minutes at room temperature, followed by rinsing three-times with phosphate-buffered saline. The slides were incubated with a horseradish peroxidase/Fab polymer conjugate (EnVision™ Detection systems Peroxidase/DAB, Rabbit/Mouse [K5007 HRP; DaKo]) for another 30 minutes. After rinsing, the chromogen was developed with 3,3′-diamino-benzidine tetrahydrochloride as the substrate and hematoxylin as the counterstain. The staining intensity in the cancerous tissue was examined by two pathologists (Kun-Tu Yeh and Tzu-Cheng Su). The scoring criteria were similar to those in [25]. Specifically, the stain intensity was graded as negative-0, indeterminate-±, weak positive-1 +, moderate positive-2 + or strong positive-3 +. For staining intensity negative-0, there is no expression of the detected protein; for indeterminate-±, the staining is weak and its percentage cannot be accurately counted; for weak positive-1 +, there is less than 5% expression of the detected protein in the cancer cells; for moderate positive-2 +, there is focal expression in 5% to 20% of the cancer cells; for moderate positive-3 +, there is diffuse expression over 20% of the cancer cells. The staining intensity in the cancerous tissue was categorized as overexpression or underexpression through comparison with noncancerous colonic mucosa.

Fisher’s Exact Test

Fisher’s exact test was conducted to determine independence between IHC expression levels of each predicted SL pair and clinic-pathological features. IHC expression levels of each protein were dichotomized into two classes, abnormal and normal, according to the criteria in Table 2. Further, clinical factors of 171 CRC patients were also dichotomized as follows; survival time < 3 years versus ≥ 3 years; tumor size > mean versus tumor size ≤ mean; tumor grade poorly-differentiated cells versus well and moderately differentiated cells; metastasis yes versus no; lymph node metastasis N1-2 versus N0; and stage III and IV versus I and II.

Table 2

List of IHC-stained proteins and the percentages of 171 patients with abnormal IHC.

No.	Protein name	Criterion of abnormality	% of patients
1	ABL1(C)*	≥ 2 +	56
2	AURKB(C)	≥ 1 +	94
3	BCR(C)	≥ 2 +	28
4	BRCA1(C)	< 1 +	28
5	BRCA2(C)	< 2 +	32
6	BRCA2(N)*	< 1 +	34
7	CSNK1E(C)	< 2 +	54
8	CTNNB1(N)	≥ 1 +	44
9	MSH2(C)	< 1 +	8
10	MSH2(N)	< 1 +	28
11	MYC(N)	≥ 2 +	83
12	P53(N)	≥ 1 +	78
13	PARP1(N)	≥ 1 +	87
14	POLB(N)	≥ 1 +	66
15	RAD54B(N)	≥ 1 +	82
16	RB1(N)	< 1 +	24
17	SGK2(C)	< 2 +	73
18	SKP2(C)	≥ 2 +	23
19	WNT5A(C)	≥ 1 +	60

The notation (C) and (N) represent cytoplasm and nucleus, respectively.

Cell Culture

HCT-116 cell was obtained from ATCC; RKO, and TP53−/− derivatives of HCT116 and RKO were obtained from GRCF Cell Center and Biorepository of John Hopkins University. HCT-116 cells were cultured in DMEM, and RKO cells were cultured in McCoy 5A, all supplemented with 10% FBS and antibiotics (Invitrogen). Cells were seeded at density of 105/ml one day prior to shRNA transfection, lentivirus infection and small-molecule inhibitor treatment.

RNAi Knockdown and Small-Molecule Inhibitor Experiment

Plasmids each harboring gene-specific or scrambled shRNA were ordered from the RNAi core, Genomics Research Center, Academia Sinica. shRNAs were transiently transfected into isogenic HCT-116 cells using Lipofectamine-2000 (Invitrogen) following the manufacturer’s protocol, harvested for MTT assay and RT-PCR 48 hours post-transfection. Isogenic RKO cells were infected with lentiviruses harboring the shRNAs and harvested 4 days after puromycin selection. Cells were treated with IC261 and D4476 (Sigma) (HCT116 and RKO), and 2,4-Diamino-Quinazoline (Sigma) (HCT-116) at indicated concentrations. Cells were harvested for MTT assay 48 hours after treatment. For The TRC IDs corresponding to the shRNA sequences used please see Supplementary Materials and Methods.

Real-Time PCR

Real Time PCR was performed using Roche UPL Real Time-PCR System (Roche) according to manufacturer’s protocol; mRNA levels were measured using the ΔΔCT method, with beta-actin as internal control for each sample. For primer sequences and corresponding UPL probes used please see Supplementary Materials and Methods.

MTT Assay

MTT assay was performed using manufacturer’s protocol (Invitrogen).

Mouse Model Experiment

Animal care in accordance with institutional guidelines and authority to perform in vivo work was granted by the Home Office (China Medical University). 2 x106 HCT-116 cells (TP53 and TP53) were injected subcutaneously into both flanks of groups of five 6- to 7-week old NOD.CB17-Prkdc/JNarl male mice purchased from the National Laboratory Animal Center, Taiwan and treated for eight consecutive days with either IC261 (20.5 mg/kg, dissolved in 5% Kolliphor®EL (Sigma), 5% DMSO, and 90% PBS) or DMSO 10 days after cell inoculation. Tumor volumes were measured every day and calculated as tumor width × tumor length × tumor height. Significance of the difference in tumor volume was tested as described in Supplementary Text S1.

Mantel-Haenszel Log-Rank Test

For each predicted SL pair, we generated Kaplan-Meier survival curves of the high and low risk groups of patients using the software R, where the high and low risk groups comprise patients whose paired IHC expression levels were both abnormal and others, respectively. The log-rank test was used to calculate the significance of the differences between the survival curves of the two groups.

Cox Proportional Hazard Regression

For each predicted SL pair and clinicopathological feature, we computed a P value based on univariate Cox regression model using global permutation test, which permuted survival time and censor status simultaneously. We then included IHC of CSNK1E-P53, CTNNB1-P53, MSH2(C)-RB1 and BRCA1-WNT5A, which were significantly correlated with survival, four clinicopathological features, and statuses of chemotherapy treatment and radiation oncology treatment as covariates to conduct a multivariate Cox regression by stepwise variable selection (using the software R).

Results

Initial Panel of Candidate Tumor-Dependent Gene Pairs for Colorectal Cancer

To understand the role that SL pairs play in cancer cells, we first investigated the gene expression patterns of the known SL pairs collected. The more closely functionally related a SL pair is, more frequently the genes may be simultaneously differentially expressed in cancer cells. We gathered a list of 663 pairs of SL gene pairs verified in various cancer cell lines; most of the SLs were validated by genome-wide RNAi screenings [15], [16], while the remaining SLs were verified via individual RNAi knockdown experiments using cancer cell lines [11], [23], most of which were conducted based on hypothesis generated from pathway information and/or conservation of SLs from other species.

Validated SL pairs comprise oncogenes, tumor-suppressor genes, stability genes [2] and others. In general, tumor cells depend on the overexpression of oncogenes and/or under-expression of tumor suppressor genes for survival or/and proliferation. Therefore, for each SL pair we looked for the simultaneously up-regulated pattern (abbreviated as the [up, up]) one up-regulated and one down-regulated patterns ([up, down] and [down, up]) and the simultaneously down-regulated pattern ([down, down]), where the cutoff of up- and down-regulation were 1 and -1 in log2 ratio of gene expression of cancer versus non-cancerous tissues, respectively. Specifically we computed the percentages of the aforementioned four patterns of the collected SL gene pairs using the log ratios of gene expression (see Materials and Methods). To include as many potential TD pairs as possible, we first included the gene pairs whose percentages of any of the four patterns were at least one percent. Because RNAi knockdown is easier than overexpression of a particular gene, we sorted the gene pairs by the percentages of the (up, up), (up, down), (down, up), and then (down, down) patterns to result in Table 1, which included genes highly mutated in CRC such as TP53 [26] (see Materials and Methods).

Although the over-expression of tumor-suppressor and stability gene pairs, e.g., DNA repair and checkpoint classes, in Table 1 seems surprising at first glance, it is consistent with the drastic increase in genomic instability and DNA replication caused by mutant oncogenes such as KRAS and MYC.

Correlating IHC Stainings of TD Pairs With Clinicopathological Features Identifies SL Pairs

We first filtered out FEN1, KRAS and PTEN from the candidates in Table 1 resulting in 17 genes for IHC staining. FEN1 was not chosen because no isogenic FEN1 mutant CRC cell line is available. KRAS and PTEN were omitted because we aimed to uncover novel SL pairs, and multiple SL targets of KRAS and PTEN have been identified through genome-wide RNAi screening [15], [16] and small-molecule library screening, respectively [27]. IHCs of these remaining 17 proteins at multiple cellular locations were conducted (Table 2, Supplementary Figure W2) to confirm the trends of their mRNA expression using CRC tissues dissected from 171 patients from Changhua Christian Hospital, Changhua City, Taiwan (Supplementary Table W2, see Materials and Methods).

To extend the analysis beyond the protein pairs listed in Table 1, we then considered all possible pairs (169 pairs) of proteins, from the 19 candidates in Table 2 (excluding the pairs of the same protein stained in different cellular locations). Note that this procedure allows novel SL pairs to be included in this study, which is essential for uncovering novel SL pairs. To see how the TD pairs were related to the phenotypes of CRC, we correlated IHCs of each of these 169 protein pairs with five clinicopathological features of the patients: overall survival, tumor size, histological grade, lymph node metastasis, and metastasis. Please note for simplicity, henceforth we have omitted the location attached to the protein name, e.g., CSNK1E(C) is written as CSNK1E, unless ambiguity arises.

A pair of proteins is called to have a synergistic effect if their paired abnormal IHCs are significantly correlated (P < 0.05; Fisher’s exact test) with a clinicopathological feature (e.g., overall survival), but the abnormal IHC of each single protein is not. Any protein pair passing this synergy screening was predicted to be SLs, because this synergistic effect implies that abnormal IHCs of a protein pair is associated with a cancer phenotype (e.g., a poor overall survival), but the abnormal IHC of each individual protein is not. In other words, CRC cells seem to depend on each of these protein pairs, but not each individual protein, for maintenance of viability. Thus, simultaneous mutation of each TD pair may eliminate the tumor cells, which is synthetic lethality. Note that in model organisms [6], [8], [9], coexpression and phenotype have been shown to be useful in predicting SL interactions. However, the proposed synergistic effect is novel in (1) prioritizing the gene pairs co-expressed in cancer cells and (2) correlating their IHC expression to clinicopathological features (served as phenotypes) of CRC patients.

Correlating the paired and individual IHCs of 169 protein pairs to each of the five clinicopathological features resulted in 11 predicted SL pairs (Table 3); Fisher’s exact test (see Materials and Methods) was applied to test the significance of the correlation (P < 0.05), and the false discovery rate of a test was estimated by the q-value (qvalue in R), which measured the proportion of false positives incurred when the test was claimed significant (those with P < 0.05) [24]. The complete results are shown in Supplementary Table W4.

We first checked the 11 predicted SL pairs (Table 3) against the literature (searched up to February 2014). Out of these predicted SL pairs, two pairs MSH2-POLB and CSNK1E-MYC are known SL pairs, and have been validated in cell lines [28] and via large scale RNAi screening [17], respectively. The gene expression of CSNK1E was shown to be associated with that of MYC in colon cancer. Because CSNK1E was implicated in cancer only recently [17], these findings demonstrate that our method is able to uncover SL gene pairs without being limited by pathway information or without prior knowledge of whether a target gene is involved in cancer.

Table 3 shows that correlating IHCs of the protein pairs to overall survival resulted in the largest number of predicted SL pairs, and included P53. Note that TP53 has a very large sequence mutation score in CRC (Table 3 in [26]) and is mutant in about half of all cancers. Further, the fact that CSNK1E-MYC ranked third in Table 3A was a verified SL pair suggested that the top-two pairs in Table 3A, CSNK1E-P53 and CTNNB1-TP53 may be novel SL pairs; we validated these using RNAi, small-molecule inhibitor and mouse model as follows.

A rule of thumb for prioritizing predicted SL pairs for validation: first, incorporating prior knowledge, e.g., known SL pairs or/and pathway information. Second, using biological importance, e.g., those having frequently mutated cancer genes in the cancer of interest will be on top of the list.

CSNK1E and CTNNB1 are Validated as SL to TP53 by RNAi Knockdown (CTNNB1 with TP53), Small-Molecule Inhibition and In Vivo Model (CSN1KE with TP53)

We selected two pairs of isogenic CRC cell lines, HCT-116 and RKO, varying only in TP53 status (TP53WT and TP53-/-), as in vitro models to validate that TP53- CSNK1E and TP53-CTNNB1 are SL pairs. We transiently transfected HCT-116 cells and infected RKO cells with lentiviruses with shRNAs targeting CSNK1E (CSNK1E-1 and CSNK1E-2) and CTNNB1 (CTNNB1-1 and CTNNB1-2). Knockdown of each gene markedly reduced the viability of TP53−/− cells as compared to TP53WT cells (Figure 1A-D, for details please see Supplementary Text S2). This result further implies that the proposed approach is suitable for the effective prediction of SL gene pairs. The other predicted SL pairs shown in Table 3 are also promising and could be validated in the future.

Figure 1

RNAi knockdown of CSNK1E (A and B) and CTNNB1 (C and D) selectively killed TP53−/− cells. Cell viability was measured by MTT assay. Relative mRNA expression level was measured using the comparative CT (ΔΔCT) method; beta-actin mRNA was used as an internal control. Relative viability of HCT-116 and RKO cells (TP53WT and TP53-/-) transfected with (A) CSNK1E-targeting shRNAs (#1 and #2) or (C) CTNNB1-targeting shRNAs (#1 and #2), and control shRNA. Relative mRNA level of (B) CSNK1E following transfection of CSNK1E-targeting shRNAs (#1 and #2) or (D) CTNNB1 following transfection of CTNNB1-targeting shRNAs (#1 and #2), and control shRNA. Relative viability of cells (TP53WT and TP53-/-) following treatment with IC261 in (E) HCT-116 and (F) RKO. Relative viability of cells (TP53WT and TP53-/-) following treatment with D4476 in (G) HCT-116 and (H) RKO.

To evaluate the specificity and efficacy of small-molecule inhibitors of CSNK1E and CTNNB1 in killing TP53−/− cells, we treated the aforementioned cell lines with inhibitors of CSNK1E, IC261 and D4476 (Figure 1E-H). Upon treatment with IC261, we observed a reduction in cell viability in TP53−/− HCT-116 cells (IC50 of 1.37 μmol/L as compared to 8.30 μmol/L, which was about 6-fold higher); in RKO cells, about 1.8-fold difference in IC50 was observed (0.62 μmol/L as compared to 0.34 μmol/L; Figure 1E and F). Treatment with D4476 yielded similar fold changes of IC50 of TP53 cells and TP53−/− cells (about 4.7-fold) in both cell lines (250.72 μmol/L as compared to 53.22 μmol/L for HCT116 cells, 60.79 μmol/L as compared to 12.76 μmol/L for RKO cells; Figure 1G and H). However, treatment of 2,4-Diamino-Quinazoline (inhibitor of CTNNB1) in HCT-116 cells did not exhibit selective toxicity against TP53−/− cells (IC50 of 995 μmol/L and 821 μmol/L for TP53WT cells and TP53−/− cells, respectively). We wonder whether this is due to the non-specific targeting of 2.4-Diamino-Quinazoline, which targets the CTNNB1/TCF complex. Since IC261 and D4476 were shown to inhibit microtubule polymerization [29] and ALK5 (transforming growth factor beta1 type I receptor) [30], respectively, the use of these compounds to asses the physiological role of CSNK1E needs to be cautiously interpreted.

We further evaluated the efficacy of IC261 in an in vivo preclinical model. TP5 (TP53) HCT-116 cells were injected subcutaneously into the right (left) flank of groups of five NOD/SCID mice and allowed to grow for 10 days (tumor volumes ranged from 12.50 mm3 to 41.72 mm3 and 13.02 mm3 to 48.19 mm3 for TP53 and TP53−/− tumor xenografts on the 10th day, respectively). Subjects were then treated for nine consecutive days with either IC261 or DMSO before sacrifice. Tumor volumes were recorded daily after the start of the treatment.

At the end of experiment, tumor volumes ranged from 165.60 mm3 to 487.20 mm3 and 195.86 mm3 to 725.91 mm3 for DMSO-treated TP53 and TP53−/− tumor xenografts, respectively. For IC261-treated mice, tumor volumes ranging from 50.82 mm3 to 250.88 mm3 and 25.67 mm3 to 342.90 mm3 were observed for TP53 and TP53−/− tumor xenografts, respectively.

For each day, the tumor volumes of five TP53 xenografts were adjusted to their mean of the same day. Two sample t-test was conducted using TP53−/− xenograft tumor volumes under treatment IC261 versus DMSO, calibrated according to the volumes of their TP53 counterparts, to test the hypothesis that tumor volumes of IC261-treated TP53−/− xenografts were less than their DMSO-treated counterpart based on the definition of synthetic lethality. Calibrated TP53−/− xenograft tumor volumes were 490.93 mm3 ± 323.35 mm3 and 138.03 mm3 ± 79.56mm3, for DMSO-treated and IC261-treated samples, respectively (for details, please see Supplementary Text S1).

Consistent with our in vitro data, IC261-treated HCT-116 TP53−/− tumor (n = 5) in NOD/SCID mice showed markedly reduced tumor formation (Figure 2) as compared to the corresponding treatment with DMSO (n = 5), with the last three tumor volume differences of about 2.0-, 2.4- and 3.6-fold, respectively (P value = 0.033, 0.045, 0.023, for the last three days, respectively; two-sample t-test). Additionally, our preventive model, by pre-treating TP53−/− cells with DMSO and IC261 before cell inoculation, also showed selective inhibition of TP53−/− tumor growth by IC261 (Supplementary Materials and Methods and Supplementary Figure W3).

Figure 2

IC261 treatment inhibited TP53−/− tumor growth in vivo. (A) Representative images of xenografts obtained with HCT-116 cells (image for each treatment was the same mouse) in NOD/SCID mice after treatment with either DMSO or IC261. (B) Mean tumor size over the treatment regimen with IC261 and DMSO control (n = 5 for each group, error bars indicate standard error, the notation * indicates P < 0.05).

CSNK1E-P53, CTNNB1-P53, MSH2(C)-RB1 and BRCA1-WNT5A are Correlated with Poor Prognosis

Prognosis provides patients with information to make decisions about adopting further treatments, thus it is important in clinical medicine. The results of the log-rank test in Supplementary Table W4A suggested that of the 171 CRC patients, those with abnormal IHCs of CSNK1E-P53, CTNNB1-P53, MSH2-RB1 or BRCA1-WNT5A had a significantly different survival curve than that of patients in the remaining subgroup (P = 0.005, 0.018, 0.002 and 0.041, respectively; Mantel-Haenszel log-rank test; see Materials and Methods).

Kaplan-Meier survival analysis was performed, and the estimated survival curves of patients with abnormal IHCs for each of these four pairs was indeed separated from that of the remaining subgroup (Figure 3). Interestingly, none of individual IHC of CSNK1E, P53, CTNNB1, MSH2(C), RB1, BRCA1 or WNT5A could separate the estimated survival curves of patients into two subgroups (P = 0.067, 0.119, 0.536, 0.338, 0.374, 0.118 and 0.510, respectively; log-rank test; Supplementary Figure W4).

Figure 3

Abnormal IHCs of the four protein pairs are correlated with poor prognosis Kaplan-Meier survival curves of CRC patients divided by the paired IHCs of (A) CSNK1E and P53, (B) CTNNB1and P53, (C) MSH2(C) and RB1, and (D) BRCA1 and WNT5A. The curves for the patients with paired abnormal IHCs are plotted in solid line and the curves of the remaining patients are plotted in dashed line; the symbol “-” denotes abnormal IHC of the corresponding protein.

CSNK1E-P53, CTNNB1-P53, MSH2(C)-RB1, BRCA1-WNT5A are Independent Prognosis Markers from Stage for CRC Patient Survival

Potential prognostic markers can be identified from gene, microRNA or protein expression data [31], [32], [33], [34]. We further investigated which IHC of the 169 protein pairs were correlated with CRC patient survival, because no IHC of any individual protein was significantly correlated with overall survival. A univariate Cox proportional hazard regression analysis, applied to an independent cohort of 171 CRC patients, suggested that IHC of CSNK1E↓-P53↑, MSH2(C)↓-RB1↓, CTNNB↑1-P53↑ and BRCA1↓-WNT5A↑ (the symbols “↑”and “↓”denoting abnormal IHC [overexpression and underexpression, respectively] of the corresponding protein) were relevant prognosis factors with P values equal to 0.006, 0.015, 0.018 and 0.047, respectively [32], [33] (see Materials and Methods); see Supplementary Table W5 for the complete results. Moreover, survival analysis of four clinical factors indicated that stage was correlated with survival (P = 0.034), while age, sex and grade were not (P = 0.289, 0.351, 0.466, respectively).

Subsequently, a multivariate Cox proportional hazard regression analysis, with stepwise variable selection, was applied to IHCs of the aforementioned four protein pairs, the four clinicopathological factors and the statuses of chemotherapy treatment (CT) and radiation treatment (RT) [32], [33]; age, grade, CT and RT statuses were filtered out in the first four steps.

The four protein pairs and stage were found to be significantly associated with patient survival (P = 0.001, 0.007, 0.027, 0.010 and 0.010, respectively; with hazard ratio = 7.53, 1.65, 1.50, 1.93 and 1.58, respectively; Table 4). Another multivariate Cox regression model applied to the five significant factors excluding CTNNB1-P53 resulted in about the same hazard rates for these factors as shown in Table 4 (with minor changes in the first decimal places); the result was similar when CSNK1E-P53 was deleted. Thus, the correlation between these two P53 pairs did not seem to influence the hazard ratio result.

Table 4

Overall survival of 171 CRC patients in relation to clinical factors and IHC of protein pairs.

A. Univariate Cox regression
Variable	Subset	Hazard ratio (95% CI)	P value
Age	age > 65/age < = 65	1.19 (0.86-1.66)	0.289
Gender	male/female	1.17 (0.84-1.63)	0.351
Grade	poor/well and moderate	1.27 (0.67-2.42)	0.466
Stage	III-IV/I-II	1.42 (1.03-1.97)	0.033
(MSH2(C), RB1(N))	(↓,↓)*/otherwise	4.34 (1.58-11.95)	0.015
(CSNK1E(C), P53(N))	(↓,↑)/otherwise	1.62 (1.15-2.26)	0.006
(CTNNB1(N), P53(N))	(↑,↑)/otherwise	1.50 (1.07-2.09)	0.018
(BRCA1(C), WNT5A(C))	(↓,↑)/otherwise	1.64 (1.02-2.63)	0.047
B. Multivariate Cox regression†
Variable	Subset	Hazard ratio (95% CI)	P value
Stage	III-IV/I-II	1.58 (1.11-2.24)	0.010
(MSH2(C), RB1(N))	(↓,↓)*/otherwise	7.53 (2.25-25.24)	0.001
(CSNK1E(C), P53(N))	(↓,↑)/otherwise	1.65 (1.14-2.38)	0.007
(CTNNB1(N), P53(N))	(↑,↑)/otherwise	1.50 (1.05-2.14)	0.027
(BRCA1(C), WNT5A(C))	(↓,↑)/otherwise	1.93 (1.17-3.19)	0.010

The symbols “↑”and “↓”denote abnormal IHC (overexpression and underexpression, respectively) of the corresponding protein.

Variables were selected by using stepwise selection method.

CSNK1E-P53 is an Early-Stage Prognosis Marker While the Remaining Three Protein Pairs are Markers Throughout the Stages

These four protein pairs and stage are independent prognostic factors, thus the former could serve as markers for CRC patient survival. We further studied whether these four paired markers depended on stage or not as follows. Each pair of IHC, stage and their interaction was fitted by a multivariate Cox regression model, which resulted in CSNK1E-P53 being dependent on stage. Note that this pair is an early-stage prognosis marker (both Stage I-II and Stage I-III) (P = 0.032 and 0.007, n = 86 and 141, respectively; log-rank test). For details, please refer to Supplementary Table W6 and Figure W5.

Discussion

We have developed a cost-effective approach to uncover SL gene pairs for therapeutics of CRC, which predicted SL pairs efficiently via finding synergistic correlation of IHC (protein) pairs with clinicopathological factors (phenotypes) of CRC patients. Of the 11 predicted SL pairs, MSH2-POLB and CSNK1E-MYC were known SL in human, and we further validated CSNK1E-TP53 and CTNNB1-TP53. Our method suggested that clinical factors may be used as phenotypes to predict synthetic lethality. In addition, MSH2-RB1 was previously implied as SL in a mouse study [35]. The study demonstrated that mouse with deficiency in both MSH2 and RB1 showed delayed development of lymphoma, associated with increased apoptosis of lymphoma cells, as compared to MSH2-deficient mice. With rich repertoire of public domain gene expression data (e.g. The Cancer Genome Atlas), commercially available tissue arrays with clinical data, it is conceivable that our approach can be applied to other cancers in the future.

While safety regarding inhibition of TP53 is controversial [36], our validation of CSNK1E as SL to TP53 in vitro (further supported by a previous study in mouse MEF cells [37]) and in vivo showed that this pair may be translated to preclinical study. Note that IC261, the inhibitor of CSNK1E used in both in vitro and in vivo studies also inhibited microtubule polymerization [29], which may lead to increased cancer cell death [38]. However, in our in vivo study the IC261-treated TP53−/− xenograft tumor volumes were calibrated to the volumes of their TP53 counterparts, which adjusted the effect of microtubule polymerization, and the IC261-treated tumors still shrunk more significantly than the DMSO-treated ones. This suggests that synthetic lethality between TP53 and CSNK1E plays a significant role in killing cancer cells. CSNK1E and CTNNB1 are both known components of WNT signaling pathway. Both activated WNT signaling [39] and TP53 inactivation [40] are frequently implicated in cancer tumorigenesis. Previous study showed that loss of TP53 lead to activation of WNT signaling pathway, mediated through miR-34 [41], suggesting that our result is consistent with current literature regarding synthetic lethality, which arises primarily from compensatory pathways [4], thus implying a buffering relationship between p53 and WNT signaling pathway in CRC. Since TP53 is mutated in half of all cancers, CSNK1E may be a promising druggable target for a broad population of TP53-mutant patients.

CSNK1E↓-P53↑, MSH2(C)↓-RB1↓, CTNNB↑1-P53↑ and BRCA1↓-WNT5A↑ were revealed to be prognosis markers, with the first pair for early-stage and the remaining for all stages, independent from stage for CRC, but interestingly none of ~ 20 individual IHC was correlated with overall survival. Previous studies identified prognostic markers via cellular phenotypes (such as cell differentiation) [42] and protein pairs [43], separately, through insights from pathway information, while our approach identified markers without requiring prior knowledge of molecular mechanism.

Finally, our method can be extended as follows. First, we can pair a gene of interest, preferably a frequently mutated gene in cancers, with the remaining genes in human genome to constitute the initial panel and the remaining procedures follow similarly. Second, this approach can be applied to several cancers to predict their unique/common SL pairs and prognosis markers.

Supplementary Figures

Supplementary Tables

Supplementary Dataset S1

Supplementary Text S1

Supplementary Text S2

Supplementary Materials and Methods