RBMS1 Suppresses Colon Cancer Metastasis through Targeted Stabilization of Its mRNA Regulon

Metastatic progression in colorectal cancer is accompanied by widespread gene expression reprogramming. Cancer cells often co-opt post-transcriptional regulatory mechanisms to achieve pathologic expression of gene networks that drive metastasis (1–3). Colorectal cancer is the third most commonly diagnosed cancer (4); therefore, understanding the underlying regulatory programs that drive metastatic progression in this disease is a crucial step toward improving patient outcomes. Notably, there are not many predictive computational methods aimed at the discovery of unknown regulatory networks. By relying on annotated post-transcriptional regulatory pathways, for example those mediated by miRNAs, existing methods fail to capture previously unknown regulatory interactions. To tackle this problem, we have developed a computational approach called PRADA (Prioritization of Regulatory Pathways based on Analysis of RNA Dynamic Alterations) that identifies post-transcriptional master regulators responsible for aberrant mRNA stability and gene expression in cancer cells. By applying this tool to a large compendium of gene expression profiling data from patient samples, patient-derived xenograft models, and established colon cancer cell lines, we identified a novel regulatory program involved in colorectal cancer metastasis. We find that this previously unknown regulatory pathway, which controls mRNA stability and is mediated by the RNA-binding protein RBMS1, is often silenced in highly metastatic colorectal cancer tumors. We demonstrate that RBMS1 stabilizes transcripts by binding to the last exons of target mRNAs in concert with the RNA-binding protein ELAVL1, and that, in highly metastatic colorectal cancer cells and patient tumors, RBMS1 downregulation is associated with poor clinical outcome. We identify mRNA targets of RBMS1 that are functionally relevant to metastasis and reveal that RBMS1 silencing can be accomplished through epigenetic dysregulation. This study not only describes a disease-relevant post-transcriptional regulatory mechanism that governs the stability of a sizeable regulon, but also demonstrates the value of bottom-up discovery strategies such as PRADA that do not rely solely on prior knowledge of annotated regulatory programs.

Metastasis is a complex multistep process and requires modulation of many cellular pathways and functions. As such, increased metastatic capacity often involves broad reprogramming of the gene expression landscape in cancer cells. Thus, a mechanistic dissection of cancer progression relies on the systematic identification of the underlying regulatory programs that drive pathologic cellular states. To accomplish this, we developed PRADA, which uses regulatory network predictions to identify the key RNA-binding proteins (RBP) that may act as master regulators of mRNA stability and gene expression. Conceptually, PRADA solves a regression problem to predict changes in gene expression as a function of the expression of RBPs and their known or predicted targets across the transcriptome. In this customized regression analysis, the coefficient assigned to each RBP reflects its strength as a regulator of gene expression and the direction of its effect (i.e., activator or repressor). Given our limited knowledge of post-transcriptional regulatory pathways and their role in human disease, PRADA provides an opportunity for the discovery of previously unknown disease pathways. To assess the contribution of post-transcriptional regulatory programs to colon cancer metastasis, we took advantage of a publicly available set of gene expression profiles from established colon cancer cell lines (GSE59857). We first categorized these cell lines as poorly metastatic or highly metastatic based on their liver colonization capacity in xenograft mouse models (Supplementary Fig. S1A; ref. 5). We then performed differential gene expression analysis to compare gene expression changes across the transcriptomes of the two groups. Finally, we applied PRADA to this dataset to find RBPs whose differential activity is most informative of metastasis-associated gene expression changes. PRADA identified the protein RBMS1 as the RNA-binding protein with the strongest regulatory potential in this dataset (Fig. 1A). The size and the direction of the regulatory coefficient assigned to RBMS1 implies that RBMS1 levels are strongly informative of changes in the expression of transcripts with sequence matches to RBMS1-binding sites. To confirm this, we performed motif enrichment analysis using FIRE, which uses mutual information to capture the association between the presence and absence of a given binding site and genome-wide transcriptomic measurements (6). As shown in Fig. 1B, the RBMS1 consensus binding site is strongly informative of gene expression differences between poorly and highly metastatic cells, showing a highly significant enrichment among the transcripts that have decreased expression in highly metastatic cells. To ensure that this result is not sensitive to the choice of a specific RBMS1 binding motif, we tested two additional sequence models: (i) an independently derived representation of RBMS1-binding preferences (7) and (ii) predictions based on DeepBind models (ref. 8; Supplementary Table S1). Each of these three predicted RBMS1-binding motifs gave largely identical results, and we observed a similar decrease in the expression of the putative RBMS1 regulon derived from each motif (Supplementary Fig. S1B). For consistency, we have used the DeepBind-derived RBMS1 regulon (referred to as the RBMS1 regulon hereafter) in our subsequent analyses, as DeepBind is the state-of-the-art approach for predicting protein–nucleic acid interactions. Importantly, concordant with the results reported by PRADA and the lower expression of its regulon, the expression of RBMS1 was also significantly lower in highly metastatic cells (Fig. 1C). Moreover, consistent with RBMS1 acting as a post-transcriptional regulator of these genes, we observed a significant correlation between the expression of RBMS1 and the average expression of its regulon in multiple independent datasets (Fig. 1D; Supplementary Fig. S1C and S1D).

