Integrated Metabolomics and Transcriptomics Suggest the Global Metabolic Response to 2-Aminoacrylate Stress in Salmonella enterica

In Salmonella enterica, 2-aminoacrylate (2AA) is a reactive enamine intermediate generated during a number of biochemical reactions. When the 2-iminobutanoate/2-iminopropanoate deaminase (RidA; EC: 3.5.99.10) is eliminated, 2AA accumulates and inhibits the activity of multiple pyridoxal 5’-phosphate(PLP)-dependent enzymes. In this study, untargeted proton nuclear magnetic resonance (1H NMR) metabolomics and transcriptomics data were used to uncover the global metabolic response of S. enterica to the accumulation of 2AA. The data showed that elimination of RidA perturbed folate and branched chain amino acid metabolism. Many of the resulting perturbations were consistent with the known effect of 2AA stress, while other results suggested additional potential enzyme targets of 2AA-dependent damage. The majority of transcriptional and metabolic changes appeared to be the consequence of downstream effects on the metabolic network, since they were not directly attributable to a PLP-dependent enzyme. In total, the results highlighted the complexity of changes stemming from multiple perturbations of the metabolic network, and suggested hypotheses that will be valuable in future studies of the RidA paradigm of endogenous 2AA stress.


Introduction
The metabolic network of an organism consists of a complex system of biochemical reactions that together result in the behavioral characteristics of the organism. In general, the reactions that make up the metabolic network, and their integration, are governed by chemical and physical constraints placed on the cell by its environment. Microbes in particular face rapidly changing environments and their metabolic networks must be structured to absorb and/or respond to internal and external and external insults [1][2][3]. Technological advances continue to expand and provide deeper and more detailed global snapshots of features in a cell or population. As these large datasets accumulate, the focus turns to analyzing these data as a means to glean insights about fundamental processes in biological systems. To that end, untargeted metabolomics and comparative transcriptomics data are increasingly being used as the basis for building metabolic models that describe the physiological state of a cell [4][5][6]. Examples of this approach can be found in disciplines ranging from mammalian genetic disorders and cancer biology, to antibiotic modes of action and pathogen-host interaction [7][8][9][10]. While these analyses can contribute to defining the response to various cellular perturbations [11,12], they are most significant when confirmed by rigorous genetic, biochemical, and molecular biological experimentation. In theory, the combination of global discovery approaches with focused reductionist analyses provides an optimal approach to generate insights into the fundamentals of the metabolic network structure and cellular strategies [13,14].
Overall, S. enterica has over 40 different PLP-dependent enzymes with functions in amino acid metabolism, cofactor/coenzyme biosynthesis, and the synthesis of cell structure components [26,27]. While a handful of specific 2AA targets have been defined, less is known about the extent of 2AAdependent damage or how such damage, (i.e., enzyme inhibition) impacts the metabolic network.  Overall, S. enterica has over 40 different PLP-dependent enzymes with functions in amino acid metabolism, cofactor/coenzyme biosynthesis, and the synthesis of cell structure components [26,27]. While a handful of specific 2AA targets have been defined, less is known about the extent of Metabolites 2020, 10, 12 3 of 14 2AA-dependent damage or how such damage, (i.e., enzyme inhibition) impacts the metabolic network. Further, 2AA stress is expected to elicit a broad shift in the metabolome, observed as the cumulative response to damaging a subset of PLP-dependent enzymes. A previous comparative transcriptomics study highlighted transcriptional changes in a strain lacking ridA, and led to the identification of an uncharacterized role for RidA in motility [28]. A majority of the most significant transcriptional changes (fold-change > 2) had no direct connection to a PLP-dependent enzyme, suggesting they were a downstream consequence of 2AA damage.
This study was initiated to further our understanding of the global metabolic response of S. enterica to 2AA stress. Specifically, untargeted NMR metabolomic data were sought to complement existing transcriptomics data. Our understanding of the RidA system at a physiological, genetic, and biochemical level provided a critical basis to guide the interpretation of the global data sets.

