Protein interaction and functional data indicate MTHFD2 involvement in RNA processing and translation

Cell culture

HCT-116 cell lines were obtained from the National Cancer Institute. The MTHFD2 CRISPR^−/− knockout (D2-KO) and parental wildtype (WT) HCT-116 cells were obtained from the lab of Dr. Rabinowitz [7]. Cells were grown in RPMI-1640 medium supplemented with 5% FBS (Life technologies) and Penicillin/Streptomycin (Gibco) and maintained at 37 °C and 5% CO2.

For serum starvation experiments, 200,000 cells were plated in six-well plates. The next day, the medium was removed, the cells were washed with PBS and then serum-free RPMI-1640 medium was added. Cyclohexamide was obtained from Sigma Aldrich (C4859). For the serum re-stimulation experiments, after 48 h of serum starvation, wells were washed with PBS and then fresh RPMI-1640 containing 20% FBS was added to the cells.

Real-time PCR

RNA was isolated from cells using the Qiagen RNeasy mini kit and quantified using a Nanodrop ND-1000 spectrophotometer. The RNA was reverse transcribed to cDNA using the VILO cDNA synthesis kit (Invitrogen). The qPCR reactions were performed in triplicate on a StepOne Real-Time PCR machine (Thermo Fisher Scientific) using the Fast SYBR Green master mix (Thermo Fisher Scientific). Oligos were ordered from Sigma Aldrich. Primer sequences are listed in Additional file 1. Relative amounts of mRNA were calculated using the ΔΔC_T method normalized to RPLPO mRNA as a reference.

Immunoblotting

Total protein was isolated from cells using RIPA buffer (Thermo Fisher Scientific), with the addition of × 100 Halt protease inhibitors (Thermo Fisher Scientific) just prior to the extraction. Protein from nuclear and cytosolic cellular compartments were fragmented as described previously [17]. Protein levels were quantified by the BCA assay (Pierce). Equal amounts of protein (10–20 μg) were loaded per lane of 10% SDS-PAGE gel. The following primary antibodies were incubated overnight: anti-MTHFD2 (Proteintech 12270-1-AP); anti-MTHFD1 (Proteintech 6113-1-AP); anti-MTHFD1L (Proteintech 10794-1-AP); Lamin A/C (Thermo Fisher Scientific MA3-1000); HSP60 (Abcam ab46798); Tubulin (Abcam ab6046); and COX IV (Abcam 33985). Quantification of protein bands was performed using the ImageJ software and values were normalized to tubulin for each sample.

Co-immunoprecipitation

WT and D2-KO HCT-116 cells were lysed using the moderate Pierce IP lysis buffer (Thermo Fisher Scientific) with addition of × 100 Halt protease inhibitors (Thermo Fisher Scientific). All steps were performed at 4 °C to reduce protein disassociation. The protein lysate was quantified using the BCA assay, and 1 or 2 mg of lysate was used for immunoprecipitation with no cross-linkage. Immunoprecipitation was performed using two separate anti-MTHFD2 antibodies (Proteintech 12270-1-AP (AP) and Genetex N1C3 (NC)) or an anti-rabbit immunoglobulin (IgG) (Santa Cruz). The immunoprecipitation was performed at 4 °C for 2 h using Dynabeads Protein A (Thermo Fisher Scientific) following the manufacturer’s instructions.

In-gel digestion and mass spectrometry

Immunoprecipitated proteins were eluted from Dynabeads through a combination of elution buffer and heating and separated on a 10% SDS gel. In-gel staining of proteins was performed using the SilverQuest silver staining kit (Thermo Fisher Scientific). Protein lanes were excised manually and in-gel digested using a MassPREP robotic protein-handling system (Waters, Millford, MA, USA). Gel pieces were distained twice with 100 μL of 50 mM ammonium bicarbonate (AmBic) containing 50% acetonitrile at 40 °C for 10 min. Proteins were reduced with 10 mM DTT in 100 mM Ambic for 30 min at 40 °C and alkylated with 55 mM iodoacetamide in 100 mM AmBic at 40 °C for 20 min followed with in-gel digestion with 0.3 μg trypsin (sequence grade, Promega, Madison, WI, USA) in 50 mM at 40 °C AmBic for 5 h. The tryptic peptides were extracted with 1% formic acid in 2% acetonitrile, followed by 50% acetonitrile twice. The liquid was evaporated to dryness, and the peptides were reconstituted in 10 μL of 2% acetonitrile/0.1% formic acid. Five microliters of the samples was injected onto the LC-MS/MS system (Ultimate™ 3000 RSLCnano chromatography system and Q Exactive HF Orbitrap mass spectrometer, Thermo Scientific). The peptides were separated on an Easy-C18 column, 50 cm (Thermo Fisher Scientific) at 55 °C in a 120-min gradient at a flow rate of 300 nL/min. The linear gradient was from 5 to 26% of buffer B (98% acetonitrile/0.1% formic acid) in 115 min and to 95% of buffer B in 5 min.

