Intra-Ramanome Correlation Analysis Unveils Metabolite Conversion Network from an Isogenic Population of Cells

ABSTRACT To reveal the dynamic features of cellular systems, such as the correlation among phenotypes, a time or condition series set of samples is typically required. Here, we propose intra-ramanome correlation analysis (IRCA) to achieve this goal from just one snapshot of an isogenic population, via pairwise correlation among the cells of the thousands of Raman peaks in single-cell Raman spectra (SCRS), i.e., by taking advantage of the intrinsic metabolic heterogeneity among individual cells. For example, IRCA of Chlamydomonas reinhardtii under nitrogen depletion revealed metabolite conversions at each time point plus their temporal dynamics, such as protein-to-starch conversion followed by starch-to-triacylglycerol (TAG) conversion, and conversion of membrane lipids to TAG. Such among-cell correlations in SCRS vanished when the starch-biosynthesis pathway was knocked out yet were fully restored by genetic complementation. Extension of IRCA to 64 microalgal, fungal, and bacterial ramanomes suggests the IRCA-derived metabolite conversion network as an intrinsic metabolic signature of isogenic cellular population that is reliable, species-resolved, and state-sensitive. The high-throughput, low cost, excellent scalability, and general extendibility of IRCA suggest its broad applications.

(MWAS) of metabolomic data sets from high-resolution mass spectrometry (MS), the links among target metabolites that characterize a certain disease (4)(5)(6). Therefore, strategies and methods for rapidly and comprehensively detecting and profiling such links are of interest.
To detect such links, variations in metabolism-related phenotypes are required. As a result, typically, a time series or condition series set of samples is characterized for metabolite contents first and then correlated across the samples, with the resolution of detection generally dependent on sample size (7)(8)(9). For example, in order to understand the mechanism of oil production in microalgae, the lipidomes of an industrial microalga, Nannochloropsis oceanica, were measured in triplicate by electrospray ionization MS over eight time points (3,4,6,12,24,48,72, and 96 h) under nitrogenreplete (N1) or nitrogen-depleted (N-) conditions; then the profiles were correlated across the 48 samples, which revealed several temporal patterns among the triacylglycerol (TAG) species that suggest distinct regulation mechanisms; e.g., TAGs with lower degrees of desaturation were induced at the early stages of N-, while those containing polyunsaturated fatty acids (PUFAs) increased considerably only at later stages (10).
However, within a single sample of an isogenic population of cells, at any given time or condition, variations in metabolism-related phenotypes among individual cells are inherent and universal, due to the stochasticity of gene expression in cells (11,12). Therefore, an intriguing question is, can such phenotypic variations among individual cells, instead of those among multiple samples, be exploited to predict the links among metabolism-related phenotypes? Specifically, there are three hypotheses: (i) every cell can be quantitatively profiled for its various metabolism-related phenotypes as an independent biological replicate; (ii) intercellular variations of the phenotypes can be measured rapidly and simultaneously for many cells; (iii) intercellular correlation of the phenotypes can unveil important functional features of the system.
To test these hypotheses, we employed single-cell Raman microspectrometry, which captures the in vivo chemical profiles of a cell in a rapid, label-free, and nondestructive manner (13)(14)(15). In a single-cell Raman spectrum (SCRS) (13)(14)(15), each of (or a combination of) its thousands of Raman peaks potentially represents a specific metabolism-related phenotype, such as the presence and the concentration of a metabolite synthesized by the cell (16). Therefore, just like a portrait can reveal multiple facial features from a human individual, an SCRS can unveil cellular phenotypes (i.e., functions) in a "landscape-like" manner, i.e., simultaneously revealing multiple metabolism-related phenotypes of the cell at that particular state (16). For example, in a single microalgal cell (e.g., Chlamydomonas reinhardtii and N. oceanica), starch content can be quantified by the specific Raman peaks of 478 cm 21 (C-C-C deformation) and 940 cm 21 (C-O stretching; C-O-C and C-O-H deformation; a-helix C-C backbone), while the content of triacylglycerol (TAG) can be modeled by 1,441 cm 21 (alkyl C-H 2 bend) and 2,851 cm 21 (C-H 2 ; C-H 3 asymmetric and symmetric stretches) (17,18). In addition, the degree of lipid unsaturation in a cell can be measured by I 1,658 /I 1,441 , i.e., the ratio of 1,658 cm 21 (allyl C=C stretches which are proportional to the amount of unsaturated C=C bonds) and 1,441 cm 21 (alkyl C-H 2 bends which are proportional to the amount of saturated C-C bonds) (18,19). On the other hand, based on the full spectra of SCRS, the contents of many compounds in a cell can be derived, such as protein, starch, and TAG in C. reinhardtii (18) and astaxanthin, and b-carotene and chlorophyll in Haematococcus pluviali (20), via chemometric multivariate methods including partial least square regression (PLSR) and multivariate curve resolution (MCR).
We previously proposed the "ramanome" concept (16,21). A ramanome is the collection of SCRS (one from each cell) randomly sampled from a given instance of an isogenic cellular population and thus represents a single-cell-resolution metabolic snapshot of the population (16,21). Here, by treating every SCRS in a ramanome as one biological replicate and each of its peaks or combination of peaks as a metabolismrelated phenotype, we established a formal framework called intra-ramanome correlation analysis (IRCA; Fig. 1). From just one single ramanome, by pairwise correlating all the Raman peaks of SCRS among the individual cells, IRCA unveils a network of potential metabolite conversions. Based on a series of ramanome-profiling experiments on C. reinhardtii mutants and additional microalgae, fungi, and bacteria, we showed that IRCA is a rapid, low-cost, high-throughput, landscape-like, and universally applicable method for unraveling metabolic features of cellular systems.