Metabolite Levels Are Altered in an S. enterica RidA Mutant
Comprehensive 1 H NMR analysis was performed on the cellular contents of an isogenic pair of S. enterica LT2 strains that differed only at the ridA locus (DM9404 and DM3480) after growth to stationary phase in minimal glucose medium. These data were viewed through a biological lens, specifically considering the role PLP-enzyme-dependent enzymes have in the ridA paradigm [29,30], to extract insights into the effect of 2AA stress on the metabolic network that goes beyond information gleaned from past genetic, biochemical, and transcriptomic analyses.
In total, 37 metabolites could be structurally identified from the 1 H-NMR spectral data (Tables 1  and 2). Table 1 contains data from the 23 metabolites with significantly altered concentrations between wild-type and ridA strains, using a false discovery rate (FDR) adjusted p value < 0.05. Table 2 shows the 14 metabolites detected in all samples that did not meet the statistical threshold for a significant difference between strains. The data in Table 1 emphasize that the metabolic environment in a cell is impacted by the status of ridA. The relevant metabolites represented diverse areas of the metabolic network and did not suggest a simple model for specific site(s) of 2AA-dependent perturbation. A majority of the metabolites (10/23) were related to, or the product of, amino acid biosynthetic pathways and a significant number (5/23) were components nucleotide metabolism. The former was expected based on the prevalence of PLP-dependent enzymes (targets of 2AA damage) in amino acid metabolism [22][23][24][25]31], while a direct role of PLP-dependent enzymes in nucleotide metabolism was not obvious. The 37 metabolites that could be confidently identified represented a small sample from all metabolites produced by S. enterica and thus limited the conclusions that could be made about the global metabolic state of S. enterica ridA mutants.

Transcriptome Analyses of RidA Mutant Complements Metabolomics Data
A transcriptomic dataset obtained with the same wild-type and ridA mutant strains [28] was considered in combination with the metabolomics data to build a more complete picture of the metabolic changes that reflect the cellular response to 2AA accumulation. All gene expression data discussed were extracted from Borchert and Downs, 2017 [28] and is accessible at the National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO) (GSE103146; http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103146). As reported, the transcription of 186 genes in a S. enterica ridA mutant were significantly (FDR < 0.05) elevated and the transcription of 227 genes was significantly decreased when compared to the wild-type parental strain [28]. Of the 413 genes that were differentially expressed, 113 had a greater than 2-fold difference in expression between ridA and WT. These data emphasized that the cellular response to 2AA stress generated transcriptional changes of modest magnitude. Over-representation analysis was performed on all differentially expressed genes, regardless of fold-change intensity, to determine pathways that showed significant enrichment for the differentially expressed genes in a ridA mutant (Table 3) [32,33]. These analyses confirmed that gene networks involved in flagellar biosynthesis, chemotaxis, and epithelial cell invasion were significantly down-regulated in a ridA mutant, observations that previously led to the characterization of a motility defect in ridA mutants [28]. Genes encoding ribosomal proteins were also significantly down-regulated (Table 3). Since ribosome synthesis is tightly regulated in accordance with growth rate and amino acid availability, [34][35][36] this change was assumed to reflect the general cell status of ridA mutants, i.e., perturbed amino acid metabolism and reduced growth rate [25].
Over-representation analysis showed the KEGG pathways for thiamine metabolism (stm00730), biotin metabolism (stm00780), amino acid biosynthesis (stm01230), one-carbon metabolism (stm00670), ABC transporters (stm02010), and glycine, serine and threonine metabolism (stm00260) were all significantly over-represented as up-regulated in a ridA mutant (Table 3). Since many of these KEGG pathways contain PLP-dependent enzymes, they were the focus of hypotheses to define the metabolic perturbations resulting from 2AA accumulation. In cells that lack RidA, accumulated 2AA targets PLP-dependent enzymes, and precedent with IlvE, GlyA, and Alr [21][22][23] suggests the activity of their population is decreased by 20-60%. The decreased activity is expected to collectively dampen flux in multiple metabolic nodes and cause metabolite imbalance(s). A subset of these imbalances would be detectable by an untargeted metabolomics approach. Such metabolite imbalances have the potential to trigger transcriptional regulatory responses, such as those identified in the transcriptome data.

