Metabolism-dependent secondary effect of anti-MAPK cancer therapy on DNA repair Open Access

Abstract

Amino acid bioavailability impacts mRNA translation in a codon-dependent manner. Here, we report that the anti-cancer MAPK inhibitors (MAPKi) decrease the intracellular concentration of aspartate and glutamate in melanoma cells. This coincides with the accumulation of ribosomes on codons corresponding to these amino acids and triggers the translation-dependent degradation of mRNAs encoding aspartate- and glutamate-rich proteins, involved in DNA metabolism such as DNA replication and repair. Consequently, cells that survive MAPKi degrade aspartate and glutamate likely to generate energy, which simultaneously decreases their requirement for amino acids due to the downregulation of aspartate- and glutamate-rich proteins involved in cell proliferation. Concomitantly, the downregulation of aspartate- and glutamate-rich proteins involved in DNA repair increases DNA damage loads. Thus, DNA repair defects, and therefore mutations, are at least in part a secondary effect of the metabolic adaptation of cells exposed to MAPKi.

Graphical Abstract

Open in new tabDownload slide

Introduction

Some amino acids correspond to several codons, bringing a redundancy into the genetic code, which plays a role in gene product expression level. Indeed, the nature of mRNA codons affects the dynamics of ribosomes, through the time required to ‘decode’ each codon, with consequences on translation efficiency, protein folding and mRNA co-translational degradation through the translation-dependent mRNA decay (TDD) pathway (1–7). For instance, sub-optimal codons or clusters of rare codons modulate gene product expression levels in a gene- and species-selective manner (1,3,4,7–9). At least in humans, sub-optimal and rare codons are more frequently A/T-ending codons, while optimal codons are more frequently G/C-ending codons likely because of the more unstable interactions between A/T-ending codons with their cognate anti-codons as compared to G/C-ending codons (4,10–14).

While codon usage is an intrinsic parameter of gene products with respect to their translatability and stability, there are many extrinsic parameters that modulate codon-depending effects in a cell type- and context-dependent manner. Among the extrinsic parameters, enzyme-dependent biochemical modifications of anti-codons have been reported to modulate the interactions between codons and anti-codons (15,16). For example, the cytosolic thiouridylase 2 (CTU2) enzyme that biochemically modifies the tRNA wobble uridine – thereby affecting the decoding of some A-ending codons, such as the AAA and GAA codons – is required for the efficient translation of a subset of mRNAs that promote survival and resistance to therapy of BRAF-mutated melanoma cells (17,18). In addition, variation in the expression levels of different classes of tRNAs can change the codon-dependent effect on gene product expression level in a cell type- and context-dependent manner (19,20). For instance, proliferative cells express higher levels of tRNAs corresponding to A/T-ending codons than differentiated cells, leading to the differential expression levels of gene products enriched in either A/T- or G/C-ending codons (21,22). Along the same line, the expression levels of aminoacyl transferases, which load amino acids onto tRNAs, can modulate the codon-dependent effects on gene product expression levels (23,24). An example is given by the leucyl-tRNA synthetase (LARS), which increases the selective loading of tRNA-Leu^CAG isoacceptor and thereby affects the translatability of mRNAs containing the CAG codon (25).

The bioavailability of amino acids is of particular interest as an extrinsic parameter that modulates codon-dependent effects on gene product expression levels. Indeed, the dynamics of ribosomes depend in part on the intracellular concentration of loaded-tRNAs, which itself depends on the intracellular concentration of amino acids. Thus, the translation of mRNAs requiring an amino acid whose bioavailability decreases can be impacted as shown for numerous amino acids, including non-essential ones such as glutamine (5,7,26–29). As a consequence of the link between amino acid bioavailability and translation, the synthesis of proteins with amino acid composition biases depends on the cell metabolism (23,25,28,30,31). For example, a large amount of proline can be produced only under certain metabolic conditions, which therefore determines whether some proline-rich proteins of the extracellular matrix (e.g. collagen) are produced (29,32,33). This illustrates how the cellular metabolism – through amino acid bioavailability – is coupled to the nature of cell-expressed proteins on which the cell phenotype depends.

The coupling between cell metabolism and the protein-dependent cell phenotype can be illustrated by the competition between energy production and gene product biogenesis on which depends cell proliferation because some amino acids like glutamine (Gln), glutamate (Glu) and aspartate (Asp) are at the crossroads between several metabolic pathways and the gene expression process. Indeed, the carbon skeleton of these amino acids can either be degraded and end up in the production of energy, or ‘recycled’ to synthesize other amino acids and nucleotides (and therefore gene products) (34–39). As a consequence, some amino acids may either be used by the cells to produce energy through their complete degradation or be used for the synthesis of large amount of gene products as during cell proliferation. This may explain why such amino acids play a particularly important role in cancer cells, which have a high proliferation rate that requires gene product synthesis but are in a resource-impoverished micro-environment as a consequence of cell proliferation (35,37,40–45). Accordingly, cancers cells are often addicted to certain amino acids, such as Gln in melanoma cells (37,41,44,45). The link between cell metabolism and gene expression–dependent cell phenotypes could have consequences in cancer cells exposed to anti-cancer agents, such as MAPK-inhibitors (MAPKi) used to treat melanoma, since these molecules modify the cancer cell metabolism (46–51). In other words, anticancer therapies may impact cellular phenotypes because of metabolic-dependent effects on gene product expression levels.

Here, we report that MAPKi treatment of a melanoma cell line downregulated mRNAs that encode proteins enriched for certain amino acids, including Glu and Asp, whose intracellular concentration decreased in MAPKi-treated cells. Interestingly, MAPKi-downregulated mRNAs encoded proteins involved in cell proliferation and DNA repair, two classes of proteins that are globally enriched in Glu and Asp residues. In line with this observation, MAPKi-treated cells show DNA repair defects. Our results support a model in which the metabolic-dependent effects of MAPKi therapy combined with modification of the expression level of tRNA-charging enzymes could result in secondary defects of DNA repair, which could increase the probability of genetically-adapted cancer cells to emerge after MAPKi therapy.

Materials and methods

Cell culture and persister cells

The human melanoma A375 cell line (ATCC) was cultured at 37°C and 5% CO₂ in Dulbecco's modified Eagle's medium (DMEM; Gibco) supplemented with 10% FBS, 2 mM glutamine and penicillin–streptomycin. Cells were split at 80% confluence, three times a week. MAPKi-treated cells were cultured in a medium containing 1 μM vemurafenib/cobimetinib (with 500 nM vemurafenib and 500 nM cobimetinib) (Euromedex) for 24 h. Glutamine-depleted cells were cultured in a glutamine-free medium for 24 h. Persister cells were obtained after being cultured for 72 h in a medium containing 1 μM vemurafenib/cobimetinib in three experimental batches (A–C). Batch A was harvested after three days and batches B and C were cultured in drug-free DMEM for two and nine additional days. The human M249 cells were cultured in Roswell Park Memorial Institute medium (RPMI 1640 Eurobio Scientific) supplemented with 10% FBS, 2 mM glutamine and penicillin–streptomycin.

