To decipher a catalogue of metabolic biomarkers pertinent to glutamatergic neuronal commitment, without the interference of complex genetic background diversity, we used hESC (H9) line. hESCs were maintained in an established feeder-free culture system in vitro prior to neuro-ectodermal lineage commitment (Szczerbinska et al., 2019). The stage-specific differentiation to glutamatergic neurons is comprised of four distinct stages “hESCs-to-Embryonic Bodies (EBs)-to-Rosettes-to-human Neuroprogenitors (hNPCs)-to-Neuron.” These stages were established in vitro using modified protocols (Cho et al., 2008; Boisvert et al., 2013). To attain these stages, we prepared four different combinations of media for sequential progression of neural commitment: namely, neural induction medium (NIM), neural proliferation media 1 and 2 (NPM1, NPM2), and terminal differentiation media (TDM). We used NIM to induce neuro-ectodermal lineage specification, marked by formation of EBs on 4^th day in vitro (DIV). During induction and the early stage of neural rosette formation, we used both BMP/TGF inhibitors. To block non-neuronal lineages, we used Dorsomorphin (4 μM) + SB431542 (10 μM) combination. Following induction, we used NPM1 with bFGF/EGF (20 ng/mL) and Cyclopamine (1 μM) to induce proliferation of neural rosettes. Rosettes were dissociated into a single cell suspension with Accutase (Gibco) and seeded directly onto freshly prepared poly-L-ornithine (PLO, 10 mg/mL) and Laminin (LN, 1 mg/mL)-coated culture plates in NPM1 to develop hNPCs. Dorsal patterning of hNPCs was continued using NPM2 with a specialized cocktail of growth factors IGF+ BDNF+ GDNF (all 10 ng/mL) along with Cyclopamine (1 μM) for inhibiting ventral patterning into GABAergic neurons. Committed dorsally patterned hNPCs were seeded onto new PLO/LN coated culture dishes (100 mg/mL and 10 mg/mL respectively) and were driven to maturation using TDM with IGF, BDNF, and GDNF (all 10 ng/mL), and cAMP (100 μM). To prevent clumping of mature neurons in culture, we used 2 μg of LN in the media.

To characterize stage-specific populations, we employed flow cytometry (BD Bioscience) (Yuan et al., 2011). We used cell specific markers to examine cell population purity at each stage: CD24 for EBs, ZO-1 for Rosettes, Nestin for hNPCs, and MAP2 and vGLUT1 for neurons. We used two markers for neurons to establish the glutamatergic neuron population in our cultures. Flow cytometry data were analyzed using FlowJo software.

We validated hESCs and their guided neural progeny in vitro using primary antibodies against stage-specific markers such as CD184 (Anti-CD184 Mouse, Invitrogen) and CD24 (Anti-CD24 Rat, Invitrogen) for EBs; Pax6 (Anti-Pax6, Avain, Santa Cruz) and ZO-1 (Anti-ZO-1 Rabbit, Invitrogen) for Rosettes; Nestin (Millipore) for hNPCs; VGLUT1(Anti-VGlut1 rabbit, Invitrogen), MAP2 (Anti-Map2 Mouse, Invitrogen), and PSD95 (Anti-PSD95 rabbit, Invitrogen) for neurons. To identify the stages, we plated the cells onto German cover slips (18 mm) coated with appropriate substrates (Cho et al., 2011). Approximately 5×10⁴ cells were plated per coverslip at each stage for immunofluorescent analysis. 4% paraformaldehyde (5–10 min) was used for fixation, followed by 1X phosphate-buffered saline (PBS) washes. Fixed cultures were kept for 1 h at room temperature in freshly prepared blocking solution in 1X PBS with 1% (w/v) bovine serum albumin (BSA; 0.3% (v/v) Triton X-100 and 0.05%) sodium azide and then treated with stage-specific primary antibodies prepared in the blocking buffer, overnight at 4°C. Next, secondary antibodies were applied at 1:500 dilution for 60 min at room temperature in the dark. A 1 mg/mL 4′,6-diamidino-2-phenylindole (DAPI, Life Tech, 62,248) stock solution was diluted in PBS (1:1,000) before adding to the wells to stain the nuclei. After 10 min, the solution was washed off twice with PBS before mounting the coverslips on glass slides with the anti-fading reagent (Molecular probes, P36930). The cells were visualized at 20x magnification using Zeiss inverted fluorescence microscope (Zeiss 710 microscope).

