Stochastic variation in the FOXM1 transcription program mediates replication stress tolerance

TL;DR

The study shows that variable expression of the FOXM1 gene helps cancer cells survive DNA-damaging drugs, offering new insights into drug resistance.

Contribution

The study reveals that partial suppression of FOXM1 and its targets can protect cells from DNA damage during replication stress.

Findings

37 genes, including FOXM1 targets, are differentially expressed in drug-tolerant versus sensitive cells.

Partial FOXM1 knockdown reduces DNA damage and improves recovery from replication stress.

UBE2C and MKI67, FOXM1 targets, play unexpected roles in the replication stress response.

Abstract

Oncogene‐induced replication stress (RS) is a vulnerability of cancer cells that forces reliance on the intra‐S‐phase checkpoint to ensure faithful genome duplication. Inhibitors of the intra‐S‐phase checkpoint kinases ATR and CHK1 have been developed, but resistance to these drugs remains problematic. Understanding drug tolerance mechanisms is impeded by analysis of bulk samples, which neglect tumor heterogeneity and often fail to accurately interpret cell cycle‐mediated resistance. Here, by combining intracellular immunostaining and single‐cell RNA‐sequencing, we characterized the transcriptomes of oncogenic RAS‐expressing cells with variable levels of RS when challenged with a CHK1 inhibitor combined with gemcitabine. We identified 37 genes differentially expressed between tolerant and sensitive cells, including several FOXM1 targets. While complete knockdown of FOXM1 impeded cell proliferation, partial knockdown protected cells against DNA damage, and improved recovery from drug‐induced RS. Remarkably, knockdown of individual FOXM1 target genes UBE2C and MKI67 also mitigated DNA damage, uncovering unanticipated roles for these in the replication stress response. Our results suggest that low levels of FOXM1‐dependent gene expression during S and G2 phase protects cells against excessive DNA damage during drug‐induced replication stress.

Cellular heterogeneity is a major cause of drug resistance in cancer. Segeren et al. used single‐cell transcriptomics to investigate gene expression events that correlate with sensitivity to the DNA‐damaging drugs gemcitabine and prexasertib. They show that dampened expression of transcription factor FOXM1 and its target genes protected cells against excessive DNA damage during S‐ and G2‐phases of the cell cycle.

Figures

γH2AX is a replication stress marker suitable for flow cytometry of DSP‐fixed cells. (A) Schematic overview of the technique to combine immunostaining and single‐cell RNA sequencing. Cells are fixed with DSP [dithiobis(succinimidyl propionate)], permeabilized and stained using fluorescent antibodies. Next, cells are sorted based on fluorescent intensity. After de‐crosslinking, cells are subjected to single‐cell RNA sequencing. (B) Flow cytometry data showing the intensity of potential RS markers in individual RPE‐HRASG12V cells treated for 24 h with 10 nm CHK1 inhibitor (CHK1i) + 100 nm gemcitabine or vehicle (Veh). Unstained control refers to control cells not incubated or, when applicable, incubated with the secondary antibody only. Inset zooms in on cells stained for γH2AX after treatment with CHK1i + gemcitabine to indicate heterogeneity in γH2AX level. Representative of 3 independent experiments. (C) Representative example of γH2AX immunostaining on RPE‐HRASG12V cells treated for 24 h with 10 nm CHK1i + 4 nm gemcitabine or vehicle (Veh). Representative of 3 independent experiments. (D) Immunoblot showing synergistic induction of RS by 10 nm CHK1i + 4 nm gemcitabine in RPE‐HRASG12V cells, as indicated by phospho CHK1 S345 and γH2AX. The absence of phosphorylation of CHK1 on its autophosphorylation site S296 indicates effective inhibition by CHK1i. Representative of 3 independent experiments.

High quality single‐cell RNA sequencing data of fixed cells with known level of replication stress. (A) Violin plot representing the average numbers of genes detected per cell in fresh (n = 232) and DSP [dithiobis(succinimidyl propionate)] fixed RPE‐HRASG12V (n = 237) cells. (B) Scatter plot showing the average gene expression in DSP‐fixed and fresh cells. R value indicates Pearson Correlation. The red line indicates x = y. (C) Scatter plot showing the correlation between the gene detection rate and average gene expression in fresh and DSP‐fixed cells. R value indicates Pearson Correlation coefficient between fresh and DSP‐fixed cells. (D) Scatter plot showing the correlation between the coefficient of variation and average gene expression in fresh and DSP‐fixed cells. R value indicates Pearson Correlation coefficient between fresh and DSP‐fixed cells. (E) Dimensionality reduction using tSNE (t‐distributed stochastic neighbor embedding) of DSP‐fixed cells. Cells are color coded according to γH2AX signal. (F) Feature plot in which cells on tSNE plot in E are color coded according to mAzami Green‐Geminin1‐110 signal. (G) Feature plots in which cells on tSNE plot in E are color coded according to the expression of S‐phase (CDC6 and E2F1) or G2‐phase (PLK1 and CCNB1) markers.

Single‐cell RNA sequencing of S‐phase cells sorted based on γH2AX level. (A) Flow cytometry data of RPE‐HRASG12V and control cells treated for 16 h with 10 nm CHK1i + 4 nm gemcitabine or vehicle (Veh). Sorting strategy is shown: first S‐phase cells were selected based on DAPI (4′,6‐diamidino‐2‐phenylindole) signal. For drug‐treated cells equal number of cells with no, low, medium or high levels of γH2AX were sorted. The percentages of cells in these different categories before sorting are indicated and show an increase in the cell population with high level of RS after treatment with CHK1 inhibitor (CHK1i) + gemcitabine. Respectively one and three 384‐well plates were sorted for each vehicle and drug‐treated group. (B) Dimensionality reduction using tSNE (t‐distributed stochastic neighbor embedding) of all cells (HRASG12V and control) before and after treatment with 10 nm CHK1i + 4 nm gemcitabine shows separate clusters of cells. (C) Feature plot in which cells on tSNE plot in B are color coded based on the different conditions. Inset shows an example of cells from the drug‐treated group clustering among the Veh‐treated cells. (D) Feature plot in which cells on tSNE plot in B are color coded according to γH2AX signal.

Identification of genes correlated with levels of replication stress. (A) Heatmap of genes differentially expressed and downregulated in RPE HRASG12V γH2AXlow versus γH2AXhigh control RPE cells 16 h after treatment with 10 nm CHK1 inhibitor (CHK1i + 4 nm gemcitabine). (B) Heatmap of genes differentially expressed and upregulated in γH2AXlow versus γH2AXhigh RPE HRASG12V cells 16 h after treatment with 10 nm CHK1i + 4 nm gemcitabine. (C) FOXM1 gene expression as measured by single‐cell RNA sequencing in RPE‐HRASG12V cells with low, medium, and high γH2AX levels after 16 h of treatment with 10 nm CHK1i + 4 nm gemcitabine. The upper and lower edge of the boxplot indicates the 1st and 3rd quartile, the middle line indicates the median, and the upper and lower edges of the whiskers indicate the 95% confidence interval. No significant differences were identified based on the Wilcoxon rank sum tests for multiple group comparisons. (D) Correlation between pFOXM1 and γH2AX in RPE‐HRASG12V cells treated for 16 h with 10 nm CHK1i + 4 nm gemcitabine, as measured by flow cytometry. The trendline represents the data fit to a linear model, with the shaded region around the trendline indicating the 95% confidence interval. R value indicates Pearson Correlation coefficient. (E) Correlation between DAPI (4′,6‐diamidino‐2‐phenylindole) and pFOXM1, as measured by flow cytometry in the same experiment as described under D. The trendline represents the data fit to a linear model, with the shaded region around the trendline indicating the 95% confidence interval. R value indicates Pearson Correlation coefficient. Plots in A–C show n = 188, 230, and 210 cells in the γH2AXlow, γH2AXmedium, and γH2AXhigh groups, respectively.