siRNA transfections, cell harvesting, RNA extraction and qRT-PCR

Two different siRNAs (Merck, Supplementary Table S2) that target constitutive exons with minimum predicted off-targets were designed and pooled together. Cells were reverse-transfected with lipofectamine RNAiMAX (ThermoFisher), following the manufacturer's instructions. After transfection of 6 well plates (200 000 cells/well), cells were washed twice with ice-cold PBS, scraped with 1 ml of PBS and pelleted by centrifugation (500 × g for 1 min at 4°C). Cells were then suspended in 1 ml lysis buffer (10 mM Tris–HCl pH 7.5, 5 mM MgCl₂, 100 mM KCl, 1% Triton X-100) and incubated on ice for 10 min. Cellular lysates were centrifuged (1000 × g for 10 min, 4°C), and RNA was extracted using TRI Reagent (Sigma), following the manufacturer's instructions. For RT-qPCR, 1 μg of extracted RNA was retro-transcribed using the Maxima First Strand cDNA Synthesis Kit (ThermoFischer), following the manufacturer's instructions. qPCR reactions were run in triplicate on a LightCycler 480 (Roche) in 10 μl reactions (see Supplementary Table S2 for primer sequences). The absence of siRNA off-target effects was tested by transfecting individual siRNAs; effects on mRNA targets were then compared by RT-qPCR to the effects of the mixed siRNAs (Supplementary Figure S4e and g). The amino acid intracellular concentration was measured using cytoplasmic extracts of the A375 cell line from 5 × 10⁷ cells and was performed by the AltaBioscience and Xell companies. The intracellular concentration of Glu in the M249 cell line was performed using the glutamate assay kit from Abcam (ab83389).

Gene annotation

Bioinformatics analyses were performed using the GRCh38.p13 assembly and NCBI’s annotation. Only genes with at least one coding sequence (CDS) and one start and stop codon were kept. Merged genes, and CDS that overlap several exons or an ambiguous coding frame were filtered out. A total of 19 143 coding genes and 196 652 CDS were selected. In the following analyses, the mRNAs are the concatenation of the CDS of a gene and a gene was associated with only one mRNA.

RNA-seq and QuantSeq

From 6-well plates, 10⁶ cells were used to prepare total RNA, and 500 ng of RNA were used for RNA-seq libraries that were prepared and sequenced by Novogene (rRNA depletion library preparations and sequencing was performed on Novaseq 6000 2 × 150). QuantSeq libraries were prepared using the QuantSeq 3′ mRNA-Seq Library Prep Kit (Lexogen). 500 ng of RNA was spiked-in with 1 μl of a 1:100 dilution of ERCC Spike-In Mix (ThermoFisher) prior to library preparation. QuantSeq libraries were then quantified, pooled and sent for sequencing at Novogene (Novaseq 6000 2 × 150). Only fastq files containing forward reads were used for QuantSeq analyses, while the entire pair of fastq files were used for RNA-seq analyses. Adapters were removed from raw reads and trimmed using fastp version 0.20.1 (52) with the following parameters: qualified_quality_phred 30 -l 25 –detect_adapter_for_pe) (RNA-seq) and qualified_quality_phred 20 -3 –cut_tail_window_size 10 –adapter_sequence = AGATCGGAAGAGCACACGTCTGAACTCCAGTCA –adapter_sequence = AAAAAA (QuantSeq). Reads were mapped against the human genome GRCh38.p13 with HISAT2 version 2.2.1. Reads were counted on exons with htseq-count version 0.13.5. Unexpressed genes (genes with raw read counts < 2 across all conditions) were filtered out only for RNA-seq data. Differential expression analysis was performed with the DESeq2 package version 1.34.0 using the option lfcThreshold = 0.585 (only RNA-seq data). mRNAs in the RNA-seq dataset with a DESeq2 corrected P-value ≤0.05 and an average DESeq2 normalized expression ≥10 across conditions (basemean ≥ 10) were considered to be significantly differentially regulated.

TDD experiments and analysis

TDD monitoring was performed as previously described in (53). Briefly, 2 × 10⁵ cells were plated into 6-well plates, and the culture medium was removed after 24 h and replaced with fresh medium containing 1 μM vemurafenib/cobimetinib (with 500 nM of each) or DMSO. After another 18 h, cells were treated with fresh cycloheximide (100 μg/ml) or DMSO for 5 min and then treated with tryptolide (25 μM) or DMSO. This was prepared in four identical batches (A–D): batch A was harvested immediately upon ± tryptolide treatment (T0 samples), and the batches B, C and D were harvested after 3 (T3) and 5 (T5) hours of treatment, respectively. Computation of the TDD index was performed using pre-processed reads, mapped and counted with htseq-count (see above). From the raw count tables obtained with htseq-count, the CPM was computed for each gene. Then, we searched for genes with a stable expression in MAPKi and DMSO conditions for normalization. A gene was considered as stable if (i) its CPM count was >0.2 in initial condition, (ii) its CPM count after 3 or 5 h of transcription inhibitor treatment was at least 10% greater than its CPM at initial condition and (iii) its CPM at initial condition was greater than 10% of its CPM after 3 or 5 h of transcription inhibitor treatment. Only ‘stable genes’ in all replicates were kept. Stable genes were used as normalization factors in DESeq2 package to normalize reads counts. The TDD index was next computed for an mRNA produced by a gene G using the following formula:

where |${R_0}$| is the normalized number of G reads at 0 h (initial condition), |${R_{Tci + Tli,t}}$| is the normalized number of G reads after t hour of exposition to a transcription (Tci) and a translation inhibitor (Tli) and |${R_{Tci,t}}$| is the normalized number of G reads after t hour of exposition to a Tci. The TDD index was computed, for each replicate, at T3 and T5 for cells treated with MAPKi and DMSO. The TDD index according to the condition (MAPKi or DMSO), the time (T3 or T5) and the replicate was modeled using a linear model (in R function lm) for each mRNA. With these models, a Student's test was computed to determine if the TDD index of a mRNA was significantly affected (P-value ≤ 0.05) in response to MAPKi as compared to the DMSO control condition. Finally, mRNAs with a significant TDD index variation were classified as having an increased or decreased TDD index, according to whether their average TDD index (across replicates and T3 and T5 treatment times) was higher or lower, respectively, in the MAPKi condition compared to the DMSO condition.

Ribosome profiling and analysis