To further assess the association between RBMS1 and colorectal cancer metastasis in more directly disease-relevant models, we took advantage of patient-derived xenograft (PDX) models of colon cancer metastasis to the liver. We used RNA-sequencing (RNA-seq) data from three parental colorectal cancer PDX models (CLR4, CLR27, and CLR32) and their highly liver-metastatic derivatives (CLR-LvM). The highly metastatic CLR-LvM models were derived from the parental PDXs by repeated splenic delivery and growth of cells in the liver of immunocompromised mice (NOD/SCID gamma), recapitulating the major site of colorectal cancer metastasis in humans (9). We observed that the increase in the metastatic capacity of the CLR-LvM PDX models compared with their parental PDXs was accompanied by a significant reduction in the expression of the RBMS1 regulon in the highly metastatic CLR-LvM models (Fig. 1E). More importantly, all three models also showed a significant and concomitant decrease in the expression of RBMS1 (Fig. 1F).

In addition to these cell line and PDX models of colorectal cancer liver metastasis, we also performed RNA-seq on matched primary tumors and liver metastases biopsied and frozen from two patients. As shown in Fig. 1G, in both cases the RBMS1 regulon was enriched among transcripts that were downregulated in the metastases relative to their primary tumors. Interestingly, although the gene expression changes between the primary tumors and metastases were not generally correlated (R = 0.02, P = 0.1), the RBMS1 regulon was independently downregulated in both metastatic tumors (Supplementary Fig. S1E). Consistently, RBMS1 was also strongly silenced in both metastases (Fig. 1H). Collectively, these findings in cell lines, PDXs, and clinical samples implicate RBMS1 silencing in both the downregulation of its putative regulon and in colorectal cancer metastatic progression.

The correlated expression of RBMS1 and its putative regulon defined from its binding site preferences implies that RBMS1 modulates gene expression. To further assess the potential role of RBMS1 as a post-transcriptional regulator, we performed RNA-seq of RBMS1 knockdown and control SW480 colon cancer cells. We chose the SW480 cell line for this experiment because: (i) RBMS1 expression in SW480 is among the highest in colon cancer lines tested by us and others and (ii) it is an established xenograft model of liver metastasis and is considered to be poorly metastatic (relative to other colorectal cancer lines listed in Supplementary Fig. S1A). We used two independent shRNAs to silence RBMS1 in SW480 cells and confirmed its knockdown with both qPCR and Western blotting (Supplementary Fig. S2A). A 2.5-fold reduction in RBMS1 protein expression induced major changes in gene expression (Supplementary Fig. S2B) and resulted in a significant decrease in the expression of RBMS1 regulon (Fig. 2A). Because we are focused on the role of RBMS1 as an RNA-binding protein, we reasoned that its effect on gene expression is likely through regulation of the stability of its target RNAs. To test this hypothesis, we took advantage of REMBRANDTS, a computational framework we have developed to estimate RNA stability from RNA-seq data (10). We have previously established that REMBRANDTS accurately recapitulates experimental RNA stability measurements (10). Application of this method to RNA-seq data from RBMS1 knockdown and control cells found a significant enrichment of RBMS1 targets among transcripts destabilized upon RBMS1 silencing (Fig. 2B). REMBRANDTS relies on the comparison of exonic and intronic reads to measure changes in RNA stability. As shown in Supplementary Fig. S2C, RBMS1 silencing results in lower expression of its putative regulon as measured by exonic reads; however, intronic reads, which reflect changes in pre-mRNA levels (and transcription rates), do not significantly change in response to RBMS1 depletion. To further strengthen our claim that RBMS1 functions as a post-transcriptional regulator of RNA stability rather than as a transcriptional activator, we also performed whole-genome RNA stability measurements in control and RBMS1 knockdown cells. For this, we used α-amanitin–mediated inhibition of RNA-polymerase II followed by RNA-seq at two time points (0 and 9 hours). Consistent with our analyses with REMBRANDTS, we observed a similar reduction in the stability of the putative RBMS1 regulon upon RBMS1 silencing (Fig. 2C). We used RNA-seq data from the three matched PDX models (CLR panel) to compare RNA stability in highly liver metastatic models to that of their parental PDXs (using REMBRANDTS). As shown in Fig. 2D, lower RBMS1 expression in the LvM models (Fig. 1F) accompanies a reduction in the stability of its target regulon in all three independently derived models. Together, these findings establish RBMS1 as a post-transcriptional regulator of RNA stability with broad functional consequences for the transcriptome and clear implications for colorectal cancer progression.

To this point, our analyses were performed using in silico predicted binding targets of RBMS1. To create a transcriptome-wide snapshot of RBMS1-binding sites in colon cancer cells at nucleotide resolution, we performed UV cross-linking immunoprecipitation followed by sequencing (CLIP-seq). We carried out irCLIP (11) for endogenous RBMS1 in SW480 cells and identified hundreds of high-confidence RBMS1-binding sites across the transcriptome (Fig. 3A). Importantly, we observed that 90% of the CLIP-identified RBMS1 targets overlapped with our computationally derived putative RBMS1 regulon (Supplementary Table S1). Consistently, these RBMS1-bound transcripts in our irCLIP data (RBMS1 targets hereafter) follow the same expression patterns as the computationally derived RBMS1 regulon described above, exhibiting both decreased stability and expression upon RBMS1 silencing (Fig. 3B; Supplementary Fig. S3A).

In our analysis of RBMS1-binding sites, we noted a strong enrichment of RBMS1 binding to the last exon (on or close to 3′ UTRs; Fig. 3A) of mRNAs. This is consistent with RBMS1 acting as a post-transcriptional regulator of RNA stability, and suggested that RBMS1 3′ UTR binding could result in increased RNA stability. To assess this possibility, we built a bidirectional CMV promoter that drives the expression of both GFP and mCherry. We then cloned nine CLIP-seq–derived RBMS1-binding sites, plus approximately 150 nucleotides flanking the binding sites, downstream of GFP, and asked whether there was a reduction in GFP mRNA relative to mCherry mRNA upon RBMS1 depletion. As shown in Fig. 3C, there was a decrease in GFP reporter transcript levels in almost all instances. To further ensure that this effect was RBMS1-dependent, we also tested the two reporters with the strongest reduction in transcript levels (DAG1 and CD9) in LS174T cells, where RBMS1 is endogenously silenced. In this instance, as expected, there was no response in reporter mRNA levels upon transfection of RBMS1-targeting siRNAs (Supplementary Fig. S3B).