Partial knockdown of FOXM1 improves tolerance to replication stress without affecting cell proliferation. (A) Schematic representation of the experimental design to identify and validate putative RS‐tolerance genes. (B) Relative FOXM1 expression of RPE‐HRASG12V following transfection with four individual siRNAs targeting FOXM1 (1 nm each) as measured by qPCR. Gene expression was normalized to the average of two housekeeping genes (GAPDH, 18S). Error bars indicate mean ± SEM. Significant differences determined by one‐way ANOVA with Geisser Greenhouse correction followed by Dunnett's multiple comparison test. P < 0.01, N = 3. Representative of 3 individual experiments. (C) Flow cytometry data of RPE‐HRASG12V cells which were arrested in G1‐phase after 24 h treatment with a CDK4/6i (1 μm) and individual siRNAs against FOXM1. Subsequently, cells were released in the absence or presence of 10 nm CHK1i + 4 nm gemcitabine and harvested after 14 h to enrich for S/G2‐phase cells. DAPI (4′,6‐diamidino‐2‐phenylindole) staining was used to determine cell cycle progression (top row) and γH2AX staining was used to determine the degree of replication stress (bottom row). Number in bottom right corner of bottom row plots indicates the Geminin1‐110 + cells as percentage of the total cells. Number in the top right corner of bottom row plots indicates γH2AX+ cells as a percentage of Geminin1‐110 + cells. Representative of 2 individual experiments. (D) Outgrown colonies of RPE‐HRASG12V cells transfected for 24 h with siRNAs against FOXM1 followed by treatment for 48 h with the 2 nm CHK1i and 4 nm gemcitabine. After removing drug‐containing media, colonies were allowed to grow for 5 days. Representative of 2 individual experiments. (E) Quantification of colonies presented in panel D, quantified as the IntensityPercent (which takes into account both the area covered by the cell growth as well as the pixel intensity of the covered area) and presented as IntensityPercent relative to untreated cells for each transfection condition. IntensityPercent was quantified using the ImageJ ColonyArea plug‐in. Error bars indicate mean ± SEM. Significant differences were determined by ordinary One‐way ANOVA followed by Dunnett's multiple comparison test. *P < 0.05, P < 0.01, N = 3.

Full text

Introduction

1

Oncogene‐induced replication stress (RS) is a vulnerability of cancer cells that can be exploited by anti‐cancer therapies. Seminal studies in the beginning of this century already showed that oncogenes, such as RAS, induce DNA damage in precancerous lesions [1, 2]. Further research revealed that oncogene‐induced RS underlies the elevated levels of DNA damage and that RS is present in the vast majority of human tumors. As a result, RS is proposed as an emerging hallmark of cancer [3].

RS is defined as stalling of the replication fork, which can arise due to shortage of substrates, collisions between replication and transcription machinery, or DNA lesions or secondary structures that hinder the replication machinery. Unresolved RS can progress to replication fork collapse, resulting in single‐ and double‐stranded DNA breaks. To prevent this, cells respond to RS by triggering the intra‐S‐phase checkpoint. Briefly, this checkpoint is initiated when Replication Protein A (RPA) binds to single‐stranded DNA that is exposed upon uncoupling of helicase and polymerase activity during fork stalling. This triggers recruitment and activation of ATR and its downstream kinase CHK1, which together induce a cascade of kinase activation that acts to stabilize and repair the stalled replication fork, fire dormant origins in the vicinity of the stalled fork, attenuate global DNA replication and slow down cell cycle progression. This multifaceted response ensures faithful genome duplication before mitosis [4].

In general, loss of ATR or CHK1 is lethal in cells where oncogenes are activated [5, 6, 7, 8]. On the basis of this knowledge, inhibitors against key‐players of the intra‐S‐phase checkpoint are developed and currently evaluated in clinical trials [9]. To potentiate the effect of intra‐S‐phase checkpoint ablation, ATR and CHK1 inhibitors can be combined with a low dose of chemotherapeutic drugs [10, 11]. However, drug resistance remains a major problem [12]. The limited in vivo activity of drugs which exacerbate RS suggests that cancer cells employ strategies to tolerate RS. Indeed, stabilization of the replication fork [13], increased expression of RPA [14], increased dormant origin firing [15], and modulation of cell cycle regulators [16] grant RS tolerance. The RS response is tightly linked with cell cycle progression, as multiple intra‐S‐phase checkpoint kinases play a role in curtailing proteins involved in the S‐G2 transition [16, 17]. Moreover, unbiased screening approaches uncovered that cell cycle related genes mediate resistance to intra‐S‐phase checkpoint inhibitors [18, 19, 20]. However, since these studies employed bulk sample approaches, transcriptional heterogeneity was neglected, rare resistance‐conferring events missed, and the role of cell cycle progression potentially misinterpreted. As a result, the development of novel clinical strategies based on these studies is rare.

The importance of single‐cell data in drug resistance studies is highlighted by Shaffer et al. [21], who unveil that rare cancer cells express resistance genes prior to treatment to resist therapy. In support of this notion, treatment with RS‐inducing drugs leads to a reduction in the number of transcriptionally distinct clones [22], suggesting selection pressure for cells harboring drug‐tolerant characteristics. Besides pre‐existing heterogeneity, it is becoming increasingly evident that cancer cells modulate their transcriptome upon treatment to circumvent therapy. For example, chemotherapeutic drugs may induce a transient drug‐tolerant state in a subpopulation of cells [23, 24]. It is hypothesized that this provides a time window in which permanent resistant cells can arise. Because transcriptional heterogeneity could result in resistance to RS‐inducing drugs and consequently tumor relapse, the mechanisms underlying RS tolerance warrant further investigation.

Here, we employed a strategy in which we combine immunostaining of RS markers and information on cell cycle phase with single‐cell RNA‐sequencing. This allowed us to shed light on the biological variability in response to RS. We uncovered a subset of genes with an altered expression profile in cells that maintained low levels of RS despite challenge with RS‐inducing drugs. We also identified genes that make cells more sensitive to replication stress, which included several FOXM1 target genes. Consistent with this, partial knockdown of FOXM1 mitigated DNA damage and improved cell survival following treatment with RS‐inducing drugs. These findings provide potential new avenues for development of synthetic lethality strategies and identification of biomarkers to optimize anti‐cancer therapy.

Methods

2

Key resources

2.1

Key resources are listed in Table S5.

Cell lines

2.2

The parental hTERT RPE‐1 cell line (RRID:CVCL4388) was obtained from ATCC and cultured at 37 °C, 5% CO_2_ in DMEM supplemented with 10% FBS and 1% pen/strep. Cell lines were authenticated via genetic fingerprinting and were used no longer than 20 passages upon thawing cryovials of the HRAS^G12V^‐transduced cells. Cell lines were regularly tested and confirmed mycoplasma negative.

Overexpression of HRAS^G12V^ was induced by adding 0.2 μg·mL^−1^ doxycycline to the culture medium. Gemcitabine, prexasertib and palbociclib were purchased from Selleck chemicals (Selleck Biotechnology, Cologne, Germany) and used at a final concentration of 4 nm, 10 nm and 1 μm respectively, unless stated otherwise.

RPE cell lines harboring the Tet Repressor, doxycycline‐inducible HRAS^G12V^, FUCCI4 system and fluorescent tagged truncated 53BP1 were created using the third‐generation lentiviral packaging system as previously described [25].

DSP fixation and antibody staining of single cells

2.3

