Overcoming Therapy Resistance in Colorectal Cancer: Targeting Rac1 Signaling Pathway as a Potential Therapeutic Approach

1. Introduction

Colorectal cancer (CRC) is the third most commonly diagnosed type of cancer worldwide, accounting for ten percent of all malignant diagnoses and responsible for numerous deaths, ranking third in terms of incidence, but second in terms of mortality [1] with estimations of increasing incidence for the next years [2]. Increased risk of CRC is normally associated with age, western lifestyle, diet, personal history and chronic intestinal diseases [3].

Most CRCs arise from previous adenomatous polyps. Although it remains uncertain the time it takes for an early adenoma to progress to CRC, the detection and removal of precancerous lesions before progressing to malignancy and metastasis provides a prevention strategy [1,4].

CRC treatment depends on factors such as patient’s health, tumor size and location, and the presence of metastasis, but surgical removal is the most common option [5] followed by administration of adjuvant chemotherapy when there is local-regional or distant invasion.

Since its discovery in 1957, chemotherapies based on 5-fluorouracil (5-FU, a fluorinated uracil analog) have remained the mainstay of adjuvant and palliative therapies for CRC patients [6]. 5-FU interferes with nucleosides metabolism by inhibiting the action of thymidylate synthase and missincorporating metabolites into DNA and RNA, leading to cytotoxicity and cell death [6,7,8]. Despite the clinical benefits and the extended use of 5-FU, response rates to 5-FU monotherapy are below 20% as most patients do not completely eliminate tumor cells, and tumor recurrence leads to poor outcomes [9]. To counteract this, various strategies have been developed to increase 5-FU effectiveness through the modulation of its intracellular and biochemical metabolism, mainly combining 5-FU with other cytotoxic drugs with different mechanisms of action. Combined chemotherapies FOLFOX (leucovorin+5-FU+xaliplatin), FOLFIRI (leucovorin+5-FU+irinotecan) and FOLFOXIRI (leucovorin+5-FU+oxaliplatin+irinotecan) are the most used in the clinic and have been used as the standard therapy for advanced CRC, increasing response rates up to 40-50% [10,11]. Unfortunately, despite the increased response CRC patient’s disease-free survival has not been efficiently extended [12,13], and half of CRCs are resistant to 5-FU-based therapies. Therefore, there is a need for studies focusing on characterizing resistance-mediators biological factors or identification of biomarkers to define which CRC population is most likely to respond to 5-FU-based therapies [14,15,16].

One of the main causes of failure in cancer treatment is the development of drug resistance by cancer cells. This is a very serious phenomenon since it causes the recurrence of the disease or even death [17]. Therapy resistance occurs when diseases become tolerant to pharmaceutical treatments. This concept was first considered when bacteria became resistant to certain antibiotics, but similar mechanisms have since been found in other diseases, including cancer. Some resistance mechanisms are specific to each pathology, while others, such as drug efflux can be observed at the microbial level and also in tumors, making them evolutionarily conserved [18]. Although many types of tumors are initially susceptible to chemotherapy, over time they may develop resistance to through this and other mechanisms, such as DNA mutations and metabolic changes that promote the inhibition and degradation of drugs [19]. Tumor resistance is not limited to conventional chemotherapeutic drugs, but also appears associated with the use of target therapies or biotherapies [20,21]. Consequently, numerous efforts are focused on exploring biomarkers that can predict or indicate the success of therapy.

In the case of CRC, 5-FU-based therapies increased response rates up to 40-50%, therefore strategies to improve clinical outcomes are required. To address this, a deeper understanding of the mechanisms associated with resistance or recurrence after 5-FU-based therapies is imperative.

On the other hand, Rac1 (Ras-related C3 botulinum toxin substrate 1) is a key member of the Rho GTPases family. It is well known that Rac1 is a regulator of actin-based cytoskeletal dynamics, modulating cell adhesion, morphology, and movement. Rac1 is highly expressed in different tumor types, and related to poor prognosis [22]. In tumors, it was described that Rac1 modulates cell cycle, apoptosis, proliferation, invasion, migration and angiogenesis. Rac1 also plays a key role in anti-tumor therapy and participates in immune escape mediated by the tumor microenvironment [22]. Increasingly, studies are reporting the role of Rac1 as a potential target for tumor therapy [23].

To contribute to understanding the mechanisms underlying CRC resistance to therapies, we have conducted a study integrating in silico, in vitro and in vivo approaches. First, we compared microarray gene expression data from 5-FU-treated CRC patients with and without recurrence after 5-FU monotherapy. Then, we extended the analysis to 5-FU-based therapies. Comparisons uncovered common enrichment pathways associated to chemotherapy resistance that allowed us to select drugs to overcome this phenomenon. Interestingly, inhibition of Rac1 by 1A-116 compound decreases the growth of 5-FU-resistant CRC, sensitizes cells to 5-FU therapy and prevents metastasis development, suggesting that therapies based on Rac1 inhibition could be of benefit to overcome therapy resistance.

2. Materials and Methods

2.1. Gene Expression Data Collection

