¹H NMR studies distinguish the water soluble metabolomic profiles of untransformed and RAS-transformed cells

DNA constructs and viral transduction

The retroviral pBABE KRASV12, HRASV12 and NRASQ61 DNA constructs were obtained from Addgene. Stable transduction of the pBABE empty vector and the RAS constructs into BEAS2B cells was performed as previously described (Patel et al., 2012). Transduced cells were selected in 2.5 µg/ml puromycin-containing complete culture media for a minimum period of 5–7 days (corresponding to the time taken for untransduced BEAS2B cells to die completely in selection media). Oncoprotein overexpression relative to the control cells was verified via Western blotting as previously described (Patel et al., 2012).

Cell culture

BEAS2B cells were obtained from the American Type Culture Collection. All cells were grown at 37 °C in 21% oxygen and 5% CO₂. BEAS2B cells and their derivative lines were maintained in DMEM:F12 complete base media supplemented with 10% fetal bovine serum and 100 units/ml penicillin-streptomycin. All cell culture reagents were obtained from Life Technologies. For each cell line, ten biological replicates were generated by initially seeding ten different 15 cm dishes (Nunclon) with an initial plating of approximately 1 × 10⁶ cells for the control cells, 1.5 × 10⁶ cells for the HRAS- and KRAS-transformed cells, and 2 × 10⁶ for the NRAS-transformed cells (differences in the initial seedings were used to compensate for differences in cellular growth rates so that by the end of the growth period, approximately the same number of cells for each cell line was obtained). In total, all forty plates were seeded at approximately the same time and were allowed to proliferate for a period of four days with the media changed every 48 h. After four days, the cells were trypsinized for approximately two minutes and counted using a Moxi automatic cell counter (VWR) with size parameters adjusted to exclude apoptotic cells. The average final cell counts were (1.107 ± 0.050) × 10⁷ cells per control sample, and (1.558 ± 0.291) × 10⁷, (1.486 ± 0.124) × 10⁷, and (1.613 ± 0.156) × 10⁷ cells per HRAS-, KRAS-, and NRAS-transformed sample, respectively. This corresponded to an average population doubling time of 27.68 h for the control cells and 28.43 h, 29.02 h, and 31.89 h for the HRAS-, KRAS-, and NRAS-transformed cells, respectively. After counting, the cells were pelleted at 1,500 rpm for 5 min at 4 °C with the pellets immediately snap-frozen in liquid nitrogen and stored at –80 °C.

Metabolite extraction

The extraction of hydrophilic metabolites from cell pellets was performed using previously established procedures (Gottschalk et al., 2008). Briefly, cell pellets in a 1.5 ml Eppendorf tubes were resuspended by adding 500 µl of a 2:1 (v∕v) ice-cold solution of methanol (Sigma-Aldrich) and chloroform (Sigma-Aldrich) followed by 3–5 min of vortexing and manual mixing for at least 10 min until a clear solution was obtained. Next, 250 µl of ice-cold chloroform and 250 µl of ice-cold water were each added to the sample, which was then vortexed for 5–7 min to yield a cloudy solution. The sample was sonicated at room temperature for ten minutes followed by centrifugation at 13,000 rpm for 5 min at 4 °C in order to yield three layers. The hydrophilic layer was transferred to a fresh Eppendorf tube followed by bubbling with nitrogen gas (Airgas) to remove any residual methanol. The samples were placed under a high speed vacuum concentrator at room temperature until dried, and the dried hydrophilic layer was stored in a –80 °C freezer until needed.

NMR sample preparation, acquisition, and processing

The dried hydrophilic layer was resuspended in 400 µl of deuterated PBS at pH = 7.6 that was prepared as previously reported (Sambrook, Fritsch & Maniatis, 1989). The pH of each sample was adjusted to 7.6 by the addition of either dilute HCl or NaOH as needed to ensure that each metabolite appeared at the same chemical shift in all samples. In each sample, 0.5 µl of a 0.1 M aqueous solution of DSS (Sigma Aldrich) was added for chemical shift referencing. After vortexing, each sample was transferred into a 5 mm NMR tube.