Folate Metabolism Is Perturbed in a RidA Mutant
Collectively, the transcriptional and metabolic shifts detected in a ridA mutant are consistent with a perturbation in folate metabolism. Tetrahydrofolate (THF) and its derivatives are essential cofactors used to facilitate one-carbon unit transfers in many important pathways [37]. The transcription of multiple genes involved in 5,10-methylenetetrahydrofolate (5,10-mTHF) biosynthesis (glyA, gcvTHP, tdh, kbl, pabB, and folE) were increased in a ridA mutant, as compared to wild-type (Table 3, Figure 2A). This expression profile indicates a response in cells that have been compromised in THF pools. Extensive biochemical genetic analyses support this finding, showing that accumulated 2AA from an S. enterica ridA mutant inhibits GlyA and reduces both glycine and 5,10-mTHF levels [22,25]. Interestingly, the expression of glyA and gcvTHP are induced in a ridA mutant, a response which may serve a role in counteracting the consequence of 2AA-dependent GlyA damage. The genes encoding threonine dehydrogenase (Tdh; EC: 1.1.1.103) and 2-amino-3-ketobutyrate CoA ligase (Kbl; EC: 2.3.1.29), which are involved in the conversion of threonine to glycine, are also induced in a ridA mutant. Induction of these genes could lead to the rerouting of threonine for glycine production, ultimately generating additional 5,10-mTHF via the glycine cleavage complex (GCV) (Figures 2A and 3). Transcription of pabB (1.45-fold) and folE (3.21-fold), is also elevated in a ridA mutant. Expression of these genes, which are involved in 7,8-dihydrofolate production, could increase the biosynthesis of THF for use by GlyA and the GVC complex during 5,10-mTHF production. The putative 5,10-mTHF limitation in a ridA mutant would explain the increased expression of multiple pathways that are dependent upon folate metabolism (methionine, purine, histidine, thiamine, biotin, and lipoate biosynthesis) ( Figure 2).
Further, the transcriptional repressor, MetJ, binds S-adenosylmethionine (SAM) and represses expression of a number of genes whose expression is increased in a ridA mutant (metA, metBL, metC, metE, metF, metMIQ, metR, and folE) [28,38]. Increased expression of these MetJ/SAM regulated genes suggests that SAM is limiting in the ridA background. In fact, SAM biosynthesis is reliant on both 5-methyltetrahydrofolate (5-mTHF) and purines (ATP), making this conclusion consistent with the data and interpretations above. An expected consequence of SAM limitation is the decreased production of products of pathways that rely upon SAM as a cofactor. Consistently, expression of genes involved in the SAM-dependent pathways of biotin (bioA, bioBFCD, and bisC), lipoate (lipA), and thiamine (thiCEFSGH and thiMD) biosynthesis were induced, consistent with a limitation for these critical metabolites in a ridA mutant ( Figure 2B).
Finally, 1 H NMR metabolomic analysis showed that a ridA mutant had lower endogenous levels of both formate and acetate (Table 1). During mixed acid fermentation, pyruvate-formate lyase (PflB: EC: 2.3.1.54), uses Coenzyme A (CoA) and pyruvate as substrates during the production of formate and acetyl-CoA. Acetyl-CoA is further processed to acetate. Importantly, previous characterization of an S. enterica ridA mutant has revealed that 2AA-dependent damage of GlyA, and the resulting defect in 1-carbon unit production leads to a bottleneck in CoA biosynthesis [22]. Therefore, the reduced formate and acetate content in a ridA mutant could be explained by a limitation of CoA in a ridA mutant.