Ribosome profiling samples were prepared as described in (54). Briefly, 5 × 10⁷ cells were harvested with PBS and lysis buffer supplemented with 100 μg/ml cycloheximide (Sigma) and 2 mM DTT. The lysate (500 μl) was analysed at 260 nm absorbance to estimate the total quantity of material and then treated with nucleases. For every 5 units of A260 absorbance, 6 μl of MNase (1 mg/ml; Nuclease S7, Roche) was added to the lysate, along with CaCl₂ to a final concentration of 10 mM. Lysates were then incubated at 25°C for 30 min, transferred to ice and then applied to 10–50% sucrose gradients containing cycloheximide (100 μg/ml). After ultra-centrifugation at 35000 rpm for 2h and 40 min at 4°C, gradients were fractionated using a fraction collector, and the fractions containing the digested monosome fragments (80S) were kept. Fractions were supplemented with SDS (to a final concentration of 1%) and then digested with proteinase K (Roche, final concentration of 2 μg/ml) for 45 min at 42°C. Protected RNA fragments were then purified using an acidic phenol–chloroform extraction (Fischer, BP1753I) and precipitated overnight at –20°C with 0.1 × volume of sodium acetate (3 M, pH 5.2), 1 × vol isopropanol, 1 μl GlycoBlue and 10 mM MgCl₂ (to help recover smaller nucleic acids)_. Purified RNA fragments were then 3′-end dephosphorylated using PNK and fractionated on a 10% acrylamide denaturing gel, and the smears of interest (26–32 bp) were cut from the gel and purified. Size-selected fragments were rRNA-depleted by hybridization using RNA probes (55), successively RNAse H- and DNAse-treated and then purified again with a phenol–chloroform extraction before proceeding to cDNA library preparation following the Omniprep Library preparation protocol (56). An adaptator sequence was ligated to RNA, which was then retrotranscribed with barcoded primers. The barcoded cDNAs were size-selected on a 10% denaturing acrylamide gel, purified and then circularized (CircLigase, Lucigen). Amplification with barcoded primers was performed with a few numbers of PCR cycles (5–8) and a high-fidelity polymerase (Q5, NEB). Amplified libraries were size-selected on a non-denaturing 8% acrylamide gel and purified, and their quality and concentrations were assessed using the TapeStation DNA 1000 ScreenTapes. Ribosome profiling OmniPrep libraries were sequenced by GenomEast (HiSeq 4000 1 × 50bp). After removing adapter sequences from raw reads using cutadapt version 2.1 with the parameters -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -u 13 –maximum-length = 40 –minimum-length = 20 -q 28,28 for Ribo-seq data, and the parameters -a AGATCGGAAGAG -g CTCTTCCGATCT -A AGATCGGAAGAG -G CTCTTCCGATCT for RNA-seq data (raw reads, see ‘Total RNA-seq and QuantSeq’ section). A trimming step was then performed using UrQt version 1.0.18 (57), and reads were mapped to GRCh38.p13 genes sequences using HISAT2 version 2.2.1 (58) with parameters –rna-strandness ‘F’ –norc for Ribo-Seq data. Next, alignment files were converted, using deepTools version 3.0.2 (59), to bigWig files containing a count per million (CPM) mapped reads normalized coverage at one nucleotide resolution. A peak calling step was then performed as described hereafter. For a replicate i and a gene G, a normalized coverage cNorm was computed for a test T and a control C condition using the following formula:

where coverage_ribo was obtained from Ribo-seq data, and coverage_rna from RNA-seq data. Nucleotide positions at which no coverage was detected from RNA-seq data were skipped. The difference cDiff_i between normalized coverage in T and C condition was then computed for each replicate i. Thus, for a number N of replicates, a set of coverage CDIFF per replicates was obtained as follows:

The average coverage |$cMean = \{ {cMea{n_0}, \ldots ,cMea{n_p}, \ldots ,cMea{n_{L - 1}}} \}$| at each CDS position p of a gene G of length L was computed between replicates, where cMean_p was computed as follows:

Next, the average coverage MeanCov and standard error StdMean were computed for the gene G, and a coverage threshold T was defined with |$T = MeanCov + ( {StdMean} ) \times 3$|⁠, where MeanCov and StdMean were computed using the formulas:

Each region at which cMean was above T was considered as a peak. Peaks inside CDS with an average coverage below 3 across RNA-seq replicates, were removed. Peaks defined by a region where the average RNA-seq coverage was below 3 were discarded. A score was given for each peak, beginning at a position s and ending at a position e in a gene. Only peaks position with a score above 3 in two replicates were kept. The score was computed using the formula below:

The first analysis was done by taking the MAPKi or glutamine-deprived conditions as the test condition, and the DMSO or untreated condition as the control condition, respectively. For each analysis performed with this method, another was carried out by reversing the control and test conditions.

Randomization tests were performed to test whether a set of peaks P had a codon compositional bias, or whether peptides encoded in peak regions had an amino acid bias. For this, 10 000 sets of control peaks C similar to P (same peak number and peak size) located in CDS were sampled. For each peak of P and C, the frequency of a given feature X (i.e. codon or encoded amino acid) was computed using the formula: |$Freq( X ) = Count( X )/s$|⁠, where Count(X) is the number of X in a peak and s is the total number of codons in this peak. The average frequency of X, Mean_X, was then computed for P and the 10 000 sets of C|$MEA{N_C} = \{ {Mea{n_{C1}}, \ldots ,Mea{n_{C10000}}} \}$|⁠. To calculate an empirical p-value, the number of control frequencies Mean_Ci upper or equal or lower or equal than the frequency Mean_X was determined. The smaller number between these two was then divided by the number of control peak sets (i.e. 10 000). Note that the p-value cannot be lower than 1/10 000 to avoid multiple testing caveats. The P-value was then corrected using the Benjamini–Hochberg procedure. Figure 3E was generated with codons inside coding sequences (CDS) of genes producing transcripts that have at least one ribosome peak. Codons overlapping two different CDS were discarded from the analysis. The frequency of a given codon was computed in the CDS region overlapped by a peak. Starting from the central coordinate of the peak ((end – start)/2 rounded up), the frequency of codons was then computed up to 50 windows of 10 codons with a step of 1 upstream and downstream the peak.

Compositional bias analyses

To test whether the codon content of different sets of genes was different, the frequencies of each codon in genes according to their size and set was modeled with a generalized linear model for the beta distribution with zero inflation (with R glmmTMB function of the glmmTMB package using beta_family(link = ‘logit’) parameter). Then, a Tukey's test (pairwise comparison) for the ‘set’ factor was done (with R emmeans and pairs functions of the emmeans package). Control sets of genes correspond to expressed genes having a mean DESEQ2 normalized expression greater than 10 and not being in other tested sets. The same procedure was applied to test whether the amino acid content of different sets of proteins was different. When several codons or amino acids are displayed in a figure, an additional Benjamini–Hotchberg correction is performed. The relative frequency of a feature X was computed as follows:

where Freq(X)_t is the average frequency of a codon or amino acid X in a test set of mRNAs or proteins and Freq(X)_c is the average frequency of a codon or amino acid X in the set of mRNAs expressed in A375 cells or their encoded proteins.

Functional enrichment analysis and heatmap

