Early and late effects of pharmacological ALK inhibition on the neuroblastoma transcriptome

Background Neuroblastoma is an aggressive childhood malignancy of the sympathetic nervous system. Despite multi-modal therapy, survival of high-risk patients remains disappointingly low, underscoring the need for novel treatment strategies. The discovery of ALK activating mutations opened the way to precision treatment in a subset of these patients. Previously, we investigated the transcriptional effects of pharmacological ALK inhibition on neuroblastoma cell lines, six hours after TAE684 administration, resulting in the 77-gene ALK signature, which was shown to gradually decrease from 120 minutes after TAE684 treatment, to gain deeper insight into the molecular effects of oncogenic ALK signaling. Aim Here, we further dissected the transcriptional dynamic profiles of neuroblastoma cells upon TAE684 treatment in a detailed timeframe of ten minutes up to six hours after inhibition, in order to identify additional early targets for combination treatment. Results We observed an unexpected initial upregulation of positively regulated MYCN target genes following subsequent downregulation of overall MYCN activity. In addition, we identified adrenomedullin (ADM), previously shown to be implicated in sunitinib resistance, as the earliest response gene upon ALK inhibition. Conclusions We describe the early and late effects of ALK inhibitor TAE684 treatment on the neuroblastoma transcriptome. The observed unexpected upregulation of ADM warrants further investigation in relation to putative ALK resistance in neuroblastoma patients currently undergoing ALK inhibitor treatment.