Fixation of cells was performed according to the protocol described by Gerlach et al. [26]. In short, cells were collected by trypsinization, washed with PBS after which cells were fixed with 0.5 mm dithiobis(succinimidyl propionate) (DSP) in Sodium Phosphate‐buffered Saline (pH 8.4) for 45 min at room temperature at a concentration of 1 million cells per 2.5 mL. Next, DSP was neutralized by incubating the cells with quench buffer (100 mm Tris, pH 7.5, 150 mm NaCl) for 10 min and cell clumps were removed using a 70 μm cell strainer. Cells were incubated for 30 min with BP buffer (PFBB : PBS, 1 : 1, supplemented with 0.1% Triton‐X 100) to permeabilize the cells, after which samples were incubated overnight with antibodies in BP buffer. If samples were intended to use for single‐cell RNA‐sequencing, BP buffer was supplemented with 2 U·μL^−1^ RNAsin Plus. Samples were filtered using a 40 μm cell strainer and incubated with DAPI (0.2 μg·mL^−1^) prior to loading of samples on the flow cytometer.

For antibody testing, samples were loaded on a CytoFLEX flow cytometer and analyzed using flowjo v10.0 software (FlowJo LLC, Ashland, OR, USA). Index sorting of cells for single‐cell RNA‐sequencing was performed on a BD Influx cell sorter. Antibodies and dilutions used are listed in Table S6.

Microscopy

2.4

For immunofluorescence staining, cells were seeded on coverslips. Prior to fixation of cells using 4% paraformaldehyde for 20 min, pre‐extraction with 0.2% Triton‐X 100 for 1 min on ice was performed. Next, cells were permeabilized using 0.1% Triton‐X 100 for 10 min, blocked with 5% goat serum and incubated with indicated antibodies after which coverslips were mounted on slides. Samples were analyzed on a Leica SP8 confocal microscope equipped with a 20× objective. Antibodies and dilutions used are listed in Table S6.

Immunoblotting

2.5

For immunoblotting, cells were washed twice with ice‐cold PBS and lysed in ice‐cold RIPA‐buffer (50 nm Tris/HCl pH 7.5, 1 mm EDTA, 150 mm NaCl, 0.25% deoxycholate, 1% NP‐40) supplemented with NaF (1 mm), NaV_3_O_4_ (1 mm) and protease inhibitor cocktail (11873580001, Sigma Aldrich, Merck Life Science BV, Amsterdam, Netherlands) after which samples were subjected to a standard SDS/page immunoblot. Antibodies used and dilutions are listed in Table S6.

RNAi transfections

2.6

For siRNA experiments, cells were transfected with siRNA targeting the gene of interest or a scrambled control using Lipofectamine RNAiMAX according to manufacturers' instructions (Life Technologies, 13778030, Thermo Fisher Scientific, Waltham, MA USA). ON‐TARGETplus SMARTpool siRNAs were purchased as custom cherry‐pick libraries from Dharmacon and used at a final concentration of 10 nm, while individual siRNAs were used at a final concentration of 1 nm. Efficient knockdown of intended target was confirmed by quantitative PCR 24 h after transfection.

Flow cytometry

2.7

To test γH2AX levels in siRNA‐transfected cells, cells were collected by trypsinization, washed with PBS, transferred to a 96 well plate and fixed using 4% PFA for 30 min while gently shaking. Next, cells were washed with 0.1% BSA in PBS and permeabilized using 0.2% Triton for 30 min. Cells were washed once more with 0.1% BSA in PBS prior to incubation with the fluorescent linked γH2AX antibody for 1 h at room temperature. DAPI was added to the samples at a final concentration of 2.0 μg/100 000 cells to stain DNA content. Samples were loaded on a CytoFLEX flow cytometer and analyzed using flowjo v10.0 software.

To investigate the relationship between γH2AX and phospho‐FOXM1, cells were collected, fixed, and permeabilized as described above. Staining of pFOXM1 was performed as previously described [16]. Briefly, following permeabilization in 0.2% Triton‐X, cells were washed once with 1% BSA in PBS. Cells were incubated overnight at 4C with the primary antibody (diluted 1 : 500 in 1% BSA in PBS). The following morning, cells were washed 1× with 1% BSA in PBS/0.2% Triton‐X and incubated for 1 h at room temperature with the secondary antibody (diluted 1 : 800 in 1% BSA in PBS). Cells were washed again using 1% BSA in PBS/0.2% Triton‐X. DAPI was added to the samples at a final concentration of 2.0 μg/100 000 cells to stain DNA content. Samples were loaded on a CytoFLEX flow cytometer and analyzed using flowjo v10.0 software.

Clonogenic survival assays

2.8

Cells were seeded at low density (200 cells per 12 well plate) to assess colony formation. Following 24 h transfection with 1 nm siRNAs targeting FOXM1, cells were treated with 2 nm prexasertib and 4 nm gemcitabine. After 48 h exposure to the drugs, media was replaced with drug‐free media and remaining cells were allowed to grow out to form colonies. The imagej colonyarea plug‐in was used to quantify the area of the well covered by colonies as previously described [27].

Quantitative PCR

2.9

RNA isolation, reverse transcription, and quantitative PCR were performed according to manufacturers' instructions using the QIAGEN RNeasy kit, Thermo Fisher cDNA synthesis kit and Bio‐RAD SYBR Green Master mix, respectively. Quantification of relative gene transcript levels was performed using the ΔΔC t method for multiple‐reference gene correction using GAPDH or β‐Actin and RPS18 as reference genes. Primers used in this manuscript are listed in Table S7.

Single‐cell RNA‐sequencing

2.10

For single‐cell RNA‐sequencing single cells were collected in 384‐well plates containing barcoded CEL‐seq2 primers and nucleotides using index sorting and stored at −80 °C until further processing.

De‐crosslinking of the DSP‐fixed cells was performed by addition of 0.1 m DTT to the reverse transcription mix (10 mm DTT final concentration), which is part of the regular reverse transcription mix for unfixed cells. SORT‐seq sequencing and read alignment were performed by Single Cell Discoveries (Utrecht, The Netherlands) using their pipeline based on CEL‐Seq2 [28, 29]. Briefly, samples were barcoded with CEL‐seq2 barcodes and UMI during reverse transcription and pooled after second strand synthesis. The resulting cDNA was amplified with an overnight in vitro transcription reaction. From this amplified RNA, sequencing libraries were prepared with Illumina TruSeq small RNA primers, which were paired‐end sequencing on the Illumina NextSeq500 platform. Read 1 was used to identify the Illumina library index and CEL‐seq2 sample barcode. Read 2 was aligned to the human genome (hg38) transcriptome using the burrows–wheeler aligner v0.7.17. Reads that mapped equally well to multiple locations were discarded. Mapping and generation of count tables were done using the MapAndGo2 script. Downstream processing and analysis were performed in rstudio (Version 1.4.1106) and r (Version 4.0.5) using the seurat package (Version 3.2.3) [30]. Cells were filtered and selected for downstream analysis when the following parameters were met: number of detected genes > 500 and < 7500, Unique Molecular Identifier (UMI) counts > 1000 and < 50 000, and the percentage of mitochondrial counts and ERCC RNA spike‐ins below 10. Next, raw counts were normalized, and variance stabilized using the SCTransform method [31]. Subsequently, dimension reduction was performed by principal component analysis. Identified clusters were visualized with t‐Distributed Stochastic Neighbor Embedding (tSNE). Differentially expressed genes were identified using the Seurat FindAllMarkers() function with one non default argument, min.pct = 0.25 requiring a greater fraction of cells within a cluster to have expression. After this the results were filtered at a Bonferroni‐adjusted significance level of P < 0.05. Expression correlation between the differentially expressed genes was determined using Pearson correlation. All sequencing data generated in this study are available on the Gene Expression Omnibus under accession numbers GSE256134 and GSE250285.

Quantification and statistical analysis

2.11

Flow cytometry, immunoblot, and quantitative PCR experiments were performed three times unless indicated otherwise. Details on sample size and statistical methods employed are described in the figure legends. *P < 0.05, **P < 0.01, ***P < 0.001 unless indicated otherwise.

