Transcriptome-guided metabolic network analysis reveals rearrangements of carbon flux distribution in Neisseria gonorrhoeae during neutrophil co-culture

ABSTRACT The ability of bacterial pathogens to metabolically adapt to the environmental conditions of their hosts is critical to both colonization and invasive disease. Infection with Neisseria gonorrhoeae (the gonococcus, Gc) is characterized by the influx of neutrophils [polymorphonuclear leukocytes (PMNs)], which fail to clear the bacteria and make antimicrobial products that can exacerbate tissue damage. The inability of the human host to clear Gc infection is particularly concerning in light of the emergence of strains that are resistant to all clinically recommended antibiotics. Bacterial metabolism represents a promising target for the development of new therapeutics against Gc. Here, we generated a curated genome-scale metabolic network reconstruction (GENRE) of Gc strain FA1090. This GENRE links genetic information to metabolic phenotypes and predicts Gc biomass synthesis and energy consumption. We validated this model with published data and in new results reported here. Contextualization of this model using the transcriptional profile of Gc exposed to PMNs revealed substantial rearrangements of Gc central metabolism and induction of Gc nutrient acquisition strategies for alternate carbon source use. These features enhanced the growth of Gc in the presence of neutrophils. From these results, we conclude that the metabolic interplay between Gc and PMNs helps define infection outcomes. The use of transcriptional profiling and metabolic modeling to reveal new mechanisms by which Gc persists in the presence of PMNs uncovers unique aspects of metabolism in this fastidious bacterium, which could be targeted to block infection and thereby reduce the burden of gonorrhea in the human population. IMPORTANCE The World Health Organization designated Gc as a high-priority pathogen for research and development of new antimicrobials. Bacterial metabolism is a promising target for new antimicrobials, as metabolic enzymes are widely conserved among bacterial strains and are critical for nutrient acquisition and survival within the human host. Here we used genome-scale metabolic modeling to characterize the core metabolic pathways of this fastidious bacterium and to uncover the pathways used by Gc during culture with primary human immune cells. These analyses revealed that Gc relies on different metabolic pathways during co-culture with human neutrophils than in rich media. Conditionally essential genes emerging from these analyses were validated experimentally. These results show that metabolic adaptation in the context of innate immunity is important to Gc pathogenesis. Identifying the metabolic pathways used by Gc during infection can highlight new therapeutic targets for drug-resistant gonorrhea.

fastidious bacterium is a new resource for the Neisseria and microbial metabolic modeling communities. The insights into immune-driven metabolic shifts in Gc revealed by this transcriptionally guided GENRE can inform the future development of therapeutic strategies to combat antibiotic-resistant gonorrhea.