To gain insight into how RBMS1 regulates the stability of its targets, we carried out an unbiased search for its interacting protein partners. We performed immunoprecipitation of RBMS1, along with an IgG control, from SW480 cells and analyzed the resulting protein complexes by mass spectrometry (Fig. 3D). From this data, we found that the Gene Ontology terms RNA binding, RNA processing, and negative regulation of mRNA metabolism were overrepresented among the RBMS1-interacting proteins (Supplementary Table S2). Of these proteins, we have confirmed an RNA-dependent interaction between RBMS1 and two well-characterized RBPs that predominantly bind 3′ UTRs, PABPC1 and ELAVL1 (Fig. 3E). These findings suggest that RBMS1 stabilizes its targets at least partially in concert with ELAVL1 which preferentially binds A/U-rich elements (ARE; ref. 12). In line with this hypothesis, we found a significant enrichment of poly(U) motifs in the 3′ UTRs of RBMS1 targets (Fig. 3F). As shown in Fig. 3G, 72% of RBMS1-bound 3′ UTRs are also bound by ELAVL1 in vivo [as determined by photoactivatable ribonucleoside-enhanced-CLIP (PAR-CLIP); ref. 13]. Furthermore, knockdown of ELAVL1 in SW620 colon cancer cells (of the same genetic background as SW480 cells) resulted in a significant downregulation of RBMS1 targets, as determined by RNA-seq (14; Fig. 3H). Together, these results indicate that direct binding of RBMS1 to mRNA 3′ UTRs, together with other stabilizing factors such as ELAVL1, results in increased mRNA stability.

To carry out a functional study of RBMS1 and its downstream regulon in colorectal cancer metastasis, we picked an additional cell line model to complement the SW480 line. We chose the LS174T line because RBMS1 is almost completely silenced in this line as measured by qPCR and Western blotting (Supplementary Fig. S4A and S4B). LS174T is also an established xenograft model and is considered a highly liver metastatic cell line, approximately 100 times more metastatic than the SW480 line as measured by liver colonization assays (15). We first performed RNA-seq of both SW480 and LS174T lines and, consistent with RBMS1 silencing in the LS174T cells, we observed a significant reduction in the expression of the RBMS1 regulon in this line (Supplementary Fig. S4C).

To establish a causal link between RBMS1 silencing and higher metastatic capacity, we performed liver colonization assays by splenically injecting RBMS1 knockdown and control SW480 cells and measuring metastatic burden in the livers of mice over time. Although RBMS1 silencing did not have a strong effect on in vitro cell proliferation (Supplementary Fig. S4D), we observed a significant increase in liver colonization upon silencing RBMS1, based on in vivo bioluminescence measurements and gross liver mass at the experimental endpoint (Fig. 4A; Supplementary Fig. S4E). To control for possible off-target effects of the shRNAs, we repeated the experiment with an additional RBMS1-targeting hairpin and observed a similarly significant increase in metastatic liver colonization (Supplementary Fig. S4F). We also performed a gain-of-function experiment by expressing RBMS1 in the LS174T line, where RBMS1 is endogenously silenced. Consistent with our earlier findings, we observed that expressing RBMS1 in LS174T cells results in a marked reduction in their liver colonization capacity (Fig. 4B; Supplementary Fig. S4G). It should be noted that, as expected, expressing shRNAs targeting RBMS1 in highly metastatic cells in which RBMS1 is already silenced, namely LS174T, HCT116, and COLO320, did not have an impact on their liver colonization capacity in xenograft mouse assays, further supporting the on-target effect of the shRNAs used in these experiments (Supplementary Fig. S4H). In contrast, knockdown of RBMS1 in WiDr cells, in which RBMS1 is endogenously expressed, resulted in a significant increase in liver colonization in xenograft mouse models without an overt effect on in vitro cell proliferation (Supplementary Fig. S4I and S4J).

The metastatic cascade is a complex multistep process, involving cellular processes that affect not only cancer cell growth and survival, but also cell migration and invasion. As mentioned earlier, we did not observe a significant change in in vitro cell proliferation rates upon RBMS1 knockdown. In addition, transwell invasion assays of RBMS1 knockdown and control cells did not detect a significant role for RBMS1 in cancer cell invasion (Supplementary Fig. S4K). Therefore, to gain a better understanding of the key step(s) in the metastatic cascade that are regulated by RBMS1, we performed a systematic analysis of known and predicted gene sets to identify those that may be modulated upon RBMS1 silencing. We performed this analysis using iPAGE, an information-theoretic framework we have developed for this type of analysis (ref. 16; Supplementary Fig. S4L). Downregulation of genes associated with negative regulation of epithelial–mesenchymal transition [EMT(−) signature gene set defined by various studies; refs. 17–19] was the most significant pathway identified in RBMS1 knockdown cells (Fig. 4C; Supplementary Fig. S4L and S4M). For example, the canonical EMT marker E-cadherin (CDH1) is significantly downregulated by approximately 2-fold in RBMS1 knockdown cells based on RNA-seq data (Supplementary Fig. S2B). We confirmed this by performing immunofluorescence staining for E-cadherin in RBMS1 knockdown and control cells. We observed a reduction in E-cadherin expression in RBMS1 knockdown cells and also noted the appearance of spindle-like cell morphology that is associated with EMT (Fig. 4D). Moreover, the EMT(−) signature genes were also expressed at relatively lower levels in the LS174T cell line, where RBMS1 is endogenously silenced (Supplementary Fig. S4M). Finally, analysis of The Cancer Genome Atlas Colon Adenocarcinoma (TCGA-COAD) dataset revealed significantly lower E-cadherin expression in tumor samples with low RBMS1 expression (∼20% reduction, P = 0.006) and a significant general correlation between RBMS1 and E-cadherin expression (ρ = 0.1, P = 0.03), which is indicative of the clinical relevance of this regulatory axis in colorectal cancer.