Results

3

γH2AX is a replication stress marker suitable for flow cytometry of DSP‐fixed cells

3.1

To unmask transcription‐mediated RS‐adaptation mechanisms, transcriptomic information and the level of RS in single cells needs to be combined. Therefore, we adapted a previously published strategy in which cells are reversibly fixed to allow antibody staining while preserving RNA for sequencing [26]. As summarized in Fig. 1A, cells were fixed using the chemically reversible crosslinking reagent DSP. Next, cells were stained with an antibody that recognizes an RS‐specific marker and sorted based on RS levels in 384‐well plates using FACS. Subsequently, de‐crosslinking was performed using the reducing agent DTT and cells were subjected to first strand cDNA synthesis and single‐cell RNA‐sequencing.

γH2AX is a replication stress marker suitable for flow cytometry of DSP‐fixed cells. (A) Schematic overview of the technique to combine immunostaining and single‐cell RNA sequencing. Cells are fixed with DSP [dithiobis(succinimidyl propionate)], permeabilized and stained using fluorescent antibodies. Next, cells are sorted based on fluorescent intensity. After de‐crosslinking, cells are subjected to single‐cell RNA sequencing. (B) Flow cytometry data showing the intensity of potential RS markers in individual RPE‐HRASG12V cells treated for 24 h with 10 nm CHK1 inhibitor (CHK1i) + 100 nm gemcitabine or vehicle (Veh). Unstained control refers to control cells not incubated or, when applicable, incubated with the secondary antibody only. Inset zooms in on cells stained for γH2AX after treatment with CHK1i + gemcitabine to indicate heterogeneity in γH2AX level. Representative of 3 independent experiments. (C) Representative example of γH2AX immunostaining on RPE‐HRASG12V cells treated for 24 h with 10 nm CHK1i + 4 nm gemcitabine or vehicle (Veh). Representative of 3 independent experiments. (D) Immunoblot showing synergistic induction of RS by 10 nm CHK1i + 4 nm gemcitabine in RPE‐HRASG12V cells, as indicated by phospho CHK1 S345 and γH2AX. The absence of phosphorylation of CHK1 on its autophosphorylation site S296 indicates effective inhibition by CHK1i. Representative of 3 independent experiments.

Before implementation of this technique, we first investigated which antibody against RS‐induced protein modifications is compatible with DSP fixation and analysis by flow cytometry. To induce RS, we employed the frequently used chemotherapeutic drug gemcitabine in combination with the CHK1 inhibitor (CHK1i) prexasertib [25]. In response to RS, the intra‐S‐phase checkpoint kinase ATR is activated to stabilize and repair stalled forks and delay cell cycle progression. This is mediated by a sequence of events including phosphorylation of CHK1, RPA2, KAP1, and H2AX [16, 32]. Antibodies against phosphorylated variants of these proteins are previously shown to detect increased levels of RS by flow cytometry [16, 33]. We assessed if these antibodies could detect an increase in phosphorylated CHK1, RPA2, KAP1, and H2AX in DSP‐fixed cells after treatment with CHK1i + gemcitabine.

We made use of RPE‐1 cells harboring a doxycycline‐inducible variant of oncogenic RAS (hereafter referred to as RPE‐HRAS^G12V^ for cells with doxycycline‐induced expression of HRAS^G12V^ or control for their non‐induced counterparts). The advantage of this system is that adaptation to RS can be studied in the frequently occurring oncogenic context of RAS hyperactivation, while the effect of other tumor‐specific mutations is excluded. We previously described that RPE‐HRAS^G12V^ cells show mild endogenous RS and markedly enhanced sensitivity to CHK1i + gemcitabine [25]. However, control RPE cells also show RS in presence of high doses of these drugs. Accordingly, treating RPE control cells with a high dose of CHK1i + gemcitabine resulted in 2 N‐cell cycle arrest, as seen by the accumulation of cells with low DAPI signal, indicating severe stress. CHK1i + gemcitabine also triggered an abundant increase in phosphorylated KAP1, RPA2 and H2AX (Fig. 1B). However, the tested antibody against phospho‐Serine 345 on CHK1 failed to show an increase in this flow cytometric analysis of DSP‐fixed cells, excluding it as an RS‐marker for this project (Fig. 1B). While KAP1 mediates RS‐induced DNA remodeling and RPA2 protects stalled replication forks, phosphorylated H2AX is present at collapsed replication forks [32, 34]. Since the latter is the most downstream event in the RS‐cascade and indicates severe RS, the antibody against phosphorylated H2AX Serine 139 (hereafter referred to as γH2AX) was selected as a proxy for RS‐induced DNA damage. Interestingly, flow cytometry analysis of γH2AX stained cells revealed great diversity in the signal, and presumable RS‐level, between individual cells (Fig. 1B, inset). Consistent with this, heterogenous phosphorylation of H2AX S139 in response to RS was confirmed by immunofluorescence staining (Fig. 1C). In addition, immunoblotting confirmed that substantial γH2AX was observed when the CHK1i prexasertib was combined with a low dose of gemcitabine, but not with either drug alone (Fig. 1D). Based on these observations, we concluded that the antibody against γH2AX can be used to determine the level of RS induced by CHK1i and gemcitabine at a single cell resolution in DSP‐fixed cells.

High‐quality single‐cell RNA sequencing data of fixed cells with known level of replication stress

3.2

After identification of γH2AX as an RS‐marker we directly compared fresh and DSP‐fixed cells to evaluate the extent to which fixation with DSP affects the quality of single‐cell RNA sequencing data. Since the response to RS is affected by cell cycle stage, we decided to sort only cycling cells. Furthermore, we hoped to reduce transcriptional variation due to cell cycle, which is a major confounding factor in single cell analysis [35]. To this end, we made use of the fact that our RPE‐HRAS^G12V^ cells stably expressed the Fluorescent Ubiquitination‐based Cell Cycle Indicator (FUCCI4) system [36]. We sorted RPE‐HRAS^G12V^ cells expressing Geminin_1‐110_, representing S/G2‐phase, with and without treatment with RS‐inducing drugs. Half of the cells were directly sorted (fresh), whereas the other half was first fixed with DSP and stained using the aforementioned γH2AX antibody. All samples were subjected to standard cDNA preparation, including de‐crosslinking, and RNA‐sequencing. After initial quality control (described in methods section), 232 fresh (success rate = 60.42%) and 273 DSP‐fixed (success rate = 71.09%) cells were selected for downstream analysis. The mean number of identified genes (5644 in fresh versus 5044 in DSP‐fixed cells) was comparable (Fig. 2A), as was RNA count, percentage of mitochondrial genes, and spike‐in RNAs (Fig. S1A–C). Moreover, the average gene expression and gene detection rate were not affected by DSP fixation (Fig. 2B,C). In addition, the similar coefficients of variation in the two cell populations indicate that DSP fixation does not negatively impact the ability to detect expression heterogeneity (Fig. 2D).

High quality single‐cell RNA sequencing data of fixed cells with known level of replication stress. (A) Violin plot representing the average numbers of genes detected per cell in fresh (n = 232) and DSP [dithiobis(succinimidyl propionate)] fixed RPE‐HRASG12V (n = 237) cells. (B) Scatter plot showing the average gene expression in DSP‐fixed and fresh cells. R value indicates Pearson Correlation. The red line indicates x = y. (C) Scatter plot showing the correlation between the gene detection rate and average gene expression in fresh and DSP‐fixed cells. R value indicates Pearson Correlation coefficient between fresh and DSP‐fixed cells. (D) Scatter plot showing the correlation between the coefficient of variation and average gene expression in fresh and DSP‐fixed cells. R value indicates Pearson Correlation coefficient between fresh and DSP‐fixed cells. (E) Dimensionality reduction using tSNE (t‐distributed stochastic neighbor embedding) of DSP‐fixed cells. Cells are color coded according to γH2AX signal. (F) Feature plot in which cells on tSNE plot in E are color coded according to mAzami Green‐Geminin1‐110 signal. (G) Feature plots in which cells on tSNE plot in E are color coded according to the expression of S‐phase (CDC6 and E2F1) or G2‐phase (PLK1 and CCNB1) markers.