Gene expression profiles of datasets GSE81653, GSE39582, and GSE72970 were downloaded from the NCBI Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds) using the GEOquery package [24,25]. These datasets contain gene expression profiling of tumor clinical samples from patients after the exposure to 5-FU alone and 5-FU based combined chemotherapy. Sample inclusion criteria for the analysis were: 1) not proceeding from studies of patients with familial hereditary polyposis, 2) including clinical information detailing the type of chemotherapy provided, and 3) presenting clinical information about tumor recurrence status. The expression values downloaded were normalized using RMA (Robust Multichip Average) [26] method from the affy or oligo packages depending on the chip model. Probe annotation was performed using the manufacturer-supplied annotation package: “hugene20sttranscriptcluster.db” [27] for the GSE81653 series and “hgu133plus2.db” [28] for the GSE39582 and GSE72970 series (a general flowchart of this work is presented in Figure 1).

If publicly available, information on clinicopathologic characteristics of patients was retrieved to explore their association with recurrent condition using Fisher's exact test, a significant association was considered for p-value < 0.05.

2.2. Identification of Differentially Expressed Genes in Patients Treated with 5-FU Monotherapy

Data analysis was carried out using the facilities of the CCT-Rosario Computational Center, member of the High Performance Computing National System (SNCAD, MincyT- Argentina) where data were introduced into the R/Bioconductor environment. Differentially expressed genes (DEGs) were identified in GSE39582 and GSE81653 datasets independently. Each dataset was divided in two groups: 1) patients treated with 5-FU monotherapy with tumoral recurrence after surgical tumor resection and 2) patients treated in the same way but without tumor recurrence. Detection platforms and sample sizes are shown in Table S1. A more detailed description of the selected samples is described in Table S2.

Four methods were used to obtain DEGs groups; two exploratory methods: fold change (“FC”) and expression in unusual proportion; and two non-parametric algorithms: the Significance Analysis of Microarrays (SAM) [29] and RankProd (RP) [30]. For DEGs selection, we used as cutoff criterion |log₂fold change| > 1 for FC method. For UR, any observation far from the mean by more than two standard deviations was considered atypical. For the non-parametric hypothesis tests, genes with |log2fold change| > 1 and FDR ≤ 0.01 were considered.

2.3. Enrichment Analysis

Reactome signaling pathway enrichment analysis for DEGs was performed using an R package “ReactomePA” [31]. ReactomePA uses the hypergeometric model to assess whether the number of selected genes associated with a Reactome pathway is larger than expected. Pathways with an FDR < 0.05 were considered significantly enriched and were visualized using the dotplot tool. We generated an enrichment map and pre-clustered network using the “enchmap” and “cnetplot” tools to visualize relationships between pathways, and highlight genes related to the more significant terms.

2.4. Identification of Gene Expression Profile-Reversing Compounds