A genome-scale network reconstruction of Neisseria gonorrhoeae metabo lism
We generated iNgo_557, a genome-scale metabolic network reconstruction of Gc strain FA1090, the type strain of Gc which is widely used and highly annotated. A published reconstruction of N. meningitidis M58 (Nmb_iTM560) served as the starting point (23) (Fig. 1A). Nmb_iTM560 was based on the highly annotated iAF1260 recon struction for Escherichia coli and was built using the Biochemical, Genetic and Genomic (BiGG) knowledge base framework (26). We identified homologous genes between N. meningitidis M58 (AE002098.2) and N. gonorrhoeae FA1090 (AE004969.1) using a homology matrix-based workflow for generating high-quality multi-strain genome-scale metabolic models (27). Gc and N. meningitidis were found to share significant homol ogy across large stretches of the genome, particularly for metabolic genes: of the 560 genes, 1,519 reactions, and 1,297 metabolites originally present in Nmb_iTM560, 494 genes, 1,223 reactions, and 1,189 metabolites were preserved in iNgo_557 based on homology (Data Set S1). Orphan reactions from Nmb_iTM560 with no corresponding gene were included in the initial Gc reconstruction and de-orphaned or removed where possible during manual curation. The format was updated to SBML Level 3, the most up-to-date community standard (28). Gene, reaction, and metabolite annotations were updated from KEGG, PATRIC, UniProt, MetaNetX, MetaCyc, PubMLST, and BiGG databases wherever possible (26,(29)(30)(31)(32)(33)(34).
Characteristics of the original Nmb_iTM560 model were conserved, including the presence of a periplasmic compartment, simplified cytochrome respiration pathways, iron acquisition pathways from ferric iron and host proteins, and a biomass equation that reflects the neisserial cell composition, including DNA, RNA, and protein produc tion, growth-associated ATP, and ATP maintenance requirements. Targets for manual curation of the automated reconstruction in iNgo_557 included complete resolution of mass and charge balance inconsistencies, the resolution of import and export loops, removal of carbohydrate import through the phosphotransferase system (which is not functional in pathogenic Neisseria), addition of amino acid catabolism pathways, curation of lipooligosaccharide synthesis for Gc and its addition to the biomass equation, modification of the biomass composition for Gc where appropriate, and simplification of lipid biosynthesis (Data Set S1 and S2). Additionally, in Nmb_iTM560, catalytic cofactors such as biotin, thiamine pyrophosphate, pyridoxal-5-phosphate, iron, zinc, manganese, NAD, and flavin adenine dinucleotide were included as consumed reactants in reactions to reflect biological requirements for biomass production (35). While useful, the presence of these artificially consumed cofactors causes all cofactor-containing reactions in the model to be imbalanced and impedes the accurate assessment of reaction stoichio metries in the reconstruction using automated modeling tools. Furthermore, a recent reassessment of biomass equations in GENREs found that the inclusion of universally essential cofactors in the biomass equation improved gene essentiality predictions (36). These catalytic cofactors were removed in iNgo_557 as reactants from each reaction and added instead to the biomass equation.
This homology-based reconstruction process was incapable of identifying Gc-specific genes that were not present in Nmb_iTM560 (27). Therefore, to expand the metabolic coverage of iNgo_557 for metabolic pathways that are unique to Gc, genes and their corresponding reactions/metabolites were added from an automated reconstruction in the BiGG namespace of Gc FA1090 that was generated using CarveMe (37). Each of the unique genes identified by CarveMe was manually evaluated. Of the 508 genes with metabolic functions predicted by CarveMe, 388 were already present in the model. CarveMe identified an additional 39 genes and corresponding reactions that were supported by manual evaluation of the literature, and they were subsequently included in iNgo_557 ( Fig. 1A; Data Set S1) (37). The remaining 81 genes identified by CarveMe did not have sufficient evidence to support the assigned metabolic function and were not added (Data Set S1).
iNgo_557 is available as an Excel file in Data Set S3 and on Github, including as an SBML file. A comparison of the overall functions captured by iNgo_557 compared to Nmb_iTM560 and CarveMe automated models, as assessed by KEGG reaction catego ries, is presented in Fig. 1B. The overall quality of the reconstruction was assessed using MEMOTE (38). MEMOTE scores serve as a metric for benchmarking consistency in reaction stoichiometries, annotations, and adherence to community standards. Each category in MEMOTE received a score between 82% and 99% except for gene annotation (56.46%), which reflects the limited information in current gene annotation repositories for Gc. The cumulative MEMOTE score of iNgo_557 was 91%, representing an improve ment over the CarveMe automated model (23%) (Fig. 1C).

Validation and utility of predictions in iNgo_557 with experimental pheno types
In silico predictions of biomass flux and amino acid supplementation for iNgo_557 were performed and compared to experimental data to validate the model. First, the compositions of three media used for Gc growth were determined: Gonococcal Base Liquid (GCBL), Morse's Defined Media (MDM), and Roswell Park Memorial Institute 1640 media (RPMI). The metabolites present in each media were assigned to corresponding model exchanges in equivalent amounts and deemed "equally scaled" media (Data Set S4). These simulated media were used to compute biomass flux and consequent predictions of Gc doubling time. Doubling time predictions made with iNgo_557 were then compared to experimental values by conducting growth curves of FA1090 Gc in each of these media ( Fig. 2A and B). The bacterial doubling times predicted by iNgo_557 for equally scaled media were within 13, 15, and 34 minutes of experimentally determined values in GCBL, MDM, and RPMI, respectively ( Fig. 2C and D). All predicted doubling times were faster than what was measured experimentally, which is consistent with the structuring of metabolic network models to predict optimal growth (39).
As shown in Fig. 2A, growth on RPMI was the slowest experimentally, reflecting the limited nutrient content in this media relative to MDM and GCBL (Data Set S4). Specifically, metabolite concentrations in RPMI are ~2-to 10-fold less than the concentrations in MDM and GCBL. For example, glucose is found in MDM and GCBL at 27.8 and 22.2 mM, respectively, but in RPMI at 11.1 mM (Data Set S4). Based on these differences, these three media were "molarity scaled" for simulation in iNgo_557, in which exchanges were set to be equal to the molarity of each respective metabolite in the media to better represent relative availability, as has been done previously (40) (Fig. 2C and D). While using molarity-scaled media for the substrate uptake constraints did not change growth predictions for Gc in MDM and in GCBL, the predicted doubling time of Gc in RPMI was substantially slowed, from 30 to 146 minutes ( Fig. 2C and D).
To identify the substrate(s) that were limiting for Gc growth in RPMI compared with MDM or GCBL, we used iNgo_557 to predict the doubling time of Gc in a revised formulation of RPMI where the constraint for uptake of each component was individually increased to fivefold the original medium (Table S1). In the simulation of growth in equally scaled RPMI, only glucose and serine were individually predicted to increase Gc growth when supplemented at 5× the original constraint ( Fig. 3A; Table S2). In the simulation of molarity-scaled RPMI, serine, asparagine, proline, aspartate, glutamate, and glycine were individually predicted to increase Gc growth when supplemented at 5× the original constraint ( Fig. 3A; Table S2). We tested these predictions experimentally. Addition of 5× glucose, serine, or asparagine to RPMI, which was predicted to increase Gc Research Article mSystems growth the most, significantly increased the growth rate of Gc compared with unmodi fied RPMI (Fig. 3B). Addition of proline, aspartate, glutamate, or glycine, which had the smallest predicted increase in growth, did not significantly enhance growth compared to unmodified RPMI. The experimental growth rate of Gc was not significantly affected when RPMI was supplemented with 5× threonine or valine, which were not predicted to increase growth (Fig. 3B). These findings demonstrate that iNgo_557 can help predict potential nutrients that stimulate Gc growth. Gc is reported to be capable of growing on only three carbon sources: glucose, lactate, and pyruvate (3). While iNgo_557 successfully predicted the growth of Gc on glucose, lactate, and pyruvate as carbon sources in molarity-scaled MDM (27, 43, 28 minutes of doubling times), it also predicted slow growth of Gc in their absence (247 minutes of doubling time). This slow growth is due to the predicted catabolism of amino acids by Gc, for which Gc encodes the required enzymes. However, there was no significant growth of Gc in MDM that did not have glucose, lactate, or pyruvate added (Fig. S3A), even with additional amino acid supplementation (Fig. S3B). Metabolomics and C 13 metabolic flux analysis recently conducted on Gc strain MS11 demonstrated that Gc consumes a variety of amino acids when grown in a glucose-containing medium that is a derivative of MDM (41). To further explore Gc usage of amino acids, we compared the simulated growth of Gc in molarity-scaled MDM (Fig. S4) to the C 13 distributions reported in reference (41). Our in silico predictions of metabolic flux were consistent with the reported C 13 distributions: Gc exhibited a bipartite metabolism in which carbon from glucose was metabolized to acetyl-CoA, and amino acids were independently metabolized via the tricarboxylic acid cycle (TCA) cycle. Thus, a prediction made by our model, that Gc can catabolize amino acids, is supported by experimental evidence. To account for Gc usage of amino acids in the presence of its known carbon sources, amino acid catabolic pathways were left intact in iNgo_557.
We compared gene essentiality predictions yielded by iNgo_557 on GCB with a published data set comprising essential genes, which were identified through the growth of strain MS11 transposon insertion mutants (transposon sequencing, TnSeq) on GC agar (Data Set S5) (42). In iNgo_557, a gene was predicted to be essential if less than 10% of the optimal biomass of the wild type (WT) could be produced by a mutant in single-gene deletion simulations. Predicted gene essentiality had an agreement of 73% with a Mathews correlation coefficient (MCC) of 0.43 (Data Set S5). Genes identified as essential by both iNgo_557 and TnSeq included those related to lipooligosaccharide Growth was monitored by enumeration of CFU/mL, and doubling time was calculated using GrowthCurver. Results are from n = 3 biological replicates. Bars represent the mean. Error bars represent SEM. Dotted line indicates doubling time in unmodified RPMI (black bar). Metabolites predicted to increase growth are in gray bars; control metabolites predicted to not increase growth are in hatched bars. *P < 0.05 by one-tailed t-test is relative to unmodified RPMI.
Research Article mSystems and peptidoglycan biosynthesis, purine metabolism, and pyruvate metabolism. Of the genes identified as non-essential by iNgo_557 but essential by TnSeq, many encoded participants in pyrimidine metabolism, oxidative phosphorylation, and glycolysis. One such gene, encoding pyruvate kinase (pyk), was predicted by iNgo_557 as non-essential, but essential by TnSeq. We verified that pyk could be deleted from Gc and that the resulting null mutant could grow in GCBL containing glucose as the sole carbon source, albeit slower than the WT parent or when pyruvate was provided (Fig. S5). While Tnseq has been used for model validation in the past, a major caveat of using it is that genes determined to be essential by TnSeq may only appear so in a competitive setting when mixed with a library of other transposon mutants (43,44). Differences in media composition and bacterial strain backgrounds can also contribute. For these reasons, we did not use TnSeq gene essentiality data for further curation of iNgo_557. Taken together, there was good concordance between Gc growth predicted using iNgo_557 and experimental results generated in this study and from previously published work. We conclude that iNgo_557 can be used to predict growth dynamics of Gc in various in vitro conditions, which can be further validated experimentally.

Transcriptome-guided modeling of Gc metabolism during co-culture with primary human neutrophils predicts a shift in the pyruvate axis
GENREs serve as a tool for scaffolding complex metabolic information in human-inter pretable formats. One such application is the integration of transcriptional data with GENREs to develop a comprehensive picture of bacterial metabolism in complex and uncharacterized environments (45). Given that Gc is a human-specific pathogen, we sought to use the reconstruction to predict metabolic phenotypes that are consistent with Gc growth in the context of human neutrophils (PMNs), the predominant immune cell that is recruited during infection. To investigate how Gc metabolism shifts in response to co-culture with PMNs, transcriptomic data from Gc co-incubated with PMNs for 1 hour were integrated with iNgo_557 to generate contextualized models that offer insight into the metabolic state of Gc during infection.
To accomplish this goal, we applied the RIPTiDe (Reaction Inclusion by Parsimony and Transcript Distribution) algorithm, which uses RNA sequencing (RNA-seq) data to identify the most cost-effective usage of metabolism while also reflecting the organism's transcriptional investment. In brief, the model is iteratively constrained to a fraction of the maximal possible biomass. The algorithm then assigns linear coefficients to weight model reactions based on transcript abundance. Reactions corresponding to highly abundant transcripts are favored to be retained, and low abundance transcripts are favored to be removed (i.e., "pruned") from the model. Consistency of the model with the transcriptome is assessed by Spearman correlation coefficient. The fraction with the highest correlation is selected as the context-specific model (Fig. S6). RIPTiDe has been used successfully with models of P. aeruginosa and C. difficile to uncover metabolic contributors to virulence in the context of mucin degradation, biofilm formation, murine infection models, and co-culture with other microbes (17, 19,46). We reasoned this approach would generate context-specific models of the metabolism of Gc when grown with and without PMN co-culture and would identify those reactions that are likely to be differentially active in each condition.
The transcriptome data set we used was from RNA-seq of a constitutively opacity protein (Opa)-deficient isolate of strain FA1090 Gc, which was cultured in RPMI + 10% fetal bovine serum (FBS) for 1 hour. Gc was cultured in the presence or absence of primary human PMNs that were adherent and treated with the chemokine interleukin-8 to reflect the activated state of immune cells during infection (47). This Opa-negative isolate of Gc survives significantly better when exposed to PMNs than Opa-expressing strains that are rapidly internalized and killed by PMNs (48). RIPTiDe generated two context-specific models of Gc metabolism: one for Gc in a medium without PMNs and one for Gc with PMNs (Fig. S6). For each of the two models, flux samples were generated to assess all possible metabolic profiles in the two environmental contexts.
Flux samples generated with the models significantly correlated with the transcript abundances derived from RNA-seq for each condition (r = 0.242, P < 0.001 for Gc without PMNs, and r = 0.263, P < 0.001 for Gc with PMNs) (Fig. S6). These r values are consistent with those obtained from models of other organisms to which RIPTiDe was applied (17), indicating that the context-specific metabolic profiles predicted with RIPTiDe align with experimental data.
Biomass flux was significantly increased in the contextualized model of Gc co-cul tured with PMNs, compared with Gc cultured without PMNs (Fig. 4A), suggesting an overall stimulation of Gc metabolism in the presence of PMNs. Flux distributions for each model were then compared using non-metric multidimensional scaling (NMDS) of consensus reactions shared between both models to broadly identify metabolic growth patterns used by the context-specific models. NMDS revealed that the sampled flux distribution for Gc co-cultured with PMNs overlapped with, but was distinct from, the sampled flux distribution for Gc cultured without PMNs (Fig. 4B). This result reflects that the media used for growth is consistent between the two models, but co-culture with PMNs caused a shift in metabolic pathways used for growth.
We further analyzed the contextualized models to better understand the shifts in metabolism that resulted in the distinctions observed in the NMDS analysis. Reactions unique to each model (non-consensus reactions) were identified, and the absolute median activity for each reaction was determined for 1,000 simulations to examine the contribution of each reaction to biomass production (Fig. 4C). From this analysis, we identified a set of 19 reactions that were unique to Gc co-cultured with PMNs and eight reactions unique to Gc cultured without PMNs. Several reactions involved in metabolite import and catabolism were unique to Gc co-cultured with PMNs, suggesting that there are changes to the metabolites available to Gc in this condition, possibly due to the competition with or excretion by PMNs. Specifically, pyruvate and D-lactate exchange reactions were unique to Gc co-cultured with PMNs (Fig. 4C). Further, the model predicts that Gc imports these metabolites (Fig. S7), suggesting bacterial use of these alternative carbon sources in the presence of PMNs. This observation aligns with extensive evidence that PMNs secrete lactate as a byproduct of oxidative metabolism, which stimulates Gc growth (4,51). Similarly, Gc co-cultured with PMNs was also predicted to uniquely carry flux through nitrogen metabolism, in particular the import of nitrite, nitrite oxidoreduc tase, nitric oxide reductase, and nitrous oxide export ( Fig. 4C; Fig. S7). These findings align with the reported production of nitric oxide (NO) via inducible nitric oxide synthase in stimulated PMNs (52). Although NO is used by phagocytes to directly kill pathogens, Gc can exploit this aspect of inflammation by detoxifying NO to nitrous oxide or using nitrite and nitric oxide as terminal electron acceptors during anerobic growth (53). Together, these observations support the hypothesis that neutrophil byproducts mediate remodeling of Gc metabolism.
We next assessed reactions that were shared between both models of Gc cultured without PMNs and Gc co-cultured with PMNs but carried different levels of flux. From this, we identified reactions that most strongly discriminated between the metabolic activity of the two models. This analysis employed a supervised machine learning approach with Random Forest, a categorization algorithm that can segregate flux samples based on the contextualized models (Fig. 4D) (50). In brief, we trained a Random Forest classifier to predict reactions that were critical in defining whether the contex tualized model was associated with culturing with or without PMNs. Specifically, the Random Forest input was reaction flux distributions for consensus reactions between the contextualized models from the flux samples (n = 500). The classifier had an out-of-bag accuracy above 98%, indicating that membership for each contextualized model can robustly be predicted from reaction content within the two contextualized models. We then assessed mean decrease accuracy (MDA) to identify reactions that, when removed from the model, most affected the categorization predictions of the Random Forest. Gc grown in the presence and absence of PMNs was particularly distinguished by flux out of the pyruvate node, through acetate synthesis. Specifically, acetate exchange, acetate transport, acetate kinase, and acetate phosphotransacetylase were identified as reactions that impacted the categorization capabilities of the Random Forest (MDA ~13%) (Fig. 4D). Acetate production is a prominent feature of bacterial Research Article mSystems overflow metabolism, in which ATP is generated from the production of acetate from acetyl-CoA via the phosphate acetyl-transferase and acetate kinase (PTA-AckA) pathway rather than shuttled into carbon backbones for biomass (54). Visualization of flux balance analysis ( Fig. 5A and B) demonstrated a predicted increase in acetate flux in co-culture with PMNs, consistent with increased carbon flux from the addition of alternative carbon sources, such as lactate and pyruvate.
Conditionally essential genes were predicted by conducting essential gene calcula tions in each model, then comparing them (Table 1; Data Set S7). Twelve genes were predicted to be essential only when Gc was cultured without PMNs, and two genes were predicted to be essential only when Gc was co-cultured with PMNs. Of the 12 genes predicted to be essential only when Gc was cultured without PMNs, seven are within a single pathway exiting the pyruvate synthesis node (Table 1): pyruvate kinase (pyk), portions of the pyruvate dehydrogenase complex (dldH), phosphate acetyltransferase (pta), citrate synthase (gltA), aconitase (acnB), and isocitrate dehydrogenase (idh) were all predicted to be essential only for Gc cultured without PMNs.
To test the prediction that pyruvate synthesis genes were essential for Gc in rich growth medium but dispensable for Gc in the presence of PMNs, we generated a null mutant in pyruvate kinase (Δpyk), the first enzyme in this pathway. As expected, Δpyk had a growth defect in MDM, GCBL, and RPMI containing glucose as the sole carbon source, while the WT parent grew in these media ( Fig. 6A; Fig. S5A and S8). Also as  6B through E). We then measured the growth of WT and Δpyk Gc in the conditions used to collect the PMN transcriptomics data. In RPMI + 10% FBS, the Δpyk mutant stopped growing after 3 hours, and by 24 hours its viability had declined to <1% of the inoculum (Fig. 6F, H, and I). In contrast, when cultured in the presence of PMNs, Δpyk Gc grew significantly better than Gc in the absence of PMNs, and in fact increased in viability over 24 hours (Fig. 6G through I). WT Gc grew over this time whether or not PMNs were present (Fig. 6F through I). These results suggest that Gc co-cultured with PMNs has a decreased need for flux through glycolysis and instead imply that Gc has access to alternative carbon sources such as lactate and pyruvate, which support its growth in the presence of PMNs independently of the glycolytic pathway.