After confirming that we can obtain high‐quality single cell RNA sequencing data from DSP‐fixed cells, we aimed to identify gene expression programs that mediate the low level of RS in a subset of cells, which could potentially mediate drug resistance. However, when we analyzed cells treated with Chk1i + gemcitabine, the γH2AX positive and negative cells do not clearly cluster apart in the tSNE plot shown in Fig. 2E. Thus, not the level of RS, but other factors account for the clustering within the DSP‐fixed cell population. To assess if cell cycle status can explain the clustering, we plotted the protein level of the FUCCI4 cell cycle marker Geminin_1‐110_. The protein level of Geminin_1‐110_, which gradually increases during S‐phase progression, correlated well with the different cell clusters, suggesting that transcriptional events underlying S and G2 phase account for clustering (Fig. 2F). Accordingly, the expression of early S‐phase (CDC6 and E2F1) and late G2/M‐phase (PLK1 and CCNB1) markers showed that high levels of RS are predominantly present in cells in late S or G2 phase (Fig. 2G). Thus, position within S‐phase can be a major confounding factor when evaluating the transcriptomic response to RS.

Identification of putative genes that confer replication stress tolerance

3.3

To reduce the variation in level of RS caused by cell cycle status, we more stringently selected cells solely in mid‐S‐phase based on the DNA content using DAPI (Fig. 3A). Subsequently, we selected S‐phase cells negative for γH2AX, and S‐phase cells with low, medium or high levels of γH2AX using flow cytometry before and 16 h after treatment with RS‐inducing drugs. As seen previously (Fig. 1C), treatment with CHK1i + gemcitabine increased the level of γH2AX in control and RPE‐HRAS^G12V^ cells, but several cells still maintained low levels of γH2AX (Fig. 3A). Hence, to allow identification of mechanisms that facilitate resistance to RS‐inducing drugs, we collected cells with no, low, intermediate, or high levels of γH2AX staining for single‐cell RNA‐sequencing.

Single‐cell RNA sequencing of S‐phase cells sorted based on γH2AX level. (A) Flow cytometry data of RPE‐HRASG12V and control cells treated for 16 h with 10 nm CHK1i + 4 nm gemcitabine or vehicle (Veh). Sorting strategy is shown: first S‐phase cells were selected based on DAPI (4′,6‐diamidino‐2‐phenylindole) signal. For drug‐treated cells equal number of cells with no, low, medium or high levels of γH2AX were sorted. The percentages of cells in these different categories before sorting are indicated and show an increase in the cell population with high level of RS after treatment with CHK1 inhibitor (CHK1i) + gemcitabine. Respectively one and three 384‐well plates were sorted for each vehicle and drug‐treated group. (B) Dimensionality reduction using tSNE (t‐distributed stochastic neighbor embedding) of all cells (HRASG12V and control) before and after treatment with 10 nm CHK1i + 4 nm gemcitabine shows separate clusters of cells. (C) Feature plot in which cells on tSNE plot in B are color coded based on the different conditions. Inset shows an example of cells from the drug‐treated group clustering among the Veh‐treated cells. (D) Feature plot in which cells on tSNE plot in B are color coded according to γH2AX signal.

In an attempt to identify the influence of oncogenic RAS on transcriptional mechanisms of resistance, we similarly treated RPE‐HRAS^G12V^ and control RPE cells with CHK1i + gemcitabine and sorted cells with different levels of RS. We subjected these cells to single‐cell RNA sequencing and selected cells with more than 1000 unique RNA counts and expressing at least 500 genes for downstream analysis (exact parameters stated in Methods section). Next, we performed principal component analysis and tSNE visualization. This revealed eight distinct clusters of cells that correlated well with the different experimental conditions (Fig. 3B,C). Interestingly, rare cells treated with CHK1i + gemcitabine are located within the untreated cell cluster (Fig. 3C, example shown in see inset), potentially representing non‐damaged, RS‐tolerant cells. Moreover, CHK1i + gemcitabine treated cells in cluster 0 and 4 display lower levels of RS compared to cells in cluster 2 and 3 (Fig. 3B–D).

We suspect that cells able to withstand DNA damage during replication stress represent cells within a tumor that could survive treatment with RS‐inducing drugs. In an attempt to rule out the aforementioned confounding influence of position within S‐phase, we compared the DAPI signals, indicative of S‐phase progression, in cells with different levels of γH2AX signal. The DAPI signals were significantly reduced in γH2AX^negative^ cells compared to the other groups in both RPE control and RPE‐HRAS^G12V^ cells (Fig. S2A). Because we suspected that the absence of RS in γH2AX^negative^ cells could be attributed to their earlier position in S‐phase, we chose to excluded these groups from further analysis. The other groups all showed similar DAPI intensities, although γH2AX^low^ RPE‐HRAS^G12V^ cells showed a modest but statistically significant reduction compared to their γH2AX^high^ counterparts (Fig. S2A).

Differential expression analysis revealed 19 genes that were significantly downregulated in γH2AX^low^ RPE‐HRAS^G12V^ cells compared to γH2AX^high^ cells, suggesting that elevated levels of these genes are associated with sensitivity to RS‐inducing drugs (Fig. 4A and Table S1). 18 genes had significantly higher expression in γH2AX^low^ cells and thus correlated with RS tolerance (Fig. 4B and Table S1).

Identification of genes correlated with levels of replication stress. (A) Heatmap of genes differentially expressed and downregulated in RPE HRASG12V γH2AXlow versus γH2AXhigh control RPE cells 16 h after treatment with 10 nm CHK1 inhibitor (CHK1i + 4 nm gemcitabine). (B) Heatmap of genes differentially expressed and upregulated in γH2AXlow versus γH2AXhigh RPE HRASG12V cells 16 h after treatment with 10 nm CHK1i + 4 nm gemcitabine. (C) FOXM1 gene expression as measured by single‐cell RNA sequencing in RPE‐HRASG12V cells with low, medium, and high γH2AX levels after 16 h of treatment with 10 nm CHK1i + 4 nm gemcitabine. The upper and lower edge of the boxplot indicates the 1st and 3rd quartile, the middle line indicates the median, and the upper and lower edges of the whiskers indicate the 95% confidence interval. No significant differences were identified based on the Wilcoxon rank sum tests for multiple group comparisons. (D) Correlation between pFOXM1 and γH2AX in RPE‐HRASG12V cells treated for 16 h with 10 nm CHK1i + 4 nm gemcitabine, as measured by flow cytometry. The trendline represents the data fit to a linear model, with the shaded region around the trendline indicating the 95% confidence interval. R value indicates Pearson Correlation coefficient. (E) Correlation between DAPI (4′,6‐diamidino‐2‐phenylindole) and pFOXM1, as measured by flow cytometry in the same experiment as described under D. The trendline represents the data fit to a linear model, with the shaded region around the trendline indicating the 95% confidence interval. R value indicates Pearson Correlation coefficient. Plots in A–C show n = 188, 230, and 210 cells in the γH2AXlow, γH2AXmedium, and γH2AXhigh groups, respectively.