The ¹H NMR spectra were acquired on a 500 MHz Bruker Avance spectrometer (operating at 500.13 MHz for ¹H observation) equipped with a 5 mm TCI 500S2 H-C/N-D-05 Z cryoprobe head at 298 K. Each sample was tuned and matched, reshimmed, and the 90°pulse length was recalibrated (90°pulse lengths ranged between 10 and 12 µs). The ¹H NMR spectra were acquired using a standard Bruker 1D NOESY pulse program with water presaturation and spoiler gradients applied during the relaxation delay, “noesygppr1d.2” as shown in Fig. 1A. This pulse sequence provides good solvent suppression without rolling baselines (Beckonert et al., 2007; Mckay, 2011). The following experimental parameters were used in all measurements: sweep width of 10.33 kHz, 65 K acquisition points, a 2 s recycle delay during which a 93 Hz water presaturation pulse was applied, τ_mix = 101.2 ms, θ = 7.5°–9°, and 256 scans were acquired for each sample. Half-sine shaped pulsed field gradients of duration 1 s with maximum gradient strengths of G1 = 24 G/cm and G2 = –23.7 G/cm were used in Fig. 1A along with a 200 µs gradient stabilization delay placed after each gradient pulse. After acquisition, all FIDs were imported into the Chenomx NMR Suite Profiler (version 7.6., Chenomx Inc., Edmonton, Canada). The data were Fourier transformed after multiplication by an exponential window function with a line broadening of 0.5 Hz, and the spectra were manually phase corrected and baseline adjusted using a cubic-spine function. From the initial set of ten biological replicates for each cell line, only 8 of the control, 7 of the HRAS, 9 of the KRAS, and all 10 of the NRAS samples provided measureable NMR signal from resonances other than the solvent peak. Therefore, the results presented in this work represent data obtained from those N_S = 8 biological replicates of the control cells, and those N_S = 7, N_S = 9, and N_S = 10 biological replicates of the HRAS-, KRAS-, and NRAS-transformed cells.

The Chenomx NMR Suite Profiler was used to identify metabolites by fitting compound signatures from the provided NMR spectral library. In total, 37 metabolites were identified by NMR. The effective NMR metabolite concentration in each sample, S_metabolite, was calculated using the Chenomx NMR Suite Profiler by determining the heights of the compound signatures that best fit the sample spectra with the effective concentration of the internal DSS standard being set to S_DSS = 0.1248 mM, which was the actual DSS concentration in each sample. The table of identified metabolites and their signals was then exported and saved in an Excel worksheet.

Statistical analysis

The “effective” NMR cellular content for metabolite α (moles/cell) taken from the sth biological replicate of a given cell type, ${\tilde{C}}_{\alpha ,s}^{Celltype}$, was calculated by multiplying S_α by the NMR sample volume (400.5 µl) and by dividing by the number of cells used to make up each NMR sample. ${\tilde{C}}_{\alpha ,s}^{Celltype}$ is related to the actual cellular content for metabolite α, $C_{\alpha ,s}^{Celltype}$, by the relationship (1)${\displaystyle {\tilde{C}}_{\alpha ,s}^{Celltype}={\chi }_s^{Celltype}{f}_{\alpha }{C}_{\alpha ,s}^{Celltype}}$ where ${\chi }_s^{Celltype}$ and f_α are dimensionless proportionality factors. The sample- and cell type-independent factor f_α is taken to depend only upon the experimental NMR acquisition parameters (such as recycle delays, mixing times, magnetic field strength, etc.) and metabolite α’s spin topology and relaxation properties. This contribution can in principle be found by applying the 1D NOESY sequence in Fig. 1A to prepared standards. The sample- and cell type-dependent factor ${\chi }_s^{Celltype}$ in Eq. (1) is due to the overall metabolite extraction efficiency, which can vary from sample to sample and depends quite sensitively on cell handling (Duarte et al., 2009) and the particular metabolic quenching and extraction method employed in the study.