To identify genes that are regulated by RBMS1 and act as suppressors of metastasis in colorectal cancer, we used an integrated analytic approach, which relies on mining relevant datasets to prioritize target genes based on their direct interaction with RBMS1 (3′ UTR binding), their RBMS1-dependent magnitude of change in expression, the robustness of their response to RBMS1 silencing, as well as their lower expression in metastatic clinical samples. Using this approach, we prioritized two target genes, AKAP12 and SDCBP, that satisfy these criteria: (i) destabilized and downregulated in highly metastatic colorectal cancer cells and PDXs, (ii) destabilized and downregulated in RBMS1 knockdown cell lines, (iii) directly bound by RBMS1 based on our irCLIP data (Fig. 3A), and (iv) downregulated in liver metastases relative to primary colon cancers in a publicly available dataset (20). Consistently, silencing RBMS1 in SW480 cells resulted in the downregulation and destabilization of AKAP12 and SDCBP mRNAs (Fig. 5A and B). Conversely, overexpression of RBMS1 in LS174T cells resulted in upregulation and stabilization of these targets (Fig. 5A and B). In addition, in line with the concerted action of RBMS1 and ELAVL1, both AKAP12 and SDCBP 3′ UTRs are bound by ELAVL1 in vivo (13). These observations establish AKAP12 and SDCBP as direct targets of RBMS1 that are post-transcriptionally regulated by this RBP.

To test the possible role of AKAP12 and SDCBP in driving colorectal cancer metastatic progression, we silenced them individually in SW480 cells using CRISPR interference (CRISPRi; ref. 21) and performed liver colonization assays. As shown in Supplementary Fig. S5A, silencing these genes did not have a significant impact on in vitro proliferation of SW480 cells. However, reduced expression of these genes resulted in increased metastatic liver colonization in mice (Fig. 5C). To confirm that AKAP12, which, of the two targets tested, elicited a stronger phenotype when silenced, indeed functions downstream of RBMS1, we performed an in vivo epistasis experiment by generating cells with knockdown of both AKAP12 and RBMS1 in SW480 cells. As shown in Supplementary Fig. S5B, unlike in Fig. 5C, further silencing of AKAP12 failed to increase metastatic liver colonization in RBMS1 knockdown cells, in which AKAP12 is already downregulated. To evaluate the contribution of AKAP12 to the phenotype of RBMS1 silencing, we compared AKAP12 knockdown and double knockdown cells. As shown in Supplementary Fig. S5C, we observed a slight additional increase (albeit not statistically significant) in metastatic liver colonization upon RBMS1 silencing in AKAP12 knockdown cells. This suggests that AKAP12 is a major effector downstream of RBMS1, and that other downstream targets, such as SDCBP, may also play minor roles in this process. Importantly, and consistent with AKAP12 acting downstream of RBMS1, gene expression profiling of AKAP12 knockdown and control cells revealed a significant induction of an EMT signature (Fig. 5D). Notably, this included the upregulation of canonical EMT transcriptional repressors SNAI1 and SNAI2 (Supplementary Fig. S5D).

As an orthogonal approach, to assess whether the effects of RBMS1 silencing on the transcriptome are reflected in the proteome, we used tandem mass tag labeling and tandem mass spectrometry to compare the global proteomes of SW480 cells with shRNA-mediated RBMS1 knockdown and control cells (Supplementary Fig. S5E; Supplementary Table S2). As shown in Supplementary Fig. S5F, we observed a positive and significant correlation between changes in mRNA abundance (as determined by RNA-seq) and changes in protein abundance. We found that RBMS1 targets, as determined by irCLIP, were enriched among the downregulated proteins in the RBMS1 knockdown cells (Supplementary Fig. S5G). Most notably, AKAP12, which we have identified as the downstream effector of RBMS1 in EMT regulation, was significantly downregulated at the protein level as well, mirroring our transcriptomic results (Supplementary Fig. S5H). Taken together, our findings demonstrate that the RBMS1–AKAP12 regulatory axis acts as a suppressor of EMT and liver metastasis in models of colorectal cancer progression.

To further assess the clinical relevance of this previously unknown regulatory pathway, we performed a variety of measurements in clinical samples as well as analysis of clinical data to evaluate RBMS1 activity in colorectal cancer metastasis. First, we performed qRT-PCR for RBMS1 in two independent clinical cohorts, one cohort stratifying patients based on their tumor stage (n = 96) and another comparing samples from normal mucosa, primary colorectal cancer tumors, and liver metastases (n = 91). As shown in Fig. 6A and B, we observed a significant reduction in RBMS1 expression as the disease progresses, with stage IV clinical samples showing the lowest RBMS1 expression levels. We also carried out survival analysis in a large clinical cohort with publicly available expression data and clinical outcomes (22). We observed a significant association between RBMS1 silencing and reduced relapse-free survival as well as overall survival in patients with colon cancer (Fig. 6C). A multivariate Cox proportional-hazards model revealed that this association with RBMS1 silencing remains statistically significant even when controlled for other known clinical factors that may contribute to patient survival (Supplementary Fig. S6A). This observation emphasizes the relevance of RBMS1 silencing as an effective readout for risk stratification of patients based on samples collected from their primary tumors.