We used data-dependent acquisition with a survey scan range of 300 to 1650 m/z, at a resolution of 120,000 m/z and selected up to 16 most abundant features with a charge state ≥2 for HCD fragmentation at a normalized collision energy of 26 and a resolution of 30,000 at m/z 200. To limit repeated sequencing, dynamic exclusion of sequenced peptides was set to 90 s. Thresholds for ion injection time and ion target values were set to 250 ms and 5 × 106 for the survey scans, and 120 ms and 2 × 105 for the MS/MS scans. Data were acquired using the Xcalibur software (Thermo Fisher Scientific). The spectra were analyzed using the Mascot search engine v. 2.4 (Matrix Science Ltd., UK) using carbamidomethylation (C) as fixed and deamidation (NQ), oxidation (M) as variable modification. The precursor mass tolerance was set to 10 ppm while MS/MS tolerance was 0.02 Da allowing 2 missed cleavages. The SwissProt database was used with taxonomy limited to Homo sapiens.

In total, we performed 12 Co-IP experiments (for HCT-116 WT, 4 using AP antibody, 2 using NC antibody, and 2 using IgG antibody; for HCT-116 D2-KO, 2 using AP antibody and 2 using NC antibody) followed by MS analysis. A protein was considered as high confidence if at least two unique peptides were detected in 50% of IP-MS samples for each of the MTHFD2 antibodies used and was not detected in any of the IP-MS samples from D2-KO cells, nor the IgG isotype control IP-MS samples.

Protein-protein interaction network

A protein-protein interaction network was extracted from the STRING database v10.5 [21] using the set of MTHFD2 interacting proteins identified here. Only experimental data (e.g., from co-purifications and yeast two hybrids) imported in STRING from primary sources were used for constructing the network. Threshold for connections between proteins was set at medium score (more than 0.400) as calculated in STRING, and gene ontology (GO) enrichment was calculated.

Half-life analysis

For half-life analysis, we used publicly available data on global protein half-life in NIH3T3 cells [22], HeLa, and C2C12 cells [23]. The list of 1652 proteins involved in metabolism were obtained from Recon2 [24] model of human metabolism, of which 596 were measured in HeLa, 575 in C2C12, and 605 in NIH3T3.

Co-expression analysis