Gene ontology (GO) enrichment analysis was performed using DAVID Ontology (60). An annotation file containing GO terms and a gene association file (that links proteins to their most specific GO terms) were downloaded from http://geneontology.org/. A homemade tool was developed to extract all proteins associated with GO:0051301 (cell division) and GO:000628 (DNA repair). Proteins associated with child terms of these GO terms were also treated as belonging to these terms. Only child terms linked to their parents with the qualifiers ‘involved_in’, ‘located_in’, ‘is_active_in’ or ‘part_of’ were considered. In addition, proteins associated with the Uniprot keywords KW-0131 (cell cycle) and KW-0227 (DNA damage) were downloaded from https://www.uniprot.org/keywords/. Only reviewed human proteins were kept. The average frequency of each amino acid in these lists of proteins was calculated using FasterDB (http://fasterdb.ens-lyon.fr/faster/home.pl). Overall, 10 000 sets of proteins were randomly sampled for each list of proteins, and the average frequency of each amino acid for each set was computed. Finally, an empirical p-value was computed for an amino acid X in a given protein list P as:

where k is the number of controls sets with an average frequency of X higher or equal to P, and l is the number of controls sets with an average frequency of X lower or equal to P. For each list of proteins, the P-values were corrected using the Benjamini–Hochberg procedure and then transformed using the following formula: |$T = 1 - {P_{adj}} \times s$|⁠, where |${P_{adj}}$|is the corrected p-value and s= 1 if k > l; otherwise,s = − 1.

Transcriptional mutagenesis

