Chronic in vivo insulin analogue treatment

Previously, we have reported the effect of insulin analogues on tumor development [7]. Here, for the chronic exposure experiment, 200 (40 mice per treatment), 8-week-old female p53^R270H/+WAPCre mice were obtained from an in-house breeding project. The point mutation in the tumor suppressor p53 gene corresponds to the Li Fraumeni cancer syndrome mutational hotspot (R273H) in humans. Every other day these mice have been injected (subcutaneously) with either vehicle, insulin, glargine, X10, or IGF1 until tumor development. Once the tumors reached a size of 1 cm³ and 24 h after their last injection the mice were sacrificed. Tumors and other tissues were isolated. For this study, ¼ of the tumor was stored in RNALater (Ambion, Austin, Texas) at 4 °C for RNA isolation. A miRNA isolation kit (Macherey Nagel, Germany) was used to isolate and purify small and large RNA molecules in one fraction. Tumor latency time is defined as the time (in weeks) for the tumor to form, from the start of the experiment to the first time the tumor was palpated.

Single insulin analogue treatment: an animal experiment

To determine the short-term effects of insulin analogue treatment on mammary gland gene expression, a single insulin analogue exposure experiment was performed with 40 (4 mice per treatment/time point) female, 8-week-old inbred FVB/NRj mice (obtained from Janvier, rodent research models, France). This specific mouse strain was used as it is the closest relation to the p53^R270H/+WAPCre mouse strain. Mice received a single subcutaneous injection with either vehicle, insulin, glargine, X10, or IGF1. The mice were sacrificed 1 or 6 h after the injection, blood was collected (mini collect, Greiner/Omnilabo), and MGs were stored in RNALater (Ambion, USA) at 4 °C for RNA isolation. For this, a Nucleospin RNA isolation kit (Macherey Nagel, Germany) was used. For further technical details, please refer to our previous publication [7].

In vitro stimulation experiments

Next, an in vitro experiment was performed to reveal the transcriptomic effects of insulin analogue exposure on a human breast cancer cell line. MCF7 human breast cancer cells with an overexpression of the IGF1R and a stable knockdown of the IR (MCF7 IGF1R) were seeded at 60% confluence and starved for 2 days in 5% charcoal/dextran-treated FBS (CDFBS, Hyclone, USA) containing RPMI 1640 (Gibco, USA) medium. Cells were stimulated (with 10 nM compound) for 1 or 6 h after which RNA was isolated using NucleoSpin® miRNA isolation kit (Machery Nagel, Düren, Germany). Stimulations included: insulin NPH (Insuman Basal, Sanofi Aventis), insulin glargine (Lantus, Sanofi Aventis), insulin X10 (AspB10, Novo Nordisk), and IGF1 (Increlex, Ipsen). For more technical details about cell line generation and characterization or growth factor stimulation, please refer to our previous publication [7].

Next-generation sequencing and gene expression analysis

For the chronic insulin analogue exposure experiment, the quality and integrity of the RNA samples were analyzed using the bioanalyzer with an RNA nanochip. The Ion Total RNA-Seq kit was used to process the samples. Samples were Poly-A selected prior to library preparation. This library preparation included the cDNA synthesis and purification steps with the Ion Total RNA-Seq kit v2 (Life Technologies, UK) according to the manufacturer’s instructions. The Ion PI Template OT2 200 Kit v3 and Ion Sequencing 200 kit v3 (both Life Technologies) were used according to the manufacturer’s instructions for sequencing libraries on the PI chip. Sequence runs were performed on the Ion Proton Sequencer (at ServiceXS, Leiden). PI chip analysis, base calling, and quality checks were performed using the Torrent Server Suite. On average, 40 million reads per sample were sequenced with an average read length of 100 base pairs. No additional trimming or filtering of reads was performed before processing. Reads were aligned to mouse genome build GRCm38-Ensembl using Tophat2 (Version 2.0.10). Reads which could not be aligned using Tophat2 were aligned in an additional step, using Bowtie2 (Version 2-2.10) in the local, very sensitive mode. Tophat2 and Bowtie2 aligned reads were merged into a single .bam file for each sample before further analysis.