The various ${\tilde{C}}_{\alpha ,s}^{Celltype}$ were used to calculate the “effective” NMR fraction of metabolite α in each sample, $x_{\alpha }^{Celltype}$, as follows: (2)${\displaystyle x_{\alpha }^{Celltype}=\frac{{\tilde{C}}_{\alpha }^{Celltype}}{\sum _{j=1}^{37}{\tilde{C}}_j^{Celltype}}=\frac{f_{\alpha }{C}_{\alpha }^{Celltype}}{\sum _{j=1}^{37}{f}_j{C}_j^{Celltype}}.}$ As defined in Eq. (2), $x_{\alpha }^{Celltype}$ is dimensionless and independent of the number of cells in a given biological replicate that were used to make the sample. More importantly, $x_{\alpha }^{Celltype}$ is independent of the sample-dependent fluctuation factor, ${\chi }_s^{Celltype}$ in Eq. (1). The total intensity normalization in Eq. (2) is analogous to that used in spectral binning analysis commonly employed in NMR metabolomic studies. Furthermore, if the various f_α are identical for each metabolite, i.e., f_α = f for all metabolites, then $x_{\alpha }^{Celltype}$ in Eq. (2) is simply the mole fraction of metabolite α for a given cell type (in general, this is not the case, and f_α ≠ f_β).

An ANOVA test, implemented using the MATLAB function “anova1” available in MATLAB’s Statistics toolbox, was first used to test the hypotheses that $〈x_{\alpha }^{HRAS}〉=〈x_{\alpha }^{KRAS}〉=〈x_{\alpha }^{NRAS}〉=〈x_{\alpha }^{Control}〉$ for each metabolite α, where $〈x^{Celltype}〉$ represents the average value of x for a given cell type. The BY algorithm (Benjamini & Yekutieli, 2001) implemented in MATLAB (Groppe, 2010) with the false discovery rate set to 0.01 was then applied to the p-values from the ANOVA analysis to determine those metabolites where $〈x_{\alpha }〉$ significantly differed (adjusted p-values with p ≤ 0.01) between at least two cell types. For those metabolites identified by the ANOVA test, further post-hoc/multiple comparison testing using the BY algorithm was performed to identify which pair(s) of cell types $〈x_{\alpha }〉$ significantly differed (adjusted p-values with p ≤ 0.01, which are given in Table S1). Finally, a PCA of the various $x_{\alpha }^{Celltype}$ was performed using the “pca” command in the Statistics toolbox in MATLAB, which by default, centers the data before performing the PCA.

Due to the similar average values of the “effective” NMR glutamate content, $〈x_{glutamate}〉$, observed in both the control and RAS-transformed cells (Fig. S1) and the relatively large glutamate signals observed in all cells lines (only the lactate signals were larger on average), an alternative to the total intensity normalization scheme used in Eq. (2) was also investigated whereby the metabolite signals were normalized by the observed glutamate signal in each sample. In this case, the glutamate normalized signal for metabolite α is given by: (3)${\displaystyle {\xi }_{\alpha ,glutamate}^{Celltype}=\frac{{\tilde{C}}_{\alpha }^{Celltype}}{{\tilde{C}}_{glutamate}^{Celltype}}=\frac{x_{\alpha }^{Celltype}}{x_{glutamate}^{Celltype}}=\frac{f_{\alpha }}{f_{glutamate}}\times \frac{C_{\alpha }^{Celltype}}{C_{glutamate}^{Celltype}}}$ ${\xi }_{\alpha ,glutamate}^{Celltype}$ is directly proportional to the ratio of the actualcellular metabolite α to glutamate content, and, like $x_{\alpha }^{Celltype}$ in Eq. (2), it is also independent of the sample-dependent fluctuation factor, ${\chi }_s^{Celltype}$. A PCA of the various ${\xi }_{\alpha ,glutamate}^{Celltype}$ was also performed, where in this case only 36 metabolites were considered in the analysis since ${\xi }_{glutamate,glutamate}^{Celltype}=1$ in each sample by definition (Eq. (3)).

One advantage of using glutamate normalization in Eq. (3) compared with using total metabolite normalization in Eq. (2) is that the ratio of ξ_α,glutamate between different cell types is independent of the f_α factors and depends only on the actual cellular metabolite contents: (4)${\displaystyle \frac{{\xi }_{\alpha ,glutamate}^{Celltype1}}{{\xi }_{\alpha ,glutamate}^{Celltype2}}=\left(\frac{C_{\alpha }^{Celltype1}}{C_{glutamate}^{Celltype1}}\right)/\left(\frac{C_{\alpha }^{Celltype2}}{C_{glutamate}^{Celltype2}}\right).}$