Our results indicate that AKAP12 acts as a suppressor of metastasis downstream of RBMS1 and therefore is expected to show a similar association with metastatic disease. To test this hypothesis, we performed qRT-PCR measurements in the clinical cohorts mentioned above, and we observed a significant reduction in AKAP12 expression as a function of disease progression (Fig. 6D and E). Consistent with RBMS1 acting as a regulator of AKAP12, we observed a highly positive and significant correlation between the expression of these two genes (Supplementary Fig. S6B). Interestingly, our analyses indicate that the identified RBMS1 targets provide a robust gene signature that, similar to RBMS1, is highly informative of clinical outcomes. For this analysis, we defined an RBMS1 target score as an aggregate measure of expression for the RBMS1 80-gene signature set [average normalized expression of target mRNAs that (i) interact in vivo with RBMS1 on or within 200 nucleotides of their 3′ UTRs and (ii) are downregulated upon RBMS1 knockdown; this is referred to as the RBMS1 80-gene signature hereafter (Supplementary Table S1)]. We then stratified patients based on this score and we observed that lower expression of the RBMS1 80-gene signature is associated with lower relapse-free and overall survival in patients with colon cancer (Fig. 6F). Consistently, the RBMS1 signature score remained a significant covariate in a multivariate Cox proportional-hazards model of relapse-free survival (Supplementary Fig. S6C). Moreover, we have not observed significant association between the RBMS1 signature score and validated clinicopathologic features or other molecular markers (such as microsatellite instability status, left- or right-sidedness, or colorectal cancer consensus molecular subtypes). Importantly, similar to our initial observation with the RBMS1 regulon (Fig. 1D), the RBMS1 signature score was significantly correlated with RBMS1 expression in a large TCGA pan-cancer dataset (Fig. 6G). We also observed a high correlation between the RBMS1 signature score and ELAVL1 expression (Fig. 6H), consistent with RBMS1 and ELAVL1 acting in concert to promote mRNA stabilization. The association between the RBMS1 regulon and colorectal cancer metastasis is not limited to a single dataset. Consistent with our initial observations in matched samples (Fig. 1G), comparison of gene expression changes in liver metastases relative to colorectal cancer primary tumors in a publicly available dataset (23) revealed a significant reduction in the expression of the RBMS1 regulon (Supplementary Fig. S6D and S6E). Furthermore, analysis of publicly available RNA-seq data using REMBRANDTS to infer changes in mRNA stability from a set of matched primary and metastatic samples (24) also confirmed (i) RBMS1 silencing in multiple liver metastases and (ii) a concordant reduction in the stability of the RBMS1 regulon (Supplementary Fig. S6F). Together, these results establish the clinical relevance of our findings and the importance of RBMS1 silencing both as a prognostic factor and as a suppressor of liver metastasis in patients.

Our findings described above establish a previously unknown regulatory pathway driven by the RNA-binding protein RBMS1 that plays a functional role in colorectal cancer metastasis to liver. Regulatory pathways can be expanded by uncovering the upstream pathways that influence their activity. We have developed computational and data-analytic tools designed to integrate publicly available datasets from a variety of biological sources to identify such upstream regulatory mechanisms (16). Using these tools, we sought to extend the RBMS1 regulatory pathway by exploring the mechanisms for RBMS1 silencing observed in the highly metastatic LS174T colon cancer cells. First, we found RBMS1 expression to be strongly associated with promoter acetylation, but did not identify any specific transcription factors associated with RBMS1 expression (Fig. 7A). Consistent with this observation, analysis of the connectivity map dataset (25), which reports the impact of hundreds of small-molecule treatments on gene expression, identified the histone deacetylase inhibitor trichostatin A (TSA) as an activator of RBMS1 expression and of the RBMS1 regulon (Supplementary Fig. S7A and S7B). As RBMS1 is endogenously silenced in LS174T cells and expressed in the SW480 line, we tested the relative acetylation levels of the RBMS1 promoter in these two lines by performing H3K27Ac chromatin immunoprecipitation (ChIP)–qPCR. These data showed relatively low levels of H3K27 acetylation at the RBMS1 promoter in the LS174T line compared with the SW480 line, consistent with the differences in the level of RBMS1 in the two lines (Fig. 7B). TSA simultaneously inhibits multiple histone deacetylases (HDAC); however, we noted that among this group, HDAC1 is the most significantly upregulated in highly metastatic cells generally, and in LS174T cells compared with SW480 cells in particular (Fig. 7C and D; Supplementary Fig. S7C). Moreover, we observed that RBMS1 and HDAC1 expression are negatively correlated in colon cancer clinical samples (Supplementary Fig. S7D). Furthermore, we observed that silencing HDAC1 increases RBMS1 expression (Fig. 7E; Supplementary Fig. S7E and S7F). Increased expression of HDAC1, and other class I HDAC genes, has been reported as a strong predictor of survival in patients with colon cancer (26), and modulations in HDAC activity result in widespread gene expression reprogramming, including repression of tumor suppressor genes (27). Our observations here indicate that HDAC1-mediated transcriptional repression may be one possible mechanism for RBMS1 silencing in colon cancer cells.