Proton NMR spectra (1H-NMR) were generated for each cell type using the 800 MHz NMR spectrometer with cryoprobe (Avance; Bruker, Inc.) to assess specific metabolome profiles in differentiating cells. NMR was acquired on 1,000,000 cells for each population. Cells were washed with PBS and spinned down at 1,000 g for 1–2 min at 4°C to remove the remnant culture media (Ma et al., 2011). The pellet was resuspended in 450 μL of PBS, pH 7.0, in D₂O and 50 μL of 50 mM internal standard 2,2-dimethyl-2-silapentane-5-sulfonate (DSS) in D₂O. To examine the cell metabolome, samples were sonicated and kept at 4°C if data acquisition was within 3–4 h or at −80°C if data acquisition was later. Sonication breaks up the cells and allows the metabolites to exit into the PBS, increasing their mobility in the NMR field (Singh et al., 2023; Mielko et al., 2021; Leibel et al., 2022). All experiments were done in triplicate within the same culture (i.e., on the same day) to ensure reproducibility, and with 3–8 biological replicates for specificity. 1H-NMR 1D spectra were acquired for 30 min at 25°C. Before Fourier transform, FIDs were line-broadened to 1.0 Hz with an exponential weighting function, phase and baseline corrected for distortions and referenced to DSS chemical shift. To identify the J-coupled peaks that belong to the same molecule, 2D NMR was acquired for 24 h. 2D, 1H-1H NMR spectra (Total Correlation Spectroscopy (TOCSY)) was collected at 25°C. The pulse sequence dipsi2esgpph (homonuclear Hartman-Hahn transfer) was selected. Excitation sculpting with gradients was applied for water suppression. The spectral width was 13.95 ppm and the number of scans was 24. Spectra were pre-processed using TopSpin (Bruker, Inc.) and uploaded to Chenomx, Inc. software (Chenomx Inc., Edmonton, Canada) for peak identification. Chenomx has an 800 MHz library consisting of 380 metabolites, which is manually queried for each sample. We cross-referenced the spectral signature of each metabolite to the Biological Magnetic Resonance Bank (BMRB) public database.

(BMRB) public database. Upon complete query, a list of small molecules in the sample and their estimated concentration (based on the integrated area under the peaks) was derived (n = 90) (Supplementary Table S1). The identities of molecules were cross validated using 2D NMR data (Arnold et al., 2015).

As copious biological network data continues to be gathered in various omics domains, standardized methods which process complex network information into interpretable and visualizable representations can help to turn diverse datasets into novel insights. Toward this goal, we have used the combination of PCA analysis and GLASSO to analyze the relative contribution of all features in the NMR datasets and infer their relative networks. PCA is a statistical method that reduces the dimensionality of a dataset by calculating linear combinations of features which maximally capture the variance of the dataset (Gewers et al., 2022). The few dimensions that explain most of the variance (the principal components) are frequently isolated and used to cluster the data in a visualizable lower-dimensional space. Analyzing the relative contribution of each feature towards each PCA dimension can suggest a relationship between the features and samples which score highly along that dimension. When run over the full dataset, with samples from each of the five neural cell populations, PCA clearly separates the cell types in an unsupervised manner, but due to the five-dimensional classification of cell types, it is unclear how to associate even the most outlying features with specific cell types. By instead running PCA multiple times over each sequential cell-type pair, we produce a binary classification setup that allows the association of individual features (metabolites) with one of the two cell types. However, the PCA scores of most features cluster tightly around 0, and only the most outlying of the metabolites can be confidently associated.

In turn, GLASSO applies L1-regularization to the covariance matrix to produce a sparse covariance matrix. Producing a sparse matrix is especially valuable in higher-dimensional settings, like those of metabolomics, where a naive covariance matrix would yield too many connections to be easily interpretable. We used GLASSO on the full dataset (agnostic of cell type) to generate a sparse covariance network that reduces the noise of metabolome data, yielding a visually interpretable network of feature correlations revealing metabolome structure (Friedman et al., 2008; Mazumder and Hastie, 2012). This method captures both positive and negative correlations, while leaving some features with no edges removed from the graph.

Overlaying the GLASSO network with PCA scores allows the PCA scores of individual features to be placed within the context of the feature covariance relationships.

We use the PCA and GLASSO implementations from scikit-learn,¹ found in the sklearn.covariance packages, respectively. All plotting is done using matplotlib,² and data I/O makes use of pandas.³ All code is run with python version 3.7.12. The class ‘pca_bundle’ runs PCA and allows the easy visualization of PCA-related graphs, showing loadings, scores, and explained variance. All plotting functionality leverages the matplotlib package.