Gene expression was quantified using HTSeq-Count (Version 0.6.1) using the default options. Differential gene expression was analyzed for compound versus vehicle treatment and was performed using DESeq2 (Version1.2.10). For this analysis, genes with a read count of <50 reads across samples (average of <1 read per sample) were filtered out before the analysis. For the estimation of individual exon expression analysis, a RPKM table was generated with the read counts normalized to library size and gene length (using DEXSeq version 1.8.0). For the mutational profiling, the reads (unfiltered and untrimmed) were aligned to mouse genome build GRCm38.73-DNA primary assembly using TMAP within the Torrent Suite version 4.0.2. Variant calling was performed using the Torrent Suite Variant Caller version 4.0-r72612 with the settings tuned for the detection of somatic mutations at a low stringency level. The reference genome used was the same as that used for read alignment—GRCm38.73. SnpEff version 3.6c was used to filter and annotate the mutations. The list of mutations was filtered to include only exonic mutations with a quality score higher than 250. Several mutations are found in the exact same position in all tumor samples, probably strain-specific single nucleotide polymorphisms (SNPs). Known SNPs in coding regions for the mouse strain FVB/NJ (the strain most closely related to the ^{p53R270H+/–}WAPCre) were downloaded from the Mouse Genome Informatics database and these mutations were discarded from the list. Mouse homologs of the list of human tumor driver genes [10] were used to define the most clinically relevant mutations.

Phenotypic prediction based on transcriptomic data

To predict the phenotypic characteristics of the treatments of different insulin analogues using their transcriptomics, we constructed a support vector machine (SVM) classifier. We followed a procedure similar to that used in [11] to identify genes whose expression is significantly associated with cancer cell migration and proliferation across 52 breast cancer cell line data. We then applied them to the orthologs of the mouse transcriptomics data to predict the migratory and proliferative potential of 50 mouse mammary tumor samples using LIBSVM [11].

Genome-scale metabolic modeling analysis

We estimated the metabolic fluxes that are most consistent with the transcriptomics data using a computational framework called iMAT [12]. iMAT integrates the transcriptomics, as ‘soft’ constraints, by ternary partitioning the expression to lowly (–1), mediocrely (0), and highly (1) expressed genes. iMAT then attempts to collect the metabolic states that best correspond to these cues, which constructs a mixed integer linear programming (MILP) problem. We applied iMAT to the human genome-scale metabolic model Recon1 [13] with a standard medium condition (RPMI) in a condition-specific manner for the five different insulin analogues treatments. The iMAT predicts (i) the biomass production rate and (ii) the metabolic flux rates of individual metabolic reactions. The biomass production rate is the rate at which the biomass precursors are generated with appropriate proportions, and it is incorporated in the metabolic network model as a putative reaction that takes cell and energetic requirements needed to produce biomass as input [14]. With the metabolic fluxes predicted by iMAT, we performed a pathway enrichment analysis of differentially activated metabolic reactions in X10/IGF1-treated cells to the remainder conditions (insulin/glargine/vehicle). We selected the metabolic pathways whose reactions are significantly enriched in the up-/downregulated group using a hypergeometric test followed by multiple hypotheses correction with the false discovery rate (FDR) 0.05. The predicted biomass production rate does not involve standard deviation because we focused on the metabolic states where the biomass production rate is optimized (thus single-valued).

Statistical analysis

Graphpad Prism version 5.01 software was used for the statistical analysis. All standard error bars in the graphs represent standard deviations. Unpaired two tailed t tests were performed to calculate significance. Multi-experiment viewer (MeV version 4.8.1) was used for the hierarchical clustering analysis.

Genetic profiling of IGF1R signaling-mediated mammary gland tumor formation reveals cellular processes associated with decreased tumor latency

Previously, we chronically treated p53^R270H/+WAPCre mice with insulin and insulin analogues to study the effect of chronic IGF1R stimulation on tumorigenesis (Fig. 1a). We found that MG tumor latency time in X10 and IGF1 (compounds with a high affinity towards the IGF1R) treatment groups was significantly decreased (Fig. 1b). Tumors in mice treated with glargine (an insulin analogue that induces only a mild IGF1R activation in the MG) also developed earlier compared to vehicle- and insulin-treated animals, but this trend was not significant [7]. To further investigate this, we evaluated receptor gene expression levels. We found that gene expression levels of the insulin receptor (both the A and B isoform; ira, irb, respectively) and insulin-like growth factor 1 receptor (igf1r) are significantly upregulated in pre-neoplastic MG tissue of old (50 weeks on average) p53^R270H/+WAPCre mice compared to healthy 8-week-old MGs (Fig. 1c). Also, ira and igf1r levels were slightly upregulated in MG tumors. Interestingly, we found a significant and sevenfold downregulation of the B-isoform of the insulin receptor in MG tumor tissue compared to normal MG tissue. This effect was not treatment specific (Fig. 1c right). Levels of igf1r were only decreased after chronic treatment with IGF1. These data are indicative for an involvement of INSR/IGF1R signaling in tumor development and/or progression. However, the different receptor distributions cannot explain the differences in tumor latency time, since this was not a treatment-related effect.