Complex human diseases, including cancer, often accompany broad reprogramming of gene expression. In these cases, a comprehensive understanding of the disease state requires not only the identification of the differentially expressed genes, but also understanding the underlying regulatory pathways that explain the dysregulated expression patterns. A number of approaches have been developed by us (16, 28) and others (6, 29) to tackle this problem. These methods formalize the association of expression and/or activity of master regulators with those of their regulons. Because direct and indirect associations are challenging to disentangle, knowledge of binding preference (in vitro or in vivo) is also used to better define putative RNA–RBP interactions. Here, we have introduced PRADA to facilitate and formalize the discovery of post-transcriptional regulators that are involved in normal cell physiology and disease. RNA-binding proteins fall into families with highly similar binding preferences (7) and therefore have similar putative regulons. As a result, direct modeling of RBP–target interactions results in unstable predictions (i.e., the model fails to converge on similar outcomes for repeated runs) where the link between a given RBP and its putative regulon is not clear. PRADA solves this problem by integrating the knowledge of changes in the expression of RBPs (used as a proxy for their activities) directly into the analysis. This approach provides a one-step and often stable solution that effectively reveals the RNA-binding proteins whose differential expression is highly informative of changes in gene expression. As RBPs continue to emerge as key regulators of RNA dynamics with crucial roles in human disease, systematic methods like PRADA provide a suitable approach for studying this class of regulators. In this study, we used this approach to discover a previously unknown regulatory network that counteracts colon cancer progression through RBMS1-mediated RNA stabilization.

We demonstrate that RBMS1 silencing accompanies colorectal cancer progression in cell line and PDX models as well as clinical samples, resulting in destabilization of RBMS1-bound transcripts (Fig. 7F). In line with this, we found several well-known RNA-stabilizing factors interacting with RBMS1. Importantly, we found that RBMS1-bound transcripts contain ELAVL1 motifs and binding sites in their 3′ UTRs, and that ELAVL1 downregulation resulted in the decrease of RBMS1 target abundance, consistent with a model where RBMS1 and ELAVL1 have an RNA-dependent functional interaction that promotes mRNA stability. Furthermore, we observed that ELAVL1 expression is highly correlated with the expression of the RBMS1 80-gene signature set in the TCGA pan-cancer dataset. ELAVL1 has a well-documented role in stabilizing mRNA by a mechanism that is incompletely characterized. It is suggested that ELAVL1 competes for 3′ UTR ARE binding with RNA-destabilizing factors, such as TTP or BFR1. It is also proposed that ELAVL1 acts in counteracting miRNA-mediated decay and translational shutdown, thereby indirectly stabilizing mRNAs (12). ELAVL1 upregulation has been broadly associated with oncogenesis and cancer progression (30), suggesting a functional interplay between RBMS1 and ELAVL1.

Metastasis often hijacks developmental pathways to reprogram gene expression. EMT is such a pathway with a central role in cancer progression, including in colorectal cancer (31). Our results show that RBMS1 downregulation results in the suppression of negative EMT regulators, thereby promoting the mesenchymal features of the tumor cells. One of the direct targets of RBMS1, AKAP12, participates in the protein kinase A and C signaling cascades. The association between AKAP12 silencing and colon cancer recurrence and prognosis has been described previously (32). Interestingly, AKAP12 has been described as a negative regulator of SNAI1 (33), a transcription factor that represses epithelial markers during EMT. In line with this, we found that AKAP12 knockdown caused a global shift toward an EMT transcriptional program, including the upregulation of SNAI1/2. One of the main targets of SNAI1 is the epithelial marker E-cadherin, which is in agreement with our observation that E-cadherin is downregulated in RBMS1-deficient cells and clinical samples. This suggests that the RBMS1–AKAP12–SNAI1/2–E-cadherin axis is a potential route to EMT onset during colorectal cancer metastasis.

The WNT–β-catenin signaling pathway has also been implicated in both EMT and colorectal cancer metastasis. β-catenin, a central factor in the canonical WNT signaling cascade, has a dynamic role: it is known to interact with E-cadherin at the cell membrane, with TCF/LEF transcription factors in the nucleus to activate the downstream transcriptional program, and with APC, a member of the destruction complex that targets β-catenin for proteasomal degradation in the absence of pathway activation (34). APC is a known mutational hotspot in colorectal cancer, and the SW480 cells used in this study harbor an inactivating APC mutation, indicating that the WNT pathway is misregulated in this cell line (35). We have shown that RBMS1 silencing results in the downregulation of E-cadherin and have detected upregulated β-catenin levels in RBMS1 knockdown cells (Supplementary Fig. S7G). Interestingly, one of the direct RBMS1 targets validated in this study is syndecan binding protein (SDCBP), which interacts with both Frizzled family proteins (canonical WNT receptors) and syndecans (WNT coreceptors; ref. 36). As our results suggest multiple crossover points, future studies will be required to elucidate the role of RBMS1 in WNT signaling cascades and their promotion of EMT.

Mirroring our observations from the colorectal cancer TCGA dataset, we found the expression of RBMS1 and its signature gene set were highly corelated in the TCGA pan-cancer dataset. However, when examined individually, RBMS1 and its signature gene set expression were informative for specific tumor types. RBMS1 silencing was associated with poor outcome in the colorectal cancer cohort, as discussed above (Fig. 6F); in contrast, different outcomes were observed in other cancer types; for example, in the thyroid cancer cohort, patients with tumors with relatively low expression of RBMS1 had better outcomes (Supplementary Fig. S7H and S7I). We observed a similar trend in other pan-cancer cohorts, notably in a dataset of metastatic tumors (MET500; ref. 37). Although RBMS1 and its signature expression remained highly correlated, we found reduced RBMS1 levels specifically in colorectal, prostate, and gall bladder cancer metastases. Interestingly, when stratified by the metastasis site, RBMS1 expression was significantly downregulated in liver and bone marrow (independent of the primary tumor type), reflecting the principal organotropisms of metastatic colorectal cancer (Supplementary Fig. S7J and S7K). It is thus highly plausible that RBMS1 represents a regulatory node particularly suited for disease progression in the liver, potentially via activation of downstream signaling pathways as discussed above. This finding has implications for other cancer types that can metastasize to the liver.