A computational drug repositioning analysis was performed using CLUE (Connectivity map Linked User Environment: https://clue.io/query; Broad Institute, Cambridge, MA; data version 1.1.1.2; software version 1.1.1.42), a cloud-based analysis platform that catalogs 473,647 expression signatures of human cell lines treated with 25,200 perturbagens [32]. CLUE computes a CMap score (tau) that measures the similarity of a queried gene set of up to 150 upregulated and 150 downregulated genes to existing drug-matched reference gene sets, from least similar or inverse (‒100) to most similar (100). From an initial list of 426 DEGs, we selected 150 upregulated genes by screening for genes detected as DEGs by three or more methods in both data sets, as well as those associated with pathways with the highest enrichment score in the recurrent phenotype. Since the number of down-regulated DEGs was lower than 150, all of them were used. The reversing drugs of interest were selected using a negative threshold of CMap score.

2.5. Integrative Meta-Analysis

To explore whether the enriched pathways between patients with and without tumor recurrence after 5-FU treatment were maintained when patients treated with other types of chemotherapy were added to the expression matrix, we performed a meta-analysis. This technique is particularly useful for combining several datasets from the same disease when these are limited in size, therefore improving their statistical power. This time, instead of analyzing the data series independently, a large matrix was constructed merging patient’s gene expression data from GSE81653, GSE39582, and GSE72970 series. In addition to the patients treated with 5-FU already used in the previous analysis, others treated with FOLFIRI and FOLFOX were included. Data merging was performed using the empirical Bayes methods (ComBat, from Sva package) [33] to reduce confounding factors due to non-biological variations between the studies. Data about platforms and sample sizes are shown in Table S3. A more detailed description of the selected samples added to the analysis is found in Table S4.

2.6. Feature Selection

Feature selection is described as a process in which a subset of relevant features is selected from a larger data set. These features can be used for several purposes like model construction, differential analysis, enrichment exploration, etc. To select the most explanatory genes among the phenotypes from the combined gene expression matrix, we used the sigFeature package [34]. This package employs a combination of Vector Support Machine and Recursive Feature Elimination (SVM-RFE) algorithms to produce a ranked list of genes [35]. The SVM-RFE is a backward feature elimination technique that iteratively removes features based on SVM classifier weights. In each iteration, an SVM model is built based on the current features subset “F”, and the weight of each feature in “F” is calculated. The features are then ranked based on weight and the bottom-ranked features are removed from F until it is empty. The top-ranked features that are discarded in the last iteration are considered the most informative between the phenotypes. The number of features retained in the analysis depends on their particular future use. In this work, we selected the top 5000 more informative features to improve their potential to predict the biological signature of the data set.

2.7. Gene Set Enrichment Analysis (GSEA)

GSEA is a computational method used to determine whether an a priori defined set of genes shows concordant and statistically significant differences between two biological states [36]. We employed the java GSEA Desktop Application v4.2.3 to perform the GSEA analysis on 5000 more informative features selected. The gene sets available from The Molecular Signatures Database 3.0 (MSigDB) were employed [37]; genesets composed of less than 15 or more than 500 genes were excluded. The phenotype label was set as “recurrent” vs “non-recurrent”. The t-statistic mean of the genes was computed in each gene set using a permutation type test with 1000 replications. Up-regulated gene sets were defined by a normalized enrichment score (NES) > 0 and down-regulated by NES < 0. Gene sets with an FDR-P value ≤ 0.05 were chosen as significantly enriched.

2.8. Detection of Motifs and Transcriptional Factors

For the discovery of transcription factor binding sites (motifs) in the promoters of co-regulated genes, we used the Cytoscape plug-in iRegulon [38,39]. A collection of 9713 position weight matrices (PWMs) was applied to analyze 10 kb centered around the transcription start site. DNA logos corresponding to each motif were extracted and the main transcriptional factors binding to them and their sets of direct targets (metatargetnoma) were screened. Cut-off criteria used in the analysis were enrichment score threshold = 5, ROC threshold for AUC calculation = 0.03, rank threshold = 5000, minimal identity between orthologous genes = 0.05, FDR > 0.001, and normalized enrichment (NES) > 3.

2.9. Predictor Genes Detection

Important predictors genes from DEGs list and “reactome_signaling_by_rho_gtpases” gene set (https://www.gsea-msigdb.org/gsea/msigdb/cards/REACTOME_SIGNALING_BY_RHO_GTPASES) were selected via least absolute shrinkage and selection operator (LASSO) logistic regression [40] in individuals GEO series and in the merged matrix using the glmnet R package. To carry out the analysis, patients from each dataset were divided into training and validation sets in an 80-20 ratio, respectively. Optimal values for the penalty parameter λ were determined through 10-fold cross-validations. For model construction, we employed the lambda value giving minimal mean cross-validated error (lambda.min). The LASSO coefficient for recurrence predictors genes for each gene list was extracted. Receiver operating characteristic (ROC) curves were plotted to validate the prediction efficiency of the model using the ROCR and pROC R packages [41,42]. We employed survival package [43] to perform a univariate Cox regression to assess the effect of the change in expression of individual predictor gene on the survival of patients in each dataset; genes with a LogRank < 0.01 were retained.

2.10. PPI Network Construction and Hub Gene Identification

The PPI network was constructed using STRING public online database [44] and exported to Cytoscape (version 3.9.1) for viewing and analysis [39]. Cytoscape plug-in CytoHubba (version 0.1) [45] was employed to retrieve the top 20 hub genes based on maximal clique centrality (MCC) algorithm. To expand the network, main interactor genes were rescued using the network expansion function of STRING plug-in. Top 20 interactors with a selectivity of 0.5 were included in the network

2.11. Drugs

For each drug, a 10 mM stock solution was prepared and stored at −20°C. Ivermectin (Parafarm) was dissolved in 50% DMSO solution, the concentration of DMSO in the final dilutions did not exceed 0.03%. Amitriptyline (Parafarm) was dissolved in sterile distilled water. Commercial 20mg/ml irinotecan solution (Kemex) and 50mg/ml 5-FU solution (Fada Pharma) were diluted in sterile distilled water to reach final stocks concentration. 1A-116 [46] was kindly provided by Dr. Georgina Cardama (Laboratorio de Oncología Molecular, Universidad Nacional de Quilmes). 1A-116 stock was prepared by dissolving the drug in acidic water (sterile distilled water brought to pH 1-2 using a 100 mM HCl solution). Once dissolved, pH of the final solution was adjusted to 5.5-6 with 100 mM NaOH and filtered. Unless indicated otherwise, doses used of drugs were 5µM 5-FU; 20µM amitriptyline; 30µM irinotecan, 15µM ivermectin and 20µM 1A-116.

2.1. Cell Culture

5-FU resistant cell lines CT26^5FUR, HCT116^5FUR or HT29^5FUR were produced in our laboratory as described before [47] from CT26 (chemically induced BALB/c mice-derived colorectal carcinoma), HCT116 (male human colon adenocarcinoma) and HT29 (female human colon adenocarcinoma) cells respectively. Cells were cultured in DMEM (HCT116 and HT29) or RPMI (CT26) media, supplemented with 10% fetal bovine serum (FBS; Natocor, Argentina), penicillin (10 μg/ml), streptomycin (100 μg/ml) and L-glutamine (2 mM, DMEM medium). Cells were maintained at 37°C in a 5% CO₂ atmosphere and routinely tested for mycoplasma.

2.13. Cells Immunostaining Techniques

7× 10³ cells were cultured in coverslips for 24 h, treated for the 36 h, fixed in 4% paraformaldehyde and permeabilized with a 0.5% PBS-Triton and blocked with BSA (2%, Sigma). Cells were first incubated α-Tubulin (Sigma, 1:1000 dilution) and Beta-Catenin (Santa Cruz, sc-59737; 1/50 dilution) or E-Cadherin (67A4) (Santa Cruz, sc-21791; 1/50 dilution), and then with Alexa Fluor 488-conjugated secondary antibody (Invitrogen, 1:500 dilution), and counterstained with DAPI (Sigma, 1:10000 dilution) and phalloidin-rhodamine conjugate (Invitrogen, 1:2000 dilution) to observe cell nucleus and actin skeleton, respectively. Coverslips were fixed to a slide with Mowiol and observed under a Nikon Ti2E fluorescent microscope or and Zeiss LSM880 confocal microscope. Cell and nuclear areas were measured with the ImageJ software.

2.14. Viability Assays

For in vitro cell viability assays 1x10⁴ CT26^5FUR, HCT116^5FUR or HT29^5FUR cells were seeded on 96-well plates and incubated for 36 h at 37°C at a 5% CO₂ atmosphere with increasing doses of the individual drugs and combinations of selected doses with 5 µM 5-FU. Controls were treated in the same way with corresponding vehicle to each drug. After the incubation time, treatments were removed and an MTT assay (Sigma Aldrich) was performed as described before [47]. Viability was expressed as the percentage of control untreated samples. The concentration of drugs that decreased cell proliferation by 50% (IC50) as compared to controls was calculated with GraphPad Prism v8.0 (GraphPad Sofware, La Jolla, CA, USA). When indicated, viable cells were counted using Trypan Blue (0,4%; Sigma), and viability was expressed as the percentage of control untreated samples.

2.15. Immunoblotting

CRC cells were seeded in 60 mm plates until they achieved 80% confluence. Cells were lysed with RIPA buffer containing protease inhibitor cocktail (Roche Diagnostics, Mannheim, Germany) and then detached from the plates to collect protein extracts. Protein levels were quantified using the Lowry assay, and 40 µg of protein per sample was denatured for 5 minutes at 95 °C in SDS-PAGE sample buffer before being loaded onto gels. Samples were resolved on 8% SDS–polyacrylamide gels and transferred to PVDF membranes (Amersham Hybond P, GE Healthcare Life Sciences) for 90 minutes with constant current at 4 °C. Membranes were blocked with 1% (w/v) BSA in TBS-Tween (50 mM Tris, 150 mM NaCl, 0.05% Tween, pH 7.5) for 60 minutes at room temperature and then incubated with the specified primary antibody overnight at 4 °C. The primary antibodies used included β-catenin (BD Biosciences, 1:1000 dilution) and E-cadherin (BD Biosciences, 1:1000 dilution). After three washes with TBS-Tween to remove the primary antibodies, membranes were incubated with the appropriate peroxidase-conjugated secondary antibody (BioRad, 1:5000 dilution) for one hour at room temperature. Detection was carried out using chemiluminescence (Bio-Lumina; Kalium Technologies, Argentina) and imaged with a Licor C-Digit Blot Scanner (LI-COR Biosciences) according to the manufacturer's instructions. Quantification was performed by densitometry using Image J software. Cropped images are displayed in the main figures, with full-length membranes shown in the supplementary information.

2.16. GTP-Rac1 Pull-Down Assay

GST-Pak1 pull-down experiment to determine Rac1 activity in cells was performed as indicated [48]. Briefly, exponentially growing cells were harvested and resuspended in lysis buffer containing 20 mM Tris- HCl [pH 7.5], 150 mM NaCl, 5 mM MgCl₂, 0.5% Triton X-100, 10 mM beta-glycerophosphate, 1 mM DTT, Cømplete (Roche), and 10 µg of GST fusion protein containing the Pak1 Rac1 binding domain (GST-Pak1 RBD, bacterially-expressed). After incubations for 10 min on ice, cell lysates were pre-cleared by centrifugation at 14,000 rpm for 10 min at 4°C and supernatants incubated with glutathione-Sepharose beads (GE Healthcare Life Biosciences) for 1h at 4°C under gentle rotation. After extensive washes in lysis buffer, protein complexes were released by boiling in SDS-PAGE sample buffer, separated electrophoretically, transferred onto nitrocellulose filters, and analyzed by immunoblotting using an anti-Rac1 antibody (BD Biosciences, 1:1000). GTP-Rac1 levels were quantified with the ImageJ analysis software using as normalizing control the total levels of Rac1 found in each cell lysate.

2.17. Animal Studies

8 weeks-old BALB/c female mice were obtained from the School of Veterinary Sciences at the National University of La Plata and treated in accordance with the Canadian Council on Animal Care and ARRIVE guidelines. Animals were maintained in the CIPReB facilities (Centro de Investigación y Producción de Reactivos Biológicos, Medicine School, National University of Rosario). 1x10⁶ CT26^5FUR viable cells were resuspended in PBS (100 μl) and injected subcutaneously into the right flank of each animal. For all experiments, mice were distributed and treated as follows: control, consist in a daily intraperitoneal (ip) administration of vehicle (1% absolute alcohol solution); ivermectin, ip administration of 2 mg/kgBW/day (dissolved in absolute alcohol and then diluted in water to the final concentration); 1A-116, ip administration of 5 mg/kg BW/day in sterile water; 5-FU, ip administration 20 mg/kgBW/week in sterile water; iver + 5-FU, ip administration of ivermectin and 5-FU treatments; 1A-116+5-FU, ip administration of 1A-116 and 5-FU treatments.

Animals were periodically weight and checked for changes in skin, fur, eyes, secretions, excretions and autonomic activity (lacrimation, pilo-erection, unusual respiratory pattern, movement, etc).

Tumor volumes were calculated as V = 0.4ab², where a is the measurement of the tumor along its longest axis and b its shortest. When any of the groups reached the ethically permitted tumor volume, the animals were euthanized. Lungs, spleens and tumors were removed, fixed and stained with hematoxylin-eosin for histological evaluation and detection and counting of micrometastases in an Olympus BX40 microscope.

For intrasplenic inoculation of cells, mice were anesthetized by intraperitoneal injection of acepromazine/ketamine/midazolam (50mg/kg, 100mg/kg and 50mg/kg, respectively). A small incision was made to access the spleen and allow injection of 1x10⁶ CT26^5FUR viable cells resuspended in PBS (100 μl). Two days after surgery, animals were distributed in groups as described and treatments were initiated. Three weeks after injection, animals were sacrificed, for collection and weighting of spleens and livers.

For subcutaneous tumor development, 2 independent experimental rounds were performed. For the first round N = 4/group; for the second round, N=6/group. For intrasplenic injection the number of animals used was N=6/group.

2.18. Statistical Analysis

Statistical analyzes were carried out using the GraphPad Prism 8.0 software (GraphPad Software, Inc., La Jolla, CA, USA). Single comparisons between two groups were performed with the Student’s-test, whereas for multiple comparisons the ANOVA followed by Tukey’s multiple comparisons post-test was employed. Correlation between the two variables was assessed with Spearman’s rank correlation coefficient (r). For differential gene expression analysis, statistical significance was tested with the Student’s t-test followed by a False Discovery Rate (FRD) correction with Benjamini–Hochberg procedure. Survival analysis was implemented according to the Kaplan–Meier analysis and log-rank test. Overall survival (OS) was defined as the time between the date of surgery and the date of death or the date of the last follow-up. In all cases, p-values less than 0.05 were considered statistically significant and were marked with an asterisk as follows: *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001.

3. Results

3.1. Identification of Differentially Expressed Genes and Key Pathways in Recurrent CRC after 5-FU Monotherapy

To identify DEGs and key pathways associated with 5-FU resistance, raw data of the datasets GSE39582 and GSE81653 were downloaded from GEO database. Each dataset was normalized using the RMA method (Figure S1A,B). Dataset GSE39582 contained a total of 585 samples, 82 of which met the requirements described in Section 2.1. Dataset GSE81653 contained 593 samples, 192 of which met the requirements. The selected patients from each study were divided into patients with and without tumor relapse groups. Before performing the analysis, probes with intensity values close to chip background were filtered and discarded (the proportion of probes removed is shown in Figure S1C,D).

In order to explore the putative association between recurrence and clinical parameters, correlation was assessed using Fisher's exact test. The results indicated that 5-FU resistance was not associated to age, TNM stage (Tumor, Node, Metastasis staging system), or the main mutations described for CRC (KRAS, TP53 and BRAF; Table S5). As no correlation was noted, we performed enrichment analysis of gene sets associated to recurrent and non-recurrent phenotypes for both datasets. Through this approach, we observed that coordinate expression of genes grouped in categories linked to cell adhesion and migration was associated to a recurrent phenotype, while a good response to 5-FU monotherapy correlated with immune system and complement activation (Figure 2A,B and Tables S6-S7).

In parallel, we aimed to identify DEGs associated to the 5-FU-resistant phenotype. To select DEGs in each dataset independently we employed four methods, two exploratory (fold change [FC] and unusual ratio [UR]) and two non-parametric hypothesis tests (RankProd [RP] and the Significance Analysis of Microarrays [SAM]). Cut-off criterions for DEGs selection were described in Materials and Methods. DEGs lists were compared to extract common genes detected by at least two methods in both studies. The number of genes obtained by each method in each dataset was represented by a Venn diagram (Figure 2C,E). A total of 388 upregulated and 39 downregulated common DEGs were selected (Table S8). In order to detect pathways enriched in recurrent phenotype, we performed an over-representation analysis (ORA). The results of this analysis indicated that the main upregulated pathways were associated with signaling mediated by Rho GTPase proteins and their effectors (Figure 2D), while among the downregulated pathways were found those related to cytokine receptors, receptors binding to peptide ligands and the immune response (Figure 2F).

Transcription factors are key regulators of biological processes that function by binding to genes regulatory regions. Each transcription factor recognizes a collection of DNA sequences or binding sites that can be represented as motifs. Motif characterization is important for understanding the regulatory functions of transcription factors shaping gene regulatory networks. To go further in characterizing expression changes associated to resistance, we look for common transcriptional factors binding sites present at the promotor regions of upregulated genes, finding that the most relevant factors include Serum Response Factor (SRF), Myocyte Enhancer Factor 2C and 2A (MEF2C, MEF2A), Msh Homeobox 2 (MSX2), and common motives were found to FOXO1, HMGA2, HMGA1, FOXA2. FOXA1. SOX10, FOS, JUN, HLTF and RUNX3 (Figure 3A). Genes specifically regulated by each transcription factor were extracted to perform a functional enrichment analysis (P < 0.05), through a word cloud graph, which indicates that the most relevant cellular processes modulated are related to signaling, GTPases and E-cadherin (Figure 3B).

To visualize the most relevant relationships between pathways and functions associated to the 5-FU-resistant phenotype, we constructed an enrichment map that indicates over-representation of two particular clusters: Rho GTPases and TGFbeta signaling (Figure 3C,D). 5-FU-resitance associated DEGs are highly connected genes by 20 hub genes, as visualize in Figure 3E.

3.2. Determinants of Resistance to 5-FU-Based Therapies

To go further in the search for determinants to resistance in CRC, we decided to extend our approach to databases of patients receiving therapies based on 5-FU, including FOLFOX and FOLFIRI. This approach allowed us to increase our samples to 567 to perform more robust analyses. As in the case of 5-FU monotherapy, no evident correlation was observed between resistance and available clinical parameters (Table S9). Unbiased GSEAs associated with the transcriptome of recurrent phenotype following 5-FU-based chemotherapies using Reactome, GO, and KEGG databases indicated a positive enrichment for Rho GTPases activate PKNs and Formation of beta-catenin:TCF transactivating complex, and negative enrichment for immunoregulatory interactions between lymphoid and non-lymphoid cell and complement cascade (Figure 4A, Table S10). Moreover, after selecting the 5000 most explanatory genes in the recurrent phenotype after 5-FU-based chemotherapies by the Recursive Feature Elimination algorithm (RFE), and a chemical and genetic perturbations database (MSigDB) we found a positive enrichment of gene sets associated with resistance to other therapies generally used to treat tumors other than CRC such as gefitinib, tamoxifen, gemcitabine and radiation (Figure 4B, Table S11), indicating that some of the identified pathways could be orchestrating general resistance mechanisms. Surprisingly, gene sets more present in this 5000 subset were related to neutrophil degranulation and Rho GTPases, among which Rac1 seems to be the most relevant (Figure 4C). Consistently, through pull-down experiments in CRC cells, we observed a significant increase in Rac1 activity associated with 5-FU resistance (Figure 4D). As it was the case for 5-FU monotherapy, we found binding motifs for SRF, JUN and FOS transcription factors present at the promotor regions of upregulated genes (Figure 4E) associated to 5-FU-based therapies resistance.

3.3. Selection of Drugs to Overcome 5-FU Resistance in CRC

In order to overcome 5-FU resistance, we search for compounds with the ability to reverse the expression of genes and pathways identified in our previous analysis by LINCS. This tool allows entering 150 upregulated and 150 downregulated genes. For 5-FU monotherapy, we used all common downregulated genes as they were less than 150. To reduce the number of upregulated genes to 150, we selected DEGs detected by three or more statistical methods (in both databases) and those with the higher association (with a lower P value) with overrepresented pathways in functional enrichment analysis (Table S12). The platform detected similar and opposite expression profiles in nine cell lines, including HT29.

Through this approach we obtained a list of compounds that could potentially reverse 5-FU resistant phenotype (Table S13). Some of these compounds are currently used in the clinic combined to 5-FU, such as irinotecan [49,50] or to treat other types of recurrent cancers such as topotecan [51,52]; but we also found other compounds such as a norepinephrine inhibitor (amitriptyline), MEK inhibitors, a COX inhibitor, ivermectin and the Rho-associated kinase inhibitor (Rockout) among others. Accordingly, cell perturbations that could reverse the expression profile associated to 5-FU resistance determined by the same platform include inhibition of Raf, MEK, PKC or Src and loss of function of Rho GTPases activating proteins (Table S14).

3.4. Ivermectin and the Rac1 Inhibitor 1A-116 Restore 5-FU Sensitivity to 5-FU Resistant Cells

As drug repositioning provides cheaper, effective, and safe drugs with fewer side effects and fastens drug development, we focused on irinotecan, amitriptyline and ivermectin to confirm by cell studies their potential to sensitize 5-FU-resistant CRC cells. In parallel, as our studies highlighted a prominent role of the Rho GTPase Rac1 on 5-FU-resistance modulation, we decided to explore the effects of 1A-116, a novel Rac1 inhibitor previously reported to be promising to treat other cancer types [46,53,54,55]. To characterize the effect of these drugs on resistance, we used three 5-FU resistant cells lines generated as described before [47]. IC50 values for each selected drug were calculated on sensitive and resistant cell lines (Table S15).

To evaluate the potential effect of 1A-116 and selected repurposing drugs on reversing resistance of CRC cells, we treated resistant cells with 5-FU in the presence of each drug. For all the assays, we selected doses for each drug with a minimal statistical effect on cell viability as evaluated by MTT-based assays (Figure S2 and Table S15). We found that incubation of cells with irinotecan, ivermectin and 1A-116 improved responses to 5-FU in all cell lines tested (Figure 5A), while amitriptyline only had a partial effect on HT29 and CT26 resistant cells. For 1A-116 treatment, we confirmed the inhibition of Rac1 activation by this compound and its effect on reversing 5-FU resistance by counting viable cells after treatment (Figure S3A,B).

In vitro generation of 5-FU-resistance is associated to cell morphological changes as it was previously described for several cell lines including CRC [56,57,58,59]. Previous reports also described that resistance to 5-FU in CRC cells promotes the loss of epithelial markers [58,59]. To confirm the link between resistance acquisition and EMT, we analyzed the expression of epithelial markers in CRC cells, validating the loss of E-Cadherin and beta-catenin expression in resistance cells (Figure 5B,C). To assess morphological changes associated to 5-FU resistance we performed nuclear, cytoskeletal and microtubular staining of control and 5-FU-resistant cells (Figure S3D,E). In order to quantify changes in morphology, we measured nuclear and cell areas for each condition, noting a significant increase in both parameters for resistant cells (Figure 5E-G). Indeed, the relationship between nucleus/cytoplasm decreased in resistant cells and cytoskeletal architecture changes became evident. Interestingly, treatment with the Rac1 inhibitor 1A-116 was enough to reverse the expression of epithelial markers (Figure 5B,D) and morphological changes back to control-like values (Figure 5E-G and S3B,C) suggesting a critical role for Rac1 activation on triggering events leading to 5-FU-resistance in CRC.

3.5. Rac1 Inhibitor 1A-116 Reduces the Growth of CRC Resistant cells, Sensitizes Them to 5-FU and Prevents Metastasis Development

Once confirmed by in vitro assays the possibility to reverse resistance to 5-FU we moved to an in vivo assay with CT26^5FUR cells. For this experiment, we selected one drug in repositioning, ivermectin, and the Rac1 inhibitor 1A-116. Tumor growth kinetics indicated that 5-FU or ivermectin alone had not statistical effect on resistant CRC cells growing (Figure 6A,B). Surprisingly, Rac1 inhibitor 1A-116 administration statistically affected the growth evolution of tumor cells. Interestingly, both ivermectin and 1A-116 were able to induce sensitivity to 5-FU in resistant CT26^5FUR cells (Figure 6A,B and S4A,B). No signs of toxicity were associated with treatments as evaluated by general animal behavior and weight (Figure 6C and S4C) and the measurement of metabolic parameters (Figure S4D).

At the end of the experiment, we collected tumor, spleen, liver and lungs for histological observation. We did not observe histological differences between tumors by H&E staining (Figure S5), but it was evident that animals treated with 1A-116 or combination of 1A-116 with 5-FU did not develop splenomegaly (Figure 6D and Table 1). Microscopic visualization of histological sections indicated that combined treatment with 1A-116 and 5-FU reduced liver and lung metastasis development in mice (Figure 6E,F and Table 1) suggesting a role for Rac1 in CRC tumor growth, resistance and metastatic dissemination. To further confirm the antimetastatic potential of Rac1 inhibition on 5-FU resistant CRC we performed an intrasplenic cells injection, noting that 1A-116 treatment drastically reduced spleen tumor formation and resistant cells dissemination to the liver (Figure 6G-J).

3.6. Generation of Prognostic Signatures to Predict 5-FU-Based Therapies Resistance

Finally, we used our lists of genes associated to resistance to obtain by LASSO regression analysis a minimal group of predictive genes that could anticipate responses to both 5-FU monotherapy and 5-FU-based therapies (Figure 7A and Figure S6) and evaluated their predictive efficiency by generation of ROC curves (Figure 7B). Unfortunately, there were not many common genes between the different predictive lists, but all of them displayed an acceptable predictive efficiency.

As it was shown that Rho GTPases were important markers associated to resistance, we selected by LASSO the most important Rho GTPases-related genes to predict recurrence in the 5-FU monotherapy datasets and combined treatments merged matrix, and then analyzed their recurrence predictive efficiency for each signature by generation of ROC curves (Figure 7C). ROC curve analysis demonstrated the good performance of the established model, with similar predictive efficiency to those based on DEGs. Interestingly, some of these predictive genes were associated to overall CRC patients’ survival (Figure 7D).

To understand the relationship between genes predicting resistance, we constructed a protein-protein interaction network. Using the STRING database, we identified the main interacting proteins shared among these genes (Figure 7E) and extracted the top ten hub genes from this network (Figure 7F). Interestingly, top ten genes are GTPases Rac1, RhoA, and RhoB, their GEFs (Guanine nucleotide Exchange Factors Vav1, Vav2, Vav3 and Tiam1), their effector Pak1, phosphoinositide-3-kinase regulatory subunit (PIK3R2) and beta-catenin (CTNNB1). Taken together, these data indicate a prominent role for Rho GTPases, particularly Rac1, in modulating resistance to 5-FU-based therapies in CRC.

4. Discussion

CRC is the third most commonly diagnosed type of cancer worldwide. Early diagnosis, adenomas removal during screenings and improved CRC treatments have reduced morbidity and mortality rates. However, CRC incidence is increasing in middle- and low-income countries, and early-onset CRC is also emerging, positioning CRC as a growing global public health challenge [60].

Since their discovery, 5-FU-based chemotherapies have been commonly used to treat CRC. However, resistance to treatment has greatly affected 5-FU clinical use. As CRC is a heterogeneous disease characterized by different genetic scenarios [61], a CRC molecular comprehensive characterization could shed light on the understanding of resistance mechanisms and improve therapies.

In this work, we addressed the identification of genes and pathways associated with development of resistance to 5-FU and 5-FU-based treatments in CRC. Comparison of gene expression datasets revealed the enrichment of pathways previously reported to be associated with resistance such as “TGFbeta signaling” and “SMAD2/3 activation by TGFbeta” [62,63,64]. TGFbeta is a cytokine involved both in physiological and pathological processes. Canonical signaling is mediated by SMAD transcription factors, but TGFbeta “non-canonical” pathway involves activation of proteins such as MAPK, PI3K/AKT and small GTPases as Rac1 and RhoA [65]. Categories and genes related to cell stemness [66,67], EGFR signaling [68,69], apoptosis [70], DNA repair [71], and interaction between tumor cells and the immune system [72] were also present in our analysis.

When we looked for potential compounds reversing 5-FU resistance-associated gene expression, we identified molecules belonging to different categories. Some of them were involved in modulation of hormonal messages (dopamine and serotonin receptor agonists and antagonists), selective estrogen receptor modulation (this category includes tamoxifen and raloxifene, widely used in breast cancer) and adrenergic receptor inhibition such as phentolamine. Additionally, we found anti-inflammatory compounds, mainly COX inhibitors. We focused on irinotecan, amitriptyline and ivermectin, three drugs under repositioning and selected ivermectin to our in vivo approach. As predicted by our analysis, ivermectin was able to re-sensitize CRC cells to 5-FU monotherapy, but no significant metastasis prevention was associated to this treatment. Noteworthy, the antitumor mechanism of action proposed for ivermectin is associated to P21(Rac1)-Activated Kinase 1 (PAK1) function [73,74,75].

Interestingly, the most represented pathways in all datasets were directly or indirectly related to Rho GTPases activation. Indeed, activation of Rho GTPases may be under control of the Serum Response Factor as many upregulated genes in recurrent phenotype contained serum response elements in their promoter sequences. Many studies have identified SRF as a central agent in the development of multiple types of cancer, which has classified it as a potential biomarker and therapeutic target, especially for cancers with a poor prognosis. SRF controls the expression of cytoskeleton, morphogenesis and cell migration genes, but also SRF-MRTF complex activity also responds to Rho GTPase-induced actin changes, thereby coupling cytoskeletal gene expression to cytoskeletal dynamics [76,77,78]. Accordingly, it was recently described that active Rac1 modulates SRF/MRTF, which initiates a switch to a mesenchymal-like state characterized by therapy resistance in melanoma [79].

Rho-GTPases regulate a variety of important cellular activities, such as cytoskeletal remodeling, cell adhesion, cell movement, vesicle transport, angiogenesis, and cell cycle regulation [76,80,81,82,83]. Rho-GTPases are generally described as "molecular switches" because they fluctuate between their active conformation attached to GTP and the inactive conformation bound to GDP. The activation of the “molecular switch” is controlled by guanine nucleotide exchange factors (GEFs), which stimulate the release of GDP bound to the inactive form and promote combination with GTP [84]. The inactive state of Rho-GTPase is maintained by inhibitory molecules of the guanine nucleotide dissociation (GDIs) and GTPase activity activating proteins (GAPs) [84]. Through changes in Rho-GTPases protein levels, their activity or their effector proteins, abnormal signaling could contribute to different steps of cancer progression, including proliferation, survival, invasion and metastasis [22,85].

Rac1, RhoA and Cdc42 are the three classical members of Rho-GTPase family, being Rac1 the one that has received increased attention [82]. Rac1 is widely expressed in tissues, and is considered a regulatory factor related to cell movement and invasion. Rac1 is highly expressed and overactivated in many tumor types and it had lately been related to resistance to therapy in several reports [79,86,87,88,89,90].

Based on our data and the literature, Rac1 could be engaging cellular mechanisms leading to 5-FU-based resistance in CRC. To confirm our hypothesis, we treated 5-FU-resistant CRC cells with Rac1 inhibitor 1A-116 noting that doses not compromising viability of CRC cells were enough to overt morphological changes associated to resistance, and re-sensitize resistant cells to 5-FU. Moreover, Rac1 inhibition restored the expression of epithelial markers to 5-FU-resistant cells, previously characterized in a mesenchymal-like state associated with therapy resistance. Administration of Rac1 inhibitor to mice bearing CRC tumors reduced tumor growth, sensitized resistant tumors to 5-FU monotherapy, and decreased metastasis development. Interestingly, tumors derived from animals treated with Rac1 inhibitor showed an increased number of immune cell infiltrates (data not shown), suggesting that Rac1 inhibition could be acting through different mechanisms, including immune escape mediated by tumor microenvironment [91].

Additionally, our data suggest that Rac1 inhibition could be an important strategy to overcome resistance to therapy in different cancer types. Indeed, Rac1 modulated pathways could be playing essential roles in developing resistance to different therapeutic approaches such as endocrine therapies, targeted therapies and radiotherapy, as our enrichment analysis indicate common profiles between resistance to 5-FU-based therapies and gemcitabine and gefitinib in non-small cell lung cancer [92,93], dasatinib for breast, lung, and ovarian tumors [94], tamoxifen in estrogen receptor-positive breast cancer [95], and postradiation tumor escape of CRC [96].

Altogether our data point to Rac1 as a potential target to overcome therapy-resistance for CRC and other types of tumors, and suggest that Rac1 inhibitor 1A-116 could represent a good therapeutic agent to overt CRC resistance for 5-FU-based therapies.

5. Conclusions

CRC is the one of the third most diagnosed type of cancer worldwide, and the second in terms of mortality. Therapies to treat CRC are often associated with the development of resistance, so initially responding tumors became resistant to treatment. Therapies based on 5-FU have been used in the clinics since 1950, but almost half of patients develop therapy-resistance.

Our findings indicate that inhibition of Rac1 activation by the compound 1A-116 reverses 5-FU resistance in CRC cells, decreases tumor growth and prevents metastasis development, suggesting a therapeutic application of 1A-116 for the treatment of therapy-resistant CRC. Further studies are needed to fully understand Rac1's role in CRC progression and therapy resistance.