Mapped reads files (see ‘RNA-seq and QuantSeq’ section) were recovered and duplicated reads were removed using the program MarkDuplicates from picard toolkit version 2.18.11 (Picard Toolkit 2019. Broad Institute, GitHub Repository. https://broadinstitute.github.io/picard/) with the parameters VALIDATION_STRINGENCY = LENIENT REMOVE_DUPLICATES = true. Then, the number of mapped reads within each files were recovered using idxstats commands of samtools v1.11. Each bam file was sub-sampled, using samtools v1.11, to have approximately the same number of mapped reads as in the smallest bam file. The command mpileup of the program bcftools v1.16 was used to produce bcf files using the following options: -I -d 10000 -O b -a AD. The mpileup was only performed on human exonic regions. SNP and unchanging nucleotide positions were next called using the command call from bcftools and the parameters -A -V indels -m -O b. The resulting positions were filtered by depth and quality with the command filter from bcftools and the parameters -i ‘QUAL>=10 && DP>=700’ -O b. Finally, these bcf files were again filtered using a homemade Python script to keep only positions that have at least 700 nucleotides of coverage depth and an alternative allele frequency lower than 5%. A transcriptional SNP was identified by REF > ALT, where REF is the nucleotide found in the reference genome at a particular position, and ALT is the nucleotide found on mapped reads at this position with REF ≠ ALT. The number X of nucleotide positions with a coverage depth greater than 700 and containing a SNP REF > ALT was recovered. This number X was then divided by the total number of nucleotide REF with a coverage greater than 700 to obtain a proportion of sites REF with a SNP REF > ALT. Then, we tested whether the proportions of the same SNP across different conditions were different by using a logistic regression. We used the same procedure to test whether the proportion of SNP REF>* are different between conditions. A SNP REF>* corresponds to any SNP located on a given nucleotide REF with a coverage above 700 in the genome. The relative SNP frequency of a given condition compared to DMSO-treated cells was computed using the same formula as defined in ‘Compositional bias analyses section’.

Mutagenesis reporter experiments

A stable clonal A375 cell line expressing GFP and a mutated and non-fluorescent version of mCherry (CherryOFF) was obtained from retro-viral particles, prepared from the pQC-CherryOFF-GFP plasmid according to Birnbaum et al. (61). pQC-CherryOFF-GFP was a gift from Fangliang Zhang (Addgene plasmid #129101; http://n2t.net/addgene:129101;RRID:Addgene_129101). Between 5 × 10⁶ and 2 × 10⁷ genetically-modified cells were treated for 72 h with DMSO (CTRL), 1 mM ENU (+ENU, N3385, Merck), or 1 mM ENU + 1 μM vemurafenib/cobimetinib (+ENU + MAPKi). Alternatively, cells were treated with DMSO (CTRL), 1 mM ENU (+ENU), or 1 μM vemurafenib/cobimetinib (MAPKi) for 72 h before to be cultured in a drug-free medium for 6 days before to be treated with 1 mM ENU (+MAPKi + ENU) for 72 h.

Immunofluorescence

Cells were seeded into 6-well plates with coverslips and treated the following day with DMSO or 500 nM vemurafenib/500 nM cobimetinib (Euromedex) for 48 h. At the end of the treatment, cells were fixed with 4% PFA at room temperature for 20 min and then washed five times with PBS. Cells were permeabilized in 0.5% PBS-Triton, washed three times with PBS for 5 min, and then blocked in PBS containing 0.1% Tween® 20 Detergent (Sigma) and 5% bovine serum albumin (BSA) for 1 h at room temperature. Primary antibodies were diluted in blocking solution, and anti-gH2AX antibody (1 :750) (Abcam, ab26350) was incubated at room temperature for 1 h, while anti-RAD51 (1:400; Abcam, ab176459) and anti-RAD18 (1:400; Proteintech, 18333–1-AP) were incubated at room temperature for 2.5 h. Cells were then washed five times using PBS with 0.1% Tween® 20, and Alexa Fluor™ 594 (Invitrogen, A21203) or Alexa Fluor™ 488 (Invitrogen, A21206) secondary antibodies (1:1000) were incubated for 1 h at room temperature. Cells were then washed five times with PBS with 0.1% Tween® 20, and coverslips were mounted on glass slides with 20 μl of ProLong™ Diamond with DAPI (Thermo Fisher, P36962) mounting medium. Slides were imaged using Leica 3D upright deconvolution microscope with a CoolSNAP HQ camera, at the PICT-IBiSA Imaging Facility in Orsay, and analysed using a semi-automatic macro on Fiji software. Nuclei containing ≥10 distinct foci were defined as foci-positive, and the percentage of positive nuclei was calculated as [(number of foci positive nuclei)/(number of nuclei scored)] × 100. A minimum of 100 nuclei per sample were scored, and the data shown were collected from at least three biological replicates. The statistical significance of the experimental data was determined using paired T-test in GraphPad (Prism10).

PARPi sensitivity

Cells (4 × 10³) were seeded into 96-well plates. The following day, the media was replaced with 100 μl of olaparib-containing media, at final concentrations of 0, 10, 20 and 30 μM, together with 500 nM vemurafenib and 500 nM cobimetinib (Euromedex), in technical triplicates. Control cells were treated with olaparib-containing media (at 0, 10, 20 and 30 μM final concentration) together with DMSO. At the end of the treatment, cell viability was measured by a WST-1 ATP-based assay (Roche, France) after 48 h of treatment. Briefly, WST-1 reagent was added to each well (10 μL per 100 μL of medium) and incubated at 37 °C for 1.5 h. The plates were then read at 450 nm on a TriStar2 LB 942 Multimode Microplate Reader (Berthold Technologies, Germany). Cell sensitivity is represented relative to 0 μM olaparib-treated cells. Statistical significance of the experimental data was determined using two-way ANOVA in GraphPad (Prism10).

Results

Compositional biases of MAPKi-regulated gene products

We performed RNA-sequencing after culturing the A375 melanoma cell line for 24h in the absence or presence of a combination of BRAF- and MEK-inhibitors (hereinafter termed MAPK inhibitors [MAPKi]). The expression levels of 2010 or 1719 mRNAs were significantly decreased (downregulated mRNAs) or increased (upregulated mRNAs), respectively, in MAPKi-treated cells as compared to control cells (Figure 1A and Supplementary Table S1). Interestingly, 813 and 753 MAPKi-downregulated mRNAs encoded proteins associated with the GO terms ‘nucleoplasm’ and ‘cytoplasm’, respectively, and 626 MAPKi-upregulated mRNAs encoded proteins associated with the GO term ‘Integral component of membrane’ (Figure 1A). In agreement with the fact that nucleoplasmic and cytoplasmic proteins are typically hydrophilic soluble proteins, while membrane proteins tend to be hydrophobic proteins, we noticed that the hydrophobicity index of proteins encoded by upregulated-mRNAs was higher than the hydrophobicity index of proteins encoded by downregulated-mRNAs (Supplementary Figure S1a). This observation raised the possibility that proteins encoded by down- or up-regulated mRNAs had different amino acid composition biases.

Accordingly, hydrophilic residues, such as lysine (Lys), Glu, Asp and asparagine (Asn), were enriched in proteins encoded by MAPKi-downregulated mRNAs, while hydrophobic amino acids, like tryptophan (Trp), cysteine (Cys), phenylalanine (Phe), and leucine (Leu), were enriched in proteins encoded by MAPKi-upregulated mRNA (Figure 1B, C). Thus, proteins encoded by downregulated mRNAs contained on average 30%, 12%, 8%, and 7% more Lys, Glu, Asp and Asn residues, respectively, than proteins encoded by upregulated mRNAs, while the latter contained on average 24%, 11%, 10% and 9% more Trp, Cys, Phe and Leu residues, respectively. In addition, a larger part of MAPKi-downregulated mRNAs encoded for proteins with a higher frequency of Lys, Glu, Asp and/or Asn residues, while a larger part of MAPKi-upregulated mRNAs encoded for proteins with a higher frequency of Trp, Cys, Phe and/or Leu residues (Supplementary Figure S1b).

We next analysed the codon content of MAPKi-regulated mRNAs. Down- and up-regulated mRNAs were enriched for different sets of codons (Figure 1D). Indeed, A/T-ending codons were enriched in MAPKi-downregulated mRNAs, while G/C-ending codons were enriched in MAPKi-upregulated mRNAs. We also noticed a selective enrichment of a subset of synonymous codons since, for example, only the GAA (but not the GAG), the GAT (but not the GAC), and the AAT (but not the AAC) codons–corresponding to Glu, Asp and Asn, respectively–were enriched in downregulated mRNAs but reduced in upregulated mRNAs (Figure 1D, E and Supplementary Figure S1c).

In summary, mRNAs that were down- or up-regulated by MAPKi treatment contained different codon compositional biases and encoded protein sets with different amino acid compositional biases.

Compositional biases of TDD-regulated mRNAs in response to MAPKi

We next tested the possibility that MAPKi treatment could affect mRNA stability in a translation-dependent manner. For this, we first measured the TDD index of mRNAs by comparing the expression level of mRNAs in cells treated or not with MAPKi, at initial condition or 3 and 5 h after inhibition of transcription alone or after inhibition of both transcription and translation (Figure 2A). MAPKi treatment increased the TDD index of 1390 mRNAs and decreased the TDD index of 183 mRNAs (Figure 2A and Supplementary Table S1). Although not a formal proof, this observation suggested that MAPKi could affect mRNA expression levels by modulating mRNA stability in a translation-dependent manner.

Remarkably, the Lys, Glu, Asp and Asn residues were enriched in proteins encoded by mRNAs with an increased TDD index in response to MAPKi, while other amino acids, including tryptophan (Trp), glycine (Gly), alanine (Ala), proline (Pro) and leucine (Leu), were enriched in proteins encoded by mRNAs with a decreased TDD index in response to MAPKi treatment (Figure 2B and Supplementary Figure S2a). In addition, the majority of codons that were enriched in mRNAs with an increased TDD index in response to MAPKi was reduced in mRNAs whose TDD index decreased (Figure 2C). Furthermore, the majority of codons that was enriched in mRNAs whose TDD index increased in response to MAPKi corresponded to A/T-ending codons, while most codons enriched in mRNAs whose TDD index decreased in response to MAPKi were G/C-ending codons (Figure 2C and Supplementary Figure S2b). This observation agreed with a recent report showing that the MAPK pathway modulates codon optimality of A/T-ending codons (8).

We noticed that MAPKi treatment increased the TDD index of a large number of mRNAs when compared to the number of mRNAs whose TDD decreased (i.e. 1390 versus 183, Figure 2A). In addition, a large number of MAPKi-downregulated mRNAs had a significantly increased TDD index in response to MAPKi, while the TDD index of MAPKi-upregulated mRNAs was either increased or decreased (Figure 2D). This suggested that the TDD increase could contribute to the downregulation of a large subset of mRNAs in response to MAPKi. On the contrary, TDD does not appear to explain MAPKi-induced mRNA upregulation. Based on these considerations, we decided to focus our analyses on MAPKi-downregulated mRNAs whose TDD index increased in response to MAPKi.

Proteins encoded by MAPKi-downregulated mRNAs whose TDD increased had a higher frequency of Lys, Glu, Asp and Asn as compared to proteins encoded by MAPKi-downregulated mRNAs whose TDD was not affected by MAPKi (Figure 2E, upper panel). Likewise, MAPKi-downregulated mRNAs whose TDD increased had a higher frequency of AAA (Lys), GAA (Glu), GAT (Asp) and AAT (Asn) codons, as compared to MAPKi-downregulated mRNAs with a non-affected TDD (Figure 2E, lower panel).

Collectively, these observations suggested that MAPKi decrease the expression level of a subset of compositionally biased mRNAs by triggering their co-translational degradation through TDD.

Compositional biases of MAPKi-regulated ribosomal peaks

To analyse the dynamics of ribosomes on mRNAs in response to MAPKi, we performed ribosome profiling experiments and then computed mRNA ribosomal peaks (Figure 3A and Supplementary Figure S3a). mRNA regions in which the local density of ribosomes was higher in control cells than in treated cells are referred to as CC peaks, and mRNA regions in which the local density of ribosomes was higher in treated cells than in control cells, are referred to as TC peaks. In total, 1281 CC peaks were detected in 984 mRNAs in control cells, and 1974 TC peaks were detected in 1509 mRNAs in treated cells (Figure 3A and Supplementary Table S1).

We next analysed the amino acid- and codon-composition of CC peaks and TC peaks. Some amino acids were enriched in TC peaks but reduced in CC peaks and conversely some amino acids were enriched in CC peaks but reduced in TC peaks (Figure 3B and Supplementary Figure S3b). For example, TC peaks contained more frequently at least one Lys, Asp, Glu and/or Asn residue than CC peaks, and the latter contained more frequently at least one Ala and/or Gly residue (Figure 3C). Furthermore, some codons were more frequent in TC peaks than in CC peaks, while other were more frequent in CC peaks (Figure 3D and Supplementary Figure S3b). Interestingly, most codons enriched in TC peaks were A/T-ending codons, while most codons enriched in CC peaks were G/C-ending codons (Figure 3D). This observation suggested that ribosomes could spend more time on A/T-ending codons in MAPKi-treated cells compared to control cells in agreement with a recent report showing that the MAPK pathway modulates codon optimality of A/T-ending codons (8). In addition, a higher enrichment of TC peaks in the A/T-ending codons, like GAA (Glu), GAT (Asp), AAT (Asn) and AAA (Lys) was observed in contrast to the corresponding G/C-ending codons (Figure 3E and Supplementary Figure S3c). Finally, in agreement with a relationship between local density of ribosomes and translation-dependent mRNA decay, we observed that the TDD index of most MAPKi-downregulated mRNAs containing TC peaks increased in response to MAPKi (Supplementary Figure S3d). Of note as well, ∼30% of MAPKi-downregulated mRNAs contained MAPKi-induced TC peaks and/or had their TDD increased in response to MAPKi (Supplementary Figure S3e).

To summarize, codons encoding amino acids such as Lys, Glu, Asp and Asn were enriched in (i) MAPKi-downregulated mRNAs, (ii) mRNAs whose TDD increased in response to MAPKi and (iii) MAPKi-induced ribosomal peaks (Figure 3F). Furthermore, only the A/T-ending codons (AAA, GAA, GAT, AAT) corresponding to these amino acids were enriched at the expense of the corresponding G/C-ending codons, with the exception of the AAG (Lys) codon (Figure 3G). These data support a model in which MAPKi treatment affects the dynamics of ribosomes when going through A/T-ending codons corresponding to certain amino acids (e.g. Glu and Asp), which could trigger a selective-mRNA TDD-dependant degradation.

Amino acid bioavailability and codon-dependent selective effects

As the decrease in the intracellular concentration of certain amino acids can induce ribosome pauses and as the MAPK pathway in melanoma cells can affect amino acid metabolism (see Introduction), we next measured the intracellular concentration of Lys, Glu, Asp, Asn, Gln and Arg in the absence or presence of MAPKi. While the intracellular concentration of Lys, Gln and Arg was not affected, the intracellular concentrations of Glu, Asp and Asn decreased in response to MAPKi (Figure 4A). Of note, similar results were obtained using another BRAF-mutated melanoma cell line (the M249 cell line), in which MAPKi led also to a decrease in Glu intracellular concentration (Figure 4B). This result raised the possibility that MAPKi could have a selective effect on mRNAs whose translation requires a relatively large amount of specific amino acids like Glu and Asp.

To challenge this possibility, we cultured melanoma cells in the absence of Asp and/or Glu. However, we did not observe any significant effect on neither Asp- and Glu-intracellular concentration, nor on cell viability (Supplementary Figure S4a and b). One possible explanation is that Glu and Asp are produced from Gln provided by the growth medium. Supporting this possibility, Gln depletion from the growth medium decreased the intracellular concentration of Gln as well as Glu and Asp (Figure 4C).

Since Gln deprivation somehow mimics the decrease of Glu and Asp intracellular concentration as observed in response to MAPKi (compared Figure 4A and C), we analyzed the effect of Gln deprivation on gene expression and ribosome profile (Supplementary Table S1, Supplementary Figure S4c). The decrease in the intracellular concentration of Gln, Glu and Asp that was induced by Gln deprivation was associated with an enrichment of some codons corresponding to these amino acids in downregulated mRNAs and in ribosome peaks induced by Gln depletion (Figure 4D). However, in contrast to what we observed in MAPKi-treated cells, the decrease in the intracellular concentration of Glu and Asp that was induced by Gln depletion was associated with an enrichment of the G/C-ending codons (i.e. GAG and GAC) and not the A/T-ending codons (i.e. GAA and GAT) (Figure 4D). These results suggested that the decrease in amino acid bioavailability is not sufficient to explain synonymous codon-selective effects in agreement with previous reports (see Introduction).

Since aminoacyl-tRNA synthetases contribute to codon-selective effects, we inspected our RNA-seq datasets and we found that the expression level of several aminoacyl-tRNA synthetases varied in response to MAPKi treatment or in response to Gln depletion. Among these, we focused on the EPRS1 aminoacyl-tRNA synthetase, which loads Glu onto the corresponding tRNAs, because the expression level of EPRS1 was repressed in response to MAPKi but increased in response to Gln depletion, as validated by RT-qPCR (Supplementary Figure S4d). To test the potential role of EPRS1 on codon-selective effects, A375 cells were cultured in the absence or presence of Gln and in the absence or presence of EPRS1 (Supplementary Table S1). EPRS1 depletion abolished the codon-selective effect of Gln depletion with respect to the enrichment of the GAG codon to the advantage of the GAA codon suggesting that EPRS1 could at least in part participate to the codon-selective effect observed after Gln depletion (Figure 4E).

We also observed in our datasets that some enzymes that can modulate codon-selective effects by biochemically modifying tRNAs (see Introduction) were differentially expressed in MAPKi-treated cells as compared to control cells. Among these, CTU2 caught our attention for three reasons: (i) CTU2 expression was decreased in MAPKi-treated cells (Supplementary Figure S4f); (ii) an important role of CTU2 in melanoma cells has already been reported (18) and (iii) CTU2, which modifies uracil on position 34 of some tRNAs, modulates the interactions between anti-codons and some A-ending codons, in particular the AAA (Lys) and GAA (Glu) codons (18). Since we observed an enrichment of these codons in MAPKi-downregulated mRNAs, TDD-induced mRNAs, and MAPKi-induced ribosomal peaks (Figure 3G), we tested whether the CTU2 depletion could mimic the selective MAPKi-effect. Supporting such a possibility, we first observed that CTU2 depletion resulted in the downregulation of a subset of mRNAs that were also downregulated in response to MAPKi (Figure 4F and Supplementary Table S1). In addition, the AAA (Lys) and GAA (Glu) codons were enriched in mRNAs whose TDD was increased in CTU2 depleted cells, as expected (18) and as observed in MAPKi treated cells (compare Figure 4G to Figure 3G).

In sum, our observations support a model in which the decreases in Glu and Asp bioavailability in MAPKi-treated cells resulted in an increase in the local density of ribosomes going through mRNA regions that require these amino acids to be translated (Figure 3), which could result in translation-dependent mRNA degradation (Figure 2), and consequently in the decrease of the expression levels of mRNAs encoding for compositionally biased proteins (Figure 1). However, selective effects of synonymous codons–notably A/T-ending codons corresponding to Asp and Glu–probably depend on several parameters, such as the expression of aminoacyl transferases or tRNA-modifying enzymes (Figure 4, see Discussion).

Protein amino acid composition biases and cellular functions

Since MAPKi treatment triggers TDD of a subset of mRNAs (Figure 2) and since TDD is likely to be dynamic and reversible, we next wondered whether the MAPKi-dependent mRNA downregulation persists after MAPKi removal, i.e. in the so-called persister cell population. To address this question, we performed RNA-seq on persister cells (Figure 5A and Supplementary Table S1). A large number of TDD-downregulated mRNAs that we identified in cells exposed to MAPKi (Figure 2) was still significantly downregulated in persister cells (Figure 5A). In addition, the MAPKi-dependent decrease in Glu-, Asp- and Asn-intracellular concentration observed in MAPKi-exposed cells was also observed in persister cells (compare Figures 4A and 5B). Furthermore, proteins encoded by MAPKi-downregulated mRNAs that were still downregulated in persister cells were enriched in Lys, Glu, Asp and Asn (Figure 5C, left panel) and these mRNAs were enriched in the AAA (Lys), GAA (Glu), GAT (Asp) and AAT (Asn) codons (Figure 5C, right panel). In sum, at least some composition biases observed in mRNAs that were downregulated in the presence of MAPKi were also observed in persister cells after MAPKi removal.

We next questioned the biological functions of the proteins encoded by the different mRNA populations. Some terms, like cell division, cell cycle, DNA repair and DNA damage, were enriched among the biological functions associated with MAPKi-downregulated gene products (Figure 5D). Very interestingly, we noticed that the proportion of gene products associated with the proliferation and DNA repair biological functions increased among MAPKi-downregulated gene products whose TDD increased in response to MAPKi and that were still downregulated in persister cells (Figure 5D, E). For example, ∼11% of MAPKi-downregulated gene products were involved in cell proliferation and this proportion reached ∼33% in MAPKi-downregulated gene products with an increased TDD and that were still downregulated in persister cells (Figure 5D, E). Consequently, >40% of the MAPKi-downregulated gene products with an increased TDD that were still downregulated in persister cells were involved in DNA metabolism (i.e. proliferation and/or DNA repair; Figure 5D, E). This result was validated by RT-qPCR and Western blot analysis since the expression level of genes involved in cell proliferation was higher in control cells (CC) compared to MAPKi-treated cells (TC) and persister cells (PC) (Figure 5F). Of note, the expression level of pro-proliferative genes increased after nine days of MAPKi-removal (RC, Figure 5F), a time at which persister cells gave rise to a proliferative cell population similar to the initial one. Of note, similar results were obtained using the M249 melanoma cell line, in which MAPKi led also to a decrease in the mRNA and protein levels of a subset of tested genes (Figure 5G).

Because of these observations, we analyzed the composition biases of mRNAs encoding proteins involved in proliferation and DNA repair. Proteins involved in proliferation and DNA repair were enriched in a subset of amino acids, including Lys, Glu, Asp and Asn, as compared to the human proteome (Figure 6A, B). In addition, mRNAs encoding proteins involved in proliferation and DNA repair were enriched in the AAA (Lys), GAA (Glu), GAT (Asp) and AAT (Asn) codons (Figure 6B, left panel). Worth noting, the AAA (Lys), GAA (Glu), GAT (Asp) and AAT (Asn) codons were more enriched in MAPKi-downregulated gene products involved in proliferation and/or DNA repair when compared to the other MAPKi-downregulated gene products (Figure 6C).

Since the expression level of mRNAs encoding proteins involved in DNA repair decreased in MAPKi-treated cells and in persister cells, as validated by RT-qPCR and Western blot analysis, while γ-H2AX–a marker of DNA damage–was significantly increased in MAPKi-exposed cells (Figure 6D and E), one hypothesis is that MAPKi-exposed cells may have a higher probability of accumulating DNA damage that could increase the probability of genetic mutations to appear in descendant cells. However, measuring the genetic mutational rate of MAPKi-exposed cells is challenging, as genetic mutations can only be quantified after several rounds of replication cycles, while MAPKi represses cell proliferation. In order to circumvent this difficulty, we used two complementary approaches.

First, we decided to look for nucleotide variations within mRNAs. Indeed, genetic mutations are notably the consequence of nucleotide mismatches during replication that themselves are the consequence of nucleotide chemical modifications (i.e. DNA damage such as nucleotide oxidation). In this setting, nucleotide chemical modifications can also lead to nucleotide mismatches during transcription, giving rise to the so-called transcriptional mutations (62–68). In other words, unrepaired DNA damage–such as nucleotide chemical modifications–can lead to transcriptional mutations in neo-synthetized mRNAs. By comparing the transcriptome of MAPKi-exposed and persister cells to control cells at the nucleotide level, we observed that the relative frequency of transcriptional mutations affecting each of the four nucleotides increased both in MAPKi-treated cells (TC) and in persister cells (PC) (Figure 7A, upper panel). For example, there was between 5% to 20% more transcriptional mutations changing a G nucleotide into A, C or T in MAPKi-treated cells (TC) or persister cells (PC) when compared to control cells (Figure 7A, lower panel).

The second approach we used was to generate a cellular clone expressing a reporter gene that contains a TGA stop codon preventing the synthesis of the mCherry protein, that can only be expressed if nucleotide biochemical modifications – induced for example by mutagenic agents, such as ENU – result in nucleotide mismatches that change the TGA stop codon into the TGG codon coding for tryptophan (61). As expected, the number of mCherry-positive cells was increased after ENU treatment regardless of the filters used in cytometry to count positive cells (Figure 7B, upper panel). Importantly, the ENU effect was increased when cells were simultaneously exposed to MAPKi (Figure 7B, upper panel) or when cells were first exposed to MAPKi before to be exposed to the ENU (Figure 7B, lower panel). Collectively, our results point to a link between the downregulation of DNA repair genes in MAPKi-exposed cells and a higher rate of nucleotide mismatches.

In line with the MAPKi-dependent increase of γH2AX protein level (Figure 6D), we observed a MAPKi-dependent increase in the number of cells containing at least 10 γH2AX foci, which also coincided with a decrease in the number of cells containing at least 10 foci of RAD51 or RAD18, two proteins involved in homologous recombination (Figure 7C and D) (69–71). This suggested that MAPKi decreased the expression level of proteins involved in DNA repair (Figure 6), which in turn impaired homologous recombination-dependent DNA repair. In agreement with this possibility, MAPKi decreased the number of cells that survived in response to the PARP-inhibitor olaparib, a common assay used to study homologous recombination defects (72) (Figure 7E). Collectively, our results point to a link between the metabolic-dependent downregulation of DNA repair genes in MAPKi-exposed cells and MAPKi-dependent defects in homologous recombination-dependent DNA repair.

Discussion

The compositional biases of MAPKi-downregulated gene products that we observed (Figure 1) could be explained at least in two ways. First, from a gene centric point of view, the MAPK pathway could have evolved to repress the transcription of genes that are involved in DNA metabolism such as DNA replication. As gene products involved in DNA metabolism bear compositional biases (Figure 6A, B; see below), then the MAPKi-downregulated gene products would bear these function-related compositional biases. The second possible explanation, which we term a metabolic centric point of view, is that the MAPKi-dependent decrease of the bioavailability of some amino acids (e.g. Asp and Glu) results in the translation-dependent expression level decrease of Asp- and Glu-enriched gene products. Since proteins involved in DNA metabolism are enriched in Asp and Glu, then MAPKi treatment induces the expression level decrease of gene products involved in DNA metabolism in a metabolism-depending manner. As discussed below, our results and those from previous publications support this metabolic centric point of view without excluding the well-established effects of MAPKi on the transcriptional activity of genes involved in cell proliferation.

The BRAF^V600E mutation-dependent hyper-activation of the MAPK pathway in melanoma cells triggers not only cell proliferation but also the cellular addiction to some non-essential amino acids, such as Gln. Indeed, in mutated melanoma cells, the carbon skeleton of Gln fuels both non-oxidative energetic metabolism and amino acid and nucleotide biosynthetic pathways on which growing cell depends (35,45,51,73–76). Accordingly, growth medium depletion of Gln induced Glu- and Asp-intracellular concentration decrease, while Glu- and/or Asp-depletion did not result in their intracellular concentration decrease (Figure 4C and Supplementary Figure S4b), probably because Glu and Asp can be generated from Gln from the growth medium (45,76). Importantly, the use of Gln, Glu and Asp in the oxidative phosphorylation (OXPHOS) pathway, which produces energy from the complete degradation of their carbon skeleton, is reactivated by MAPKi (46–51,73). This may explain the observed Glu- and Asp-intracellular concentration decrease in MAPKi-exposed cells (Figures 4A and 5B). This together with the observed enrichment of these amino acids in (i) MAPKi-downregulated mRNAs (Figure 1B), (ii) MAPKi-induced TDD mRNAs (Figure 2B) and (iii) MAPKi-induced ribosomal peaks (Figure 3B) support a model where MAPKi-dependent effects on the cellular metabolism impacts the translation-dependent mRNA expression level through the bioavailability of amino acids.

The precise codon composition biases of MAPKi-downregulated gene products likely depend on several parameters (Figure 4, and see Introduction) that could be regulated by different mechanisms (including the MAPK-dependent phosphorylation cascade). Notably, the observed amino acid composition biases of MAPKi-downregulated gene products are particularly interesting, as they may explain how cells ‘coordinate’ their metabolic activity and their phenotype. Indeed, we observed that gene products involved in DNA metabolism (e.g. DNA replication) are enriched in charged amino acids such as Glu and Asp, consistent with the fact that these proteins are enriched in hydrophilic and charged amino acids that play a role their DNA binding- and enzymatic-activities (77–79) (80). Since charged-amino acids like Glu and Asp are at the crossroad between energetic and gene product synthesis pathways, cells that survive MAPKi treatment could correspond to cells that consume these amino acids through the OXPHOS pathway, thereby decreasing the biogenesis of Glu- and Asp-rich proteins involved in cell proliferation and in turn reducing the energetic cellular needs. Accordingly, persister cells (i.e. non-genetically modified cells that survive anti-MAPK therapy) (i) have a lower intracellular concentration of Glu and Asp as compared to the initial cell population (Figure 5B), in agreement with the OXPHOS pathway re-activation reported in these cells and (ii) express a low level of gene products involved in DNA metabolism (Figure 5D–F), in agreement with their reported slow growing rate (46–48,81–86).

While favoring survival versus proliferation, the MAPKi-dependent decrease in the intracellular concentration of Glu and Asp could have a secondary effect because of the expression level decrease of gene products involved in DNA repair (Figures 5D, E and 6D), that share the same composition biases with genes involved in cell proliferation (Figure 6A–C). The expression level decrease of gene products involved in DNA repair may seem to be of no consequence in non-dividing cells since genetic mutations can only occur during replication. However, the fact that unrepaired DNA damage can induce base-pairing mismatches during transcription – leading to transcriptional mutations (62–68) – could explain the observed larger number of nucleotides variations in the transcriptome of cells exposed to MAPKi compared to control cells (Figure 7A). Interestingly, it has been proposed that transcriptional mutations could be a ‘pre-selection step’ toward the emergence of genetically-modified and –adapted cells (62–68). Indeed, if unrepaired DNA damage leads to the synthesis of mutated gene products that contribute to the survival of a cell, this cell may have a higher probability of generating genetically-modified and -adapted descendant cells because the same unrepaired DNA damage could trigger a genetic mutation. Although this model is speculative with respect to the data we provided, persister cells have been proposed to be a reservoir of genetically modified and therapy-resistant cells (82).

In conclusion, we propose that MAPKi-induced metabolic changes result in the bioavailability decrease of amino acids such as Glu and Asp, which contributes to the expression level decrease of proteins enriched in these amino acids, including proteins involved in proliferation and DNA repair (Figure 7F). The coupling between metabolism and gene expression could, as a side effect, results in the accumulation of DNA damage – owing to the expression level decrease of DNA repair enzymes – leading first in transcriptional mutations and then in genetic mutations, increasing therefore the probability of genetically-mutated and -adapted clones to emerge in response to MAPKi.

Data availability

The raw NGS datasets from this study were deposited on Gene Expression Omnibus (GSE232604; GSE232709; GSE233122; GSE232551; GSE233120; GSE233121).

Supplementary data

Supplementary Data are available at NAR Cancer Online.

Acknowledgements

We gratefully acknowledge the support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. We thank the members of the LBMC Biocomputing Hub and V.A. Raker for manuscript editing. E.P.R is supported by the European Research Council (ERC-StG-LS6-805500) under the European Union's Horizon 2020 research and innovation programs.

Author contributions: F.A., L.G., E.V. and L.F. performed experiments; N.F. and A.L. performed bioinformatics analyses; E.P.R., D.A. and S.V. conceived the project and supervised the experimental work and the bioinformatics analyses; F.A., N.F., E.P.R. and D.A. wrote the manuscript with help and comments from all co-authors; and all authors read and approved the final manuscript.

Funding

INCa [PLBIO19-257]; Conseil Régional Auvergne Rhône-Alpes [METACOMBO]; European Research Council (ERC-StG-LS6-805500) under the European Union's Horizon 2020 research and innovation programs.

Ethics approval and consent to participate

N/A. The cell lines were authenticated by ATCC.

Conflict of interest statement. None declared.

References

Author notes