DISCUSSION
Over the last 20 years, genome-scale metabolic modeling has become a powerful tool for context-specific interrogation of complex biological networks. In this study, we developed a highly curated genome-scale metabolic network reconstruction, titled iNgo_557, for Gc strain FA1090. This model predicts the use of glucose, lactate, and pyruvate as carbon sources for Gc, and an increase in growth when selected amino acids are supplemented in a cell culture medium containing one of these carbon sources (3). iNgo_557 was contextualized using transcriptomics data that we recently generated to identify shifts in Gc metabolism that occur in response to co-culture with PMNs. These results represent the first use of genome-scale metabolic modeling in Gc for the discovery of metabolic contributors to virulence. Through the linkage of gene, reaction, and metabolite information, iNgo_557 facilitates rapid and convenient manipulation of metabolic parameters to identify contributors toward Gc pathogenesis that are otherwise complicated, time-consuming, or laborious to replicate in vitro. Independently, GENREs can be used to simulate well-defined environmental contexts, such as growth in laboratory media. iNgo_557 accurately reflects experimental growth phenotypes in simulated media and can be used to predict Gc growth phenotypes following distinct manipulations to these media. We demonstrated one such use: identification of growth-limiting nutrients in RPMI. Other applications include nutrient drop-out experiments, aerobic and anaerobic growth, and gene essentiality studies. While the predictions generated by our model were consistent with experimental results, incorrect predictions are also informative, revealing points of obscurity in our understanding of Gc metabolism. For example, iNgo_557 predicted the growth of Gc in MDM in the absence of a dedicated carbon source (glucose, lactate, pyruvate). On further interrogation, the predicted growth of Gc on MDM without a carbon source was due to the consumption of serine and alanine as carbon sources. Although Gc encodes the genes necessary to catabolize these amino acids (ALATA_L/ NGO_1047 and SERD_L/NGO_1773 and NGO_0444) and does catabolize these amino acids in vitro (41), it is unable to use amino acids as a sole carbon source (Fig. S3). A potential explanation for the inability of Gc to use amino acids as a sole carbon source includes Error bars represent SEM. n = 3-4 biological replicates. Significance was determined by two-way ANOVA with Holm-Sidak correction for multiple comparisons, *P < 0.05.