Our work in the LS174T cell line model indicated that increased activity of histone deacetylases, such as HDAC1, could be responsible for RBMS1 silencing. HDAC inhibitors have been used in the clinic to treat metastatic colorectal cancer. A number of studies have previously noted increased expression of HDAC1/2 in colon cancer (27), and in cases where RBMS1 silencing also occurs, HDAC inhibitors could be useful, as they may also restore RBMS1 expression. Although HDAC inhibitors are not specific in their effects, and their effects are not solely mediated through RBMS1, we speculate that in some cases their impact on RBMS1 may lead to better clinical outcomes. Identifying such cases, however, will require a deeper understanding of the mechanisms through which RBMS1 is silenced and the degree to which its expression can be effectively controlled.

Given our limited knowledge of the pathways and processes that regulate the RNA life cycle in the cell, analytic tools that mine quantitative measurements of mRNA dynamics to identify key regulatory interactions can provide an effective avenue for identifying previously unknown molecular mechanisms with critical functions in health and disease. However, these computational strategies must be paired with rigorous experimentation to functionally validate and characterize the putative regulons and their regulators. Using one such framework, we have established a novel functional association between RBMS1 silencing and increased liver metastatic capacity in colon cancer.

PRADA is a customized variation on lasso regression (least absolute shrinkage and selection operator). Therefore, PRADA can simultaneously perform feature selection and apply a custom penalty function as part of its regularization step. The goal of PRADA is to identify RNA-binding proteins whose differential expression explains changes in the expression of their targets that is observed in the data. We first generated an RBP–RNA interaction matrix based on binding preferences of RNA-binding proteins, as reported previously (7, 38). For this, we scanned mRNA sequences for matches to RBP-derived regular expression (as previously described; ref. 6). Given this interaction matrix, our goal is to identify RBPs whose change in expression is predictive of global changes in gene expression of their targets. In other words, we aim to solve the following: |$\Delta Exp( g ) = \sum\nolimits_i {{\alpha _i} \cdot {t_{i,g}} \cdot \Delta Exp( {RB{P_i}} )} + {c_g}$|⁠, where |$\Delta Exp( g )$| stands for changes in the expression of gene g based on dataset of interest, |${t_{i,g}}$| is a binary variable, 1 if RBP_i is predicted to bind g and 0 otherwise, |$\Delta Exp( {RB{P_i}} )$| marks changes in the expression of RBPs themselves as a proxy for their differential activity. The resulting |${\alpha _i}$| represents the strength of regulatory interactions. To ensure that the most informative RBPs are selected, we introduced a customized penalty term that extends the lasso regression framework: |$\min \frac{1}{{2n}}[ {\| {\alpha X - Exp} \| + \lambda \sum\nolimits_i {\frac{{| {{\alpha _i}} |}}{{| {\Delta Exp( {RB{P_i}} )} |}}} } ]$|Here, the coefficients are penalized by |${| {\Delta Exp( {RB{P_i}} )} |^{ - 1}})$|⁠, which is |$1/| {\hat{\beta }} |$| estimate of the linear model used to compare gene expression between the two groups. Other variations of this penalty can also be used; for example, standard effect size can account for confidence in |$\hat{\beta }$| as well (⁠|$SE/| {\hat{\beta }} |$|⁠). This custom penalty term ensures that RBPs whose activity does not change are not selected by the model. This also stabilizes the resulting model, which would otherwise be a major issue as RBPs that belong to the same family often have very similar binding preferences resulting in correlated features in the interaction matrix. After running this regression analysis, the RBPs with the largest assigned coefficients are prioritized for further study. In this study, we used PRADA to compare poorly and highly metastatic colon cancer lines, which revealed RBMS1 as the strongest candidate. The computational and experimental data flow used in this study is summarized in Supplementary Material S1. The PRADA benchmarking summary is available in Supplementary Material S2. All code and notebooks are available on Github at www.github.com/goodarzilab/PRADA.

SW480, LS174T, WiDr, HCT116, and COLO320 cell lines were obtained from ATCC. 293LTV cells were from Cell Biolabs. All cell lines were routinely screened for Mycoplasma infection by PCR (once a month in general and before every animal experiment in particular). SW480 and HCT116 cells were cultured in McCoy's 5A supplemented with 10% FBS, sodium pyruvate, and L-glutamine. LS174T and 293LTV cells were cultured in DMEM supplemented with 10% FBS, COLO320 cells were cultured in RPMI-1640 supplemented with 10% FBS, and WiDr cells were grown in Eagle Minimum Essential Medium supplemented with 10% FBS. All commonly used cell lines were routinely authenticated using STR profiling at UC Berkeley Sequencing Facility.

All animal studies were performed according to a protocol approved by UCSF Institutional Animal Care and Use Committee (AN179718). NOD/SCID gamma male mice (The Jackson Laboratory), aged 8 to 10 weeks, were used in all experiments. For splenic (portal circulation) injections, cells were injected directly into the spleen followed by splenectomy (250,000 for SW480 and LS174T lines and 100,000 for WiDr cells). In vivo bioluminescence was measured by retro-orbital injection of luciferin (Perkin Elmer) followed by imaging with an IVIS instrument (Perkin Elmer). For ex vivo liver imaging, mice were injected with luciferin prior to liver extraction, and the liver was then imaged and weighed after rinsing with PBS.