To gain more insight in to overall differences between various insulin analogues in tumor development and progression, we performed RNAseq NGS of 50 tumors (10 tumors per treatment group). We reasoned that transcriptome analysis would shed light on the cause for the differences in tumor latency time by chronic insulin analogue treatment. From the 40 tumors per treatment group, 10 epithelial to mesenchymal transition (EMT)-like tumors were selected per treatment condition for the transcriptome analysis; these 10 EMT tumors had a tumor latency time that was closest to the median latency time of that specific entire treatment group. The selected tumors have been indicated with colored dots in the Kaplan-Mayer curve (Fig. 1b). In Fig. 1d, a correlation analysis shows the sample-to-sample distance. All transcriptomics data were used to define the correlation within and between the different treatment groups. In spite of the relatively high variation in this in vivo experiment, there was a higher correlation within each treatment group than between the different treatment groups. This suggests the existence of a treatment-specific response under each treatment condition. To define these treatment-specific responses, we focused on four particular hallmarks of cancer (Fig. 1e) that might have been involved in the decreased tumor latency, either by tumor initiation or progression. These pathways were selected in close relation to the signaling pathways downstream of IR and IGF1R, and included: i) sustained growth signaling; ii) tissue invasion and metastasis; iii) deregulated cellular metabolism; and iv) genetic instability and mutations.

Chronic IGF1R activation of the MG results in tumors with a ‘sustained growth signaling’ signature

Proliferative signaling is normally a highly regulated process. In a tumor, the replicative cell homeostasis is deregulated causing sustained growth signaling [9]. Since mice were chronically exposed to insulin analogues that activate (to a lesser or greater extent) the IGF1R, we postulated that the decreased tumor latency time (after X10/IGF1 stimulation) would be a direct result of an upregulated IGF1R signaling pathway. We incorporated mouse orthologs of all human genes (~60 genes) that are directly involved in “INSR and IGF1R signaling” according to Ingenuity Pathway Analysis (IPA) (Fig. 2a). A hierarchical clustering analysis was performed for all 50 MG tumors on the relative expression levels of these genes (Fig. 2b). Two clusters could be defined. In cluster 1, tumors from low mitogenic treatments are enriched, whereas cluster 2 predominantly consists of the X10/IGF1 (high mitogenic compound) treated tumors. These two different clusters could not be linked to specific signaling pathways (e.g., PI3K, MAPK, JAK signaling cascade).

To predict the proliferative potential of all the individual tumors, we next built an SVM classifier using transcriptomic signatures associated with proliferation (Fig. 2ci). No significant effect was observed; however, a trend for the tumors in the IGF1 treatment groups to have an increased proliferative potential could be seen. As a positive control for the SVM simulation, similarly we ran the SVM simulation on our earlier micro-array data of the MCF7-IGF1R cells stimulated with the different insulin analogues. We observed that 1 h after glargine, X10, and IGF1 stimulation a strong and significant proliferative genetic signature could be detected. These results are in line with the proliferative potential as determined with functional in vitro assays using this cell line and insulin analogue exposures [5], which makes us believe that the transcriptomic profiles and used SVM simulation is predictive for the biological effect.

An effort was made to determine the growth rate of the individual tumors by calculating the tumor volume from time of detection until the time of sacrifice (data not shown). The EMT tumors in general were highly proliferative resulting in a very small timeframe (a couple of days to weeks) in which tumor volume could be assessed. Unfortunately, this resulted in just a few data points per tumor. Knowing that tumor growth is very likely not linear during tumor development we did not feel comfortable sharing these results. The proliferative potential was also assessed with ki67 staining on several tumors, but as expected there appeared to be a very small window and all tested tumors had a ki67 score >45% indicative of being highly proliferative (data not shown). Next, we attempted to isolate cells from each tumor; we anticipated that monitoring the growth rate of each mouse mammary gland tumor-derived cell line would resemble the growth rate of its originating tumor. Although it appeared to be rather easy to generate a cell line from each individual tumor, it was quite evident that just a small subset of cells would grow out from these tumors. We feared that this clonal expansion would give a skewed image of true tumor growth, and therefore we decided not to include these data.

