Oxidized phospholipids regulate amino acid metabolism through MTHFD2 to facilitate nucleotide release in endothelial cells

Oxidized phospholipids (oxPAPC) induce endothelial dysfunction and atherosclerosis. Here we show that oxPAPC induce a gene network regulating serine-glycine metabolism with the mitochondrial methylenetetrahydrofolate dehydrogenase/cyclohydrolase (MTHFD2) as a causal regulator using integrative network modeling and Bayesian network analysis in human aortic endothelial cells. The cluster is activated in human plaque material and by atherogenic lipoproteins isolated from plasma of patients with coronary artery disease (CAD). Single nucleotide polymorphisms (SNPs) within the MTHFD2-controlled cluster associate with CAD. The MTHFD2-controlled cluster redirects metabolism to glycine synthesis to replenish purine nucleotides. Since endothelial cells secrete purines in response to oxPAPC, the MTHFD2-controlled response maintains endothelial ATP. Accordingly, MTHFD2-dependent glycine synthesis is a prerequisite for angiogenesis. Thus, we propose that endothelial cells undergo MTHFD2-mediated reprogramming toward serine-glycine and mitochondrial one-carbon metabolism to compensate for the loss of ATP in response to oxPAPC during atherosclerosis.

T he endothelium maintains the vascular homeostasis and limits atherosclerosis development 1 . Activation of the endothelium leads to a vicious cycle of reactive oxygen species (ROS) formation, inflammation and recruitment and activation of monocytes. ROS generated by inflammatory cells and the endothelium promote the oxidation of lipids contained in lipoproteins that further promotes endothelial activation 2 . The process of lipid oxidation has been well studied for polyunsaturated fatty acid side chains of phospholipids in membranes and lipoproteins. Oxidized 1-palmitoyl-2-arachidonoyl-sn-glycero-3-phosphocholine (oxPAPC) is a mixture of oxidized phospholipids commonly found in oxidized low density lipoproteins (oxLDL). OxPAPC accumulates in atherosclerotic lesions and at sites of chronic inflammation 3 .
OxPAPC disturbs endothelial cells by activating both pro-and anti-inflammatory pathways 4 . Exposure to oxPAPC alters the endothelial transcriptome 5 and the complex changes are a consequence of various regulatory circuits acting in concert. If analyzed at a single time point, the endothelial responses to oxPAPC are overwhelmingly complex. To address this, we went beyond traditional transcriptome analyses of the endothelial response to oxPAPC and applied a systems level network approach.
Network models have been used for the identification of genes and gene clusters causal to cardiovascular pathology 6 , yet a comprehensive network model for endothelial cells is lacking. Here, we constructed molecular causal network models, based on a Bayesian method 7 , for human aortic endothelial cells (HAEC). Compared to other methods, Bayesian networks have several advantages for modeling complex biological processes. They model a causal predictive component for regulatory relationships as well as reveal key molecular drivers and they provide flexible platforms to incorporate various -omics data as prior knowledge 8 .
Here, we apply these techniques to study endothelial cell metabolic reprogramming in response to oxPAPC. We propose that oxPAPC targets amino acid metabolism governed by MTHFD2 to replenish endothelial purine pools. Based on this, a key role of MTHFD2 in angiogenesis as well as human atherosclerosis and CAD development could be inferred.