The ratio in Eq. (4) is equivalent to the relative fold change in the ratio of the actual cellular metabolite α to glutamate content between different cell types. In those instances where a significant difference in $〈{\xi }_{metabolite,glutamate}〉$ between at least two cell types was identified by ANOVA and post-hoc/multiple comparison testing, quantitative confidence intervals for the ratio in Eq. (4) were calculated using Fieller’s method for unpaired data (Motulsky, 1995). In this case, the (100 − α)% confidence range for ${\xi }_{metab.,glut}^{Celltype1}∕{\xi }_{metab.,glut.}^{Celltype2}$ in Eq. (4), which is denoted by ${ϵ}_{\left(100-\alpha \right)\text{\%}}\left(\frac{{\xi }_{metab.,glut.}^{Celltype1}}{{\xi }_{metab.,glut.}^{Celltype2}}\right)$, is given by (Motulsky, 1995): (5)${\displaystyle {ϵ}_{\left(100-\alpha \right)\text{\%}}\left(\frac{{\xi }_{metab.,glut.}^{Celltype1}}{{\xi }_{metab.,glut.}^{Celltype2}}\right)=\frac{{\mu }_1}{\left(1-g\right){\mu }_2}\left(1\pm t_{1-\frac{\alpha }{2},N_1+N_2-2}\sqrt{\left(1-g\right){\left(\frac{{\sigma }_1}{{\mu }_1\sqrt{N_1}}\right)}^2+{\left(\frac{{\sigma }_2}{{\mu }_2\sqrt{N_2}}\right)}^2}\right)}$ where N₁ and N₂ are the number of biological replicates of cell types 1 and 2, respectively, μ₁ and σ₁ are the average and standard deviations for ${\xi }_{metabolite,glutamate}^{Celltype1}$, respectively, μ₂ and σ₂ are the average and standard deviations for ${\xi }_{metabolite,glutamate}^{Celltype2}$, respectively, $g={\left(\frac{{{\sigma }_{2}t}_{1-}\frac{\alpha }{2},N_1+N_2-2}{{\mu }_2\sqrt{N_2}}\right)}^2$, and $t_{1-\frac{\alpha }{2},N_1+N_2-2}$ is the $\left(\mathrm{1-}\frac{\mathrm{\alpha }}{\mathrm{2}}\right)$th quantile of the t-distribution with N₁ + N₂ − 2 degrees of freedom. If g ≥ 1, the relative fold change in cellular metabolite to glutamate content cannot be calculated using Eq. (5). The MATLAB files and commands used in the statistical analysis of the metabolomics data is given as File S1.

Generation of isogenic cell lines for the study

The results of immunoblotting total protein lysates from the four cell types against the various RAS isoforms are shown in Fig. 1B. Western blotting with antibodies against HRAS (first lane), KRAS (second lane) and NRAS (third lane) confirmed that the cells expressed the appropriate RAS isoforms. Uniformity of loading was also confirmed by immunoblotting against GAPDH, a housekeeping gene, as shown in the bottom lane of Fig. 1B. Previous characterizations of these cell lines have also confirmed that the introduction of the RAS oncogene confers soft agar colony growth in these cells, which is indicative of oncogenic transformation (Rai et al., 2011; Patel et al., 2015; Giribaldi et al., 2015).

NMR-based characterization and PCA of metabolomic profiles

Representative spectra taken from a single biological replicate of the control and HRAS-, KRAS-, and NRAS-transformed cells are shown in Fig. 1C, where the spectra were normalized so that the DSS resonance at δ = 0 ppm had the same intensity in all spectra for display purposes only. Certain key metabolites are labeled using the codes, 1–21, given in Table 1.

The loadings of PC1 (score of 75.1%) and PC2 (score of 16.3%) from a PCA of $x_{metabolite}^{Celltype}$ are shown in Fig. 2A, where some of the components of both PC1 and PC2 are labeled using the codes given in Table 1. In Fig. 2B, a score plot of PC1 vs. PC2, with the corresponding 99% confidence ellipses (Hoover, 1984) drawn for convenience, shows non-overlapping grouping between the control and the RAS-transformed cells. Similar results were also observed when performing a PCA of ${\xi }_{metabolite,glutamate}^{Celltype}$ as shown in Figs. 2C and 2D.