INTRODUCTION
Neuroblastoma (NB) is a childhood cancer arising from malignant transformation of cells from the sympathoadrenergic lineage of the developing neural crest [1] and is characterized by a large clinical heterogeneity, including patients with tumours that spontaneously regress as well as patients with metastasis and refractory disease despite intensive therapy regimens [2]. For these high-risk neuroblastoma cases, a better understanding of the genes and pathways that are involved in disease development and progression is currently fuelling the identification of new molecular targets for therapy [3].
Overall, given the low frequency of mutations in neuroblastoma, options for targeted therapy are relatively limited. However, the recurrent somatic mutations of
Despite the promising therapeutic potential of ALK inhibitors, an important drawback is the rapid occurrence of resistance due to escape mechanisms including secondary mutations or activation of other kinase signaling pathways [17,[20][21][22][23][24][25][26][27][28][29][30]. A better understanding of the downstream ALK signaling cascades can inspire the development of more rationally designed combinatorial treatment approaches. In a previous study, ALK downstream signaling was characterized by in depth transcriptome analysis of neuroblastoma cells treated for six hours with the ALK inhibitor TAE684, resulting in the 77-gene ALK signature and the identification of a negative MAPK feedback loop and of RET as ALK activated target gene [31]. Moreover, we generated transcriptome data between 10 minutes and 6 hours after pharmacological ALK inhibition to show that our established 77-gene signature is gradually decreasing starting from 2 hours after the treatment. Here, we further explore this dynamic transcriptome data generated upon treatment with the ALK inhibitor TAE684 in order to look in more detail into the dynamics of early transcriptional responses following ALK inhibition to improve our understanding on downstream ALK signaling and identify novel targets for combination therapy. To this end, we performed differential gene expression analysis at different time points following TAE684 treatment of the neuroblastoma cell line CLB-GA with an activating ALK R1275Q mutation.

MAPK, PI 3 K, RET and MYC(N) signaling pathways are downregulated starting from 1 to 2 hours after TAE684 treatment
In order to gain more insights into the early dynamics of transcriptional changes after pharmacological inhibition of ALK signaling in neuroblastoma cells, the neuroblastoma cell line CLB-GA, which harbours an ALK R1275Q mutation, was treated with the ALK inhibitor TAE684 and DMSO as control. While TAE684 did not go into clinical trials, TAE684 is a valid ALK-inhibitor as we previously showed that the transcriptional responses for this drug and the more clinically relevant ones, such as crizotinib, LDK378 and X396, are very similar [31]. RNA expression profiling was performed on cells harvested 10 and 30 minutes, 1, 2, 4 and 6 hours after treatment (Supplementary Figure 1).
Our previously established 77-gene ALK signature [31] significantly decreases from 2 hours after ALK inhibitor treatment, as shown earlier by Lambertz et al. [31]. In this study, we evaluated the dynamic regulation of known ALK downstream pathways, i.e. MAPK, PI 3 K, RET and MYC/MYCN signaling pathways [4][5][6][7][31][32][33][34], by calculating signature activity scores for drugs specifically blocking these pathways. These drug signature scores summarize the transcriptional response of the components of a given signaling pathway upon pharmacological pathway inhibition. Our data confirm ALK regulation of these pathways and furthermore show that MAPK, PI 3 K and RET pathway activity is decreased as early as 1 to 2 hours after the start of TAE684 treatment, as the corresponding inhibitor scores become activated at those time points ( [32,34,[37][38][39].
To identify other dynamically regulated pathways upon ALK inhibition, functional enrichment analysis was performed on a set of k-means clusters of genes with similar transcriptional responses over the timeframe. After k-means clustering with 9 centers, as defined by the elbow method, we combined clusters with similar expression patterns over time yielding a total of 6 different clusters ( Taken together, this analysis provides insights into the dynamics of inactivation of ALK driven downstream pathways in ALK mutant neuroblastoma cells upon ALK inhibition and validates the applicability of the presented dataset to further scrutinize early temporal effects of ALK inhibition.

MYC(N) signaling is upregulated immediately following ALK inhibition
To functionally scrutinize the transcriptional changes upon ALK inhibitor treatment, Gene Set Enrichment Analysis (GSEA [42]) was performed using the 'Hallmarks v5.0' genesets from the Molecular Signatures Database (MSigDB) (software.broadinstitute. org/gsea/msigdb) at the different time points. This analysis further validates that the signaling pathways KRAS (MAPK) and PI 3 K-mTOR are significantly repressed ( Figure 3A). Remarkably, one of the MYC genesets (HALLMARK_MYC_TARGETS_1) is significantly enriched for genes upregulated 30 minutes after treatment and significantly enriched for genes downregulated 6 hours after pharmacological ALK blockade ( Figure 3A, indicated in red). Moreover, we observed a significant increase in the Fredlund MYC signature [36] 2 hours after ALK inhibition (p = 0.0038) ( Figure 1F, Supplementary Figure 2F). Therefore, we also performed GSEA on other MYC(N) activity signatures and confirmed for 4 other genesets an initial upregulation (at 10 or 30 minutes after ALK inhibitor treatment) followed by a downregulation (from 2 hours after ALK inhibitor treatment) ( Figure  3B). The observation of this initial upregulation of MYCN activity scores is intriguing and needs further investigation.

Adrenomedullin (ADM) is the earliest upregulated gene upon pharmacological ALK inhibition of NB cells
In order to identify the set of early and late genes that significantly change transcriptionally, differential gene expression analysis comparing TAE684 treated cells with DMSO treated controls was performed at each time point. This analysis revealed that no significant transcriptional changes occur as early as 10 or 30 minutes after treatment (using adjusted p-value (False Discovery Rate (FDR)) < 0.05), while starting from 2 hours after treatment an increasing number of significantly, differentially expressed genes are identified (Supplementary Table 1). Our observations are in line with these from Lai et al. [43], who have shown that transcriptional events occur between 2 and 4 hours after inhibition of a receptor tyrosine kinase.
Interestingly, one gene showed significant differential upregulation 1 hour post-TAE684 treatment, i.e. adrenomedullin (ADM) ( Figure 4A; Supplementary Figure 6A). Interestingly, using second-and thirdgeneration ALK inhibitors [31], the upregulation of ADM upon ALK inhibition was confirmed in CLB-GA and Kelly (Supplementary Figure 6B, 6C). Moreover, we could validate ADM upregulation 6 hours after treatment with TAE684 in more NB cell lines with an ALK amplification, an ALK R1275Q or ALK F1174L mutation ( Figure  4B). Furthermore, the early upregulation of ADM was validated in the cell line NB-1 upon pharmacological ALK inhibition with TAE684 ( Figure 4C), as the increase was already detectable 1 hour after the treatment. In addition, ADM expression levels are also significantly induced in the CLB-GA cell line, following PI 3 K/mTOR inhibitor BEZ-235, MEK antagonist trametinib and RET repressor vandetanib treatment [31] (Supplementary Figure 6D). Interestingly, ADM upregulation was also observed in the Anaplastic Large Cell Lymphoma (ALCL) cell line TS treated with three different ALK inhibitors for 6 hours [44] ( Figure 4D).
Furthermore, the observed ADM upregulation is intriguing and has a potentially important translational consequence. In addition to ADM activation upon hypoxia [45][46][47][48][49][50][51], ADM was also shown to act as growth factor, prevent apoptosis-mediated cell death, increase tumour cell motility and metastasis and induce angiogenesis in various cancer types [47-49, 51, 53-59]. Of further importance, it has been reported that hypoxia, which can activate ADM [45][46][47][48][49][50][51], induces resistance to ALK inhibitors in non-small cell lung cancer [60] and that ADM expression levels are increased in tyrosine kinase inhibitor sunitinib resistant renal cell carcinoma [61]. In view of these observations, we evaluated whether combining an adrenomedullin receptor antagonist (ADM22-52) [61] together with TAE684 would sensitize NB cells to ALK inhibition. No additional effects on decrease in cell viability were observed (data not shown). This observation could be due to the fact that elevated ADM levels have no effect on sensitivity for ALK inhibitors in neuroblastoma cells. Alternatively, given that the cell lines tested are very sensitive to the inhibitor with immediate strong and massive apoptotic effects, under these conditions elevated ADM levels may have no immediate effects on cell viability, while in a later stage in cells rendered resistant to ALK inhibition, combined exposure to both ALK inhibition and the ADM inhibitor might have effect.

CLB-GA time series dataset
The CLB-GA time series dataset has been described before [31]. The neuroblastoma CLB-GA cell www.impactjournals.com/oncotarget  Oncotarget 106827 www.impactjournals.com/oncotarget line (ALK R1275Q ) was cultured in RPMI-1640 medium (Invitrogen), supplemented with fetal bovine serum (10%), kanamycin (100 μg/ml) penicillin/streptomycin (100 IU/ml), L-glutamine (2 mM) and HEPES (25 mM) (Life Technologies), and maintained at 37°C in a 5% C02humidified environment. At different time points (10, 30 minutes, 1, 2, 4 and 6 hours) after treatment with 0.32 μM of the small molecule ALK inhibitor NVP-TAE684 (hereafter TAE684) (SelleckChem, S1108) or DMSO (VWR) as solvent control, cells were harvested and RNA quality was analysed using Experion (Bio-Rad) prior to gene-expression profiling. Samples were labelled and hybridized to the Sureprint G3 human GE 8x60K microarrays (Agilent Technologies), according to the manufacturer's guidelines and starting from 200 ng RNA. The data was normalized with the vsn method in R version 3.3.3 (packages vsn and limma). The R package BioMart was used to annotate gene names to their corresponding probes. Probes detecting at least a two-fold higher expression than the negative control probes of the array in at least 60% of the DMSO or TAE samples, were selected as background correction. Data can be accessed through ArrayExpress (accession number E-MTAB-3206) [31].

Signature score analysis
To establish MAPK, PI 3 K/mTOR and RET inhibitor signatures, we used published gene expression profiling data (ArrayExpress accession number E-MTAB-3206) [31]. This dataset contains expression profiling data of the CLB-GA cell line 6 hours after treatment with MAPK inhibitor trametinib, the dual PI 3 K/mTOR repressor BEZ-235 or RET antagonist vandetanib and DMSO as control. Using the limma R-package, differential expression analysis was performed comparing the DMSO-control and the inhibitor treated samples. The established signatures consist of the differentially expressed genes with adjusted p-value (False Discovery Rate (FDR)) < 0.05. In addition, a published MYCN and a published MYC signature were used to establish the MYC(N) pathway activity [35,36]. Signatures score analysis was conducted using a rank-scoring algorithm as described previously [31,36].

Gene set enrichment analysis
Gene set enrichment analysis (GSEA [42]) was performed using the MSigDB 'Hallmarks v5.0' gene sets (software.broadinstitute.org/gsea/msigdb) and an in house compiled gene set catalogue containing all MYC target genesets from this hallmark catalogue as well as publically available MYC(N) activity or target signatures [35,36,[62][63][64][65][66][67]. The genesets, showing positively or negatively enrichment at minimum one time point and with a FDR < 0.1 are plotted in a heatmap. Data was plotted as the mean of the ratio of the TAE684 treated sample and his corresponding control sample at every time point.

Cluster analysis and pathway enrichment analysis in the clusters
Using the Bayesian Information Criterion (BIC) and the elbow method, 9 was defined as the optimal number of clusters in the dataset (Supplementary Figure 4A, 4B). K-means clustering was performed on the expression level ratios of TAE684/DMSO samples for each time point in R (with k = 9). However, visual inspection of the expression patterns of the genes in the clusters showed that some clusters have a similar pattern over time (clusters 1 and 2, 4 and 7, 6 and 9) (Supplementary Figure 4B) and were therefore merged, ending with 6 final clusters (Supplementary Figure 5A). For every cluster, the average response of the genes was plotted as a line (red line). The 4 clusters that are showing a clear dynamic pattern over time, were screened for MSigDB 'c5 Gene Ontology (GO) Biological Process Ontology (BP) v6.0' (software. broadinstitute.org/gsea/msigdb), MSigDB 'c6 Oncogenic Signatures v5.0' and overrepresented Reactome pathways [40,41] using the R-packages clusterProfiler and Reactome Pathway Analysis [68,69].

Differential expression analysis
At every time point, differential expression analysis was performed using the R package limma, comparing the DMSO-control and the TAE684-treated sample. The duplicate experiments, independently generated at different time points, were set as blocking factor. Significantly differentially expressed genes refer to those with an adjusted p-value (FDR) < 0.001. Results are shown in Supplementary Table 2.

Analysis of ADM expression levels after ALK inhibition using qPCR
ADM expression levels were measured in CLB-GA (ALK R1275Q ) and Kelly (ALK F1174L ) cells that were treated for 6 hours with DMSO as control solvent and 0.09 μM of ALK inhibitor, PF06463922 acetate (Sigma-Aldrich, PZ0271) (dissolved in sterile DMSO and diluted to the final concentration in culture medium) and in NB-1 (ALK amp ) cells that were treated for 1 and 6 hours with DMSO as control solvent and 0.32 μM of TAE684 (dissolved in sterile DMSO and diluted to the final concentration in culture medium).
RNA extraction, cDNA synthesis and RT-qPCR of these samples was performed as described earlier [31]. The Cq-values for ADM expression were normalized with www.impactjournals.com/oncotarget data of at least two reference genes (TBP, YWHAZ, B2M, HPRT1 and SDHA) (primer sequences: Supplementary Table 2) using qBasePlus software (Biogazelle).

Statistical analyses
Statistical significance was calculated with GraphPad Prism7 by unpaired one-way ANOVA with Bonferroni correction when comparing more than two unmatched groups, while (un)-paired t-test was chosen when comparing two groups.

CONCLUSIONS
In summary, dynamic expression profiling following ALK inhibition of ALK mutated neuroblastoma cells revealed (1) unexpected early ADM upregulation with potential implications for design of more effective ALK targeted therapy, (2) an initial increase of MYC(N) activity immediately after ALK inhibition and (3) confirmed inhibition after 1 to 2 hours of the ALK downstream MAPK, PI 3 K-AKT, RET and MYC(N) pathways.

ACKNOWLEDGMENTS
We thank Jeroen Schacht, Fanny De Vloed, Jolien Van Laere and Els De Smet for their outstanding technical assistance.

CONFLICTS OF INTEREST
The authors declare no conflicts of interest. and an Emmanuel van der Schueren grant ('Kom op tegen Kanker'). K.D. is supported by Ghent University (BOF; BOF16/PDO/043) and R.C. by a pre-doctoral fellowship of the Research Foundation -Flanders (FWO). The authors would further like to thank the following funding agencies: the Belgian Foundation against Cancer (project 2014-175) to F.S., Ghent University (BOF10/ GOA/019, BOF16/GOA/23) to F.S., the Belgian Program of Interuniversity Poles of Attraction (IUAP Phase VII -P7/03) to F.S., the Fund for Scientific Research Flanders (Research projects G053012N, G050712N, G051516N to F.S) and 'Stichting Villa Joep' to F.S.