RBMS1 irCLIP was performed as described previously (11).

Unless otherwise specified below, RNA-seq libraries were prepared using RNA that had been rRNA depleted using Ribo-Zero Gold (Illumina) followed by ScriptSeq-v2 (Illumina), and sequenced on an Illumina HiSeq4000 at UCSF Center for Advanced Technologies. Matched primary and liver metastases were sequenced using SENSE RNA-seq Library Preparation Kit (Lexogen). AKAP12 knockdown cells were profiled using QuantSeq 3′ mRNA-Seq library prep kit fwd (Lexogen).

To measure RNA stability, SW480 RBMS1 knockdown and control cells were treated with 10 mg/mL α-amanitin (final concentration in the medium). After 9 hours, total RNA was harvested from the cells using the Norgen Cytoplasmic and Nuclear RNA Purification Kit per the manufacturer's protocol. RNA-seq libraries were then prepared and log-fold changes in RNA stability were measured by comparing log(t9/t0) in control and RBMS1 knockdown cells.

For clinical samples, 96 samples across all stages of the disease were acquired from Origene (HCRT104 and HCRT105), 14 normal, 8 stage I, 25 stage II, 32 stage III, and 17 stage IV. HPRT was used as endogenous control, and relative RBMS1 and AKAP12 expression levels were, respectively, measured across the samples (using the primers listed above). We also extracted RNA from another 100 normal mucosa, primary tumors, and liver metastases. Roughly 90 of these samples yielded sufficient RNA for qPCR, which was carried out as described previously.

Patients profiled and documented in GSE39582 were first stratified into two groups based on their RBMS1 expression: silenced (bottom 5%) and expressed (otherwise). Five percent was selected as the cutoff because it is close to two SDs away from the average RBMS1 expression across all samples. A multivariate Cox Proportional Hazard model (R package survival) was then used to evaluate 5-year disease-free survival. Univariate survival analyses, both disease-free and overall, were also carried out with this stratification. For the target regulon, the RBMS1 signature score was calculated as the median expression of target genes in each sample (normalized gene expression data). Disease-free survival was defined as the time from diagnosis to relapse or death due to any cause. We performed Cox Proportional Hazards multivariable modeling using the RBMS1 score as a continuous variable. Variables added to the model included stage, microsatellite status, KRAS and BRAF mutation status, and Consensus Molecular Subtypes.

The sequencing data generated as part of this study are accessible through Gene Expression Omnibus (GSE147749). All mass spectrometry raw files and their associated MaxQuant output data were deposited on ProteomeXchange Consortium via the PRIDE partner repository (PXD019478 and PXD019479).

B. Hänisch reports funding provided by the German National Merit Foundation and the Biomedical Education Program. R. Dienstmann reports personal fees and non-financial support from Roche, grants from Merck and Pierre Fabre, and personal fees from Boehringer Ingelheim, Ipsen, Amgen, Servier, Sanofi, and Merck Sharp & Dohme outside the submitted work. H. Goodarzi reports grants from the NIH and grants from AACR during the conduct of the study. No potential conflicts of interest were disclosed by the other authors.

J. Yu: Investigation, writing-original draft. A. Navickas: Conceptualization, funding acquisition, validation, investigation, visualization, methodology, writing-review and editing. H. Asgharian: Data curation, software, formal analysis, writing-review and editing. B. Culbertson: validation, investigation. L. Fish: Investigation, methodology, writing-original draft, writing-review and editing. K. Garcia: Validation, investigation. J.P. Olegario: Investigation. M. Dermit: Data curation, investigation. M. Dodel: Data curation, investigation. B. Hänisch: Investigation. Y. Luo: Data curation, formal analysis, investigation. E.M. Weinberg: Validation, investigation. R. Dienstmann: Conceptualization, formal analysis, investigation. R.S. Warren: Resources, validation. F.K. Mardakheh: Data curation, supervision, funding acquisition, validation, investigation, visualization, methodology. H. Goodarzi: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing-original draft, project administration, writing-review and editing.

We acknowledge the UCSF Center for Advanced Technology (CAT) for high-throughput sequencing and other genomic analyses. We thank Byron Hann and the Preclinical Therapeutics core as well as the Laboratory Animal Resource Center (LARC) at UCSF. We are also grateful for the genomic data contributed by the TCGA Research Network, including donors and researchers. We acknowledge support from our colleagues at the Helen Diller Family Comprehensive Cancer Center and Rockefeller University. We specifically acknowledge Dr. Sohail Tavazoie for providing access to data from PDX models and for reading an earlier version of this manuscript. We also acknowledge Jonathan Weissman and Luke Gilbert for CRISPRi constructs. We acknowledge Dr. David Erle for providing the backbone construct used for the reporter assay. This work was supported by grants from NIH (R00CA194077 and R01CA24098) and 2017 AACR-Takeda Oncology NextGen Grant for Transformative Cancer Research (17-20-38-GOOD; to H. Goodarzi). J. Yu was supported by an NSF Graduate Research Fellowship. A. Navickas was supported by a DoD PRCRP Horizon Award (DoD: W81XWH-19-1-0594). H. Asgharian and L. Fish were supported by the NIH training grants F32GM133118 and T32CA108462-15, respectively. M. Dermit, M. Dodel, and F.K. Mardhakheh were supported by MRC Career Development Award MR/P009417/1.

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