We then used Enrichr to investigate whether specific transcriptional programs were differentially activated in cells with high vs low levels of γH2AX (Table S3). In two separate databases, FOXM1 was returned among the top 5 transcription factors associated with genes downregulated in γH2AX^low^ cells. FOXM1 itself was not identified as a differentially expressed gene between cells with high vs low levels of γH2AX (Fig. 4C). This could mean that cell‐to‐cell variation in FOXM1 occurred at the protein rather than transcription level. To test this hypothesis, we co‐stained HRAS^G12V^ cells with γH2AX and phospho‐T600 FOXM1. FOXM1 phosphorylation is necessary for its transcriptional activity [37]. Consistent with the observed increase in FOXM1 target genes, we found that level of active FOXM1 protein correlates with the degree of DNA damage as measured by γH2AX in S‐phase cells (Fig. 4D). We did not observe a correlation between active FOXM1 and DAPI in cells gated in S‐phase, suggesting that the correlation between active FOXM1 and γH2AX is not confounded by cell cycle position (Fig. 4E).

Next, we evaluated if the genes differentially expressed in γH2AX^high^ versus γH2AX^low^ cells are coexpressed (Fig. S2B). We focused our analysis only on genes that showed a correlation coefficient greater than or equal to 0.4 with at least on other gene (genes exceeding this threshold are highlighted). Among genes downregulated in γH2AX^low^ cells, ANLN, HMGB2, CENPE, MKI67 and UBE2C are coexpressed, which is expected as they are all regulated by the FOXM1 transcription factor. We also observed a cluster of coexpressed genes containing SERPINE1, F3, CTGF, CYR61, ANKRD1, KRTAP2‐3, UGCG, and AMOTL, but our Enrichr analysis did not return any transcriptional pathways linked to these genes. Among the genes upregulated in γH2AX^low^ cells, only LGALS1 and MYL6 were correlated. Minimal correlation was observed in the rest of the genes, indicating that the majority of these genes are regulated independently of each other.

Notably, several FOXM1 target genes were also found to be downregulated in the γH2AX^low^ RPE control cells that lack expression of oncogenic RAS (Fig. S2C, Table S2). We performed a similar analysis using Enrichr with this gene set and again identified FOXM1 as a putative transcriptional program that is downregulated in the γH2AX^low^ cells (Table S4). This suggests that reduced activation of this transcriptional program in cells with decreased γH2AX levels is a general phenomenon and not necessarily linked to oncogene expression.

Altogether, these data indicate that a subset of oncogenic RAS‐expressing cells is protected from RS upon treatment with RS‐inducing drugs and that these cells transcriptionally diverge from drug‐sensitive cells, with many differentially expressed genes targeted by the transcription factor FOXM1.

Validation of putative RS tolerance mechanisms

3.4

Next, we assessed if the aforementioned genes that were differentially expressed in γH2AX^low^ versus γH2AX^high^ RPE‐HRAS^G12V^ cells could be functionally responsible for RS sensitivity and RS tolerance. To this end we knocked down these genes individually prior to treatment with CHK1i + gemcitabine and analyzed if this affected RS. We hypothesized that knocking down sensitizing genes would result in a decrease in replication stress while knocking down tolerizing genes would result in an increase in RS upon treatment with CHK1i + gemcitabine (Fig. 5A). We excluded P53 target genes (MDM2, SERPINE1, SNX5, FOSL1, MT‐CO3) as the key role of P53 in the RS response is well‐established [3]. Moreover, individual FOXM1 target genes (CENPE, UBE2C, HMGB2, ANLN, MKI67) were initially excluded from further analysis and replaced by knockdown of FOXM1 itself to address the role of this entire transcription program in the RS response [38, 39].

*Partial knockdown of FOXM1 improves tolerance to replication stress without affecting cell proliferation. (A) Schematic representation of the experimental design to identify and validate putative RS‐tolerance genes. (B) Relative FOXM1 expression of RPE‐HRASG12V following transfection with four individual siRNAs targeting FOXM1 (1 nm each) as measured by qPCR. Gene expression was normalized to the average of two housekeeping genes (GAPDH, 18S). Error bars indicate mean ± SEM. Significant differences determined by one‐way ANOVA with Geisser Greenhouse correction followed by Dunnett's multiple comparison test. **P < 0.01, N = 3. Representative of 3 individual experiments. (C) Flow cytometry data of RPE‐HRASG12V cells which were arrested in G1‐phase after 24 h treatment with a CDK4/6i (1 μm) and individual siRNAs against FOXM1. Subsequently, cells were released in the absence or presence of 10 nm CHK1i + 4 nm gemcitabine and harvested after 14 h to enrich for S/G2‐phase cells. DAPI (4′,6‐diamidino‐2‐phenylindole) staining was used to determine cell cycle progression (top row) and γH2AX staining was used to determine the degree of replication stress (bottom row). Number in bottom right corner of bottom row plots indicates the Geminin1‐110

• cells as percentage of the total cells. Number in the top right corner of bottom row plots indicates γH2AX+ cells as a percentage of Geminin1‐110
• cells. Representative of 2 individual experiments. (D) Outgrown colonies of RPE‐HRASG12V cells transfected for 24 h with siRNAs against FOXM1 followed by treatment for 48 h with the 2 nm CHK1i and 4 nm gemcitabine. After removing drug‐containing media, colonies were allowed to grow for 5 days. Representative of 2 individual experiments. (E) Quantification of colonies presented in panel D, quantified as the IntensityPercent (which takes into account both the area covered by the cell growth as well as the pixel intensity of the covered area) and presented as IntensityPercent relative to untreated cells for each transfection condition. IntensityPercent was quantified using the ImageJ ColonyArea plug‐in. Error bars indicate mean ± SEM. Significant differences were determined by ordinary One‐way ANOVA followed by Dunnett's multiple comparison test. *P < 0.05, *P < 0.01, N = 3.

First, we confirmed that small interfering RNAs (siRNAs) targeting the putative RS‐sensitizing mechanisms efficiently depleted their target gene (Fig. S3A). For initial siRNA knockdowns, we used siRNA Smartpools, which consist of four unique siRNAs targeting the same gene. Since RS‐inducing drugs only affect replicating cells, we sought to enrich the cell population for cells in late S/G2 phase of the cell cycle at the time of analysis. To this end, we depleted the gene of interest using siRNA and arrested all cells in G1‐phase using the CDK4/6 inhibitor palbociclib. Subsequently, we released the cells from palbociclib in the presence of CHK1i + gemcitabine and evaluated the level of γH2AX 14 h after release, when most cells were in late S/G2 phase (Fig. S3B). In addition to enriching for late S/G2 phase cells, this approach ensured that all cells start DNA replication in the presence of CHK1i + gemcitabine. To exclude bias from cells that failed to enter S‐phase after palbociclib release when evaluating the level of RS, we calculated the γH2AX‐positive cells as percentage of the Geminin_1‐110_ positive (i.e., S/G2‐phase) cells.

We first evaluated if depletion of genes upregulated in γH2AX^low^ cells could resensitize cells to RS (Fig. 5A). While depletion of most putative RS‐tolerance genes using Smartpool siRNAs did not affect RS, depletion of MYL6, PRDX5 and ARL4C increased RS induced by CHK1i + gemcitabine (Fig. S3C). However, when we then used individual siRNAs to deconvolute the effects of the Smartpools, none of the three targets showed a consistent RS‐sensitizing effect, suggesting off‐target effects of the individual siRNAs (data not shown).

We then shifted our focus to genes that may make cells more sensitive to RS. AMOTL2, CTGF and CORO1C were downregulated in γH2AX^low^ cells and knockdown reduced the fraction of cells with severe RS (Fig. S3C). This suggests that these genes sensitize cells to RS‐inducing drugs. Similarly, knockdown of FOXM1 bolstered resistance to RS induced by CHK1i + gemcitabine (Fig. S3C). Knockdown of CTGF and FOXM1 had the greatest effect on yH2AX levels and also appeared to alter cell cycle progression, with fewer cells in S and G2 phase. We therefore tested individual siRNAs against these two genes in the subsequent deconvolution step, but only FOXM1 knockdown showed consistent phenotypes with individual siRNAs.