Chronic IGF1R activation induces tumors with a more mesenchymal phenotype and a higher migration potential

Next, we assessed the hallmark “tissue invasion and metastasis”. We calculated the average number of tumors per mouse per treatment group. Interestingly, we found a significantly higher number of tumors in the IGF1 group (and a non-significant increase for the chronic X10 and glargine treatment groups) compared to the vehicle-treated mice (Additional file 1). This effect can be both due to an increased incidence of multiple primary tumors as well as an increased formation of metastasis in the IGF1 group. To assess the underlying genetic pathways that might have induced metastasis formation we incorporated all genes that are directly associated with the EMT or the mesenchymal to epithelial transition (MET) according to IPA. The gene expression of the mouse orthologs of these genes is presented in an unsupervised hierarchical clustering (Fig. 3a). Two clear gene groups were observed: group A containing genes associated with epithelial phenotype (Elf5, Serpinb5, Cdh1, Grhl2, Elf3, and Wnt4), and group B containing twelve genes associated with mesenchymal cells (Pdgfrb, Six1, Snai2, Tcf4, Zeb1, Klf8, Wnt5a, Snail, Vim, Twist1, and Cdh2). From this hierarchical clustering, three separate treatment clusters could be discriminated. Cluster 1 consists of tumors that show high expression of epithelial- and low expression of mesenchymal-associated genes; three out of six tumors originate from the vehicle-treated animals. Cluster 2 shows moderate expression levels of epithelial- and mesenchymal-associated genes; interestingly the majority of these tumors were from the glargine treatment condition. Cluster 3 consists of tumors that have a low expression of epithelial and high expression of mesenchymal markers. IGF1 treatments were not present at all in cluster 1 and 2, while IGF1 and X10 treatment groups were over-represented in cluster 3. This suggests that IGF1R-mediated signaling is a driver of tumors with an EMT phenotype. We substantiated the phenotypes for the three different clusters (Fig. 3b). The tumors in cluster 1 predominantly expressed epithelial cells that contain many E-cadherin-positive cells. Tumors of the cluster 2 and 3 phenotype were characterized as EMT tumors predominantly consisting of mesenchymal cells, lack of E-cadherin staining, and clear smooth muscle actin (SMA) staining (data not shown) [7].

To link the transcriptomics data further to the phenotype we performed a SVM simulation to derive predictions for the cell migratory potential of each tumor (Fig. 3ci). Overall, no significant difference in the migration potential was observed for the different treatment groups; potentially the intrinsic high variation of in vivo experiments and smaller group sizes can explain lack of power in this analysis. To evaluate whether IGF1 signaling itself is a strong activator in vitro, we also applied our SVM classifier to our previous MCF7-IGF1R transcriptomics data exposed to the different insulin analogues (Fig. 3cii). IGF1, X10, and glargine caused a significant increase in the migration potential. Altogether, these data suggest that chronic exposure to X10, and especially IGF1, induces tumors with a more mesenchymal phenotype that are possibly more aggressive in terms of their migratory potential.

Mice receiving a chronic X10 and IGF1 treatment develop tumors with a higher Warburg potential

Through the activation of the INSR, insulin and insulin analogues directly affect cell metabolism. Therefore, we wondered whether the chronic treatment with insulin analogues induce tumors with a deregulated cellular metabolism or a higher Warburg potential. To test this we used all genes directly involved in “glycolysis” and in “oxidative phosphorylation” according to IPA. An unsupervised hierarchical clustering of the expression profiles of each tumor revealed two very clear gene groups (Fig. 4a). Strikingly, the first gene group consists of eleven genes that are all directly involved in glycolysis; all the genes in the second cluster are involved in oxidative phosphorylation. Four clear tumor clusters could be detected: cluster 1 with high expression of glycolysis genes and low expression of oxidative phosphorylation genes; cluster 2 with moderate glycolysis; cluster 3 with little glycolysis and moderate oxidative phosphorylation; and cluster 4 with high oxidative phosphorylation and little glycolysis. For the insulin group almost all tumors were related to cluster 3 (little glycolysis and moderate oxidative phosphorylation: 90%). For glargine-related tumors, 70% showed enrichment for genes regulating oxidative phosphorylation. IGF1 and X10 tumors showed increased expression of genes involved in glycolysis (60% and 50%, respectively) compared to the vehicle control. This suggests that chronic insulin signaling drives tumors in a more oxidative phosphorylation mode.

