Integrated time-course omics analysis distinguishes immediate therapeutic response from acquired resistance

Background Targeted therapies specifically act by blocking the activity of genes that are critical for tumorigenesis. However, most cancers acquire resistance and long-term disease control is rarely observed. Knowing the timing of molecular changes responsible for the development of acquired resistance can enable optimization of alterations to patients’ treatments. Clinically, acquired therapeutic resistance can only be studied at a single time point in resistant tumors. To determine the dynamics of these molecular changes, we obtained high throughput omics data weekly during the development of cetuximab resistance in a head and neck cancer model. Results An unsupervised algorithm, CoGAPS, quantified the evolving transcriptional and epigenetic changes. Further applying a PatternMarker statistic to the results from CoGAPS enabled novel heatmap-based visualization of the dynamics in these time-course omics data. We demonstrate that transcriptional changes resulted from immediate therapeutic response and resistance whereas epigenetic alterations only occurred with resistance. Integrated analysis demonstrated delayed onset of changes in DNA methylation relative to transcription, suggesting that resistance was stabilized epigenetically. Conclusions Genes with epigenetic alterations associated with resistance that had concordant expression changes were hypothesized to stabilize resistance. These genes include FGFR1, which was associated with EGFR inhibitor resistance previously. Thus, integrated omics analysis distinguishes the timing of molecular drivers of resistance. Our findings are a relevant step into the better understanding of the time course progression of changes resulting in acquired resistance to targeted therapies. This is an important contribution to the development of alternative treatment strategies that would introduce new drugs before the resistant phenotype develops.

4 models and therapeutic agents are fundamental tools to overcome issues related to sample availability and 93 serial time point data analysis.

94
Even with advances in experimental sampling approaches, time-course high throughput data alone is 95 insufficient to determine molecular drivers of therapeutic resistance. A novel serial, multi-platform genomics 96 analysis is essential to untangle specific and targetable signaling changes that drive cetuximab resistance in 97 HNSCC. Current supervised bioinformatics algorithms that find time-course patterns in genomic data adjust 98 linear models to correlate molecular profiles with known temporal patterns [12][13][14][15]. However, these 99 algorithms cannot quantify the rate of genomics alterations relative to that of the input phenotype. Other

108
In this study, we used an in vitro HNSCC cell line model to induce resistance and measure the molecular 109 changes using multiple high throughput assays while the resistant phenotype developed. We measured DNA  6 Experimental protocol to establish time-course during acquisition of cetuximab resistance in SCC25

145
The HNSCC cell line SCC25 (intrinsically sensitive to cetuximab) was treated with 100nM cetuximab every 146 three days for 11 weeks (generations G1 to G11). On the eighth day, cells were harvested. Sixty thousand 147 cells were replated for another week of treatment with cetuximab and the remaining cells were separately 148 collected for: (1) RNA isolation (gene expression analysis); (2) DNA isolation (DNA methylation analysis); (3) 149 proliferation assay and (4) storage for future use. All steps were repeated for a total of 11 weeks. In parallel 150 with the cetuximab treated cells, we generated controls that received the same correspondent volume of PBS 151 (phosphate buffered saline). Cells were plated in several replicates each time at the same initial density. The    Resistance to cetuximab was induced in an independent passage of SCC25 cells. After resistance was 171 confirmed, single cells were isolated and grown separately to generate the isogenic resistant single cell 172 clones (CTXR). In total, 11 CTXR clones were maintained in culture without addition of cetuximab. With the 173 exception of one clone (CTXR6), all CTXR clones presented substantial survival advantage compared to the 174 parental SCC25, as reported by Cheng et al. [30].

175
Proliferation assay was performed to confirm cetuximab resistance in the CTXR clones compared to the 176 parental SCC25. A total of 1000 cells were seeded in 96-well plates in quadruplicate for each condition. PBS 177 or cetuximab (10nM, 100nM or 1000nM) was added after 24 and 72 hours and cells were maintained in 178 culture for 7 days. AlamarBlue reagent (Invitrogen, Carlsbad, CA) at a 10% final concentration was incubated 179 for 2 hours and fluorescence was measured according to the manufacturer's recommendations (545nm 180 excitation, 590nm emission). Resistance in the CTXR clones was confirmed when the proliferation rates were 181 higher than in the PBS treated SCC25 cells.    Unless otherwise specified, all genomics analyses were performed in R and code for these analyses is 221 available from https://sourceforge.net/projects/scc25timecourse.

222
The following filtering criterion for genes from the profiling of the time course data from generations of 223 cetuximab treated cells was used. Genes from RNA-seq data were selected if they had log fold change 224 greater than 1 between any two time points of the same condition and less than 2 between the replicate

265
Sources and analysis of human tumor genomics data to G11). Conversely, proliferation of the CTX generations progressively increased over each week ( Fig. 1).

287
Relative to the untreated controls, the growth of the treated cells was initially (CTX-G1) inhibited until CTX-

288
G3. Starting at CTX-G4, the cells became resistant to the anti-proliferative effects of cetuximab and gained 289 stable growth advantages compared to the untreated controls.

12
Comparison of proliferation rates between generations of CTX treated cells relative to generations of cells 291 treated PBS enabled us to conclude that cell growth advantages arise from chronic cetuximab treatment and 292 were associated with resistance rather than prolonged cell culturing. We mirrored the changes in proliferation 293 rates with clinical responses seen in HNSCC tumors treated with cetuximab ( Fig. 1, top panel). Specifically,

294
we inferred that the decreased growth rates in CTX-G1 to -G3 represented initial stages of treatment with a 295 decrease in tumor size. Then, the switch from decreased to increased growth rates during CTX-G3 to -G4