We observed varying levels of FOXM1 knockdown with the four siRNAs against FOXM1 present in the siRNA Smartpool, with two siRNAs accomplishing a near‐complete knockdown (< 90%) and two accomplishing a modest knockdown (50–60%) of the gene (Fig. 5B). The siRNAs that produced stronger FOXM1 knockdown also slowed cell cycle progression, as seen by a higher proportion of cells still in S‐phase and a less distinct G2 peak at 14 h following release from G1 arrest (Fig. 5C, top row). Despite the differential degrees of FOXM1 expression and cell cycle progression, all four knockdowns resulted in decreased levels of γH2AX relative to cells transfected with scrambled siRNA after 14 h of CHK1i + gemcitabine treatment in synchronized RPE‐HRAS^G12V^ cells (Fig. 5C, bottom row; quantified in Fig. S4A).

Next, we investigated the effect of knocking down the FOXM1 target genes that we originally identified in our single cell RNA sequencing data. Our strategy for investigating the target genes was the same as in the original screen: first assess the effect of the siRNA SmartPools, followed by assessment of individual siRNAs if the SmartPools show an effect (Fig. 5A). As in the screen, we first synchronized the cells in G1 with the CDK4/6 inhibitor palbociclib, and then released cells into media containing RS‐inducing drugs. Cells were collected 14 h later. Initially we saw a protective effect of all the target genes (Fig. S5A). All cells transfected with the siRNA Smartpools showed efficient knockdown of their target genes (Fig. S5B). When we performed the deconvolution step, only UBE2C and MKI67 showed a consistent effect of all 4 siRNAs (Fig. S5C) despite all siRNAs effectively reducing gene expression (Fig. S5D and data not shown). Interestingly, these knockdowns also had no effect on the amount of geminin‐positive cells, suggesting that they do not strongly influence cell cycle progression through S and G2 phase.