We ran PCA over the mean- and stddev-normalized dataset to produce 10 components, of which we use the first three. When applied to the entire dataset, GLASSO is run with an alpha selected by cross-validation via the sklearn.covariance.GraphLassoCV implementation (alpha = 0.56, 250 iterations), using the default values for the cross-validation procedure. When applying GLASSO to each cell-transition pair, we set the alpha to 0.9, as the truncated datasets for each pair are small and thus unstable under cross-validation.

The MetaboLINK code is publicly available.⁴ Additionally, we make publicly available an interactive python notebook through Google Colab, which walks through the application of PCA and GLASSO as well as their various visualizations. We hope this will facilitate the easy use of MetaboLINK for novel applications.

The metabolic maps and pathways derived from MetaboLINK analysis were manually generated by combining existing KEGG metabolic pathways (KEGG, n.d.) and cross-matching HMDB database (Human Metabolome Database, n.d.). MetaboAnalyst software (Xia and Wishart, 2016)⁵ has also been used for correlation analysis, Partial Least Squares—Discriminant Analysis (PLS-DA) and Pattern Hunter studies. Correlation analysis was used to visualize the overall correlations between different features. With PLS-DA, we calculated the Variable Importance in Projection (VIP) that is a weighted sum of squares of the PLS loadings considering the amount of explained Y-variation in each dimension. Furthermore, Pattern Hunter functionality in MetaboAnalyst was applied to enable the identification of metabolites that follow a predefined pattern of concentration changes. We applied this module to identify patterns for all metabolites across different developmental stages. The pattern to be tested needs to be specified as a series of numbers separated by “-” where each number describes the expected relative concentration change at the corresponding sampling point. Thus, a pattern designated as “1–2–3-4-5” would search for metabolites with linearly increasing values across the five corresponding sampling points (hESC-EB-Rosettes-hNPC-Neurons), whereas a pattern defined as “2–1–1-1-1” and “1–1–1-1-2” would describe a specific increase and decrease in just hESC and neurons, respectively.

To integrate multivariate data analysis, such as PCA, with the statistical network visualization conferred by GLASSO, we designed a new algorithm combining the two analytical approaches, PCA-GLASSO (Figure 1A). The algorithm processes high-dimensional datasets, including longitudinal metabolome data. To test the algorithm, we used hESCs that underwent induction and differentiation into four distinct neural lineages: EBs, rosettes, hNPCs, and neurons (Figure 1B). To ensure cell type specificity prior to NMR metabolomics profiling, we performed flow cytometry to examine cell type purity at each stage of differentiation (Figure 1B). For EBs and rosettes, we used surface markers such as CD24 and ZO-1, respectively. For neuroprogenitors, we used nestin, an intracellular marker, and for neurons, we used both MAP2 and vGLUT1 to ensure we had glutamatergic neurons in the culture (Figure 1B). Based on the flow cytometry data, we sorted EBs and plated sorted cells to induce their differentiation into rosettes. No other cell types were sorted, but at different stages we dissociated and plated onto new coverslips/wells to move forward along the differentiation cascade to neurons. Based on flow cytometry data, EB and hNPCs were least pure, thus containing a mixture of cells, and rosettes and neurons were mostly pure. In addition to flow cytometry, at each stage we performed immunocytochemistry to determine the expression of cell type-specific markers (Figure 1B).

At each stage of differentiation, we collected 1 million cells and performed NMR. The biologically different hESC-derived neural progenies had very different NMR spectra (Supplementary Figure S1), pointing at the unique metabolome at each lineage. Chenomx analysis (Chenomx, Inc.), a widely used software to identify metabolites in the NMR spectra, detected a total of 90 metabolites among the cell types. Interestingly, we could separate each cell stage based on the abundance of metabolites (Figure 2A), even though our EB and hNPC populations were not completely pure (Figure 1B). Globally, in hESCs we identified high abundance of glutamate-derived metabolites while glutathione, ATP, and NAD+ were more abundant in the neuronal cluster (Figure 2B). We observed direct linear reduction across different neural stages for ADP, AMP, adenosine, and carnitine (Figure 2B, Supplementary Figure S2A), with high concordance with nicotinamide-related metabolites (Supplementary Figure S2B). Among all the metabolites identified, amino acids represented 25% of the metabolites, while bioenergetics, peptides, and carbohydrates accounted for 18.2, 18.2, and 11.4%, respectively (Figure 2C).