Non-overlapping groupings at the 99% confidence level between all cell lines were found by plotting $x_{lactate}^{Celltype}$ vs. $x_{phosphocholine}^{Celltype}$ as shown in Fig. 3A, which were mainly due to differences in the phosphocholine levels between cell types (Fig. 4). Likewise, non-overlapping groupings between all cell lines, this time at a slightly lower confidence level of 97.5%, were also found by plotting ${\xi }_{lactate,glutamate}^{Celltype}$ vs. ${\xi }_{phosphocholine,glutamate}^{Celltype}$ as shown in Fig. 3B.

Of the 37 metabolites identified by NMR, an ANOVA analysis indicated that 18 metabolites had a significant (adjusted p-values with p ≤ 0.01) difference in their “effective” NMR metabolite fraction (x_α in Eq. (2)) between at least two of the four cell types. Box plots of $x_{metabolite}^{Celltype}$ for these 18 metabolites are given in Fig. 4. Of these 18 metabolites, post-hoc/multiple comparison testing using the BY algorithm found that $x_{metabolite}^{Celltype}$ for 17 metabolites was significantly different between the control cells and at least one of the RAS-transformed cell types whereas the “effective” NMR cellular content for 6 metabolites significantly differed in at least two of the three RAS-transformed cell lines (adjusted p-values ≤ 0.01, which are given in Table S1). It should be noted that while the ANOVA analysis indicated that $〈x_{valine}〉$ was unequal between at least two of the four cell types, post-hoc/multiple comparison testing could not identify any significant difference in $〈x_{valine}〉$ between cell lines.

Similarly, an ANOVA analysis of the glutamate normalized metabolite content indicated that 16 metabolites had a significant (adjusted p-values with p ≤ 0.01) difference in ${\xi }_{metabolite,glutamate}^{Celltype}$ (Eq. (3)) between at least two of the four cell types, and box plots of ${\xi }_{metabolite,glutamate}^{Celltype}$ for those 16 metabolites are shown in Fig. 5. Post-hoc/multiple comparison testing indicated that ${\xi }_{metabolite,glutamate}^{Celltype}$ significantly differed between the control cells and at least one of the RAS-transformed cell types for 13 metabolites whereas ${\xi }_{metabolite,glutamate}^{Celltype}$ significantly differed between at least two of the three RAS-transformed cell lines for only 6 metabolites (adjusted p −values ≤ 0.01, which are given in Table S2). It should also be noted that while the ANOVA analysis indicated that $〈{\xi }_{CP,glutamate}〉$ and $〈{\xi }_{glutamine,glutamate}〉$ were unequal between at least two of the four cell types, post-hoc/multiple comparison testing could not identify any significant difference in either $〈{\xi }_{CP,glutamate}〉$ or $〈{\xi }_{glutamine,glutamate}〉$ between the cell lines. In those instances where a significant difference in $〈{\xi }_{metabolite,glutamate}〉$ between two cell types was identified by ANOVA and post-hoc/multiple comparison testing, the 99% confidence intervals for the relative fold change in the actual cellular metabolite to glutamate content between those cell types, ${ϵ}_{99\text{\%}}\left(\frac{{\xi }_{metab.,glut.}^{Celltype1}}{{\xi }_{metab.,glut.}^{Celltype2}}\right)$, were calculated using Fieller’s method (Eq. (5)) and are given in Table 2. However, even though significant differences between the RAS-transformed and control cells for $〈{\xi }_{fumarate,glutamate}〉$ and between the HRAS-transformed and control cells for $〈{\xi }_{tyrosine,glutamate}〉$ were observed (Fig. 5 and Table S2), the 99% CIs for the relative fold change in the cellular fumarate to glutamate content and the cellular tyrosine to glutamate content could not be calculated due to the small signals and large scatter observed for both tyrosine and fumarate in the control cells (which gave g > 1 in Eq. (5)).