Finally, we performed clonogenic survival assays in to evaluate long‐term fates of cells with varying FOXM1 gene expression levels during treatment with CHK1i + gemcitabine (Fig. 5D,E). The transient effect of siRNA knockdown allowed us to focus on the consequence of reduced FOXM1 activity at the time of treatment (mimicking stochastic variation in expression levels that may occur in cancer cells) rather than a permanent change in expression levels. We quantified colony outgrowth using the IntensityPercent value calculate by the ImageJ ColonyArea plugin, which takes into account both the area covered by cells as well as the density of the cells within a colony [27]. The relative number of viable cells in the drug‐treated group relative to the vehicle‐treated group for each condition is shown in Fig. 5E and the absolute IntensityPercent values are shown in Fig. S4B. In untreated cells, the greater knockdown (siRNA#3 and siRNA#4) resulted in fewer colonies, which is expected since knockdown of FOXM1 slows cell proliferation (Fig. S4B). The partial knockdown had little to no effect on colonies in untreated cells, consistent with the minimal effect on cell cycle progression (Fig. S4B). In the scrambled condition, treatment with CHK1i + gemcitabine nearly suppressed the outgrowth of colonies once the drugs were removed. However, in the partial knockdowns, several dense colonies were able to recover even following drug exposure, suggesting that the partial knockdown of FOXM1 protected cells from RS‐induced cell death without limiting cell proliferation.

To summarize, these data indicate that partial knockdown of FOXM1 or its transcriptional targets UBE2C and MKI67 protects against drug‐induced RS while still allowing cells to proliferate. This suggests that cancer cells may use similar tuning of the FOXM1 expression program to resist the effects of RS‐inducing drugs without fully halting growth.

Discussion

4

In this study, we used single cell RNA sequencing to investigate transcriptional heterogeneity associated with differential responses to RS. We employed a novel technique wherein cells are reversibly fixed, allowing for cell sorting based on intracellular staining, and subsequent single cell RNA sequencing on sorted cells. By comparing cells with low levels of RS (as measured by γH2AX following treatment with RS‐inducing drugs) to cells with high levels of RS, we found that a moderate reduction in gene expression of downstream targets of transcription factor FOXM1 protects cells against RS induced by CHK1i + gemcitabine without significantly delaying cell proliferation.

Compared to popular screening approaches such as CRISPR or RNA interference libraries, our approach has caveats. First, a differentially expressed gene in RS‐tolerant versus RS‐sensitive cells can be either a cause or a consequence of the observed phenotype. Our observation that P53 target genes were downregulated in RS‐tolerant cells illustrates this issue (Fig. 3E, Fig. S2C). The far majority of our potential hits did not show consistent phenotypes when knocked down with siRNA, suggesting that they do not play a direct role in regulating RS. Second, genes may act in concert: the effect of the variation in expression of an individual gene can be minimal, while up‐ or down‐regulation of multiple genes can have a tolerizing effect. Further studies altering the levels of multiple genes at once would be necessary to test this hypothesis. Nevertheless, our single‐cell transcriptomics approach has the advantage over the aforementioned screening approaches that it has the potential to detect the effects of stochastic heterogenous transcription events within the physiological range.

Recent studies show that high FOXM1 activity facilitates unscheduled mitotic entry to cause RS‐induced mitotic catastrophe [16, 40, 41, 42]. In these contexts, FOXM1 likely primes cells to enter mitosis by inducing the transcription of a large set of mitotic genes, including CCNB1 and CDK1, which allows for sufficient accumulation and activation of cyclin B‐CDK1 complexes to enter mitosis [43, 44]. Curtailing CDK1 activation reduces sensitivity to RS response inhibitors because cells are less prone to enter mitosis prematurely with DNA damage sustained during RS. This can be reversed by inhibiting key regulators of CDK1 activity, such as WEE1, thus reactivating CDK1 [20].

In the present study, we observe that reducing FOXM1 expression improves cell viability (as seen by colony outgrowth) and decreases the level of DNA damage incurred in response to RS‐inducing drugs. While these responses are consistent with aforementioned published results, we assert that this mechanism is distinct from the prevention of CHK1i‐induced premature mitosis previously described [16], as we previously showed using live cell imaging that RPE‐HRAS^G12V^ cells do not respond to CHK1 inhibition by attempting a premature mitosis [25]. Instead these cells accumulated DNA damage during S‐phase and underwent a cell cycle exit from G2. One possible explanation of our observations could be that high FOXM1 expression contributes to DNA damage accumulation by triggering excessive origin firing during CHK1 inhibition. An important function of CHK1 is to mitigate DNA damage by reducing firing of late origins under conditions of replication stress [9]. Cyclin A2‐CDK1 complexes mediate origin firing, and CCNA2 is a FOXM1 target [45], thus FOXM1‐induced CCNA2 expression could exacerbate the increase in late origin firing permitted by CHK1 inhibition. Consistent with this, a recent study using ATR‐deficient B cells showed that RS triggered by loss of ATR could be reversed by suppressing origin firing, which was accomplished through partial inhibition of CDC7 and CDK1 activity [46].

Although origin firing could play a role, the impact of high FOXM1 expression seems multifactorial. We saw that knockdown of the FOXM1 targets MKI67 and UBE2C also resulted in a decrease in accumulation of DNA damage after CHK1 inhibition. This is remarkable because both these genes have described roles in mitosis [47, 48, 49], but their potential link to replication stress is less clear. A recent study showed that Ki67 interacts directly with replication machinery and that loss of Ki67 can result in unloading of origins [50]. In the context of untreated cells, this resulted in DNA damage, as measured by comet assay. However, in our experiment (Fig. S5C), origin unloading due to loss of Ki67 may counteract the excessive origin firing that occurs due to CHK1 inhibition. As a member of the APC/C^Cdh1^ complex, UBE2C is important for degrading multiple substrates involved in mitosis and its degradation of substrates in G1 is necessary for permitting entry into S‐phase [47, 51, 52]. However, as APC/C^Cdh1^ is inactive during S‐ and G2‐phase, our results could mean that UBE2C has a different role in coordination of DNA replication. Thus, the role of UBE2C in replication stress requires further exploration.

Conclusions

5

Our work describes a model in which transcriptomic variability of the transcription factor FOXM1 endows a subset of cells within a population of genetically identical cells tolerance to drug‐induced replication stress. While it may be hard to therapeutically target FOXM1 to improve efficacy of intra‐S‐phase checkpoint inhibitors, overexpression of the FOXM1 program can potentially serve as a biomarker since amplification of the FOXM1 gene occur relatively frequently in multiple cancers [53], although single‐cell analysis would need to reveal the relative heterogeneity within tumors. Furthermore, our findings support the idea that decelerated S‐phase progression could counteract CHK1 inhibitors, which also suggests that pharmacologically accelerating cell cycle progression may work to sensitize cells to this class of drugs. An excellent example of this is inhibiting WEE1, the kinase responsible for preventing CDK1/2 activation or its relative PKMYT1, which inhibits CDK1. WEE1 and PKMYT1 inhibitors force cell proliferation in the presence of RS and—at least in part—overcome resistance to intra‐S‐phase checkpoint inhibitors [20, 42, 54].

Conflict of interest

The authors declare no conflict of interest.

Author contributions

H.A.S. conceived and performed experiments, analyzed data, and wrote the manuscript. K.A.W. conceived and performed experiments, analyzed data, and wrote the manuscript. E.A.v.L performed experiments and analyzed data. F.M.R. analyzed single‐cell RNA‐sequencing data. B.W. conceived and oversaw the study and wrote the manuscript.

Peer review

The peer review history for this article is available at https://www.webofscience.com/api/gateway/wos/peer‐review/10.1002/1878‐0261.13819.

Supporting information

Fig. S1. Related to Fig. 2. (A) Violin plot showing the number of unique RNA (UMI) counts per cell in DSP‐fixed and fresh cells RPE‐HRAS^G12V^ cells. (B) Violin plot showing the percentage of RNA counts mapping to mitochondrial genes as percentage of total counts detected in DSP‐fixed and fresh cells RPE‐HRAS^G12V^ cells. C Violin plot showing ERCC spike‐in RNA counts as percentage of total counts detected in DSP‐fixed and fresh RPE‐HRAS^G12V^ cells.

Fig. S2. Related to Fig. 4. (A) The DAPI fluorescence intensity in individual FACS‐sorted cells shown in a violin plot. Horizontal lines represent median DAPI levels in each group. Wilcoxon rank sum tests for multiple group comparisons were performed on control cells and HRAS^G12V^ cells separately. ***P < 2E‐16 γH2AX‐negative versus all other groups, *P < 0.005 γH2AX^low^ versus γH2AX^medium^ and γH2AX^high^. (B) Correlation matrix displaying the correlations in normalized transcript counts between all the differentially expressed genes in RPE HRAS^G12V^ γH2AX^high^ versus γH2AX^low^ cells. Highlighted genes show a correlation coefficient greater than 0.4 with at least 1 other gene. (C) Heatmap of genes differentially expressed and downregulated in γH2AX^low^ versus γH2AX^high^ control RPE cells after treatment with 10 nM CHK1i + 4 nM gemcitabine. The heatmaps represent normalized transcript counts in single‐cell RNA‐sequencing analysis. (D) Heatmap of genes differentially expressed and upregulated in γH2AX^low^ versus γH2AX^high^ control RPE cells after treatment with 10 nM CHK1i + 4 nM gemcitabine. The heatmaps represent normalized transcript counts in single‐cell RNA‐sequencing analysis. Plots in C and D represent n = 141, 194, and 214 cells in the γH2AX^low^, γH2AX^medium^, and γH2AX^high^ groups, respectively.

Fig. S3. Related to Fig. 5A. (A) Quantitative PCR of the expression of potential RS‐tolerance conferring genes in RPE‐HRAS^G12V^ cells treated with scrambled siRNA or siRNA targeting the gene of interest. Gene expression was normalized to the average of two housekeeping genes (GAPDH, 18S). Bar represents mean ± s.e.m. (B) Flow cytometry data of RPE‐HRAS^G12V^ cells unsynchronized, arrested in G1‐phase after 24 h treatment with a CDK4/6i and at indicated hours after release in the presence and absence of CHK1i + gemcitabine to enrich for S/G2‐phase cells. DAPI staining was used to determine cell cycle progression (top row). The relationship between Geminin_1‐110_ and DAPI is shown in the bottom row. Representative of 2 independent experiments. (C) Flow cytometry data of RPE‐HRAS^G12V^ cells with the indicated genes depleted by Smartpools of four individual siRNAs. DAPI staining was used to determine cell cycle progression (top row) and γH2AX staining was used to determine the degree of replication stress (bottom row). Number in bottom right corner of bottom row plots indicates the Geminin_1‐110_ ^+^ cells as percentage of the total cells. Number in the top right corner of bottom row plots indicates γH2AX^+^ cells as a percentage of Geminin_1‐110_ ^+^ cells. Representative of 2 independent experiments.

Fig. S4. Related to Fig. 5C and E. (A) Quantification of percent of γH2AX‐positive cells in geminin‐positive cells (representative of cells in S/G2 phase) from two individual experiments. Error bars indicate mean +/− SEM. *P < 0.01. (B) Raw IntensityPercent values calculated using the ImageJ ColonyArea plug‐in. Error bars indicate mean +/− SEM. Significant differences were determined by ordinary One‐way ANOVA followed by Dunnett's multiple comparison test. *P < 0.05, **P < 0.01, N = 3.

Fig. S5. Related to Fig. 5. (A) DNA damage, measured by γH2AX flow cytometry in S and G2 RPE‐HRAS^G12V^ cells treated with indicated FOXM1 target genes depleted by Smartpools of four individual siRNAs. DAPI staining was used to determine cell cycle progression (top row) and γH2AX staining was used to determine the degree of replication stress (bottom row). Number in bottom right corner of bottom row plots indicates the Geminin_1‐110_ ^+^ cells as percentage of the total cells. Number in the top right corner of bottom row plots indicates γH2AX^+^ cells as a percentage of Geminin_1‐110_ ^+^ cells. (B) Quantitative PCR of the expression of FOXM1 target genes in RPE‐HRAS^G12V^ cells treated with scrambled siRNA or siRNA Smartpools targeting the gene of interest. Gene expression was normalized to the average of 2 reference genes (GAPDH, 18S). Bars represent mean ± s.e.m. (C) DNA damage, measured by γH2AX flow cytometry in G2 RPE‐HRAS^G12V^ cells treated with four individual siRNAs targeting the FOXM1 target genes MKI67 and UBE2C. DAPI staining was used to determine cell cycle progression (top row) and γH2AX staining was used to determine the degree of replication stress (bottom row). Number in bottom right corner of bottom row plots indicates the Geminin_1‐110_ ^+^ cells as percentage of the total cells. Number in the top right corner of bottom row plots indicates γH2AX^+^ cells as a percentage of Geminin_1‐110_ ^+^ cells. (D) Quantitative PCR of the expression of the FOXM1 target genes MKI67 and UBE2C in RPE‐HRAS^G12V^ cells treated with scrambled siRNA or 4 different individual siRNAs targeting the gene of interest. Gene expression was normalized to the average of 2 reference genes (GAPDH, 18S). Bars represent mean ± s.e.m.

Table S1. Differentially expressed genes in γH2AX^low^ versus γH2AX^high^ cells with HRAS^G12V^ overexpression. Table S2. Differentially expressed genes in γH2AX^low^ versus γH2AX^high^ control cells. Table S3. Enrichr results of downregulated genes in γH2AX^low^ HRAS^G12V^ cells. Table S4. Enrichr results of downregulated genes in γH2AX^low^ Control cells. Table S5. Key resources. Table S6. Antibodies for immunoblots and immunofluorescence staining. Table S7. qPCR primers.