To identify and rank differentiating features, we used PCA on detected metabolites from the developing neural cell lineages. This was followed by GLASSO analysis to generate a sparse feature covariance matrix and visualize feature networks. The PCA confirmed the distinct cell populations by clustering each cell lineage separately (Figure 2D). GLASSO, in turn, revealed a distinct network of sample features identifying clusters and sub-clusters among the samples analyzed (Figure 2E). While GLASSO alone revealed the structure of the metabolomics network, it did not provide insights into the relationships between individual metabolites and cell types. Therefore, we combined PCA with GLASSO to assess the relative contributions of all features in our longitudinal metabolomics dataset and to infer complex, interrelated networks among features (metabolites) (Figures 1A, 3). Subsequently, the metabolite nodes in the GLASSO network were colored based on PCA loading scores associated with cell-type transitions (see Figures 2D, 3 top panel). Each node represents a chemical compound, while the edges (lines) signify the set of reactions linking these metabolites. The network is depicted using solid lines indicating positive correlations and dashed lines representing negative correlations; the thickness of these lines reflects the strength of the interactions. Transition-specific PCA scores and the GLASSO network for each of the four cell-type transitions are overlaid on the full GLASSO network, shown in grey background (Figures 3A–D). The transition-specific networks are color coded per cell type, and the most significant metabolites are labeled for clarity. Node sizes indicate the magnitude of concentration changes for each metabolite during the transition, while the color intensity reflects the magnitude of the PCA scores. The color itself denotes the specific cell type with which the metabolite is associated. This visualization yields a graphic with a straightforward interpretation: similarly colored metabolites that cluster within the network are highlighted as functionally related sub-networks. This comprehensive visualization seamlessly integrates the overall network structure with individual transition scores and variations in metabolite concentrations. By doing so, it effectively highlights which metabolites and their associated subnetworks are implicated in each cell transition. The overarching network structure, presented in a neutral grey, establishes a reference framework that allows for easy comparison against the transition-specific networks, which are distinctly color-coded to represent their unique characteristics. This color distinction not only facilitates identification of the different transitions but also illustrates the specific cell types linked to each metabolite as well as functionally related sub-networks, enhancing interpretability and clarity.

By integrating insights from PCA with the metabolite network structure derived from GLASSO into a unified algorithm, MetaboLINK, this approach facilitates a network-level analysis of transitions between paired cell types. Notably, we identify key metabolites such as glutamate, 4-aminobutyrate, and 2-phosphoglycerate as components of the predicted subnetwork involved in the transition from hESCs to EBs (Figure 3A). These findings highlight the central role these metabolites may play in facilitating the transition process from hESCs to EBs, and potentially uncovering underlying metabolic shifts and differential cell-specific metabolic demand.

In the subsequent section, we focus on validating the biological significance of the networks highlighted by the MetaboLINK strategy, further exploring the functional implications of these key metabolites and their interactions within the metabolic landscape. This validation will enhance our understanding of how these networks contribute to cellular transitions and their broader biological contexts.

To navigate and interpret the intricate metabolic landscape of the neural cell types identified by the MetaboLINK for a deeper understanding of neural development and function, we highlight several biologically significant pathways that have emerged from the data.

First, MetaboLINK analysis of the NMR-identified metabolome of neural cells found higher abundance of glutamate and GABA in hESCs compared to other committed lineages (Figure 4A and Supplementary Figure S5). Glutamate and its related metabolite GABA are essential for brain function and neural development through several key mechanisms. Glutamate is a precursor in glutathione synthesis that contributes to maintaining cellular redox balance by influencing glutathione levels. Glutathione maintains neural redox balance, influencing reactive oxygen species (ROS) levels and progenitor cell proliferation (Renault et al., 2009). Additionally, glutamate serves as an energy source; it is metabolized to α-ketoglutarate by glutamate dehydrogenase, participating in the citric acid cycle and thereby aiding ATP synthesis, which is crucial for neuronal cell fate and mitochondrial function (Homem et al., 2014; Tohyama et al., 2016; Watkins, 2000). Mutations in genes associated with glutamate metabolism play a role in developmental disorders such as developmental and epileptic encephalopathy (MIM: 618328, 618,412), underscoring the critical role of maintaining proper glutamate levels for normal brain function (Neuray et al., 2020). The increased glutamate signaling in hESCs may result from enhanced uptake or reduced degradation. Undifferentiated hESCs exhibit highly efficient glutamine-to-glutamate conversion regulated by highly expressed mitochondrial glutaminase (GLS2) (Marsboom et al., 2016; Hensley et al., 2013). The efficient conversion of glutamine to glutamate by GLS2 in hESCs declines upon differentiation (Marsboom et al., 2016), and this may contribute to the highest glutamate level observed in hESCs.