296
represented stable disease without tumor outgrowth. Finally, the higher proliferation in cetuximab-treated 297 cells starting at CTX-G4 represented rapid outgrowth after acquired resistance.

298
Because higher proliferation in treated than untreated cells started at CTX-G4, this was the timepoint at which

308
The RNA-seq data hierarchical unsupervised clustering separated genes with expression changes in treated 309 and untreated generations ( Fig. 2A). Clustering analysis of samples (Supplemental

312
Expression changes at these distinct stages were shared between numerous genes. Confounding by 313 changes resulting from immediate therapeutic response made identification of resistance-specific gene 314 expression changes impossible with clustering.

13
Similar separation of stages of cetuximab response were observed in clustering analysis of gene signatures 316 previously described in HNSCC and non-small cell lung cancer cell line models resistant to cetuximab or 317 gefitinib (anti-EGFR small molecule), respectively [37,39] (Supplemental Fig. 4). For these genes, changes 318 during early stages of resistance clustered for CTX-G4 to CTX-G6 as distinct from later stages for CTX G7-319 11. Nevertheless, these gene signatures also clustered samples with gene expression changes at early 320 stages (CTX-G1 to G3) as distinct from samples from PBS treated generations. However, these analyses 321 were insufficient to quantify the relative dynamics of genes associated with immediate response to therapy or 322 subsequent acquired resistance.  Fig. 2). We applied the PatternMarker statistic to define genes that were uniquely associated 345 with each of these patterns. We excluded the technical, flat pattern to focus with genes with expression 346 changes. By design, genes selected with the PatternMarker statistic are selected to not be multiply regulated 347 regulated. Therefore, limiting the heatmap to these genes enabled visualization of the dynamics of gene 348 expression changes in our time-course dataset (Fig. 2B). The relative magnitude of CoGAPS pattern weights 349 for each sample quantified the dynamics of gene expression changes (Fig. 2C).

350
Similar to the separation seen with clustering (Supplemental

359
However, the transcriptional changes in this pattern were not associated with acquired resistance to 360 cetuximab, and even decreased modestly as resistance developed. Further, enrichment by transcription 361 factor AP-2alpha targets (TFAP2A; one-sided p-value of 0.05) confirmed previous work indicating that 362 transcription by AP-2alpha is induced as an early feedback response to EGFR inhibition [44]. Based upon 363 these findings, we concluded that pattern 1 was associated with immediate response to cetuximab although it 364 includes genes that were also associated with cetuximab resistance in previous studies.

365
The second CoGAPS expression pattern quantified divergence of the cetuximab treated cells from controls at 366 generation CTX-G4 (expression pattern 2, Fig. 2B and Fig. 2C, middle) which was the time point that 367 15 cetuximab treated cells presented significant and stable growth advantage over PBS controls (Fig. 1).

368
Therefore, expression pattern 2 obtained gene expression signatures associated consistently with the 369 development of cetuximab resistance. Gene set statistics of transcription factor targets of EGFR on CoGAPS 370 gene weights were significantly down-regulated in this acquired resistance pattern (Supplemental Fig. 6).

371
One striking exception was c-Myc, which trended with acquired resistance (p-value of 0.06), consistent with to immediate therapeutic response from those specific to acquired resistance.

398
Changes in DNA methylation inferred with CoGAPS were associated with resistance to cetuximab, 399 but not the immediate response to treatment observed in gene expression

400
To determine the timing of the methylation changes associated with acquired resistance, we also measured 401 DNA methylation in each cetuximab generation of SCC25 cells and PBS controls (Fig. 3A).  Fig. 7A). We 410 observed that the gene expression changes associated with acquired resistance occurred gradually and were 411 evident in early generations (Fig. 4, top). The DNA methylation was consistent in cetuximab treatment and 412 control PBS in DNA methylation patterns 2 and 3 during early generations. Then, rapid accumulation in DNA 413 methylation changes started after generations CTX-G4 and CTX-G5 in both patterns 2 and 3 (Fig. 4, bottom),

414
concurrent with the onset of the observed growth advantage over the PBS control (Fig. 1). Changes in DNA 415 methylation were delayed relative to those of gene expression in acquired cetuximab resistance (Fig. 4  The gene signatures from the CoGAPS resistance patterns for expression and methylation had low 421 correlation (Supplemental Fig. 7B)  cetuximab resistance (methylation patterns 2 and 3). We then applied correlation analysis to determine genes 432 that were epigenetically silenced. Specifically, we performed correlation analysis between DNA methylation 433 and gene expression for each of the DNA methylation PatternMarker genes (Fig. 5). FGFR1 was among 434 these genes. This finding was consistent with previous studies that associate differential expression of  Fig. 8 and 9, respectively) and cellular 446 morphology (Supplemental Fig. 10). Fig. 5A and 5B demonstrate that higher heterogeneity among single cell clones was also observed in the epigenetically regulated PatternMarker genes from the CoGAPS 448 analysis that are shown in Figure 4D. These results suggest that different mechanisms of resistance may 449 arise in the same HNSCC cell line.

450
We hypothesized that epigenetically regulated genes shared along the time course patterns and resistant 451 single-cell clones may implicate common mechanisms acquired during evolution of the stable resistance 452 phenotype. To test this hypothesis, we also performed correlation analysis for each of the epigenetically 453 regulated genes in our resistant set (Fig. 5) in the resistant clones and parental cell lines. Nine of the 454 epigenetically regulated PatternMarker genes also had significantly anti-correlated gene expression and DNA 455 methylation in the stable cetuximab resistant clones (Supplemental Fig. 11). Of these, only FGFR1 was 456 demethylated and reexpressed in a cetuximab resistant clone relative to the parental SCC25 cell line (Fig. 6).