NMR-based identification of metabolite differences among the transformed and control cell lines reflect RAS-driven physiologic alterations

Due to the Warburg effect (Warburg, 1956), it is well known that oncogenic transformed cells undergo aerobic glycolysis as opposed to oxidative phosphorylation (Dang, 2012). To establish the validity of our data against known metabolic changes, we assessed how differences in the NMR signals from lactate and alanine, two major byproducts of glycolytic metabolism (DeBerardinis et al., 2007), varied among the different cell lines. The lactate NMR signal was the largest NMR signal observed in all cell lines (Fig. 4), and lactate was also the largest component to PC1 in the PCA analyses of both x_metabolite and ξ_{metabolite,glutamate} in Fig. 2. In fact, $x_{lactate}^{HRAS}$ and $x_{lactate}^{KRAS}$ were found to be statistically larger than $x_{lactate}^{Control}$ (Fig. 4 and Table S1). The cellular alanine to glutamate content was found to be elevated between 50%–75% in KRAS- and NRAS-transformed cells relative to the control cells (Table 2), which is consistent with the reported phenotype of increased aerobic glycolysis in oncogenic RAS-transformed cells (Hahn & Weinberg, 2002), although we should point out that our study provides only a steady-state snapshot of the metabolic profile.

The cellular UDP-X (i.e., UDP-glucose, UDP-galactose, and/or UDP-glucourinate), which are important molecules in glucose metabolism and in the formation of cellular polysaccharides (Berg, Tymoczko & Stryer, 2002), to glutamate content was elevated between a factor of 1.69–2.94 in the RAS-transformed cells lines relative to the control cells (Table 2). Likewise, the cellular N-acetylcysteine, a thiolic antioxidant (Oikawa et al., 1999), to glutamate content was also elevated in all RAS-transformed cells relative to control cells (Fig. 5) with statistically significant differences occurring for the KRAS- and NRAS-transformed cells, where $\frac{C_{NAC}}{C_{glutamate}}$ was 4.11 and 5.34 times larger relative to control cells, respectively (Table 2). The elevated levels of N-acetylcysteine in RAS-transformed cells is a significant finding given that RAS-transformed cells are known to exhibit elevated redox protective mechanisms (Young et al., 2004).

The metabolomic signatures of two cellular osmolytes, taurine and myo-inositol, also showed significant differences between the RAS-transformed and control cells. The cellular taurine to glutamate content and the cellular myo-inositol to glutamate content were between 50–66% and 32–45% smaller in all RAS-transformed cells relative to the control cells, respectively (Table 2). As osmolytes regulate the apoptotic cell death pathway (Lang et al., 2005), the functional relevance of the lower values of $\frac{C_{myo-inositol}}{C_{glutamate}}$ and $\frac{C_{taurine}}{C_{glutamate}}$ observed in RAS-transformed cells may be related to their relative resistance to stress-induced programmed cell death.

An unexpected result from our study was the cellular phosphocholine to glutamate levels. Choline metabolism is an important component in lipid biogenesis (Glunde, Bhujwalla & Ronen, 2011). The cellular phosphocholine to glutamate content in our study was between 60–82% smaller in the RAS-transformed cells relative to the control cells, and statistically significant differences among the RAS-transformed cells were also observed (Fig. 5 and Table S2). Similarly, the cellular choline to glutamate content was around 1.78 times larger in the KRAS- and NRAS-transformed cells relative to the control cells (Table 2). Although there are reports indicating phosphocholine levels correlate with elevated malignancy (Aboagye & Bhujwalla, 1999; Ronen et al., 2001), exceptions in the published literature suggest that this conclusion may be specific to the RAS isoform and cell type being studied (Eliyahu, Kreizman & Degani, 2007).

We have demonstrated in this work that ¹H NMR can be used to identify unique metabolomic signatures between BEAS-2B immortalized lung epithelial cells and those transformed with the isoforms of the RAS oncogene as well as among the three RAS isoforms. Collectively, our results suggest that measuring cellular metabolomic profiles can help in distinguishing between normal and RAS-transformed cells along with potentially distinguishing among cancer cells expressing different RAS isoforms. In the future, these results may aid in the development of potential screening technology to determine particular cancer treatment regimens.