Calculation of gene co-expression was performed using a clustering-based method, as described previously [17]. Briefly, we analyzed 25,485 genes represented on a variety of microarray platforms for co-expression with MTHFD2 across 8097 human, mouse and rat data sets obtained from the Gene Expression Omnibus (GEO). From the resulting matrix of coexpression scores x_gd, we computed the overall coexpression score w_g for each gene g by simply summing over all data sets d. We selected the 50 data sets showing highest coexpression for the D2PPI genes by scoring each data set by the weighted sum \( \sum \limits_g{w}_g{x}_{gd} \), where g ranges over the 29 genes in the D2PPI gene set; and similarly for the ATF4 gene set (Fig. 3b). Example data sets from GEO (Fig. 3c–f) were analyzed as deposited on GEO, with no further normalization. In cases where multiple probes against a single gene were present on the arrays, the probe with highest mean signal across all samples was used. Enrichment analysis was done using the GSEA-P statistic [25], and an enrichment p value was calculated by gene permutation (10,000 permutations). An independent analysis was made using the SEEK tool (http://seek.princeton.edu) with default settings [26].

Connectivity Map (CMap) analysis

Expression data from the Connectivity Map (CMap) project in Level 5 (signature) GCTX format was obtained from GEO (accession GSE92742), and data for 99 experiments with five distinct shRNA hairpins against MTHFD2 across 12,328 measured genes was extracted using R v.3.3.3 and the cmapR package (https://github.com/cmap/cmapR). Since the three hairpins TRCN0000036550, TRCN0000036551, and TRCN0000036553 exhibited strongest knockdown on MTHFD2 itself, these were used for subsequent analysis. Data from the ASC and NPC cell lines were discarded due to few replicates. The z-scores for the remaining 9 cell lines were averaged into a final z-score vector, which was used for enrichment analysis using the GSEA-P statistic [25]. The genesets examined were downloaded from ConsensusPathDB collection (http://cpdb.molgen.mpg.de/). An enrichment p value was calculated by gene permutation (10,000 permutations).

Connectivity scores between hairpins targeting three one-carbon enzymes included in CMap (MTHFD2, MTHFD1, and SHMT1) and hairpins targeting 3799 examined genes were analyzed using the CMap online tools (https://clue.io/cmap). Connectivity scores of more than 90 or less than −90 were considered as significantly enriched. Clustering analysis of connectivity scores was performed using ClustVis [27] with average linkage and correlation metric. Ribosomal genes were selected for analysis based on the KEGG ribosome gene set (entry hsa03010).

RNAi and CRISPR screening data analysis

Global RNAi and CRISPR screening data were performed as part of the Achilles project in 501 and 342 cell lines respectively [28, 29]. We also examined a second CRISPR study conducted in 14 AML cell lines [30]. To identify genes with similar dependencies across examined cell lines in the RNAi and CRISPR screens, we calculated Pearson’s correlation for MTHFD2 against all examined genes in both datasets.

Protein interaction partners of MTHFD2

Since a non-metabolic function of MTHFD2 would likely involve direct physical contact with other proteins, we decided to first perform co-immunoprecipitation (Co-IP) of MTHFD2-interacting proteins followed by mass spectrometry (MS) to identify its binding partners. To avoid potential artifacts due to unspecific antibody binding, we chose to perform Co-IP experiments in wildtype (WT) and CRISPR MTHFD2 knockout (D2-KO) HCT-116 cells [7], using two distinct MTHFD2 antibodies, 12270-1AP (AP; Proteintech) and N3C3 (NC; Genetex) as well as an isotype control (IgG) antibody (Fig. 2a). We first confirmed by immunoblotting and IP that the MTHFD2 protein was absent from D2-KO cells (Additional file 2). Next, we performed Co-IP followed by SDS-PAGE and silver-staining to verify the efficient immunoprecipitation of MTHFD2. This analysis revealed a strong band at the expected size of MTHFD2, as well as a large number of other immunoprecipitated protein bands (Fig. 2b). In D2-KO cells, the MTHFD2 band was not detected, but other protein bands were still observable, presumably due to non-specific binding of the antibodies used (Additional file 2). Detected peptide fragments in all examined Co-IP lysates are listed in Additional file 3. Reassuringly, MTHFD2 itself was reliably identified from WT cells using both antibodies and strongly decreased in the D2-KO lysates (Fig. 2c).

Gene Ontology (GO) analysis of these proteins confirmed a highly significant enrichment of RNA binding proteins (24/29 proteins; q < 10⁻²⁰) (Fig. 2d; Additional file 4) and that a large fraction of these proteins were also nuclear (Fig. 2e). To determine if the D2PPI proteins belonged to any known complexes, we extracted a protein-protein interaction (PPI) network over these 29 proteins from the STRING database [21]. We obtained 47 physical interactions, which is two times more than expected by chance (p < 10⁻⁷); hence, the D2PPI proteins are likely part of previously known protein complexes, including the small ribosomal subunit proteins (four RPS proteins), RNA processing proteins (several hnRNP family members; SF3B3; ILF2), protein chaperones (several heat shock proteins; NPM1), and histone/DNA repair proteins (three histone proteins; XRCC6; RPN1) (Fig. 2f). Hence, the proteins shown here to directly interact with MTHFD2 indicate a role for the protein in metabolism and translation of RNA.

Co-expression analysis of MTHFD2 and interacting proteins

To systematically assess if the D2PPI proteins indeed share a common function with MTHFD2, we next asked whether their mRNAs are also frequently co-expressed with MTHFD2. We have previously shown that frequent co-expression across diverse experimental conditions is highly predictive of closely related biological function [32, 33]. Using a previously described method [17], we scored 25,845 genes for co-expression with MTHFD2 in 8067 human, rat, and mouse microarray datasets (Additional file 5). Remarkably, D2PPI gene set was clearly enriched for coexpression with MTHFD2 (p = 0.003, permutation test), with 13 of 29 of genes within the 95th percentile of coexpression scores (Fig. 3a). An exception was the histone proteins H3F3A, HIST1H3A, HIST3H3, and HIST2H3A that did not exhibit co-expression with MTHFD2; however, these proteins were not well-represented on the microarrays used, so this result may reflect low power to detect coexpression. An independent analysis of 3356 human gene expression datasets using the SEEK gene co-expression analysis tool [26] gave similar results (Additional file 5). To investigate factors that might drive the observed coexpression of the D2PPI gene set, we identified specific gene expression data sets where this set was coexpressed (Fig. 3b, Additional file 6). Among the top scoring data sets were several tumor gene expression studies, where expression varied between subtypes; for example, in one leukemia study, highest expression was found in T-ALL and lowest in AML subtypes, suggesting specific disease contexts where the D2PPI proteins may be relevant (Fig. 3c). We also noted a concerted induction of the D2PPI proteins in stimulated T cells (Fig. 3d), in line with previous observations of MTHFD2 induction in this setting [4], and also suppression of D2PPI genes in HCT-116 cancer cells in response to treatment with a CDK inhibitor (Fig. 3e). These data indicate that MTHFD2 and the D2PPI genes are frequently coexpressed and responsive to mitogenic stimuli and anti-proliferative drugs.

As previously shown [17], a set of ATF4-responsive mRNAs involved in amino acid metabolism and aminoacyl-tRNA synthesis, including the mitochondrial folate-coupled enzymes SHMT2 and MTHFD1L, was also strongly co-expressed with MTHFD2. However, we found that these genes co-express with MTHFD2 in different conditions than the MTHFD2-interacting proteins (Fig. 3b, Additional file 6), suggesting that the underlying mechanism is different. For example, the mitochondrial enzymes were induced in contexts of amino acid starvation, possibly reflecting the involvement of this pathway in amino acid synthesis, while the D2PPI genes were not (Fig. 3f). In summary, the D2PPI genes are co-regulated with MTHFD2 in specific biological conditions, providing independent evidence that they share a biological function with MTHFD2, which appears to be distinct from ATF4-driven induction of MTHFD2.

MTHFD2 is a rapidly regulated protein with a short half-life

Given the above finding that MTHFD2 and the D2PPI gene set is acutely regulated by pro- and anti-proliferative stimuli, we examined in more detail the dynamics of MTHFD2 regulation. Generally, metabolic enzymes are long-lived proteins, while proteins involved in RNA metabolism, cell cycle, and signaling have shorter half-lives [22]. In a reanalysis of global protein half-life data [22, 23], we noted that the half-life of the MTHFD2 protein was consistently within 17–22 h in HeLa, C2C12, and NIH3T3 cells, which was markedly lower than the other folate one-carbon enzymes (Fig. 4a) and within the 5–10th percentile of all detected enzymes (as defined by Recon2 [24]; 596–605 enzymes were detected across three cell lines) (Fig. 4b). We confirmed this short half-life of MTHFD2 by treatment with the cyclohexamide inhibitor of protein synthesis, which resulted in an approximate halving of protein levels after 24 h (Fig. 4c). Therefore, MTHFD2 is unusually short-lived for a metabolic enzyme in general, and one-carbon enzymes in particular.

Considering the short half-life of MTHFD2, we next examined the regulation of the MTHFD2 protein compared to other enzymes in the mitochondrial folate pathway. In HCT-116 cells, after 48 h serum starvation, MTHFD2 and SHMT2 mRNA was significantly decreased, while SHMT1, MTHFD1, and MTHFD1L was unaffected (Fig. 4d). At the protein level, MTHFD2 and MTHFD1L decreased noticeably at the 24 h time point, while MTHFD1 did not (Fig. 4e), consistent with the shorter half-life of MTHFD2. In serum restimulated HCT-116 cells, MTHFD2 protein was induced within 3 h (Fig. 4f), consistent with our previous report [17]. Moreover, induction of the protein was observed in both the nuclear and cytoplasmic compartments (Fig. 4h), and also on the mRNA level (Fig. 4g). Our data show therefore that enzymes of the one-carbon pathway differ in their response to serum depletion and stimulation, with MTHFD2 being particularly responsive.

Transcriptional response to MTHFD2 knockdown in cancer cells

To study the effects of loss of MTHFD2, and how these relate to its identified interacting proteins, we analyzed transcriptomics data from MTHFD2 shRNA knockdown cells generated by the Connectivity Map (CMap) project [34]. Interestingly, the D2PPI gene set was clearly decreased (ES -0.64, p < 10⁻⁴) in MTHFD2 knockdown cells (Fig. 5a) and was among the most negatively enriched among examined collection of 3581 gene sets (Fig. 5b). Hence, loss of MTHFD2 affects expression of the D2PPI genes, consistent with a shared function.

The CMap data also allows assessing functional relationship between genes by scoring the similarity of transcriptional response of cancer cells to individual gene knockdown. In this regard, knockdown of MTHFD2 was highly similar (score > 90/100) to four out of the 14 D2PPI genes represented in CMap: HSPA8, HSPA9, HSPD1, and RPS3A. In contrast, with this criterion, other folate-metabolizing enzymes present in the CMap data set (MTHFD1, SHMT2, GLDC, AMT, TYMS, MTR, MTRR, FTCD, and DHFR) were not similar to MTHFD2, with the exception of MTHFD1 (score = 91.0, rank 414/3798), possibly reflecting the fact that loss of MTHFD1 also gives severe growth phenotypes [7]. In particular, SHMT2 was not closely related (score = 18.74 of 100, rank 2790/3798), consistent with the observations that loss of MTHFD2 is often detrimental to cells [4], while SHMT2 is not [2], and suggesting that at least some effects of MTHFD2 suppression are unrelated to mitochondrial folate metabolism.

Finally, we investigated the sets of genes that in CMap were highly similar to MTHFD2. We noted that, across all 3798 genes for which shRNA data was available in the CMap data set, MTHFD2 was highly similar (absolute score > 90) to 510 genes, compared to only 266 for SHMT2 and 169 for MTHFD1, indicating that loss MTHFD2 induces a response that more commonly occurs with gene knockdown. In addition, the set of 510 genes similar to MTHFD2 in this respect were clearly different from those of SHMT2 and MTHFD1 (Fig. 5c), again consistent with MTHFD2 having a distinct function. In particular, knockdown of MTHFD2, but not SHMT2 or MTHFD1, was highly similar to knockdown of a group of ribosomal proteins (Fig. 5d). Therefore, MTHFD2 is unique in that suppression of this protein causes cellular transcriptomic responses that are highly similar to those observed following the targeting of ribosomal proteins.

Comparison of RNAi and CRISPR suppression of MTHFD2

At first, the observation that MTHFD2 knockdown is similar to knockdown of ribosomal proteins, but not the upstream enzyme SHMT2, seemed counter-intuitive. A possible explanation is that shRNA knockdown cells still contain some residual protein, which may be sufficient to retain activity of a metabolic enzyme, but disrupts a more sensitive non-metabolic function such as control of translation. Indeed, there is evidence that a small amount of residual enzymatic activity can be sufficient to maintain necessary metabolic output, as seen for example with SHMT2 and formylmethionyl-tRNA pools [35]. To investigate this hypothesis, we compared the effects of shRNA knockdown to effects of CRISPR knockout, using genome-wide shRNA and CRISPR screening data from the Achilles project, measuring growth phenotypes across 501 and 341 cell lines, respectively. Here, genes with shared functions are expected to have similar dependency profiles across cell lines [30]. Overall, 977 genes were similar to MTHFD2 in the shRNA screening data, and 161 with CRISPR screening (Pearson’s correlation, Bonferroni-corrected p < 0.05). Interestingly, genes involved in folate metabolism (SHMT2, GART, GCSH, FPGS) were similar to MTHFD2 in the CRISPR, but not in the shRNA screens (Fig. 6a, b, Additional file 7). Moreover, in a similar analysis of an independent CRISPR dataset comprising 14 AML cell lines, SHMT2, SLC25A32, and MTHFD1L were strongly correlated with MTHFD2, while cytosolic one-carbon enzymes showed weak correlations (Fig. 6c). However, in these data sets, the D2PPI proteins frequently showed strong growth phenotypes across all cell lines and were not significantly correlated with MTHFD2 (data not shown). This is consistent with the notion that in CRISPR MTHFD2 knockouts, the metabolic disruption dominates the effects on cell proliferation, while for shRNA knockdown, other functions of MTHFD2 may play a role.