Results
Differential connectivity clusters in response to oxPAPC. To determine gene clusters and key drivers that shape the response of endothelial cells to pro-atherogenic lipids, an integrative network approach was applied ( Supplementary Fig. 1). In the first step, expression profiles of HAEC obtained from 147 heart transplant donors were reanalyzed. In the initial experiments by Romanoski et al., HAECs of this cohort were exposed to oxPAPC (40 µg ml −1 ) and vehicle control for 4 h 5 . Rather than focusing on canonical analysis of differentially expressed genes or co-expression modules we went further and identified differentially connected gene pairs 9 under the two conditions. This differential co-expression analysis is more sensitive to identify disease-or treatment-induced deregulation among interacting genes 9 . In total, 26,759 differentially connected gene pairs were significantly differentially coexpressed, among them 50.4% showed gain of connectivity (GOC), meaning enhanced co-regulation between genes with oxPAPC treatment. In all, 49.5% showed loss of connectivity (LOS). Clustered differentially connected gene pairs yielded nine significant GOS clusters (Fig. 1a) and 11 significant LOS clusters ( Fig. 1b) with co-regulations elicited by oxPAPC. Gain of connectivity cluster 6 showed the most coherent differential connectivity changes compared to all other clusters. Thus, cluster 6 contains and reflects the strongest endothelial responses to oxPAPC. To identify differentially regulated biological processes, we compared each cluster with canonical pathways for enrichment ( Fig. 1c, d, Supplementary Tables 1, 2). Gain of connectivity cluster 6 was significantly enriched for genes involved in amino acid-related biological processes (Fisher's exact test (FET) p value = 1.63E-09) (Supplementary Table 3), hereafter termed amino acid cluster. The applied differential connectivity clustering approach suggests that amino acid metabolism experiences a massive remodeling in response to oxidized phospholipids.
MTHFD2 is a key driver for the HAEC Bayesian network. To explore the mechanisms underlying deregulated amino acid metabolism in oxPAPC-treated endothelial cells, we constructed molecular Bayesian networks. This approach allows for the identification of corresponding causal regulators. For network construction, we integrated the above used expression profiles and previously published expression quantitative trait loci 10 of the same cohort. This allowed us to construct two Bayesian networks for HAEC: one for the control situation without oxPAPC (Fig. 2a, Supplementary Data 1) and one for oxPAPC exposure (Fig. 2b, Supplementary Data 2). Assessment of the control and oxPAPC Bayesian network against widely used databases showed that the estimated accuracy of these Bayesian networks is significantly higher than that of random networks ( Supplementary Fig. 2a). Also, the proposed networks showed significantly better prediction of gene−gene interactions in endothelial cells as compared to global gene networks from the Human Protein Reference Database ( Supplementary Fig. 2b, c). To identify causal regulators and the biological processes targeted by oxPAPC, we identified key drivers 7,11 in the two networks and analyzed their associated subnetworks for enriched canonical pathways (Fig. 2c, d). The identified top ranked key drivers and their associated subnetworks (>100 nodes in subnetwork) are listed in the supplement (Supplementary Tables 4,5). MTHFD2 was identified as a key driver specific for the oxPAPC Bayesian network (Fig. 2d). The subnetwork regulated by MTHFD2, termed MTHFD2 network, was significantly enriched for genes involved in amino acid metabolism (FET p value = 7.32E-08) (Supplementary Table 6) and significantly overlapped with the identified amino acid cluster (FET p value = 2.33E-28). Thus, MTHFD2 may facilitate a shift in amino acid metabolism in HAEC in response to oxidized phospholipids.
MTHFD2-associated amino acid metabolic gene cluster. The endothelial MTHFD2 network comprised 114 genes and was significantly enriched for cytosolic tRNA-aminoacylation (FET p value = 2.53E-06), membrane transport proteins of the solute carrier family (FET p value = 1.82E-04), the majority of which are involved in the import of amino acids, and transcriptional regulation in response to stress (FET p value = 2.17E-05) (Fig. 3a). Importantly, the MTHFD2 network was also enriched for genes involved in glycine-serine metabolism (FET p value = 2.53E-06). The genes within the overlap of the MTHFD2 network and this gene set category control mitochondrial folate cycle, glycine and serine metabolism and contribute to aspartate, asparagine, methionine, and cysteine interconversion (Fig. 3b).
MTHFD2 is highly induced in many cancers and is essential for mitochondrial one-carbon metabolism-the highest scoring metabolic pathway across human cancers 12 . MTHFD2 is a mitochondrial enzyme with 5,10-methylene-tetrahydrofolate (5,10-meTHF) dehydrogenase and cyclohydrolase activity. The mitochondrial folate cycle generates glycine from serine, which is catalyzed by the enzyme mitochondrial serine hydroxymethyltransferase 2 (SHMT2) (Fig. 3b). The serine consumed for this reaction is synthesized by the cytosolic phosphoglycerate dehydrogenase (PHGDH) and the phosphoserine aminotransferase 1 (PSAT1) from the glycolysis intermediate 3-phosphoglycerate 13 . During the reaction of SHMT2 the co-factor tetrahydrofolate (THF) is converted to 5,10-meTHF. The enzyme MTHFD2 uses 5,10-meTHF to produce the highly reactive 1-C donor 10-formylTHF (10-fTHF) which can be recycled to THF by MTHFD1L. Both 10-fTHF and glycine are indispensable for the de novo synthesis of purines 14 . Interestingly, the function of MTHFD2 in differentiated cells has been widely overlooked and its role in endothelial biology or in lipid signaling remains, to our knowledge, unknown. Given the important role of metabolism for endothelial cell reprogramming 15 , we further analyzed this interesting subnetwork.
MTHFD2 mediated oxPAPC-dependent changes. To better understand and validate the MTHFD2 network, we assessed the role of MTHFD2 in endothelial cells. MTHFD2 expression was decreased by siRNAs in HAEC and gene expression changes were determined with and without oxPAPC. To determine the hierarchical position of MTHFD2 within the subnetwork, we additionally silenced PSAT1, one of its direct downstream nodes (Fig. 3a). Expression of MTHFD2 was induced by oxPAPC as well as by knockdown of PSAT1 (Fig. 3c). Importantly, also the expression of the key enzymes of the serine-glycine synthesis pathway (SHMT2, PSAT1, and PHGDH) was increased by oxPAPC, and this effect was potentiated by knockdown of MTHFD2 (Fig. 3d-f). OxPAPC and MTHFD2 knockdown, both and in an additive manner, increased the expression of the transcription factor CEBPB (Fig. 3g), the aminoacyl-tRNA synthetases for glycine (GARS) and cysteine (CARS) (Fig. 3h, i) and the solute carrier transporters SLC7A5 and SLC7A1 (Fig. 3j, k). Expression of genes in the MTHFD2 network was similarly but to a lesser extent affected by PSAT1 knockdown. Importantly, MTHFD1L, which contributes to mitochondrial one-carbon cycle, but does not belong to the MTHFD2 network, was unaffected by oxPAPC and unaltered by siRNA against PSAT1 or MTHFD2 (Fig. 3l). These results validate MTHFD2 as a key regulator of the MTHFD2 network in oxPAPC-exposed HAEC. The data indicate that oxPAPC elicits a metabolic shift promoting the cellular uptake of amino acids and the de novo synthesis of glycine.
MTHFD2 regulates 18 genes related to amino acid metabolism. To address the importance of MTHFD2 for endothelial gene expression, we performed RNAseq of HAEC with and without siRNA against MTHFD2 (Fig. 4a). Gene set enrichment analysis of the MTHFD2 RNAseq signature confirmed a strong enrichment for amino acid metabolism (FET p value = 7.54E-14) (Supplementary Data 3). We then projected the MTHFD2 RNAseq signature onto the MTHFD2 network. This demonstrated that MTHFD2′s close neighbors were upregulated by MTHFD2 siRNA (Fig. 4b). Furthermore, the MTHFD2 RNAseq b a Gain of connectivity Cluster 1 2 3  4  5  6  7  8  9 Cluster 1 2 3 4 5 6 7 8 9 10 11 signature overlapped with the MTHFD2 network as compared to the whole oxPAPC Bayesian network (FET p value = 8.34E- 16). Finally, the MTHFD2 RNAseq signature also significantly overlapped with the amino acid cluster as compared to all other differential connectivity clusters (FET p value = 1.19E-17).
To identify indispensable pathway genes, an integrative analysis of all three data sets, the amino acid cluster, the MTHFD2 network and the MTHFD2 RNAseq signature, was carried out. This yielded 18 commonly shared genes, which included key enzymes of the glycine-serine-synthesis pathway, aminoacyl-tRNA synthetases, solute carrier transporters as well as the transcription factor CEBPB and the mitochondrial phosphoenolpyruvate carboxykinase 2 (PCK2) (Fig. 4c). Collectively, these data further support the role of MTHFD2 in amino acid metabolic reprogramming.
OxPAPC and MTHFD2 siRNA deplete intracellular glycine pools. We next identified the physiological function of this MTHFD2 response. Knockdown of MTHFD2 decreased intracellular glycine and increased the intracellular serine concentrations ( Fig. 4d), suggesting that the de novo mitochondrial glycine synthesis is highly active in endothelial cells and more important than cellular uptake or cytosolic glycine synthesis. Importantly, oxPAPC decreased the intracellular glycine level, even though it induced the expression of MTHFD2 (Fig. 4d).
OxPAPC elicit endothelial purine nucleotide release. Why do endothelial glycine levels decrease upon oxPAPC stimulation? Glycine is an important substrate for numerous pathways including heme synthesis, glutathione synthesis and for protein de novo synthesis but particularly important for purine synthesis. Purines, in the form of ATP and GTP, are essential for the cellular energy homeostasis and for DNA and RNA de novo synthesis. Purines are also important signaling transmitters and endothelial cells release purines as signaling autacoids in response to numerous stimuli including shear stress 19 . Through the activation of purinergic receptors, purines impact on vascular homeostasis, coagulation, inflammation and on the control of vascular tone 20 . Thus, we speculated that oxPAPC induces the release of ATP and other nucleotides and thereby depletes endothelial cells of these molecules. Indeed, exposure to oxPAPC reduced the endothelial intracellular purine but not pyrimidine nucleoside pools (Fig. 5a-d).
Moreover, the extracellular degradation products of ATP and GTP, but not of pyrimidines, were increased in response to oxPAPC exposure in the cell culture supernatant of HAEC (Fig. 5e-h). MTHFD2 contributes to purine synthesis in the form of glycine and 10-fTHF: Two carbon units are derived from serine which is converted into glycine and incorporated into the purine backbone (Fig. 5i). Two additional single carbon units are derived from two serine molecules and are incorporated into the purine backbone as 10-fTHF. Additionally, the mitochondrial glycine cleavage system can cleave glycine into CO 2 and 5,10-meTHF which both can be incorporated into the purine backbone. We therefore analyzed the flux of serine-and glycine-derived carbons in response to oxPAPC by heavy isotope tracing. With and without oxPAPC, labeled serine was detected intracellularly in serine labeled cells and labeled glycine was detected to a lesser extent intracellularly in glycine labeled cells, indicating that HAEC take up serine as well as glycine ( Supplementary Fig. 5a, b). Extracellular AMP, a degradation product of ATP, contained roughly 90% heavy carbons in heavy serine labeled as well as in heavy glycine labeled cells (Fig. 5j, k). This indicates that HAEC incorporated imported serine as well as imported glycine into the purine backbone of adenosine nucleotide derivatives which are then released from HAEC. OxPAPC-treated HAEC showed comparatively more heavy labeled extracellular AMP which was decreased after silencing of MTHFD2 (Fig. 5l). In contrast, the purine-derivatives NAD and inosine showed only about 20% incorporation of heavy serine-derived carbon units ( Supplementary Fig. 5c, d) and NAD showed only about 10% incorporation of heavy glycine derived carbons ( Supplementary  Fig. 5e). The incorporation of heavy carbons was decreased upon oxPAPC treatment. This suggests that serine-glycine metabolism in HAEC contributes to the synthesis of adenosine nucleotide derivatives for release which is enhanced by oxPAPC. We therefore hypothesized that oxPAPC activates HAEC to synthesize purine nucleotides. Conforming to this hypothesis, oxPAPC exposure increased the metabolic activity of HAEC as scored by the mitochondrial oxygen consumption rate (OCR) (Supplementary Fig. 5f).
Blockade of purine release prevents MTFHD2 induction. If purine release was indeed the mechanism driving the cellular responses to oxPAPC, blockade of ATP release should attenuate the effects of oxPAPC. Flufenamic acid (FFA) blocks ATP release 21 and inhibits nucleotide transporter (VNUT)-mediated ATP transport 22 . Indeed, FFA not only blocked the oxPAPCmediated increase in extracellular purines (Fig. 5m), but also prevented the induction of MTHFD2 and PHGDH in response to oxPAPC (Fig. 5n, o). Similarly, the VNUT inhibitors αglycyrrhetinic acid and A438079 prevented oxPAPC-mediated ATP release (Supplementary Fig. 5g). Moreover, all inhibitors partially restored sprouting after oxPAPC treatment of human umbilical vein endothelial cells (HUVEC) (Fig. 5p, q, Supplementary Fig. 5h−j). In contrast, glycine alone could not rescue oxPAPC-inhibited sprouting ( Supplementary Fig. 5k-m).
Amino acid deprivation induces the MTHFD2 network. Our observations raise the possibility that the MTHFD2 network induction in response to oxPAPC is a direct consequence of glycine depletion. To address this, we deprived HAECs of glycine and serine which induced the expression of MTHFD2 and PHGDH (Fig. 6a, b). Additionally, we elicited amino acid deprivation by asparaginase which depletes cells of asparagine 23 and determined the effect of the histidine analog L-histidinol which inhibits activation of histidine by histidyl-tRNA synthetase 24 . Indeed, both treatments induced the expression of genes in the MTHFD2 network ( Supplementary Fig. 6a, b). These findings suggest that the MTHFD2 network constitutes an amino acid response that compensates for oxPAPC-induced loss of glycine.
MTHFD2-knockdown-impaired angiogenesis is glycine dependent. In line with this conclusion, we hypothesized that supplementation of glycine should overcome the effects of MTHFD2 siRNA. Indeed, glycine but not serine or asparagine, fully prevented the expression of the tested genes belonging to the MTHFD2 network as well as ATF4 in HUVEC (Fig. 6c-f). To determine the importance of this response we measured endothelial migration in a scratch wound assay and the angiogenic capacity as determined by spheroid outgrowth in response to VEGFA. MTHFD2 siRNA impaired migration of HUVECs which was prevented by glycine supplementation (Fig. 6g, Supplementary Fig. 6c). Interestingly, also folate but not formate restored normal migration upon knockdown of MTHFD2 ( Supplementary  Fig. 6d, e). Impairment of sprouting by silencing MTHFD2 was prevented by glycine, but not by folate or formate supplementation (Fig. 6h, i, Supplementary Fig. 6f−i). Importantly, proliferation was only slightly decreased upon knockdown of MTHFD2 and this effect was rescued by glycine whereas oxPAPC inhibited proliferation in an MTHFD2-independent manner ( Supplementary Fig. 6j). To confirm the importance of MTHFD2 for angiogenesis, we tested the contribution of MTHFD2 in an ex vivo organoid mouse model of angiogenesis. Murine aortic rings with siRNA-mediated knockdown of Mthfd2 showed decreased vessel formation which was restored by glycine     (Fig. 6j, k). Furthermore, zebrafish embryos exposed to oxPAPC showed deregulated vessel formation and increased mRNA expression of mthfd2 (Fig. 6l-n). Additionally, zebrafish embryos with a morpholino-mediated knockdown of mthfd2 showed deregulated vasculature morphology ( Fig. 6o-r,  Supplementary Fig. 6k). Collectively, our findings indicate that oxPAPC induced endothelial glycine depletion; that MTHFD2 is induced to compensate for glycine deprivation; and that the MTHFD2-dependent glycine production is essential for endothelial cell function. We suggest that oxidized phospholipids induce a metabolic reprogramming response of endothelial cells which allows for the replenishment of cellular nucleotide pools that become depleted as a consequence of oxPAPC-induced nucleotide release (Fig. 6s).
MTHFD2 impacts on plasma metabolites in humans. To support the role of the MTHFD2 network and serine-glycine metabolism, we speculated that alterations of the MTHFD2 network might impact on the plasma concentrations of relevant metabolites in humans. It should be mentioned that this approach would not allow us to link the findings to any specific vascular alterations. We obtained data from a human genome-wide association study that provides the plasma abundance of 400 metabolites 25 allowing us to test whether genes in the MTHFD2 network were associated with the metabolite abundances. For this, we searched for SNPs (meta-analysis p value < 1×10 −4 ) significantly associated with human plasma metabolites. An SNP rs10174907 in MTHFD2 was significantly associated with plasma N-acetylglycine concentration, SNPs associated with N-acetylglycine concentration were enriched for the genes in the MTHFD2 network (FET p value < 0.0045) and genes of the MTHFD2 network showed associations with plasma metabolites directly or indirectly related to glycine metabolism (Table 1). In total, 60 SNPs in 28 genes in the MTHFD2 network were associated with plasma metabolite concentrations (Supplementary Data 4). This analysis suggests that the MTHFD2 network, as constructed by us for HAECs, is essential in the regulation of plasma metabolite concentrations. To address a link to vascular disease, the plasma glycine to serine ratio in a human cohort with carotid artery disease was determined. Interestingly, the plasma glycine to serine ratio was decreased in subjects with unstable atherosclerotic plaques (n = 26) and decreased by tendency in subjects with stable atherosclerotic plaques (n = 26) compared to plaque-free subjects (n = 26) (Fig. 7a).

MTHFD2 expression is increased in cardiovascular diseases.
Since atherosclerosis is driven in part by oxidized lipids, we hypothesized that the MTHFD2 network is active in human atherosclerotic samples. To identify whether the MTFHD2 network has any association with cardio-vascular disease we obtained the data set of the CARDIoGRAMplusC4D genomewide association (GWAS) study 26 . For our analysis, we searched for SNPs (p value from logistic regression <1×10 −4 ) significantly associated with CAD and identified that genes of the MTHFD2 network contained several CAD-associated SNPs ( Table 2, Supplementary Data 5). Moreover, the CAD risk loci were enriched for genes of the MTHFD2 network (FET p value < 0.0077). Collectively, the data raise the possibility that the MTHFD2 network contributes to CAD development.
To address whether the expression of the MTHFD2 network is changed in human atherosclerosis, we reanalyzed data from a study comparing carotid artery gene expression between 32 healthy and 32 plaque-carrying specimens (GSE43292) 27 . Expression changes of the MTHFD2 network genes in the carotid arteries in response to atherosclerosis were similar to those of endothelial cells in response to oxPAPC (Supplementary Fig. 7a). The majority of the MTHFD2 network genes that were induced in response to oxPAPC were also increased in atherosclerosis. Moreover, an association was also observed for oxPAPC-and atherosclerosis-mediated gene repression. In line with this observation, expression of MHTFD2 in 126 human carotid plaque samples of the BiKE study (GSE21545) 28 was positively correlated with genes of the MTHFD2 network as well as Nrf2 and ATF4 but not with MTHFD1L which does not belong to the MTHFD2 network (Fig. 7b, Supplementary Fig. 7b).
To correlate our findings with late stage atherosclerosis, we tested carotid endo-atherectomy specimens from the same human cohort as in Fig. 7a (patient characteristics, Supplementary Table 7). The tissue samples were classified as stable plaques (n = 20) or unstable plaques (n = 20). SHMT2 expression was increased in plaque material of subjects with unstable plaques and MTHFD2 expression was increased in plaque material of both subpopulations (stable and unstable plaques) compared to healthy arteries (Fig. 7c, d). Conforming to the mRNA expression pattern, MTHFD2 expression was increased in vessels of subjects with atherosclerotic disease as compared to subjects without plaque (Fig. 7e, f). Subjects with unstable plaques exhibited a significantly decreased plasma glycine to serine ratio and had a stronger increase in MTHFD2 expression than subjects with stable plaques. Unstable plaques are inflammatory activated but are also enriched in endothelial cells compared to stable plaque and thus are characterized by increased angiogenesis. Taken together, these data support the notion that MTHFD2 is expressed and upregulated in human atherosclerosis and seems to affect the amino acid glycine and serine metabolism.
A problem of these human studies, however, is that the analyses focused on a rather late time point-established vascular Direction of expression of nodes of the MTHFD2 network in RNAseq signature is indicated by the node color. Node size reflects out-degree. c Venn diagram of genes in the amino acid cluster, MTHFD2 network and MTHFD2 RNAseq signature (FDR < 0.01, Benjamini−Hochberg). Genes belonging to all three gene sets are listed. d Heatmap of amino acid profile. HAEC were treated with three different siRNAs against MTHFD2 or scrambled control and exposed to medium (1% FCS) with or without oxPAPC for 4 h. Amino acids in cell lysates were measured by mass spectrometry (n = 6-9). e Relative mRNA expression of MTHFD2 in HAEC pretreated with N-acetylcysteine (NAC) (5 mM) and glycine (500 µM) for 1 h and then exposed to medium (1% FCS) with (oxP) or without (Ct) oxPAPC or tunicamycin (10 µg ml −1 ) for 4 h (n = 5). f qRT-PCR detection of MTHFD2 in HAEC exposed to medium (1% FCS) with or without oxPAPC and treated with tBHQ (20 µM) or DMSO as control for 24 h (n = 4). g HAEC were exposed to medium (1% FCS) with or without oxPAPC and additionally treated with torin (100 nM), ML385 (10 µM) or DMSO as control (n = 4). h, i Western blot detection (h) and quantification (i) of phosphorylated S6, S6 and MTHFD2 in HAEC pretreated with rapamycin (20 nM) or DMSO as control overnight and exposed to medium (1% FCS) with or without oxPAPC for 4 h (n = 4). j HAEC were treated with scrambled control siRNA (siCtr) or siRNA against ATF4 and exposed to oxPAPC or control medium for 4 h (n = 5). Data are represented as mean ± SEM, *p ≤ 0.05 (oxP or tunicamycin vs Ct), # p ≤ 0.05 (inhibitor present vs absent or ATF4 vs Control siRNA) (ANOVA with Bonferroni post-hoc test) disease in a heterocellular context (whole organ) whereas our initial experiments were carried out in endothelial cells at an early time point. To bridge this gap, we returned to mouse models. Exposure of the murine thoracic aorta to oxPAPC in organ culture increased the expression of the network genes Mthfd2, Phgdh, Shmt2 and Slc3a2, an effect prevented by rapamycin (Fig. 7g-j). Moreover, the MTHFD2 network was induced in endothelial cells harvested from carotid arteries in the partial ligation model, an in vivo mouse model of blood flow withdrawalinduced accelerated atherosclerosis 29 . Mthfd2 and Shmt2 were significantly induced in the endothelium of partially ligated arteries compared to non-ligated arteries (Fig. 7k, l). Importantly, endothelial expression of genes within the MTHFD2 network was already induced 24 and 48 h after partial ligation in a previously published data set 29 (Supplementary Fig. 7c).
These findings may suggest that the MTHFD2 network induction is fundamental to the atherosclerotic disease at large and not only observed upon direct stimulation with oxPAPC. To test this hypothesis, we subjected ApoE−/− mice to a brief course of an atherogenic high-fat Paigen diet. Paigen diet is a specific high fat diet that contains not only triacylglycerides (15.5%) and cholesterol (1.25%) but also cholate (0.5%). When administered, high plasma lipid and cholesterol levels, as well as systemic inflammation, are induced and massively accelerate atherosclerosis development. ApoE−/− mice are particularly sensitive to this diet as these mice already exhibit high VLDL levels under basal conditions 30 . Importantly, already 4 days after the initiation of the diet, the expression of Mthfd2 was increased in the endothelium of the murine carotid artery (Fig. 7m). Encouraged by this finding, we wondered whether also atherogenic lipoproteins obtained from human subjects with vascular disease would induce the MTHFD2 network in endothelial cells. To study this, we obtained high-density lipoproteins (HDL) from healthy subjects (n = 10) or stable CAD patients (n = 10) (patient  Table 8). It was previously established that HDL from CAD patients loses its protective effects and becomes atherogenic 31 . Indeed, HDL derived from CAD patients compared to HDL from healthy subjects induced the network genes MTHFD2, PHGDH, CEBPB, and PCK2 in HAEC (Fig. 7n-q). This effect was sufficiently strong to translate into protein expression: As determined by western blot analysis, MTHFD2 protein abundance was increased in HAEC exposed to CAD patient-derived HDL (Fig. 7r, s).

Discussion
Using an integrative network modeling approach, we identified an unexpected association of atherogenic lipids and atherosclerosis with the mitochondrial one-carbon metabolism. In human endothelial cells, oxidized phospholipids activated a network in which MTHFD2 was identified as a key driver to increase the generation of glycine. Glycine was required for the formation of purines to compensate for their loss from endothelial cells in response to oxidized phospholipids. Endothelial metabolism has been recognized as an important determinant of endothelial cell function 32 , controlling angiogenesis and the adaptation of endothelial cells to different environmental conditions 33 . Oxidized phospholipids, however, have not yet been linked to endothelial metabolism. Deregulated amino acid metabolism and in particular overexpression of serineglycine-synthesizing enzymes and mitochondrial one-carbon metabolism is associated with rapid proliferation of cancer cells 34 . Thus, MTHFD2 is highly expressed in many cancers or the embryo 35 , whereas expression is low in most postmitotic cells 12,36 . We show that MTHFD2 contributes to glycine synthesis in endothelial cells in a disease-relevant manner. Plasma glycine is inversely associated with myocardial infarction risk 37 and in our present work we found an association between plasma amino acids, the MTHFD2 genotype and cardiovascular risk. Moreover, knockdown of MTHFD2 reduced the angiogenic function of cultured cells. Importantly, glycine prevented this effect and also blocked MTHFD2 network induction in response to MTHFD2 siRNA. Currently, MTHFD2 is tested as a potential drug target for cancer therapy 38 . Our data suggest that MTHFD2 inhibition might also attenuate tumor angiogenesis.
MTHFD2 contributes to NADPH production that is required for nucleic acid and lipid synthesis and to glutathione and thioredoxin recycling 39 . Glycine is also needed for glutathione synthesis. In line with this, the MTHFD2 network comprised genes responsive to amino acid starvation, ER stress, mTORC1, and Nrf2 activation. This response pattern reflects the central position of glycine in cellular metabolism.
Purines are products of glycine in the intermediate metabolism and induction of the MTHFD2 network facilitated the endothelial production of purines in response to oxPAPC. In fact, purine nucleotides are important autacoids which, if released from endothelial cells, elicit a plethora of responses. The local release of purine nucleotides from atherosclerotic plaques and the subsequent activation of cellular purinergic receptors is thought to promote the development of atherosclerosis 40 . Nucleotide release is an endothelial stress reaction and occurs in response to thrombin, shear stress, and oxidative stress induced by copperoxidized LDL 19,41,42 . The mechanism of stress-induced purine nucleotide release is not fully understood but a role of gap junctions has been suggested 41,43 . Indeed, gap junction blockers prevented ATP release in response to oxPAPC in the present study. The compounds also prevented the induction of the MTHFD2 network in response to oxPAPC suggesting that the induction of the network is a compensatory reaction to the loss of purines.
OxPAPC and knockdown of MTHFD2 induced ATF4 and several ER stress response genes like DDIT3 and sXBP1 which were also members of the MTHFD2 network. This observation is in line with the fact that ER stress is present in atherosclerotic plaques and promotes the atherosclerotic process through apoptosis 44 . Important pro-atherosclerotic lipids like 7-ketocholesterol through this mechanism fuel the atherosclerotic process 45 . ER stress is also present in regions of the arterial endothelium susceptible to atherosclerosis 46 . This and the inflammatory activation of the endothelium in the atherosclerotic process explain why we observed an activation of the MTHFD2 network in atherosclerotic plaques. Moreover, our finding that decreasing MTHFD2 expression induces endothelial cell dysfunction could serve as an explanation of why genetic variations in genes of the MTHFD2 network were associated with CAD risk.
Oxidation of lipoproteins which has been best studied for LDL occurs in the intima of the arterial wall and thus oxLDL is scarcely found in the plasma 47 . Nevertheless, we observed that the highly inflammatory and atherogenic lipid-and cholate-rich Paigen diet in ApoE−/− mice as well as HDL particles obtained from the plasma of patients with CAD increased endothelial MTHFD2 expression. Although these data link MTHFD2 to the development of atherosclerosis, they must not raise the impression that this response is mediated through oxPAPC. As we observed in the present study, MTHFD2 is a highly stresssensitive gene and induced by numerous pathways. In fact, the dysfunctional HDL from CAD patients induces endothelial oxidative stress 31 and also any high cholesterol diet increases endothelial ROS formation 48 . Thus, although the present work on oxPAPC-mediated induction of MTHFD2 was instrumental in linking signaling of oxidized lipids to amino acid metabolism in endothelial cells, the findings extend beyond this pathway: Endothelial cell activation in general appears to induce MTHFD2 which is a meaningful response given that ATP release is a stereotypic response of endothelial cells to stress.
In conclusion, by taking advantage of an integrative network modeling approach we unveiled a new reaction pattern of endothelial cells. Our work not only demonstrates the power of Bayesian network analysis to uncover novel signaling mechanisms in vascular disease, it also illustrates the enormous secretory capacity of endothelial nucleotide release in the signaling context. It appears that endothelial cells prioritize autacoid release to such an extent that they would rather become deprived of vital metabolic components than neglect their signaling function.
Gene set enrichment analysis. Genes of interest (e.g. genes within each subnetwork or differential connectivity cluster compared to whole network or all differential connectivity clusters) compared against a collection of canonical and noncanonical gene sets (MSigDB v6.0) 49 by using Fisher's exact test in R (p values correspond to Fisher's exact test). Overlaps between the MTHFD2 RNAseq signature, the amino acid cluster, and the MTHFD2 network were computed using Fisher's exact test in R.
RNA sequencing and data analysis. RNA isolated from HAEC was treated with DNase (Qiagen, Cat# 79254). Library construction (LncRNA library, Ribo-Zero), quality assessment, and sequencing (HiSeqSE50) were performed by Novogene. Differentially expressed genes were identified using the following procedure: firstly, sequencing reads were aligned to the human reference genome hg19 using TopHat 50 . Next, the mapped sequences were aligned with htseq-count to quantify the read count for each gene 51 . The edgeR package 52 was then used to identify differentially expressed genes between control and siRNA treatment. We detected 399 differentially expressed genes at FDR of 0.05 (Benjamini−Hochberg).
Differential connectivity clusters. Differentially co-expressed gene pairs and differential connectivity (DC) clusters were computed as the following: In short, Spearman correlation coefficients of each gene pair were computed and transformed into Fisher's Z-statistics 9 . Next, Q statistics were computed and 997 permutations were conducted by randomly assigning sample labels. A final cutoff of Q 0 = 22.5 that corresponds to FDR = 4.89% was used to detect DC pairs. In addition a given gene pair had to be significantly co-expressed (Spearman's correlation p value< 0.01) in either control or oxPAPC group to be called a differentially co-expressed DC gene pair. If the differentially co-expressed gene pair was significantly co-expressed in control group but not oxPAPC group it was assigned to the LOC category and if the DC gene pair was significantly co-expressed in oxPAPC but not control group it was assigned to the GOC category.
HAEC Bayesian networks. The accuracy of Bayesian networks is highly dependent on the number of samples used in the network construction process. We showed in a series of simulation studies that accuracy of reconstructed networks significantly improves when sample size increases from 100 to 1000 53 . In the same paper, we showed that integration of genetic information as structure prior can significantly improve accuracy of constructed Bayesian networks when sample size is between 100 and 300. Experimentally, we systematically validated networks constructed with genetic priors and inferred key drivers in yeast 7 , mouse 54 and human 55 with sample sizes in the above range. In this study, we constructed HAEC Bayesian networks based on 147 samples with genetic data integrated as priors.
Bayesian networks are directed acyclic graphs in which the edges of the graph are defined by conditional probabilities that characterize the distribution of states of each node given the state of its parents. The network topology defines a partitioned joint probability distribution over all nodes in a network, such that the probability distribution of states of a node depends only on the states of its parent nodes: formally, a joint probability distribution pðXÞ on a set of nodes X can be decomposed as pðXÞ ¼ Q i pðX i jPaðX i ÞÞ, where PaðX i Þ represents the parent set of X i . In our networks, each node represents transcription expression of a gene. These conditional probabilities reflect not only relationships between genes, but also the stochastic nature of these relationships, as well as noise in the data used to reconstruct the network.
Bayes formula allows us to determine the likelihood of a network model M given observed data D as a function of our prior belief that the model is correct and the probability of the observed data given the model: P(M|D)~P(D|M)*P(M). The number of possible network structures grows super-exponentially with the number of nodes, so an exhaustive search of all possible structures to find the one best supported by the data is not feasible, even for a relatively small number of nodes. We employed Monte Carlo Markov Chain (MCMC) 56 simulation to identify potentially thousands of different plausible networks, which are then combined to obtain a consensus network. Each reconstruction begins with a null network. Small random changes are then made to the network by flipping, adding, or deleting individual edges, ultimately accepting those changes that lead to an overall improvement in the fit of the network to the data. We assess whether a change improves the network model using the Bayesian Information Criterion (BIC) 57 which avoids overfitting by imposing a cost on the addition of new parameters. This is equivalent to imposing a lower prior probability PðMÞ on models with larger numbers of parameters.
Even though edges in Bayesian networks are directed, we can't infer causal relationships from the structure directly in general. For example, in a network with two nodes, X 1 and X 2 , the two models X 1 ! X 2 and X 2 ! X 1 have equal probability distributions as pðX 1 ;X 2 Þ ¼ pðX 2 jX 1 ÞpðX 1 Þ ¼ pðX 1 jX 2 ÞpðX 2 Þ. Thus, by data itself, we can't infer whether X 1 is causal to X 2 , or vice versa. In a more general case, a network with three nodes, X 1 , X 2 , and X 3 , there are multiple groups of structures that are mathematically equivalent. For example, the following three different models, M1 : X 1 ! X 2 ; X 2 ! X 3 , M2 : X 2 ! X 1 ; X 2 ! X 3 , and M3 : X 2 ! X 1 ; X 3 ! X 2 , are Markov equivalent (which means that they all encode for the same conditional independent relationships). In the above case, all three structures encode the same conditional independent relationship, X 1 6 ?X 3 X 2 j , X 1 and X 3 are independent conditioning on X 2 , and they are mathematically equal Thus, we can't infer whether X 1 is causal to X 2 or vice versa from these types of structures. However, there is a class of structures, V-shape structure (e.g. Mv : X 1 ! X 2 ; X 3 ! X 2 ), which has no Markov equivalent structure. In this case, we can infer causal relationships. There are more parameters to estimate in the Mv model than M1, M2, or M3, which means a large penalty in BIC score for the Mv model. In practice, a large sample size is needed to differentiate the Mv model from the M1, M2, or M3 models.
Incorporating genetic data as a structure prior in the HAEC Bayesian networks. In general, Bayesian networks can only be solved to Markov equivalent structures, so that it is often not possible to determine the causal direction of a link between two nodes even though Bayesian networks are directed graphs. However, the Bayesian network reconstruction algorithm can take advantage of the experimental design by incorporating genetic data to break the symmetry among nodes in the network that lead to Markov equivalent structures, thereby providing a way to infer causal directions in the network in an unambiguous fashion 58 . We modified the reconstruction algorithm to incorporate eSNP data as prior as following: genes with cis-eSNP 59 are allowed to be parent nodes of genes without cis-eSNPs, but genes without cis-eSNPs are not allowed to be parents of genes with cis-eSNPs, pðtrans À > cisÞ ¼ 0. We have shown that integrating genetic data such as cisacting eSNP or eQTLs (excluding edges into certain nodes) improves the quality of the network reconstruction by simulations 53 and by experimental validations 7,58 . We note that in applying this particular version of the Bayesian network reconstruction algorithm (incorporating genetic information as a prior); if genetic   8). e, f Western blot analysis of MTHFD2 expression (e) and quantification (f) of plaque material from human subjects with unstable atherosclerotic plaque (SP) (n = 20), stable atherosclerotic plaque (SP) (n = 20), or non-atherosclerotic artery (NP) (n = 8). g-j Relative mRNA expression of Mthfd2, Phgdh, Shmt2, and Slc3a2 in mouse thoracic aortic rings kept in organ culture and exposed to medium (1% FCS) with or without oxPAPC and rapamycin as indicated for 8 h (normalized to 18SrRNA) (n ≥ 4). k, l Relative mRNA expression of Mthfd2 and Shmt2 in the endothelium of partially ligated left carotid artery (LCA) compared to healthy right carotid artery (RCA) 48 h post ligation (normalized to 18S rRNA) (n = 3). (*p ≤ 0.05 Student's t test). m Relative mRNA expression of Mthfd2 in the endothelium of the left carotid artery of ApoE−/− mice which were fed with high fat diet (HFD) for 0, 1, or 4 days (normalized to 18S rRNA) (n = 5). n-q Relative mRNA expression of MTHFD2, PHGDH, CEBPB and PCK2 in HAEC exposed to HDL from healthy human subjects (n = 10) or human subjects with CAD (n = 10) for 4 h. (*p ≤ 0.05 Mann Whitney test). r, s Western blot analysis of MTHFD2 (r) and quantification (s) of HAEC treated as in f for 24 h (n = 10) (*p ≤ 0.05 Mann Whitney test). Data are represented as mean ± SEM, *p ≤ 0.05 (ANOVA with Newman−Keuls post-hoc test if not otherwise indicated) information is not available or is ignored, the population is simply treated as a population with random genetic perturbations.
Averaging network models. Searching optimal Bayesian network (BN) structures given a data set is an NP-hard problem. We employed an MCMC method to do local search of optimal structures as described above. As the method is stochastic, the resulting structure will be different for each run. In our process, 1000 BNs were reconstructed using different random seeds to start the stochastic reconstruction process. From the resulting set of 1000 networks generated by this process, edges that appeared in greater than 30% of the networks were used to define a consensus network. A 30% cutoff threshold for edge inclusion was based on our simulation study 53 , where a 30% cutoff yields the best tradeoff between recall rate and precision. The consensus network resulting from the averaging process may not be a BN (a directed acyclic graph). To ensure the consensus network structure is a directed acyclic graph, edges in this consensus network were removed if and only if the edge was involved in a loop and the edge was the most weakly supported of all edges making up the loop.
Identification of key drivers and associated subnetworks. For each Bayesian network, we further identified the key causal regulators by examining the number of N-hob downstream nodes (NHDN) for each gene in the directed network 53 . For a given network, let μ be the numbers of N-hop downstream nodes and d be the out degrees for all the genes. The genes with the number of N-hop downstream nodes (NHDN) greater than μ þ σðμÞ are nominated as key drivers. The key drivers with degree above d þ 2σðdÞ, where d denote the number of downstream genes, become key drivers. These criteria identified genes with number of downstream nodes and number of outlinks significantly above the corresponding average value. Based on the causal networks which we constructed, we identified 29 key drivers in control Bayesian network and 27 key drivers in oxPAPC Bayesian network. The subnetwork associated with each key driver was defined as downstream nodes (i.e. neighbors in outdegree direction) with the key drivers as the seeding point. Edges remained the same as in the complete network. Cytoscape 3.4 was used for network visualization.
Assessment of Bayesian networks. To assess the accuracy of the HAEC Bayesian networks, we compared the BNs with several widely used databases of gene networks and gene sets: (1) 37,080 interactions covering 9465 genes from Human Protein Reference Database (HPRD) database 60 , (2) 195,859 high confident interactions covering 12,427 genes from STRING database 61 , (3) 1329 canonical pathways covering 8439 genes from MsigDB databases 62 , and (4) 11,174 Gene Ontology (GO) annotation sets (sets with size ≥ 200 are excluded) covering 11,508 genes. In particular, we calculated the percentage of our inferred gene-gene connections that are in existing protein/gene network databases, or within the same pathway in gene set databases. For comparison, we generated 100 random networks for corresponding BNs by using the degree.sequence.game function in the igraph R package, and estimated the accuracy of random networks.
Additionally, we assessed the predictive power of our BNs by using gene sets closely regulated in endothelial cells. In particular, we used two independent gene sets, which include 77 gene sets related with endothelial cell from MsigDB, and 400 siRNA gene signatures in HUVEC 63 . For siRNA gene signatures, we downloaded and preprocessed microarray data as previously described 63 , and genes with a z-score >2 and <−2 were defined as significantly upregulated and downregulated genes, respectively. Based on these gene sets, we compared accuracy of our BNs to that of widely used gene networks including HRPD database by calculating the percentage of gene−gene connections that are within the same gene set.
Human carotid atheroma, Biobank-of-Karolinska-Endarterectomy, and partial carotid artery ligation data. We retrieved human atherosclerotic plaque data 27 , which included 32 paired samples of atheroma plaque and macroscopically intact tissue. The mRNA expression levels were downloaded from Gene Expression Omnibus (GEO) with the corresponding GEO accession ID GSE43292. Each platform's probe ID was mapped to the corresponding gene symbol and the expression levels were averaged over multiple probes mapped to the same gene symbol. The significance level of differentially expressed genes was calculated by Wilcoxon rank sum test.
Similarly, mRNA expression levels in 126 human carotid plaque samples from the Biobank Each (BiKE) were downloaded (GSE21545) 28 . For correlation analysis of gene expression comparing MTHFD2 and other genes within the MTHFD2 network, Pearson correlation coefficient was calculated in R by using the cor.test function.
Our previously published microarray 29 was reanalyzed for expression levels of genes within the MTHFD2 network 24, 48 h and 2 weeks after partial carotid artery ligation surgery.
Genome-wide association studies (GWAS) of plasma metabolites and coronary artery disease. We collected significant SNPs (meta-analysis p value < 1×10 −4 ) associated with one of more than 400 metabolites in human blood in a genome-wide association study (http://metabolomics.helmholtz-muenchen.de/ gwas/) 25 . Candidate SNPs were mapped to genes if their physical locations were within ±5 kb of gene bodies. We also compared candidate SNPs in CARDIo-GRAMplusC4D (Coronary ARtery DIsease Genome wide Replication and Metaanalysis (CARDIoGRAM) plus The Coronary Artery Disease (C4D) Genetics) consortium. In particular, we collected significant SNPs with p value < 1×10 −4 based on CARDIoGRAMplus4D 1000 Genome-based GWAS study, which is a meta-analysis of GWAS studies using 1000 genomes with 38 million variants 26 . We searched significant SNPs whose physical location is within ±500 kb of coding region of genes within our MTHFD2 network.
SiRNA transfection and overexpression and proliferation. For siRNA treatment, endothelial cells (80-90% confluent) were transfected with Lipofectamine RNAiMAX transfection reagent (Thermo Fisher, Cat# 13778150) according to the instructions provided by Thermo Fisher. The target-specific siRNAs for MTHFD2 and PSAT1 were obtained from Invitrogen (Stealth RNAi, Cat# 1299001) and for ATF4 from Sigma-Aldrich (Cat# NM_001675). The scrambled siRNAs were obtained from Invitrogen (Stealth RNAi, siRNA Negative Control Med GC Duplex). Plasmid overexpression in endothelial cells was performed with the Neon electroporation system (Invitrogen) with FLAG-tagged ATF4 (gift from Yihong Ye, Addgene plasmid pRK-ATF4 #26114) 64 . Cells were lysed 24 h after electroporation. For proliferation assay cells were counted 48 h after transfection with a CASY cell counter (OMNI Life Sciences).
Oxygen consumption rate. The cellular OCR was analyzed using a Seahorse 96 extracellular flux analyzer (Agilent). HAEC were plated in Seahorse 96-well cell culture plates 1×10 4 cells/well 1 day before the assay and equilibrated for 1 h in Krebs Henseleit buffer (111 mM NaCl, 4.7 mM KCl, 1.25 mM CaCl 2 , 2 mM MgSO 4 , 1.2 mM NaH 2 PO 4 ) supplemented with 11 mM L-Glucose and 2 mM L-Glutamine. Cells were treated with 2.5 µM oligomycin A, 1 µM CCCP, 1 µg ml −1 antimycin A and 1 µM rotenone as indicated. OCR was normalized to protein content of the wells.
Amino acid profiling by mass spectrometry. For amino acid profiling, HAEC were lysed in 85% Ultra LC-MS methanol (Roth, Cat# HN41.1). Lysates were centrifuged for 10 min at 17,000 × g and 50 µl of supernatants were processed using the EZ:faast LC-MS free amino acid analysis kit (Phenomenex, Cat# AL0-7500) according to the manufacturer's instructions with minor modifications. An internal standard (10 µl) was applied to all samples and to the standard curve. The internal standards included homoarginine, methionine-D 3 and homophenylalanine. Analysis of metabolites was performed by LC-MS/MS using an Agilent 1290 Infinity LC system (Agilent) coupled to a QTrap 5500 mass spectrometer (Sciex). The intensity of the measured metabolite was normalized to internal standards and protein content of cell lysate pellets. Analyst 1.6.2 and MultiQuant 3.0 (Sciex, Darmstadt, Germany) were used for data acquisition and analysis.
ATP measurement. Extracellular ATP was measured in cell supernatant centrifuged for 3 min at 5000 × g and then 10 min at 17,000 × g using CellTiter-Glo Luminescent Cell Viability Assay (Promega, Cat# G7570). Luminescence was detected using an Infinite 200Pro plate reader (Tecan) and normalized to intracellular RNA concentration.
Cell migration. A scratch wound-healing assay was performed in 24-well plates. HUVEC were scratched and cultured in EBM containing 1% FCS. Endothelial cell migration into the scratched area was monitored by live cell imaging (Zeiss TIRF System LASOS77). The distance migrated was calculated using ImageJ software.
Spheroid outgrowth assay. HUVEC spheroids were generated as described 65 . Images were acquired with an Axiovert135 microscope (Zeiss). For quantification of the cumulative sprout number, ten spheroids per condition were analyzed with the help of the AxioVision software (Zeiss). Treatments with VEGF-A 165 were performed for 16 h with a concentration of 10 ng ml −1 .
Animal experiments. All animal experiments were performed in accordance with the National Institutes of Health Guidelines on the Use of Laboratory Animals. All animal procedures were approved by the local committees or regulatory authorities. Mouse experiments were carried out in the C57BL/6 strain and only male mice were used (Jackson Laboratory, Bar Harbor). Mice were fed ad libitum with standard chow diet (unless otherwise indicated) and experiments were carried out at an age of 8-12 weeks.
Aortic ring assay. Aortic ring assay was conducted as the following: aortae from ten male wild-type 6-8-week-old mice were removed, cleaned, and embedded in a collagen type I gel (BD Biosciences, Heidelberg, Germany) in a 48-well plate containing EBM medium supplemented with murine orthologous serum (2%). For Mthfd2 knockdown, scrambled control siRNA (Ambion) and Mthfd2 siRNA (Sigma-Aldrich) in a final concentration of 100 µM were added in the growth media for 24 h by use of the Lipofectamine RNAi reagent (Invitrogen). Glycine treatment was performed in 500 µM final concentration every 48 h. Tube-like structures were allowed to develop over 5 days. Thereafter, the samples were fixed (4% paraformaldehyde) and endothelial cells were visualized using antibody against Pcam1. The knockdown was evaluated by immunostaining against Mthfd2.
High-fat diet and oxPAPC ex vivo treatment. All animal procedures were approved by the local government authority Regierungspräsidium Darmstadt. For high-fat diet male ApoE−/− mice (C57BL/6 background) were purchased from Taconis M&B A/S (Ry, Denmark, strain B6.129P2-Apoetm1Unc N6) and bred at the local facility. Paigen-type (19.6% protein, 15.8% fat, 4.3% fiber, 5% ash, 24.4% starch, 7.9% sugar; 1.25% cholesterol, 0.5% Cholate) diet was purchased from Ssniff Spezialdiäten (Soest, Germany). Animals received paigen-type diet ad libitum for 1 or 4 days. After diet feeding, endothelial-specific RNA of carotid arteries were isolated 29 . For ex vivo oxPAPC treatment C57BL/6J mice were euthanized with isoflurane. Mice were coated with 70% ethanol to decontaminate the incision area and then perfused with cold HANKS solution. The thoracic aorta was detached from the thoracic vertebrae and cleaned from adipose tissue. Vessels were cut into 1 mm segments and incubated in EBM 1% FCS (endothelial basal medium + 1% fetal calf serum). Each vessel was split into 12 × 1 mm segments, after which four segments were used for each treatment group. All vessels were incubated overnight in EBM + 1% FCS and treated the following day. RNA was isolated according to the manufacturer instructions (Bio&Sell RNA Isolation Kit).
Mouse partial carotid ligation and RNA isolation and qRT-PCR. All animal procedures were approved by the Animal Care and Use Committee at Emory University. Mice were anesthetized with 3.5% isoflurane initially and then anesthesia was maintained with 1.5−2% during the entire procedure. Partial carotid ligation surgery was performed on the left carotid artery. Briefly, the LCA bifurcation was exposed by blunt dissection and three of four caudal LCA branches (left external carotid, internal carotid, and occipital arteries) were ligated with 6-0 silk sutures, leaving the superior thyroid artery patent. The contralateral RCA was left intact as an internal control. Mice were monitored post-surgery and were allowed to recover. Following surgery, analgesic buprenorphine (0.1 mg kg −1 ) was administrated. Post procedure, d-flow in the LCA of all animals was confirmed by ultrasonography to verify the success of partial ligation surgery.
After 48 h of partial carotid ligation surgery, endothelial-enriched RNA was collected from the carotids by flushing the lumen quickly with Qiazol. Briefly, LCA and RCA were quickly flushed with 150 μl of QIAzol lysis reagent (Qiagen) using 29G insulin syringe into a microfuge tube. The eluate was then used for total intimal RNA isolation using the miRNeasy mini kit (Qiagen) following the manufacturer's protocol.
Total RNA was polyadenylated and reverse transcribed for use in a two-step qPCR setup using High-capacity cDNA synthesis kit (ABI) and using Brilliant II SYBR Green QPCR Master Mix (Stratagene) with custom-designed primers on a Real-Time PCR System (StepOne Plus, Applied Biosystems). Fold changes between