Research Article mSystems
transcriptional regulation, which prevents the use of serine as a sole carbon source in P. aeruginosa (55). Post-transcriptional regulation including allosteric modulation of metabolic enzymes (56) and substrate competition for available enzymes or transport ers (57) can also prevent the catabolism of amino acids. Our results suggest that a similar form of transcriptional or post-transcriptional regulation may also dictate Gc carbon source utilization. Regulatory features that can impact carbon source usage are missing from GENREs, an inherent limitation in genome-scale metabolic modeling. These discrepancies serve as points for further investigation and facilitate hypothesis generation.
Incorporation of additional layers of regulatory information can improve model accuracy, particularly for the modeling of complex environments such as during co-culture with other species or cell types, which is impeded by the lack of knowl edge of the metabolite environment. As a human-adapted mucosal pathogen, Gc must co-exist with a complex assortment of human microbiota, epithelial cells, and mucosal immune cells. The recruitment of PMNs and the inflammation associated with gonococcal infection further complicate an already complex metabolic environment. Transcriptomic integration with metabolic models serves to deconvolute the modeling of these complex settings through unsupervised contextualization of GENREs for a specific environment. As such, we leveraged RIPTiDe with iNgo_557 to better understand the metabolic pathways enabling Gc growth during co-culture with PMNs and to predict the behaviors of this host-associated bacterial species. Intriguingly, several genes are predicted to be essential only when Gc is cultured without PMNs, but not in the context of PMNs, including pyruvate kinase (pyk) and pathways exiting the pyruvate synthesis node. Metabolic modeling using iNgo_557 predicts that this effect is due to the bypass of pyruvate synthesis through the import of alternative carbon sources, including lactate and pyruvate, when in the presence of PMNs, which is supported by our growth data. PMNs are highly glycolytic cells that consume glucose and secrete lactate when they are activated (51). The ability of Gc to use lactate would relieve competition between Gc and PMNs for glucose during infection and enable Gc to exploit this byproduct of PMN glycolysis. In support of this possibility, lactate use has been found to be required for Gc survival from PMNs, within cervical epithelial cells, and in the female mouse genital tract (4,5,58). The increase in biomass flux predicted for models of Gc cultured with PMNs compared to Gc cultured without PMNs is further consistent with reports that Gc growth on lactate stimulates Gc metabolism (51,58). Together these results provide evidence that Gc utilizes additional alternative carbon sources, such as lactate and pyruvate, when co-cultured with PMNs to enhance its growth. Future metabolomics studies can measure concentrations of carbon sources and other nutrients at the mucosal sites Gc inhabits to better predict bacterial growth during human infection.
Regardless of the source, carbon exiting the pyruvate synthesis node can proceed in one of the two pathways in Gc: acetate production or oxidation through the TCA cycle. Acetate production through the PTA-AckA pathway is a prominent feature of Gc growth on glucose, lactate, and pyruvate (54,59). Downstream of pyruvate kinase, iNgo_557 predicted increases in Gc acetate production when in the presence of PMNs. In N. meningitidis, acetate is secreted following growth on glucose, lactate, and pyruvate, and the highest activity of the PTA-AckA pathway occurs when all three carbon sources are present, compared with glucose alone (59). Our results are consistent with this obser vation. Alternatively, glucose, lactate, and pyruvate can instead be further catabolized by the TCA cycle. In N. meningitidis, pyruvate dehydrogenase (dldH), citrate synthase (gltA), aconitase (acnB), and isocitrate dehydrogenase (idh) reaction activities were all demonstrated to be high in the presence of glucose but decreased in the presence of pyruvate (59). Consistent with the stimulation of these enzymes in the presence of glucose compared to pyruvate, iNgo_557 predicted dldH, gltA, acnB, and idh to be essential only in the absence of PMNs, in which glucose is the sole carbon source available. The alleviation of the requirement for acnB in the context of PMN co-culture is notable in light of a recent study that identified compensatory mutations within acnB that enabled the recovery of antibiotic-resistant penA mutant Gc from the mouse genital tract (60). Together our results highlight the pyruvate node as a critical pivot point in Gc metabolism, particularly in the context of an inflammatory environment created by PMNs. Overall, the predictions generated here by contextualized models of iNgo_557 reveal new insights into Gc pathogenesis, highlighting it as a viable platform for the discovery of metabolic pathways associated with virulence and antibiotic resistance.
Treatment options for Gc have become increasingly limited over the last two decades, and only a single recommended antibiotic remains for the treatment of gonorrhea (61). The development of new potential therapies is essential to avoid the threat of completely antibiotic-resistant Gc. Targeting essential bacterial metabolic pathways during infection represents a promising approach, one that was first shown decades ago in the context of sulfonamide antibiotics, which directly inhibit folate synthesis (62). Novel approaches for the treatment of antibiotic-resistant infections have included the application of metabolites to shift the metabolism of pathogens toward a less favorable state (63). There is a need for a revisitation of Gc metabolism and physiology in light of the approaching post-antibiotic era for gonorrhea (64). Technologies such as RNA-seq, forward and reverse genetic screens, and metabolic modeling can all provide insights into Gc metabolism. Here, the integration of transcriptomics with genome-scale metabolic modeling is synergistic, providing more insight into the remodeling of Gc metabolism in the context of PMN co-culture than could be discerned from each technique alone. In sum, this study highlights the opportunities afforded by genomescale metabolic modeling for the targeted identification of context-specific essential metabolic pathways that enable Gc to thrive within the human host, with further predictions and discoveries remaining to be made.