Branched-Chain Amino Acid (BCAA) Metabolism Is Altered in a RidA Mutant
Some of the most significant changes, by magnitude, in a ridA mutant involved metabolites related to branched-chain amino acid biosynthesis (Figure 3). Metabolomic analysis with 1 H NMR showed that isoleucine and leucine content were decreased in a ridA mutant, while valine, 2-isopropylmalic acid, and threonine were increased (Table 1). IlvE is a target of 2AA damage and the resulting decrease in activity could contribute to the reduced levels of leucine and isoleucine [21]. Down-regulation of the leu operon and ilvC could also play a role in reducing leucine and isoleucine content in the ridA mutant background (Figure 3). Pyruvate accumulates during exponential growth of ridA mutants as a result of CoA limitation [22], and since pyruvate accumulation can increase flux toward L-valine biosynthesis [39,40], this could lead to the increased valine and 2-isopropylmalic acid levels observed in a ridA mutant ( Table 1). As described above, the final steps of leucine and isoleucine biosynthesis involve a known target for 2AA damage (IlvE). It is possible that constriction of leucine and isoleucine generation, owing to damage of IlvE by 2AA, causes the accumulation of their precursors, 2-isopropylmalate and threonine, respectively ( Figure 3).
Metabolites 2020, 10, x 7 of 14 content in the ridA mutant background (Figure 3). Pyruvate accumulates during exponential growth of ridA mutants as a result of CoA limitation [22], and since pyruvate accumulation can increase flux toward L-valine biosynthesis [39,40], this could lead to the increased valine and 2-isopropylmalic acid levels observed in a ridA mutant ( Table 1). As described above, the final steps of leucine and isoleucine biosynthesis involve a known target for 2AA damage (IlvE). It is possible that constriction of leucine and isoleucine generation, owing to damage of IlvE by 2AA, causes the accumulation of their precursors, 2-isopropylmalate and threonine, respectively ( Figure 3).  Grey gene products showed no significant difference between wild-type and ridA cells, while grey metabolites were either not identified or not significantly altered. A bold border around a gene product designates enzymes that use PLP as a cofactor. GlyA is a known target for 2AA damage and is shown with a red border. Dashed arrows and triangles indicated general pathways which utilize the metabolite at the base of the arrow, while blue and red colored arrows denote genes involved in biosynthesis of that compound up-or down-regulated, respectively. and triangles indicated general pathways which utilize the metabolite at the base of the arrow, while blue and red colored arrows denote genes involved in biosynthesis of that compound up-or downregulated, respectively.

Many Metabolic Perturbations Consistent with Transaminase Damage by 2AA
The 1 H NMR metabolomics data indicated there was a decrease in lysine, serine, and phenylalanine, and an increase in glutamate in a ridA mutant. Synthesis of lysine, serine, and phenylalanine all rely upon the activity of PLP-dependent transaminases, which use glutamate as an amino group donor [41][42][43]. Thus, these changes could reflect damage to one or more of these transaminases by 2AA. For example, IlvE, AspC, and aromatic amino acid aminotransferase (TyrB; EC: 2.6.1.57) facilitate the synthesis of phenylalanine from phenylpyruvate [42]. Damage of IlvE [21] and AspC [24] by 2AA may explain the phenylalanine decrease in a ridA background. Similarly, phosphoserine transaminase (SerC, EC: 2.6.1.52) facilitates the conversion of 3-phosphooxypyruvate to O-phospho-L-serine during serine biosynthesis, while SerC and N-succinyldiaminopimelate aminotransferase (ArgD, EC: 2.6.1.17) catalyze N-succinyldiaminopimelate production from Nsuccinyl-L-amino-6-ketopimelate during lysine biosynthesis (Figure 3) [41,43]. The susceptibility of SerC to 2AA damage has not been confirmed, but reduced SerC activity would be expected to decrease both serine and lysine levels. Metabolites and gene products involved in branched-chain amino acid biosynthesis are shown using circles and rectangles, respectively. Blue shading represents significantly (FDR < 0.05) increased abundance of the respective metabolite/transcript in ridA cells, while red shading indicates significantly decreased abundance in ridA cells. Grey gene products showed no significant difference between wild-type and ridA cells, while grey metabolites were either not identified or not significantly altered. A bold border around a gene product denotes enzymes that use PLP as a cofactor. IlvE (red border) is a known target for 2AA damage. Dashed arrows and triangles indicated general pathways that utilize the metabolite at the base of the arrow.

Many Metabolic Perturbations Consistent with Transaminase Damage by 2AA
The 1 H NMR metabolomics data indicated there was a decrease in lysine, serine, and phenylalanine, and an increase in glutamate in a ridA mutant. Synthesis of lysine, serine, and phenylalanine all rely upon the activity of PLP-dependent transaminases, which use glutamate as an amino group donor [41][42][43]. Thus, these changes could reflect damage to one or more of these transaminases by 2AA. For example, IlvE, AspC, and aromatic amino acid aminotransferase (TyrB; EC: 2.6.1.57) facilitate the synthesis of phenylalanine from phenylpyruvate [42]. Damage of IlvE [21] and AspC [24] by 2AA may explain the phenylalanine decrease in a ridA background. Similarly, phosphoserine transaminase (SerC, EC: 2.6.1.52) facilitates the conversion of 3-phosphooxypyruvate to O-phospho-L-serine during serine biosynthesis, while SerC and N-succinyldiaminopimelate aminotransferase (ArgD, EC: 2.6.1.17) catalyze N-succinyldiaminopimelate production from N-succinyl-L-amino-6-ketopimelate during lysine biosynthesis (Figure 3) [41,43]. The susceptibility of SerC to 2AA damage has not been confirmed, but reduced SerC activity would be expected to decrease both serine and lysine levels.