Since IGF1 and X10 were mostly related to glycolysis we performed an independent bioinformatics analysis on the activity of various metabolic processes. This revealed several significantly down- and upregulated metabolic programs in IGF1/X10-derived tumors compared to the vehicle/insulin/glargine tumors (Fig. 4b). Importantly, as expected, oxidative phosphorylation was significantly downregulated in the IGF1/X10 group (P = 0.000107) while glycolysis was significantly upregulated. Again, this analysis was also performed on the microarray data of MCF7-IGF1R cells exposed to the different insulin analogues. Similar to the in vivo data, IGF1 and X10 treatment significantly downregulated the citric acid cycle as well as oxidative phosphorylation (both P < 0.0005) and upregulated glycolysis (P < 0.00001) (Fig. 4c). These data suggest that both IGF1 and X10 promote a Warburg effect in mammary gland tumors (chronic exposure) and human tumor cell lines (a single 1 h exposure).

Biomass production rate is highly upregulated in tumors of chronic X10/IGF1-treated mice

Tumor mice receiving a chronic X10 and IGF1 treatment showed increased expression of glycolysis related genes (Fig. 4) and are associated with a decreased tumor latency time (Fig. 5a). We evaluated whether this difference in metabolic capacity is related to enhanced growth. To test this hypothesis, we predicted the rate of the accumulation of biomass in all these different tumors based on the gene expression profiles of biomass-producing metabolic processes. In Fig. 5b, a bar graph of these results is presented. A highly increased mammary gland tumor biomass production rate was observed in chronic X10- and IGF1-treated tumors. This finding is in agreement with the observed increased proliferative potential in the X10/IGF1 tumors. We hypothesize that the decreased tumor latency time by X10 and IGF1 treatment was caused by enhancing the tumor development rather than interfering with the tumor initiation.

Ezh2 and Hras mutations are enriched in chronic IGF1R-activated tumors

Since IGF1 and X10 promote cell proliferation, this could lead to a manifestation of mutations and consequently modulation of tumor development and progression. To evaluate if chronic IGF1 or X10 treatment affects the number of mutations, a mutational analysis was performed on all tumors based on the 40 million 100 base pair reads for each tumor. The number of different mutations was ~3000 mutations per tumor and no significant difference between treatment groups could be detected. Furthermore, there was no correlation between the number of mutations per tumor and the tumor latency time (best fit correlation values, slope = 0.00001207, r ² = 0.0000099, P value = 0.9827) (Fig. 6a).

Next, we focused on specific clinically relevant mutations that are part of the ~140 known human tumor drivers according to Vogelstein and colleagues [10]. In the 50 mouse tumors that we sequenced, 102 of these tumor drivers were mutated. On average, each tumor had about 35 of these tumor driver mutations and again it seems that the chronic treatment did not affect the overall number of tumor driver mutations per tumor (Fig. 6b). Also, for the tumor driver mutations, there was no correlation between the number of mutations and latency time of each tumor. Furthermore, no treatment-specific effects could be detected for the average number of point mutations, the average number of frame shifts, the average number of start CODON insertions, or the average number of stop CODON insertions (Additional file 2).

We next determined the mutations that are over-represented (>50%) in the X10/IGF1 treatments (Fig. 6c). All mutations are involved in cancer development and/or progression [10, 15]. Although no specific core pathway links these individual mutations, strikingly Ezh2 and Hras were strongly over-represented in IGF1/X10-treated tumors. In a similar way we looked at specific mutations that have been over-represented in insulin analogue treatment altogether and highly under-represented in the vehicle condition. No general mechanism was identified that could explain the enrichment of these mutations in the treatment groups. However, a set of five mutations was not detected at all in the vehicle control, yet these were specifically affected in the insulin treatment group. Overall, these data suggest that the IGF1/X10 enhanced tumor formation is associated with several key candidate cancer drivers that could contribute, including Ezh2 and Hras which are known drivers in human breast cancer [16, 17].