Genome-scale metabolic reconstruction
To generate a GENRE for Gc, we used N. meningitidis M58 Nmb_iTM560 as an initial template for the automated multi-strain model, reconstruction pipeline (27). In brief, the pipeline used bidirectional best hit BLAST to identify genes with >80% homology between N. meningitidis M58 (AE002098.2) and N. gonorrhoeae FA1090 (AE004969.1) to generate a homology matrix for the two species. A secondary comparison using BLAST on nucleotide sequences was conducted to identify potential homologs with poor open reading frame (ORF) annotation. These automated calls were inspected and reassessed for each gene present in Nmb_iTM560 as indicated in Data Set S1. Using the homol ogy matrix, a draft strain-specific model was generated using COBRApy (65). Metabolic genes (and the corresponding reactions and metabolites) specific to Gc FA1090 were added to the reconstruction using CarveMe when supported by literature evidence (37). Exchange reactions that were missing for extracellular metabolites in the reconstruction were added. The model was then further manually curated to de-orphan reactions and incorporate published metabolic functions for Gc according to literature evidence where possible (Data Set S1). Final gene and reaction calls, along with decision annotations, can be found in Data Set S1. Annotation data were automatically assigned using ModelPo lisher (66). Reaction and stoichiometric inconsistencies were corrected for each reaction. All formulas were mass and charge balanced using the BiGG database, when possible, to maintain a consistent namespace (26). A list of mass and charge imbalanced reactions and their corrections are provided in Data Set S1. Additional annotations were collected and added to the annotation field dictionary for all model components from KEGG, PATRIC, UniProt, MetaNetX, MetaCyc, PubMLST, or BiGG databases (26,(29)(30)(31)(32)(33)67). The pipeline for development of the reconstruction is available in the GitHub repository associated with this study (https://github.com/aimeepotter/Gc_GENRE_2022).

Assessing reconstruction quality
Modeling assessments, including flux balance analysis, flux variability analysis, singlegene knockout analysis, were conducted using COBRApy (65). Model quality was assessed with MEMOTE using a local installation v0.13.0 (38). Gene essentiality predic tions were compared to a published data set of essential genes for growth on solid, rich media for Gc strain MS11, which was aligned to Gc FA1090 by bidirectional best-hit BLAST as above. Prediction accuracy was calculated as the number of correct predictions divided by the number of total predictions for genes present in both data sets and the MCC was calculated as in reference (68).
A protocol for defining realistic modeling constraints for in silico media was recently described, in which metabolite exchanges are scaled based on the maximum possi ble usage defined by the concentration of metabolites in millimole per liter (40). We therefore generated two in silico exchange reaction constraints for each simulated media: equally scaled, to avoid constraining the model with incorrect assumptions, and molarity-scaled, to match the maximum possible use of metabolites. The concentration of metabolites present in each media and their corresponding assignments to in silico media constraints are detailed in Data Set S4. Biomass flux and subsequent doubling times for simulated growth in GCBL, MDM, and RPMI were compared to experimental values. Predictions of Gc doubling time were calculated assuming a biomass equation scaled to 1 g dry weight of bacteria based on the following formula: Doubling Time = ln(2) * 60/(objective value) Experimental doubling times were determined using GrowthCurver implemented in R for both optical density (OD) and CFU per milliliter, with stationary phase values trimmed (69).

RIPTiDe contextualization and analysis
Transcriptomic data retrieved from the GEO database (GSE123434) for Gc cultured without and with PMNs over the course of 1 hour were mapped to the correspond ing FA1090 gene IDs using the conversion table provided in Data Set S1 of reference (47). For RIPTiDe contextualization, an unsupervised approach was used in which all exchange reaction bounds were set to ±10, except oxygen, which was set at ±20. The transcriptomic data were then integrated with the model using RIPTiDe using the maxfit_contextualize() function (minimum fraction 0.1, maximum fraction 0.8, n = 1,000) to produce contextualized models for Gc grown in the presence or absence of PMNs (70). Flux samples were gathered from consensus reactions between both contextualized models (n = 500 samples per model). Bray-Curtis-based NMDS (k = 4, trymax = 25) and permutational multivariate analysis of variance (PERMANOVA) (perm = 999) analyses were accomplished using the Vegan R package (71). Supervised machine learning was accomplished with the implementation of AUC-Random Forest also in R on the sampled flux distributions for shared reactions to determine metabolic functions that distinguish between the two models (ntree = 1,500, importance = TRUE, err.rate = TRUE, mtry = 20), and importance was reported as mean decrease in accuracy for the top 20 predictors (49).
Statistical analysis was performed in R v4.1.0. Visualizations of flux balance analysis were performed using Escher (72).
WT Gc were grown on Gonococcal Medium Base (GCB, Difco) plus Kellogg's supple ments at 37°C with 5% CO 2 (76,77). Δpyk strains were grown on GCB plus Kellogg's supplements with glucose replaced by pyruvate (36 mM) as in reference (78). For preparation of mid-logarithmic phase bacteria, Gc were grown in liquid medium (GCBL) or carbon-matched GCBL containing pyruvate (45 mM) as the sole carbon source, where appropriate, for successive rounds of dilution, and enriched for piliation, as previously described (79). Spectinomycin was used for the selection of the pyk mutation at 80 µg/mL.

Growth curves
Gc in mid-logarithmic phase were pelleted, resuspended in the indicated media, and diluted to ~5*10 7 CFU/mL in 6 mL of media in 15-mL conical tubes (Sarstedt). The bacterial suspension was incubated with rotation at 37°C. Bacterial growth was measured by OD 550 and CFU enumeration at specific time points. CFUs are presented relative to 0 hour (100%). Gc was grown in GCBL, HyClone RPMI 1640 media without glutamine (Catalog no. SH30096.FS) (Cytivia), or carbon-matched Morse's defined media (MDM) containing either glucose (27 mM), lactate (54 mM), or pyruvate (54 mM) (80). Doubling times were calculated from best-fit logistic curves generated with Growth Curver (69) for the lag and exponential phases of each growth curve for at least three experimental replicates and averaged. Significant differences for growth over time were determined by one-tailed t-test in GraphPad Prism v9.

Gc-PMN co-culture
PMNs were isolated from venous blood as previously described and used within 2 hours of isolation (79). Subjects gave informed consent in accordance with an approved protocol by the University of Virginia Institutional Review Board for Health Sciences Research (no. 13909). Synchronized Gc infection of PMNs in suspension was conducted as previously described (81). PMNs were resuspended in RPMI (Cytivia) containing 10% heat-inactivated FBS (Gibco) at 1 × 10 6 PMN/mL, and Gc was added to each tube at a multiplicity of infection of 10. One millimolar of IPTG was added where indicated to induce the expression of pyk. Six milliliters of the suspension was incubated in 15-mL conical tubes with rotation at 37°C. Bacterial CFUs were enumerated at specified time points and expressed relative to the CFU at 0 hour (100%). Data are expressed as the mean ± SEM of at least three replicate experiments. Significant differences were determined by two-way analysis of variance (ANOVA) with Holm-Sidak correction for multiple comparisons in GraphPad Prism v9.