Bacterial Strains, Chemicals and Media
Strains used in this work are derivatives of Salmonella enterica subsp. enterica serovar Typhimurium str. LT2. The ridA null mutant (DM3480) used in this study contains a MudJ1734 transposon insertion [44] disrupting the ridA locus (ridA3::MudJ1734) [45]. Strain DM9404 is isogenic to DM3480 and has a wild-type ridA locus. Minimal medium was no-carbon E medium (NCE) containing 1 mM MgSO4 [46], trace metals [47], and 11 mM D-glucose. Difco nutrient broth (NB) (8 g/L) supplemented with 5 g/L NaCl was used as rich medium. All chemicals were purchased from the Sigma-Aldrich Chemical Company (St. Louis, MO, USA).

Metabolomics Cell Preparation
Ten biologically independent cultures each of wild-type (DM9404) and ridA mutant (DM3480) strains were grown aerobically overnight in NB medium at 37 • C and used to inoculate (1% inoculum) 250 mL minimal glucose medium in 500 mL non-baffled flasks. Flasks were randomly arranged in an Innova ® 44 incubator and cultures were allowed to grow to stationary phase (16 h) shaking at 180 RPM and 37 • C. Cultures were cooled on ice for 5 min and then harvested by centrifugation at 7000× g for 10 min at 4 • C. The supernatant was decanted, pellets were resuspended in 10 mL ddH 2 O, and transferred to sterile 15 mL conical tubes in which they were pelleted at 7000× g 10 min at 4 • C. Final supernatant was decanted before pellets were frozen in liquid nitrogen and stored at −80 • C prior to cell extractions.

Metabolite Extraction
Each frozen bacterial pellet (~1 g) was thawed on ice and homogenized 3 times. For homogenization, 2.  The 1D proton spectra were collected using the pulse sequence 'noesypr1d' from the Bruker library. A mixing time of 5 ms was used, and during the acquisition, 65,536 complex datapoints were collected for the FID using 128 scans, with four additional dummy scans for equilibration and 4 s between scans. The spectral width was 20 ppm. For processing, 32,768 points were used for the spectrum, using an exponential window function of 0.3 Hz before the Fourier transform (FT). After FT and phase correction, a polynomial baseline correction of order 3 was used, and the frequency calibrated to the DSS peak (0.00 ppm).