Additionally, rosettes exhibited elevated levels of glycerol and serine, which are key components of phospholipids (Figure 4B). Glycerol serves as the backbone for phospholipid molecules, while serine contributes to the synthesis of phosphatidylserine and other essential phospholipids. The increased presence of these metabolites in rosettes highlights their critical role in membrane biosynthesis and cellular structure within these neural cell aggregates, as previously suggested (Arai et al., 2015; Sánchez-Ramírez et al., 2024).

In turn, hNPCs showed higher concentrations of branched-chain amino acids (BCAA) valine and leucine, supporting adaptive changes in BCAA metabolism for an energy shift (Figure 4C, Supplementary Figure S4) (Bifari et al., 2020; Abdi et al., 2022). BCAA metabolism is critical for learning and memory, hippocampal neurogenesis, and neuronal differentiation, underscoring its critical role in cognitive function and neural development (Abdi et al., 2022). Furthermore, hNPCs had increased content of arginine and related metabolites compared to other neural cell types (Geiger et al., 2016). Arginine induces a metabolic switch from glycolysis to oxidative phosphorylation (OXPHOS), which is highly relevant for neurogenesis (Iwata et al., 2023; Kondoh et al., 2010).

Finally, neurons exhibited higher levels of specific lipids, such as sn-glycerol-3-phosphocholine, crucial for membrane build-up and signaling compared to other cell types (Figure 4D and Supplementary Figure S6). Sn-glycerol-3-phosphocholine, an intermediate in glycerophosphocholines-derived phospholipids like phosphatidylcholine, is significantly abundant in mature neurons. Beyond their structural role in membranes, these lipids act as bioactive molecules, impacting prostaglandin pathways, peroxisome proliferator-activated receptors responses, and more generally, lipid metabolism (Furse and de Kroon, 2015). The pro-neurogenic effect of phosphatidylcholine (Magaquian et al., 2021; Kuge et al., 2014) aligns with our findings, emphasizing the necessity of precise control over mature neurons membrane composition. Furthermore, ethanolamine, a phospholipid component, was higher in neurons compared to other cell types, potentially influencing neurogenesis especially under inflammatory stress (Montaner et al., 2018; Magaquian et al., 2021) (Figure 4D). Finally, neurons exhibited elevated levels of glutathione and myo-inositol (Figure 4D and Supplementary Figure S2C), both known as ROS scavengers, supporting their role to counterbalance ROS generated during OXPHOS-ATP production (O’Brien et al., 2015).

Among the various metabolic hallmarks, a notable shift in bioenergetics was observed as cells transitioned from hESCs to mature neurons (Supplementary Figure S3). The metabolic fate of neuroprogenitors is intricately linked to a metabolic rheostat governing mitochondrial dynamics, an evolutionary conserved phenomenon from flies to humans (Bonnefont and Vanderhaeghen, 2021; Iwata and Vanderhaeghen, 2021). ATP generation pathways, such as glycolysis and mitochondrial OXPHOS, play a critical role in fulfilling the distinct energetic demands of hESCs versus the needs of differentiating progenies (Folmes et al., 2012; Becker et al., 2010; Birket et al., 2011; Zhang et al., 2014). Glycolysis, although less ATP-efficient than OXPHOS, is faster and suitable for scenarios demanding quick energy bursts, making it advantageous for immature or rapidly proliferating cells like hNPCs (Seo et al., 2018; Tohyama et al., 2016; Armstrong et al., 2010). Additionally, glycolysis supports overall cellular growth and proliferation. Conversely, mitochondria respiration through OXPHOS stands out as the most efficient pathway for ATP generation, which is the preferred mechanism for cells with high energy demands such as neurons (O’Brien et al., 2015; Lange et al., 2016; Ito and Suda, 2014; Roy et al., 2018; Iwata and Vanderhaeghen, 2021). The balance between these pathways is crucial in mitigating mitochondrial ROS production (Stacpoole et al., 2011) and meeting cell-specific energy demands across different developmental stages. We observed higher concentrations of AMP, ADP, lactate, pyruvate, and alanine in hESCs. Conversely, neurons displayed the highest levels of ATP and NAD+, indicating a transition towards OXPHOS activity (Figure 2A and Supplementary Figure S2A, S3). Increase in NAD+ (oxidized form of NADH) to NADH ratio activates a NAD + -dependent deacetylase to inhibit hNPC vs NPC self-renewal and promote their differentiation (Hisahara et al., 2008). Our data support a model in which glycolysis sustains stemness, while the transition to OXPHOS is closely linked to neuronal differentiation (Agathocleous et al., 2012; Zheng et al., 2016; Bifari et al., 2020).