RESULTS
IRCA predicts the interconversion among starch, protein, and TAG from one single instance of an isogenic microalgal population. To test the IRCA concept ( Fig. 1), we employed the unicellular microalga of C. reinhardtii CC124 as the first model. Triplicate isogenic cultures of CC124 grown under nitrogen depletion (N-) were sampled from 16 time points from 0 h to 8 days (i.e., 48 ramanomes in total; see Fig. S1; Materials and Methods). For cellular protein, starch, and TAG contents, each sample was analyzed via two strategies: (i) at the population level, via various conventional techniques which all require metabolite extraction from cell lysates (bicinchoninic acid [BCA] protein assay kit, amyloglucosidase/a-amylase method and thin-layer chromatography plus gas chromatography/mass spectrometry [TLC-GC-MS]; Fig. 2A to C; also see Text S1) and (ii) at the single-cell level, via full-SCRS-based quantitative partial least square regression (PLSR) models for metabolite contents, which are noninvasive and label-free (18) (;20 randomly selected cells per sample and thus 60 per time point; Fig. 2D to F). The PLSR model for starch content, for example, was established using the population-level measurements (i.e., via the amyloglucosidase/a-amylase method; see Text S1) and the averaged SCRS of 20 cells in the corresponding biological replicate (i.e., one ramanome). The correlation coefficients (R 2 ) of protein, starch, and TAG contents between the two strategies above are 0.9924, 0.9892, and 0.9686, respectively, confirming high accuracy of simultaneous quantification of protein, starch, and TAG for a C. reinhardtii cell via its full SCRS (18). Notably, for each phenotype, at any instance of the population, the degree of intercellular heterogeneity is high, and can vary greatly along the 8 days (see Fig. S1G to I). Specifically, for protein, both minimal and maximal single-cell contents in the population were decreasing with time (see Fig. S1G), while for starch, the trends are exactly the opposite (see Fig. S1H). As for TAG, the maximal content was increasing, while the minimal content remained at the baseline (i.e., a subpopulation of non-TAG-producing cells was always present), with the degree of within-population heterogeneity continuing to grow (i.e., the "delta" in Fig. S1I).
Notably, at the population level, significant correlation (defined by the Pearson correlation coefficient, or r ) was observed in each of three phenotype pairs among the 48 populations-protein-starch (r = -0.926; Fig. 2A), starch-TAG (r = 0.639; Fig. 2B), and protein-TAG (r = -0.848; Fig. 2C). This indicates protein-starch conversion, protein-TAG conversion, and accordant change of starch and TAG contents (22,23). Intriguingly, at the single-cell level, correlation for the 960 cells collectively reached a consistent conclusion, albeit with lower correlation coefficients (r = -0.879, 0.286, and -0.229, respectively; Fig. 2D to F). Thus, the interphenotype correlation among pooled cells from all populations can recapitulate that among the populations.
We then probed whether such interphenotype correlation can be detected via cells from just one instance of a population, i.e., at each of the 16 time points (Fig. 2G to I). Intriguingly, for protein-starch, negative correlation (NC) is phase-specific and shows a temporal trend of gradual weakening, quite strong at each of the time points only before 2 days (r # -0.6, P , 0.01; Fig. 2G) yet absent at 7 days or 8 days (Fig. 2J). Thus, it is possible that the protein-starch conversion took place only at the early phase of N-. In fact, this prediction is supported by the routing of hydrolytic release of carbon skeletons from protein to the synthesis of starch granules during the first 2 days in C. reinhardtii under N-, as previously observed (22,23).
For starch-TAG, the situation is the opposite, as the trend of NC among cells within a population intensified along time ( Fig. 2H and K). Thus, in contrast to the positive correlation (PC) between starch and TAG at the population level (r = 0.639, P , 0.01; Fig. 2B), the within-population NC between starch and TAG was prominently present at the late phase (i.e., starting from 3 days; Fig. 2K), which predicts conversion between starch and TAG inside the cell. This single-cell-based prediction, which the population- contents (G and J), between starch (x) and TAG (y) contents (H and K), and between TAG (x) and protein (y) contents (I and L) are also shown. (M to P) Phenotypic correlation between starch (x) and TAG (y) contents modeled by pair of Raman peaks of 940 cm 21 and 2,851 cm 21 . Correlation of the starch and TAG contents at the population level (M) and single-cell level (N) was compared. Phenotypic correlation between starch (x) and TAG (y) contents (O and P) at each time point were shown. In the curves of temporal dynamics, r is the Pearson correlation coefficients of two phenotypes among single cells (**, P , 0.01; *, P , 0.05), with value indicating the mean of triplicates and the error bar indicating the standard deviation. Absence of correlation (r = 0) or presence of strong correlation (r $ 0.6 or r # -0.6) is highlighted by red horizontal lines. level analysis missed, is actually supported by two observations: (i) competition between starch and lipid synthesis for their shared biosynthetic precursors intensifies temporally (24)(25)(26)(27), and (ii) early build-up of starch may serve as the carbon source for lipid synthesis at a later subsequent phase (28)(29)(30)(31)(32).
For protein-TAG, unlike their very strong among-population NC (r = -0.848, P , 0.01; Fig. 2C), no strong correlation was found for each population, although the trend of NC gradually intensified ( Fig. 2I and L). In particular, at 2 days and beyond, the protein contents of individual cells were all already depleted to a very low level, despite their relatively wide range of TAG contents, which indicates the decoupling of protein-synthetic and TAG-synthetic pathways. This distinction in the temporal withinpopulation NC pattern between protein-starch ( Fig. 2J) and protein-TAG (Fig. 2L) suggests that (i) prior to 2 days, the majority of released carbon skeletons from proteins were routed to starch rather than to lipids and (ii) when starch biosynthesis was saturated at 2 days, substantial accumulation of TAG occurred which is converted partially from proteins (23). Therefore, by correlating SCRS-derived phenotypes among cells, IRCA can reveal intermetabolite conversions from just one instance of an isogenic population.
Notably, besides the full SCRS (18), individual Raman peaks can also quantitatively model single-cell starch, protein, and TAG contents. For example, 940 cm 21 (C-O stretching, C-O-C, C-O-H deformation, and a-helix C-C backbone) and 2,851 cm 21 (C-H 2 , C-H 3 asymmetric and symmetric stretches) can model starch and TAG contents, respectively (the correlation coefficient R 2 between bulk-biomass-based conventional method and SCRS-based method being 0.884 and 0.954, respectively; see Fig. S2A and B). Moreover, correlations via just these two peaks among populations (Fig. 2M) or among cells (Fig. 2N) are consistent with those based on the full SCRS ( Fig. 2B and E). When using only these two peaks to model single-cell starch and TAG, strong within-population starch-TAG NC was absent at the early stage of N-(0 h to 1 day, 3 days) yet emerged at the later stage of N-(2 days, 4 days to 8 days; Fig. 2O and P). These results from C. reinhardtii CC124, which are consistent with the full-SCRS-derived findings described above, raise the possibility of reconstructing a network of potential correlations among metabolites by (i) treating the intensity of each peak as the content of a potential metabolite (or class of metabolites), and (ii) pairwise correlating all the ;1,600 peaks in an SCRS among the cells in a ramanome.
Genetically validating IRCA by knockout and then complementation of starch synthetic genes. To validate such IRCA-based prediction of intermetabolite conversions, we employed a genetic approach (see Fig. S2C to K; Fig. 3; also see Table S1). CC4325, a mutant directly derived from CC124 by X-ray mutagenesis, is deficient of starch due to the knockout of sta1-1 (which contributes to starch synthesis by catalyzing synthesis of the activated glycosyl donor, ADP-glucose, from Glc-1-P and ATP [33]), and thus, the starch-TAG conversion (i.e., NC between starch and TAG peaks in IRCA) FIG 3 Validation of the IRCA method via a genetic approach using a series of targeted mutants of C. reinhardtii. Correlation between starch (x) and TAG (y) contents (modeled by 940 cm 21 and 2,851 cm 21 , respectively) among cells within a ramanome for the mutant strains (Table S1)

Intra-Ramanome Correlation Analysis
® should attenuate in this mutant. As expected, IRCA of the time series CC4325 ramanomes revealed that the strong starch-TAG NC among single cells was no longer present (0 to 96 h under N-; see Fig. S2D), in contrast to CC124 ( Fig. 2O and P). To validate the specificity of IRCA in revealing the starch-TAG conversion, we further collected the time series ramanomes for CC406, a wall-less strain derived from CC124 (34). For CC406, IRCA revealed strong NC (r # -0.6) starting at 6 h under N-and remained so afterwards, suggesting starch-to-TAG conversion (see Fig. S2E). Therefore, the cell wall does not interfere with the detection of starch-TAG conversion via IRCA.
Interestingly, for CC4326 (i.e., sta1-2), another direct CC4324 derivative whose starch synthesis pathway is only partially disrupted (34), the temporal correlation pattern in IRCA is distinct from CC4333 and CC4334, which are fully starchless (under N-; Fig. 3C); instead of the latter's full ablation of starch-TAG NC throughout 96 h, CC4326 exhibited strong starch-TAG NC, yet only at the middle phase of from 12 h to 48 h, consistent with transient accumulation and conversion of a low level of starch to TAG, and then ceased the conversion upon starch depletion. Notably, under N1, which serves as a control condition (no TAG is accumulated in C. reinhardtii under abundant medium N [34]), CC4326 did not exhibit the starch-TAG NC until at the very late phase (e.g., 96 h, when the medium N was depleted by algal growth); in contrast, for CC4333 and CC4334, no starch-TAG NC was apparent throughout the 96 h under N1, consistent with their genotypes (i.e., complete disruption of the starch synthetic pathway; see Importantly, for both CC4565 and CC4566, both genetically complemented strains with fully restored starch synthetic capability (i.e., the sta6 gene was reintroduced into CC4333 by overexpression via pSL18-STA6 [34]), the temporal patterns of starch-TAG NC under N-are identical to those of CC4324, which is the direct progenitor of CC4333 and carries an intact starch synthetic pathway; they all exhibited significant NC throughout the 96 h (under N-; Fig. 3D and E). Altogether, the starch-TAG NC in IRCA that indicates starch being converted to TAG is directly correlated with starch-synthetic genotypes. Therefore, IRCA can detect and model the interconversion between metabolites from a single instance of isogenic population.
IRCNs reveal novel product-product and substrate-product links from an isogenic population. In a ramanome, the 1,581 Raman peaks in each of the many SCRS represent over one million pairwise correlations among individual cells. Such richness of information indicates an enormous number of possible between-phenotype links (and, in particular, between-metabolite conversions). For example, for the CC124 N-7-day ramanome that consists of 60 cells (which exceeds the minimal sampling depth for IRCA; see Text S1; also see Fig. S4), all pairwise among-cell correlations of peaks unveiled, in total, 60,994 strong NC "links" (r # -0.6, P , 0.05) that formulate a 1,581node network ( Fig. 4A; Materials and Methods). In such an intra-ramanome correlation network (IRCN), a node represents a Raman peak in SCRS and thus a potential phenotype (e.g., a metabolite), while an edge indicates a tentative link between two phenotypes (e.g., conversion between two metabolites).
An IRCN can reveal valuable features of the system. For example, in a 51-peak module of the CC124 N-7-day IRCN ( Fig. 4B; also see Fig. S5A), besides 478 cm 21 and 938 cm 21 for starch and 1,658 cm 21 and 2,853 cm 21 for TAG, many peaks of known or unknown assignments are present, suggesting additional between-metabolite conversions. In particular, 526 cm 21 (phosphatidylserine; PS) and 577 cm 21 (phosphatidylinositol; PI) exhibit strong NC with 968 cm 21 , 1,302 cm 21 , and 1,263 cm 21 , which are all TAG markers (18) (Fig. 4B). As both PS and PI are main components of membrane lipids, these observations suggest the conversion between membrane lipids and TAG under N-at 7 days (37,38). These findings are supported by (i) the concomitant TAG accumulation and membrane lipid degradation under stress for C. reinhardtii (39) and (ii) the stable, modest rise of cellular FAME (fatty acid methyl ester) levels yet dramatic fall of chloroplast membrane lipid level at the early stage of N-in C. reinhardtii (40). Notably, comparison of the 16 time series IRCNs of CC124 under N-revealed a deepening trend of the PS-TAG correlation and of the PI-TAG correlation (see Fig. S5B to E), Contents of starch and TAG were quantified by 940 cm 21 and 2,851 cm 21 , respectively. Protein content was modeled via the full SRCS. The degree of DU was quantified by the ratio between 1,658 cm 21 (unsaturated C=C bonds) and 1,441 cm 21 (saturated C-C bonds). The D 2 O incorporation rate was quantified by the CD ratio, i.e., the ratio between the C-D bond area (2,040 to 2,300 cm 21 ) and the area of C-D plus C-H bonds (2,040 to 2,300 cm 21 and 2,800 to 3,050 cm 21 ). Blue lines and blue arrows represent strong conversions.
Intra-Ramanome Correlation Analysis ® suggesting the increasing extent of such conversions with time. Moreover, the NC between membrane lipids and TAG occurred as early as 8 h, much earlier than that between starch and TAG ( Fig. 2O and P; 2 days), indicating that the membrane-lipid-TAG conversion precedes the starch-TAG conversion. On the other hand, in the N-7day IRCN, the starch-TAG NC is present, an indication of the starch-TAG conversion (Fig. 4B). Therefore, this local IRCN module detects multiple concomitant conversion processes that all produce TAG and reveals their dynamic features.
Notably, the IRCN can be readily expanded to include nodes representing those phenotypes that are underpinned by multiple Raman peaks. For example, the degree of unsaturation (DU), a key feature that determines the application and value of lipids, can be quantified by the ratio between the two Raman peaks of 1,658 cm 21 (unsaturated C=C bonds) and 1,441 cm 21 (saturated C-H 2 bonds) (18,41,42). Inclusion of DU into the aforementioned time-series IRCNs revealed a local module of 17 peaks with functional assignment, where DU is highly negatively correlated with two starch peaks (938 cm 21 and 1,125 cm 21 ) and one protein peak (1,007 cm 21 ) at the late phase of N-( Fig. 4C; also see Fig. S6A and B; 2 days to 8 days). This suggests starch and proteins are converted to unsaturated lipids at the late phase of N-, consistent with the concordance between (i) degradation of proteins and starch and (ii) accumulation of unsaturated lipids (23,(28)(29)(30)(31)(32).
Moreover, IRCN can reveal the link between substrate intake and metabolite production of the cell. For example, cellular intake of D 2 O resulted in substitution of C-H bonds in intracellular macromolecules by C-D bonds. Therefore, the CD ratio, i.e., the ratio between the C-D bond area (2,040 to 2,300 cm 21 ) and the area of C-D plus C-H bonds (2,040 to 2,300 cm 21 and 2,800 to 3,050 cm 21 ), can measure cellular metabolic activity (43,44). In a separate time-resolved experiment, D 2 O was fed to CC124 immediately after the removal of medium nitrogen (i.e., sampled for ramanomes in triplicates at 0 h, 12 h, 1 day, 2 days, 3 days, and 4 days under N-; see Fig. S6C; Materials and Methods). IRCA revealed that (i) specifically at 2 days, a module of 17 characteristic peaks was formed, in which the CD ratio exhibited strong NC (r # -0.6, P , 0.05) with 5 TAG peaks (1,263, 1,443, 1,656, 2,855, and 3,011 cm 21 ; Fig. 4D, also see Fig. S6C and D), suggesting that at 2 days, higher TAG-content cells exhibit lower D 2 O-assimilating activity (and visa versa). Thus, at 2 days, the algal population was at a most "diverse" metabolic state, where both low-TAG, active-"drinking" cells and high-TAG, inactivedrinking cells were abundant, while before 2 days the former dominated, and after 2 days the latter dominated. (ii) Between CD ratio and DU, no strong NC (r # -0.6, P , 0.05) was observed at any of the time points (Fig. 4D, also see Fig. S6C and D), except for weak correlation at 3 days (r = -0.50, P , 0.01), suggesting that the DU of synthesized lipids increased with reducing C. reinhardtii vitality under the nitrogendepletion stress (this finding is supported by GC-MS data [18]). (iii) In contrast to the TAG peaks, peaks of starch (the major carbon storage form of C. reinhardtii under N-) showed no correlation with CD ratio at any of the time points. Thus, starch synthesis appears to be mainly supported by endogenously derived H, while TAG synthesis requires exogenously supplied H. Therefore, IRCA is a new strategy to track the cellular destination of target substrate (in this case, the water).
Altogether, a choreography of interplay among water intake and major cellular products was revealed (Fig. 4E).  number of all possible edges for the same nodes; Fig. 5E) and the average degree (ave_Degree, average number of adjacent edges; Fig. 5F), both depicting the degree of network complexity, increased. In contrast, the average Pearson correlation coefficient (ave_PCC; Fig. 5G) became more negative (from -0.000078 at 2 h to -0.0335 at 7 days), indicating more frequent and more active conversions between the implicated metabolites.
Across all time points, the dynamics of num_Node, num_Edge, Density, ave_Degree, and ave_PCC of 16 IRCNs were significantly correlated with the starch-TAG conversion (PCC of -0.653, -0.807, -0.807, -0.807. and 0.81, respectively; P , 0.05; Fig. 5H to L). In addition, Raman peaks at the carbohydrate and lipid region were usually prominently found in the nodes with the highest degree (see Fig. S7B). Moreover, the 1,457 cm 21 region (the alkyl C-H 2 bend of saturated lipids) dominated at the later phases under N-(i.e., 4 days, 5 days, 7 days, and 8 days), consistent with the turnover between lipid classes (saturated switching to unsaturated). Therefore, topological features of an IRCN can reveal both metabolically active compounds and their interconversions.
Based on the matrix of r (r # -0.6, P , 0.05), the 16 CC124 IRCNs under N-can be classified as clusters I and II ( Fig. 5M and 6N); I consists of those time points before 24 h, which are characterized by the global features of low num_Edge and low intensity of ave_PCC, and II includes those after 24 h, with global features of a sharp increase of num_Edge and of ave_PCC. Thus, 24 h is the temporal transition point (consistent with population-level measurement; see Fig. S1), when protein stopped degradation, starch stopped accumulation, and the early build-up of starch began to contribute to other cell activities, such as serving as the carbon source for lipid synthesis.
In addition, the global heatmap of Dr (the difference between max r and min r in a ramanome series) reveals the dynamic change of the "hot spots" of metabolite conversion. For example, in the 16-time point CC124 N-process (see Fig. S7C), the Dr between 1,441 cm 21 (alkyl C-H 2 bend, CH 2 scissoring, and CH 3 bending; i.e., saturated lipids) and 1,402 cm 21 (bending modes of methyl groups; i.e., proteins) is one of the most prominent hot spots (1.18; highlighted in Fig. S7C), with a max r of 0.619 at 2 h (strong positive correlation) and min r of -0.557 at 7 days (strong NC). This suggests the temporal switch, from the synergic degradation of lipids and protein at 2 h to the protein-to-lipid conversion at 7 d, is a key dynamic feature of this process.
Furthermore, such heatmaps of r can measure the global similarity of metabolite conversion modes among ramanomes from different strains (see Fig. S7D). For example, at N-72 h, comparison among IRCNs of CC4324 (wild-type), CC4333 (starchless mutant), and CC4565 (genetically complemented strain) reveal much higher similarity with CC4324 for CC4565 than for CC4333. Specifically, CC4333 is distinct due to the absence of several regions (dotted rectangle in Fig. S7D) found in both CC4324 and CC4565, such as (i) 867 to 970 cm 21 (conversion of starch to lipids) and (ii) 1,441 to 1,658 cm 21 , and 1,441 to 1,744 cm 21 (conversions of saturated lipids to unsaturated lipids; solid rectangles in Fig. S7D). These IRCN-derived phenotypes of the three strains are consistent with their genotypes (Fig. 3).
IRCN is a new "state"-specific metabolic signature of an isogenic population for diverse organisms. A ramanome is essentially a single-cell-resolution metabolic snapshot of the cellular population (16,21). To maximally extract the metabolic features from a ramanome (Fig. 1), three signatures are proposed, metabolite profile (MP; derived via the mean SCRS of a ramanome; Fig. 6A and B), metabolite interaction (MI; derived via the correlation matrix of a ramanome consisting of significant correlations of all pairwise Raman peaks [P , 0.05]; Fig. 6C and D), and metabolite conversion (MC; derived via a correlation matrix of a ramanome that consists of strong NC of all pairwise Raman peaks [r # 0.6, P , 0.05]; Fig. 6E and F). To test their general applicability across organisms, MP, MI, and MC were derived for each of the 64 ramanomes from C. reinhardtii (wild type [WT] and the starchless mutant series), N. oceanica, Saccharomyces cerevisiae, and Escherichia coli (see Table S1; Materials and Methods). To facilitate interspecies comparison, the temporal ramanome data set of CC124 under N-along 16 time points  Table S1. Details for the Raman barcodes and local IRCNs are also provided in Table S1. Global IRCNs were derived from Raman peaks from 600 cm 21 to 1,800 cm 21 . Clusters are colored based on HCA.

Intra-Ramanome Correlation Analysis
® was regarded as a standard "protein-starch-TAG" (PST) process, onto which the entire 64 ramanomes were "projected." Similarity-based clustering of MP reveals six modes ( Fig. 6A and B). The MP, MI, and MC signatures, although inherently linked, can be highly distinct (see Fig. S8). Specifically, MP appears to show more species specificity, while MI and MC exhibit more state-specificity. For example, E. coli ramanomes are clustered with algae and yeast in MC-1; in contrast, E. coli ramanomes are solely clustered as MP-6. However, MC can be more sensitive in detecting cellular state change than MP and MI. For example, at 4 h under N-, the metabolite content of CC124 has yet to change (i.e., clustered with 0 h in MP-1; it did not switch to MP-2 until 6 h), yet the mode of metabolite conversion has already altered (i.e., clustered with 6 h to 96 h in MC-2).

DISCUSSION
Exploiting a fundamental and inherent property of all cellular systems, i.e., the heterogeneity in single-cell metabolism-related phenotypes (e.g., metabolite contents and substrate intake) among individual cells, here, we proposed and then biochemically and genetically validated the IRCA approach. IRCA is advantageous, as it unveils a comprehensive and landscape-like network of potential links among metabolic activities, in this case metabolite conversions or links between substrates and products, from a single instance of an isogenic population of living cells, rather than requiring a time or condition series of samples. This ability is of particular importance for those phenotyping experiments where spatiotemporal or condition-resolved sampling is a hurdle or a constraint.
Despite its need of just one snapshot of an isogenic population, an IRCN is information rich. From an IRCN, a large number of hypotheses on potential links among metabolism-related phenotypes can be generated simultaneously, based on the potentially millions of pairwise correlations from the Raman peaks that each potentially represents a phenotype or from the combinations of Raman peaks that model particular phenotypes, among the individual cells sampled for SCRS in a ramanome. The scope of such phenotypes is broad (16,45), including but not limited to substrate intake (46), product synthesis (17,18), and response to environmental stress (e.g., antimicrobial susceptibility [21,44]). Based on the long and expanding list of such SCRSderived phenotypes (16,45), IRCNs can be constructed and then interpreted to mine the ramanome data space for new interphenotype links without a priori hypotheses. Moreover, since Raman microspectroscopy is label-free, noninvasive, and generally applicable to any cells, the strength of IRCA also includes high throughput, low cost, excellent scalability, and broad applicability. For example, acquisition of SCRS can be automated via flow-mode Raman cytometry or sorting (47,48), suggesting the possibility of ultrafast, robust yet low-cost acquisition of ramanomes for IRCA. Therefore, IRCA presents a new dimension of "metabolism-related phenome" for cellular systems, which serves as a highly species-and state-specific signature of metabolic activity that captures not just the profile of Raman-sensitive metabolites (as well as other SCRS-derived phenotypes) but their dynamic links. This capacity of IRCA would enable a data-driven research strategy for profiling cellular metabolism.
On the other hand, the potential of IRCA is limited by its frequent inability to unambiguously assign Raman peaks to specific metabolites. For example, in this study, although ;1,600 spectral resonance features were detected in each SCRS, only interconversion of macromolecules such as lipids, protein, and starch can be revealed, because a large number of vibrational features were shared among metabolites and cannot be assigned uniquely to individual metabolites. Even though the list of models that link spectral resonance features to a metabolite (or other metabolism-related phenotypes) has been rapidly growing (16), only a small portion of Raman peaks in an SCRS can be assigned to a defined metabolism-related phenotype at present. In addition, biological hypotheses for many of the edges, i.e., phenotype pairs showing significant NC, in an IRCN remain difficult to interpret or confirm. Therefore, new experimental and computational methods, such as stable-isotope probed SCRS (to track the assimilation of target substrate) and multivariate curve resolution-alternating leastsquares (MCR) algorithms (to deconvolute macromolecular components from overlapping Raman peaks [49,50]) should be developed to establish new assignments or to improve their specificity for the Raman peaks, so as to enrich and expand the actual information content of an IRCN.
Despite these present limitations, as demonstrated here using a number of microalgal, fungal, and bacterial species as examples, IRCA can serve as a valuable tool to rapidly discover features of metabolic dynamics of cellular systems, such as metabolite conversions. The capability of revealing such features from just one snapshot of an isogenic population has profound implications for designing phenotyping experiments.

MATERIALS AND METHODS
Strains and growth conditions. For Chlamydomonas reinhardtii, a total of nine wild-type (WT), lowstarch mutant and starchless mutant strains were employed (see Table S1). Specifically, CC124, CC406, and CC4324 were WT strains. The low-starch mutant CC4325 was derived by X-ray mutagenesis from CC124. The low-starch mutant CC4326 and the starchless mutants CC4333 and CC4334 were derived from CC4324 by random integration of cassette pARG7 into the nuclear genome (34,35,51). The starchless mutant CC4333 is deficient in the catalytic (small) subunit of ADP-glucose pyrophosphorylase, which interrupts synthesis of the ADP-glucose, a substrate for starch biosynthesis. The CC4334 mutant contains a disrupted isoamylase gene; thus, the level of starch is severely attenuated, but it accumulates a soluble glycogen-like product. CC4565 and CC4566 were genetically complementary strains of CC4333, derived by complementation with plasmid pSL18-STA6 (34). All these mutants can be obtained from the Chlamydomonas Resource Center (http://www.chlamycollection.org).
The C. reinhardtii cells were inoculated into the TAP (Tris acetate phosphate) liquid medium with or without arginine (100 mg ml 21 ) supplement under one-side continuous light (approximate 150 mmol photons m 22 s 21 ) at 25°C bubbled with air to ensure mixture and to prevent settling and were then grown to the late log phase in the N-replete TAP medium. Then they were reinoculated at 1 Â 10 6 cells/ ml in parallel into the nitrogen-replete TAP medium (N1), the nitrogen-depleted TAP medium (N-; in which NH 4 Cl was omitted), or the 25% D 2 O nitrogen-depleted TAP medium (25% D 2 O N-), each in triplicates. Cultures at each of a series of time points were sampled in triplicate (see Table S1A).
For the Nannochloropsis oceanica IMET1 strain, cells were cultured in a modified f/2 liquid medium with 4 mM NO 3  Acquisition of single-cell Raman spectra from an isogenic population of cells. Raman spectra of individual cells were acquired using modified Raman microspectroscopy equipped with a confocal microscope with a Â50 PL magnifying dry objective (numerical aperture [NA] = 0.55, BX41; Olympus, UK) and a 532-nm Nd:YAG laser (Ventus, Laser Quantum Ltd., UK). The scattered photons were collected by a Newton electron multiplying charge coupled device (EMCCD) (Andor, UK) utilizing a 1,600 Â 200 array of 16-mm pixels with thermoelectric cooling down to 270°C for negligible dark current. Before measurement, each sample was washed three times and resuspended in double-distilled water (ddH 2 O) to remove the culture medium. For algal and yeast samples, cells were loaded into a capillary tube (50-mm length by 1-mm width by 0.1-mm height; Camlab, UK) for SCRS acquisition. The power out of the objective was 100 mW. For each SCRS, the signal acquisition time was 2 s for C. reinhardtii and N. oceanica and 3 s for yeast. For C. reinhardtii and N. oceanica, prior to Raman signal acquisition, the cell was quenched with a 532-nm laser until the signal of fluorescent and resonantly enhanced biomolecules was no longer detectable. For each individual cell, a background spectrum was generated as the average of four spectra acquired from the liquid around the cell. For E. coli samples, cells were loaded onto a clean CaF 2 slide and air dried before Raman measurement, and the SCRS of E. coli samples were collected as described (21).
Notably, a eukaryotic cell is generally larger than the laser focal spot obtained with a lens objective. In the case of a C. reinhardtii (;10 mm in diameter), N. oceanica (2 to ;3 mm in diameter), or S. cerevisiae (2 to ;3 mm in diameter) cell, which was held in the single-beam gradient force trap under the aqueous conditions, the cell was rolling in random orientation when undergoing SCRS acquisition. Therefore, the SCRS acquired represents the overall metabolic state of the cell.
Intra-ramanome correlation analysis. The raw SCRS was first preprocessed with LabSpec 5 (Horiba Scientific, France), including background subtraction and baseline correction by a polynomial algorithm (with degree of 7). Then SCRS were normalized, followed by IRCA, which used a customized computational pipeline for data analysis and result visualization.
Due to potential technical variation (e.g., change of data collector, batch, etc.), SCRS from different ramanomes may have distinct spectral ranges and resolution. Therefore, prior to computing the MP signature and the MI and MC networks, the 64 ramanomes were standardized in three steps: (i) for spectral range, only the "fingerprint area" (600 to 1,800 cm 21 ) was extracted; (ii) spectral resolution was simulated to 1 cm 21 , via the interpolation algorithm; (iii) spectral normalization was performed by division by its area.
For each IRCN, the correlation matrix of each ramanome was constructed by calculating the Pearson correlation coefficient (PCC; r ) of all possible pairwise combinations of Raman peaks, among all the 60 cells sampled for the ramanome. Those pairs of Raman peaks with significant correlation were considered candidates that potentially indicate links between two metabolites (P , 0.05), while those with strong negative correlations suggested potential conversions among two metabolites (r # -0.6, P , 0.05). The igraph package in R (http://www.r-project.org) was used to derive the key network properties and visualize the IRCNs (either from specific subsets of Raman peaks or from all of the Raman peaks). To probe the global features of an IRCN, key network parameters including number of nodes (num_Node), number of edges (num_Edge), number of modules (num_Module; each module is a sub-IRCN not connected with any other nodes in the IRCN), size of the largest module (size_largest_Module; number of nodes in the largest module of each IRCN), density (Density, ratio of the number of edges divided by the number of all possible edges of the same nodes) and average degree (ave_Degree, average number of adjacent edges), and average Pearson correlation coefficient (ave_PCC, sum of all significant strong negative correlations divided by all nodes), were derived. For global IRCNs, all Raman peaks were used. For simplified versions of IRCN that facilitate visualization, characteristic marker Raman peaks were used (669, 783, 814, and 1,481 cm 21 for nucleic acids; 622, 643, 758, 1,003, 1,176, 1,211, 1,246, 1,584, 1,606, and 1,619 cm 21 for proteins; 865, 940, 1,033, 1,045, and 1,127 cm 21 for carbohydrates; and 725, 971, 1,083, 1,265, 1,305, 1,441, 1,450, 1,658, and 1,742 cm 21 for lipids; see Table S1B).
To characterize, compare, and cluster ramanomes, the three signatures of MP, MI, and MC were proposed. For MP, the mean SCRS of each ramanome was used. The MI network of a ramanome was generated as follows: (i) a correlation matrix was constructed by calculating the PCC of all possible pairwise Raman peaks; (ii) significant correlations (P , 0.05) were tabulated, while those with no significant difference were counted as 0 (no correlation). The MC network of a ramanome was generated by including only those significant, strongly negative correlations (r # -0.6, P , 0.05). To measure the similarity, pairwise Euclidean distances were calculated. Then hierarchical cluster analysis (HCA) was performed with Ward's algorithm, and six clusters were produced (52,53). Principal-component analysis (PCA; the factoextra package in R) was used to visualize the MP, MI, and MC signatures.
Data availability. The data that support the findings of this study are available from the corresponding author, Jian Xu, upon reasonable request.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. TEXT S1, DOC file, 0.3 MB.