J-RES
The 2D J-RES spectra were collected using the pulse sequence 'jresgpprqf' from the Bruker library. 32,768 points were used for the FID in the direct dimension and 64 points for the indirect dimension, with a spectral width of 20 ppm for the direct dimension and 50 Hz for the indirect dimension. Eight scans were acquired for each t 1 increment, along with four dummy scans for the equilibrium and 2 s between scans. For processing, 32,768 points were used for the direct dimension of the spectrum and 512 points for the indirect dimension. An unshifted sine window function was also used in each dimension, along with backward linear prediction of 64 points in the indirect dimension to increase the sensitivity as well as the resolution [48]. A polynomial baseline correction of order 3 in the direct dimension and order 5 in the indirect dimension was applied after the FT in both dimensions. Spectra were processed using in-house MATLAB scripts (https://github.com/artedison/Edison_Lab_Shared_Metabolomics_UGA) and all raw and processed metabolomics data are available on the Metabolomics workbench (http: //www.metabolomicsworkbench.org/).

Compound Identification/Database Matching
Two-dimensional 1 H-13 C heteronuclear single quantum correlation (HSQC) and 1 H-13 C HSQC-TOCSY (HSQC-total correlation spectroscopy) experiments were used to aid in annotation of 2D J-RES metabolites and were not used for statistical analysis. A total of 37 metabolites were identified in bacterial extractions using COLMARm [49], comprehensive metabolite identification strategy using multiple two-dimensional NMR spectra of a complex mixture implemented in the COLMARm Web Server and assigned a confidence level ranging from 1 to 5, as previously described [50]. In short: (1) putatively characterized compound, (2) matched reported 1D spectra, (3) matched reported HSQC spectra, (4) matched reported HSQC and HSQC-TOCSY spectra, and (5) validated by spiking putative compound into sample. A list of these assignments can be found in Supplementary Table S1. In addition to the database matching, using 2D NMR data, we identified nicotinate mononucleotide, which was not in a spectral database. We verified this identification with a spiking experiment (Supplementary Figure S1).

Quantification of Metabolomics Data
Two-dimensional (2D) J-RES spectra were used for quantification and annotation of metabolites. A representative 2D J-RES spectra, with significantly altered metabolites annotated, is provided in Figure S2. The binned feature used for integration of each metabolite included in the statistical analysis are listed in Table S1. Univariate statistics were performed on metabolites identified using 2D J-RES data after PQN-normalization. A student's t-test with an FDR correction was used to determine metabolites significantly altered in bacterial extracts from wild-type (n = 9) and ridA mutants (n = 10; FDR-corrected p value < 0.05).

Pathway Analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that were induced or repressed in S. enterica ridA mutants were determined using the Database for Annotation, Visualization and Integrated Discovery (DAVID) pathway enrichment analysis tool [32,33]. Input was a list of the locus tags (STM numbers) from the significantly (FDR < 0.05) differentially regulated genes separated into "up-regulated" and "repressed" subsets. Of the 413 locus tags used as input, 228 (55.2%) could be placed into a KEGG pathway. The output for up-regulated and down-regulated genes and which KEGG pathway they mapped to is provided in Supplementary File 1 (Supplementary Tables S2  and S3, respectively). These tables contain values representing the total gene count of pathway, hit counts, over-representation analyses (ORA) of raw p values, gene identities for hits, and ORA Benjamini-Hochberg FDR values [51]. Significantly over-represented pathways were defined as pathways with an FDR < 0.1. Figures 2 and 3 were generated using Cytoscape 3.7.1 software [52].

Conclusions
The current study used comparative transcriptomics and untargeted 1 H NMR metabolomics analysis to uncover global metabolic consequences of removing RidA. Importantly, all S. enterica ridA mutant-associated phenotypes to date are not only eliminated by expression of wild-type ridA in trans, but exist as the consequence enamine species (2-aminoacrylate) accumulation [15,17,22,23,25,28,31,45,[53][54][55]. The transcriptional changes seen for a ridA mutant permeated multiple KEGG pathways, including amino acid, coenzyme, and folate metabolism. Complementing this finding, 1 H NMR revealed multiple metabolic changes that were involved in amino acid, mixed acid fermentation, and nucleic acid metabolism. Some of the patterns from the transcriptomics and metabolomics data were explained by extrapolating known consequences from 2AA stress. For instance, induced expression of folate-related genes can be tied to damage of GlyA, and decreased isoleucine and leucine content with damage to IlvE. Indeed, many of the shifts in threonine, serine, lysine, and phenylalanine were consistent with the 2AA-dependent damage of PLP-dependent enzymes involved in the metabolism of these metabolites.
Significantly, numerous other metabolic and transcriptional changes in a ridA mutant could not easily be attributed to the inhibition of one or more PLP-dependent enzyme(s), the known targets of 2AA. These changes are expected to be the consequence of downstream effects on the metabolic network. Importantly, while they are difficult to tie to a specific 2AA target, these effects often prove important for the physiological state of the organism and cause visible phenotypes. This study took advantage of a mechanistically understood but complex paradigm of metabolic stress. The accumulation of 2AA causes a cellular stress that is multipronged, yet subtle, since the activity of multiple enzymes is reduced by 20-60%. The data herein showed that defining changes to the metabolic network created by multiple small perturbations to the system is not straightforward with the tools currently available. The RidA system can provide a valuable template for the continued refinement of global approaches to understand metabolic network structure, since it is understood at a biochemical level and can be manipulated genetically to test hypotheses that arise.