key: cord-0031597-2mqrcez5
authors: Ambikan, Anoop T; Svensson-Akusjärvi, Sara; Krishnan, Shuba; Sperk, Maike; Nowak, Piotr; Vesterbacka, Jan; Sönnerborg, Anders; Benfeitas, Rui; Neogi, Ujjwal
title: Genome-scale metabolic models for natural and long-term drug-induced viral control in HIV infection
date: 2022-05-10
journal: Life Sci Alliance
DOI: 10.26508/lsa.202201405
sha: 4230dc25187f012807b5b41b0fd08ba68212fd6f
doc_id: 31597
cord_uid: 2mqrcez5

Genome-scale metabolic models (GSMMs) can provide novel insights into metabolic reprogramming during disease progression and therapeutic interventions. We developed a context-specific system-level GSMM of people living with HIV (PLWH) using global RNA sequencing data from PBMCs with suppressive viremia either by natural (elite controllers, PLWH(EC)) or drug-induced (PLWH(ART)) control. This GSMM was compared with HIV-negative controls (HC) to provide a comprehensive systems-level metabo-transcriptomic characterization. Transcriptomic analysis identified up-regulation of oxidative phosphorylation as a characteristic of PLWH(ART), differentiating them from PLWH(EC) with dysregulated complexes I, III, and IV. The flux balance analysis identified altered flux in several intermediates of glycolysis including pyruvate, α-ketoglutarate, and glutamate, among others, in PLWH(ART). The in vitro pharmacological inhibition of OXPHOS complexes in a latent lymphocytic cell model (J-Lat 10.6) suggested a role for complex IV in latency reversal and immunosenescence. Furthermore, inhibition of complexes I/III/IV induced apoptosis, collectively indicating their contribution to reservoir dynamics.

During HIV-1 infection, cellular metabolic activity, combined with glycolytic enzymes, regulates susceptibility to HIV-1 at the cellular level (Clerc et al, 2019; Valle-Casuso et al, 2019) . Elevated oxidative phosphorylation (OXPHOS) and glycolysis thus favor infection in CD4 + T cells (Clerc et al, 2019; Valle-Casuso et al, 2019) . CD4 + T cells up-regulate glycolysis to meet the energy-demanding turnover for HIV-1 virion production, resulting in their eventual death (Hegedus et al, 2014; Palmer et al, 2014) . After initiation of combination antiretroviral therapy (cART), virus-induced short-term metabolic changes do not restore the transient metabolic modulation caused by the infection. Altered amino acid (AA) metabolism has been reported in untreated people living with HIV-1 (PLWH) as well as within the first 2 yr after initiation of cART compared with the HIVnegative controls (Cassol et al, 2013; Peltenburg et al, 2018) . In our recent extensive multi-omics system biology studies on cohorts from India (Babu et al, 2019; Gelpi et al, 2021) , Cameroon (Gelpi et al, 2021) , and Denmark (Gelpi et al, 2021; Villumsen et al, 2022) , we mapped the in-depth metabolomic dysregulation associated with long-term treatment in PLWH. Our group (Babu et al, 2019; Gelpi et al, 2021; Villumsen et al, 2022) , and others (Mukerji et al, 2016; Rosado-Sánchez et al, 2019; Valle-Casuso et al, 2019; Meeder et al, 2021; Shytaj et al, 2021) , have highlighted how the coordinated modulation of central carbon metabolism, AA metabolism, glutaminolysis, and fatty acid biosynthesis can potentiate accentuated immune aging and cognitive decline in a subset of PLWH on therapy who have dysregulated metabolic profile.

Elite controllers (EC) are a unique group of PLWH that naturally control viral replication and maintain a low viral reservoir. Our recent study indicated that EC had a distinct lipid profile, reduced inflammation, and increased antioxidant defense which may contribute to the EC status . Moreover, the integrative proteomic and transcriptomic analysis suggested that the EC group had a unique metabolic uptake and flux profile through hypoxia-inducible factor signaling and glycolysis, which could contribute to the natural control of HIV-1 infection . A study also showed how suboptimal inhibition of glycolysis in CD4 + T cells decreased the latently infected reservoir (Valle-Casuso et al, 2019) . However, EC is heterogeneous, and one mechanism of elite control does not exist (Zhang et al, 2018; Akusjärvi et al, 2022) . Instead, PLWH on long-term successful therapy with prolonged suppressive viremia are more homogenous in their immune profile (Zhang et al, 2017) . A deep understanding of the immune profile of these groups of HIV-1-infected individuals could help to develop strategies for analytical treatment interventions (ATI) to achieve a clinically relevant ART-free HIV cure or remission (Julg et al, 2019) .

Genome-scale metabolic models (GSMMs) can provide novel insights toward understanding host-pathogen interactions and metabolic reprogramming during acute and chronic infections. When applied to PBMCs, GSMM can contribute to unraveling the mechanistic processes at the systems level . By combining contextualized PBMC-specific biological network analysis, GSMMs, and multi-omics integration, one can attain holistic and dynamic characterizations of complex rearrangements during disease progression or therapeutic interventions Zeybel et al, 2021) .

In the present study, we sought to understand and infer changes in HIV-1 infection at the system level by comparing successfully treated PLWH on prolonged cART (herein PLWH ART ) with the HIV-seropositive ECs (herein PLWH EC ) and HIV-negative controls (herein HC). Contextualized PBMC-specific GSMMs and biological networks were thus developed for PLWH ART and PLWH EC to identify the metabolic alterations during prolonged therapy. We further modulated the key pathways pharmacologically to determine their role in HIV-1 reservoir dynamics and immune senescence profile. By combining the multidimensional omics data, our study is the first to provide a comprehensive mapping of the immunometabolic dysregulations using GSMM in PLWH ART with successful long-term treatment. Furthermore, our comparative analysis with PLWH EC offers mechanistic insights into natural immune control.

The study population included three PLWH cohorts, where two groups had suppressed viremia (PLWH ART and PLWH EC , n = 19 each), and one untreated group was viremic (herein PLWH VP , n = 19). In addition, we included 19 HC. Given that extensive transient metabolic changes occur in the PLWH VP due to the acute viremic phase, we used this group to develop the cART-specific model only (see the Materials and Methods section). The clinical characteristics are given in Table  S1 . The median (IQR) duration of diagnosed HIV-1 seropositivity infection for PLWH EC was 9 (5-14) yr, and none had received treatment. In PLWH ART , the median duration of suppressive therapy was 13 (7-17) yr with no viral blips except for two individuals. At the time of sample collection, both PLWH ART and PLWH EC had undetectable plasma viral load (<20 copies/ml) and CD4 + T-cell counts (>500 cells/μl) indicative of significant immune reconstitution.

System-level PBMC-based gene expression identifies dysregulation of OXPHOS in PLWH ART To identify the system-level host response during HIV-1 infection, we performed transcriptomic profiling of total RNA isolated from PBMCs. The differential gene expression analysis was performed between all pair-wise comparisons among the four groups (adjusted P < 0.05, Supplemental Data 1). No genes were found to be dysregulated between PLWH EC and HC, whereas 949 genes were differentially expressed between PLWH ART and HC (adjusted P < 0.05). To identify whether the changes in gene expression between the groups were due to the altered cell-type proportions, we performed digital cell quantification (DCQ), estimating cell-type proportions in each group (Racle & Gfeller, 2020) . We characterized 18 immune cell populations (Figs 1A and S1). As expected, the proportion of several cell types was significantly different in PLWH VP compared with the other groups. No significant difference in the proportion of cell types was observed between PLWH EC and HC and the only difference in regulatory T cells (Tregs) was observed between PLWH EC and PLWH ART (P < 0.05). Based on the differentially expressed genes, we identified 1,037 specifically dysregulated genes in PLWH ART (see the Materials and Methods section) with an explicitly differential regulation in PLWH ART (Supplemental Data 2). Sample clustering using the cART-specific genes separated PLWH ART samples from the other groups ( Fig 1B) . One PLWH ART sample (marked by an arrow) (Fig 1B) was identified as belonging to a patient who had been classified in the past as an EC but started treatment 23 yr after HIV-1 diagnosis (date of diagnosis: 05 January, 1989, treatment initiation 11 December, 2012) due to two successive viral loads were above the detection limit (240 and 185 copies/ml, respectively). The patient maintained viral load below detection level following treatment. Gene set enrichment analysis (GSEA) using MSigDB hallmark gene sets on the ART-specific genes highlighted that the primary mechanisms related to the long-term treatment were OXPHOS (adjusted P < 0.05) and reactive oxygen species (ROS) pathways (adjusted P < 0.1), as the top significantly regulated gene sets ( Fig 1C) .

Larger viral reservoir and up-regulated OXPHOS differentiate PLWH ART from PLWH EC Next, we performed a comparative analysis between PLWH EC and PLWH ART to identify the immune signature during suppressive viremia that is naturally controlled, or cART induced. First, we performed relative reservoir quantification on total PBMC HIV-1 DNA and observed that PLWH ART had a significantly larger reservoir than PLWH EC (Fig 2A) . Furthermore, we performed differential gene expression analysis between PLWH EC and PLWH ART to identify the cART related changes during suppressive viremia. We identified 1,061 significantly dysregulated genes in PLWH ART compared with PLWH EC (adjusted P < 0.05). There were 400 genes up-regulated and 661 genes down-regulated in PLWH ART compared with PLWH EC (Fig 2B) . The dysregulated genes displayed distinct expression patterns in the two groups and hierarchical clustering, showing apparent clustering of PLWH ART and PLWH EC samples ( Fig 2C) . No other factors like age, duration of treatment, and gender showed any clustering pattern. The GSEA analysis using the MSigDB hallmark gene set showed that OXPHOS and ROS pathways were significantly enriched with most of the genes up-regulated in PLWH ART (Fig 2D) (false discovery rate [FDR] < 0.2). Pathways with most genes downregulated in PLWH ART were not statistically significant. Pathways such as mTORC1 signaling and glycolysis also appeared in the analysis, with most of the genes up-regulated in PLWH ART but without passing the significance threshold (FDR > 0.2). OXPHOS was identified as significantly altered in long-term treated patients. Therefore, we looked at OXPHOS in detail to find which complexes were most affected. Among the genes in the five complexes of OXPHOS (I to V), complexes I (34%), III (45%), and IV (45%) were primarily affected in PLWH ART compared with PLWH EC (Fig 2E) . We also checked the overlap between the ART-specific genes (n = 1,037) and the differentially regulated genes between the PLWH ART and PLWH EC (n = 1,061). We observed that 557 genes were overlapping between the two sets of genes. The gene list enrichment analysis using MSigDB hallmark gene set identified OXPHOS (adjusted P < 0.001), MYC targets V1 (adjusted P = 0.004), and ROS pathway (adjusted P = 0.04) as significant pathways. Combining all the data, up-regulation of the OXPHOS was the hallmark of PLWH ART and complexes I, III, and IV were primarily affected.

Altered flux balance in PLWH on cART is linked to OXPHOS, glycolysis, and TCA cycle Given that significant metabolic pathway-centered perturbations were found in PLWH ART , we next performed reporter metabolite analysis to identify metabolites around which most of the transcriptional changes occurred. Five metabolites, namely, superoxide, ubiquinol, ubiquinone, ferrocytochrome C, and ferricytochrome C were significantly up-regulated in PLWH ART compared with PLWH EC (adjusted P < 0.2) ( Fig 3A) . In addition, nicotinamide adenine Nodes are genes and edges represent association with pathways. Node size is relative to the mean expression of the genes among the PLWH ART . Genes overlapping between pathways and high abundance genes are labeled. dinucleotide hydrogen, S-adenosyl methionine, and S-adenosyl-Lhomocysteine were predicted to be significantly dysregulated in PLWH ART (adjusted P < 0.2) ( Fig 3A) . The overall results suggested a significant change in porphyrin, glycine, serine, and threonine metabolism, and a positive regulation in OXPHOS. The reactions involving significant reporter metabolites, catalyzed by genes in complexes I, III, and IV of OXPHOS (Fig S2) , had a distinct expression pattern in PLWH ART compared with PLWH EC . Next, we performed context (disease state)-specific GSMM and flux balance analysis (FBA) to calculate the metabolic flux in response to transcriptional changes in the PLWH ART , PLWH EC , and HC cohorts ( Fig 3B) . Contextspecific metabolic models for PLWH ART , PLWH EC , and HC having 6,179, 6,237, and 6,199 reactions and 1,799, 1,842, and 1,834 genes/ transcripts catalyzing them, respectively, were developed (available: github.com/neogilab/LongART). After excluding the reactions with same directional fluxes in all the three cohorts and reactions with insignificant flux (<10 −5 mmol/h/gDCW), 80 reactions (Supplemental Data 3) were found to be uniquely regulated in PLWH ART Figure 2 . Comparative analysis of PLWH ART and PLWH EC . (A) Relative reservoir quantification using total HIV-1 DNA in PLWH ART (n = 17) and PLWH EC (n = 14). (B) MA plot showing differential gene expression results of PLWH ART versus PLWH EC . Negative log 2 -fold change values represent down-regulation and positive values represent up-regulation in PLWH ART . Grey-colored dots denote nonsignificant genes (adjusted P > 0.05). (C) Heatmap showing the expression pattern of significantly regulated genes between PLWH ART and PLWH EC (adjusted P < 0.05). Column annotation denotes cohort, age, gender, and duration of combination antiretroviral therapy of the corresponding samples. Row and column clustering was performed using Euclidian distance. (D) Gene set enrichment analysis results using MSigDB hallmark gene set between PLWH ART versus PLWH EC . A positive enrichment score represents up-regulation and negative score represents down-regulation in PLWH ART . Statistically significant pathways are labeled and highlighted by asterisk. Bubble size is relative to the adjusted P-values of the pathways. *Indicates FDR < 0.2. (E) Schematic visualization of the five complexes of OXPHOS pathway. The heatmap shows expression pattern of genes belonging to OXPHOS pathway in PLWH ART and PLWH EC . Column annotation denotes OXPHOS pathway complexes and row annotation denotes the cohort. The bottom annotation shows the log 2 fold change values of the genes. Red color represents up-regulation and green color represents down-regulation of the gene in PLWH ART compared with PLWH EC . compared with PLWH EC and HC cohorts. These reactions belonged to the AA, nucleotide, carbohydrate, and energy metabolism pathways. There were also 33 significant transport reactions that were transporting metabolites between cell compartments. Of the energy metabolism, pathways surrounding the tricarboxylic acid (TCA) cycle, including glycolysis, glutaminolysis, and OXPHOS, were affected in PLWH ART (Fig 3C) . The OXPHOS reaction converting ADP to ATP (HMR-6916) had a positive flux in PLWH ART whereas no flux was shown in PLWH EC and HC, indicating that higher energy was required in PLWH ART . There were also cytoplasmic reactions that appeared to increase the production of α-ketoglutarate (αKG) in PLWH ART . Reactions producing fructose-6-phosphate (HMR-7749 and HMR-4489), which further feeds the reaction producing glutamate (HMR-4300), showed a positive flux in PLWH ART , whereas showing a negative or no flux in PLWH EC and HC, respectively, indicative of higher glutamate production and conversion in PLWH ART . The reaction converting glutamate and oxaloacetic acid (OAA) to αKG and aspartate (HMR-3829) also showed a flux towards αKG production in PLWH ART and the opposite direction in PLWH EC and HC. Also, the transporter reaction (HMR-6024) transporting αKG from extracellular space to the cytoplasm showed a flux in PLWH ART . The reactions mentioned above signify increased accumulation of αKG in the cytoplasm in PLWH ART that can feed the TCA cycle in the mitochondria. The reaction producing OAA from pyruvate (HMR-4143) and the reaction producing αKG from oxalosuccinate (HMR-4112) had positive flux in PLWH ART indicating activation of the TCA cycle. To further understand the metabolic rearrangements, we performed a topological analysis on metabolic networks generated for PLWH ART , PLWH EC , and HC cohorts. The metabolic networks were generated by drawing edges between reactants, products, and associated genes of the reactions found to exhibit significant (>10 −5 mmol/h/gDCW) and diverging flux among the three cohorts. Communities were identified and betweenness centrality of the nodes was calculated to rank the genes and metabolites for their influence in the network. The top five metabolites and genes in PLWH ART , PLWH EC , and HC based on node centrality measurements are shown in Fig 3D. The metabolites fructose-6-phosphate, OAA, glutamate, and pyruvate uniquely play a central role in PLWH ART indicative of a role of TCA cycle and glycolysis in differentiating PLWH ART from PLWH EC and HC. Transporter genes of the SLC16 gene family (monocarboxylate transporters) SLC16A3, SLC16A6, and SLC16A7 were central in all three groups further suggesting a role for pyruvate and lactate transport. Combining these results, it can be concluded that reactions surrounding the TCA cycle including glutaminolysis, OXPHOS, and glycolysis differentiate PLWH ART from PLWH EC and HC.

The earlier used antiretrovirals frequently induced severe adverse effects that were linked to the occurrence of oxidative stress and mitochondrial damage (Smith et al, 2017) . As we observed an upregulation of superoxide, ubiquinol, ubiquinone, ferricytochrome C, and ferrocytochrome C in PLWH ART we evaluated total cellular ROS levels in different PBMCs subpopulations from PLWH ART (n = 16), PLWH EC (n = 16), and HC (n = 18), using flow cytometry (Fig S3A) . The distribution of CD4 + and CD8 + T cells, classical (CM), intermediate (IM), and non-classical monocytic (NCM) populations are depicted in UMAP (Fig 4A) . The percentage of CD4 + T cells were decreased, whereas CD8 + T cells were increased in PLWH EC and PLWH ART compared with HC ( Fig S3B) . PLWH EC also exhibited a decreased proportion of CM compared with HC, but no other differences were identified on the monocytic subpopulations ( Fig S3C) . We did not observe any significant differences in ROS on CD4 + or CD8 + T cells (Figs 4B and S3D ). In CM, ROS was significantly higher in PLWH ART samples compared with HC (Fig 4B) . Some of the PLWH ART had higher ROS, whereas others expressed lower ROS on lymphocytic cell populations (Fig 4C) . Therefore, to determine if long-term successfully treated HIV-1 infection influenced ROS production, we divided arbitrarily the PLWH ART group into long-term ART ([>10 yr, n = 8] with a median of 19 [16] [17] [18] [19] [20] [21] [22] yr treatment) and short-term cART ([<10 yr, n = 8] with a median of 7 [6-8] yr treatment). Interestingly, levels of ROS were increased in long-term ART compared with short-term ART on CD4 + T cells and compared with short-term ART and PLWH EC on CD8 + T cells (Fig 4D) . ROS levels were not affected by cART treatment on the monocytic cell populations (Fig S3E) . These data highlights the effects of longterm cART treatment on oxidative stress and redox homeostasis in lymphocytic cell populations.

In the ex vivo part of this study, we showed how up-regulation of OXPHOS was a signature of PLWH ART that differentiated them from the PLWH EC and how long-term treatment influenced oxidative stress and redox homeostasis on lymphocytic cell populations. Therefore, we decided to study the effect of inhibiting OXPHOS complexes I-V in a lymphocytic latency cell model (J-Lat 10.6) together with the parental cell line (Jurkat) by targeting complex I (metformin), complex II (D-α-tocopheryl succinate, aTOS), complex III (antimycin), complex IV (arsenic trioxide), and complex V (oligomycin) (Fig 5A) with respect to apoptotic properties, latency reversal and cellular senescence. The drugs did not have any effect on cell viability (Figs 5B and S4A ), although inhibition of complex I, III, and IV in J-Lat 10.6 increased the levels of Annexin V (a marker of apoptosis) compared with the respective untreated control, whereas only inhibition of complex I and IV showed the same effect in Jurkat cells (Figs 5C and S4B ). This indicates the role of OXPHOS complex III in the apoptotic properties of the HIV-1 latent cell model. It was only when inhibiting complex IV a significant increase in HIV-1 reactivation was observed in J-Lat 10.6 cells (Figs 5D and S4C) . Several studies including ours have shown that PLWH ART has a potential for attenuated immune aging due to a shift in glutaminolysis, in a subset of PLWH ART who had dysregulated metabolic profiles. A recent pivotal study also indicated a role of glutaminolysis in senolysis (removal of senescence cells) as senescent cells are dependent on glutaminolysis (Johmura et al, 2021) . To prove this, we measured the senescence markers CD57, Ki-67, and PCNA using flow cytometry and DNA damage marker H2A.X (S139) by Western blot. Cell surface expression of CD57 was not altered compared with the respective control when inhibiting the OXPHOS complexes although a baseline increase in CD57 was seen in J-Lat 10.6 cell compared with Jurkat (Fig S5A-C) . The proportion of Ki-67-negative J-Lat 10.6 cells increased after inhibiting complexes II, III, and IV, whereas only inhibition of complex IV increased the proportion of Ki-67 negative Jurkat cells (Figs 5E and S5D) . A mild decrease in KI-67 negative cells was also observed in J-Lat 10.6 cells when inhibiting complex I (Figs 5E and S5D) . Inhibition of complex III increased PCNA negative J-Lat 10.6, whereas no significant differences were observed in Jurkat cells (Figs 5F and S5E) . Furthermore, phosphorylation of H2A.X (S139) increased when inhibiting complex IV and decreased when inhibiting complex I, irrespective of the cell type (Figs 5G and H and S5F) . The original blots were presented as a source file to Fig 5. Collectively, our data highlighted the potential role of pharmacological inhibition of the OXPHOS complexes with differential regulation of latency reversal, apoptotic properties, and cellular senescence in lymphocytic HIV-1 latent cell model, depending on compounds and targeted complexes.

In the present study, we combined system-level blood cell transcriptomics and developed context-specific GSMM to provide a comprehensive system-level characterization of HIV-1 infected individuals with suppressive viremia either by natural (PLWH EC ) or drug-induced (PLWH ART ) control. The transcriptomic data identified up-regulation of OXPHOS as the characteristic feature of PLWH ART , differentiating them from HIV-1 seropositive PLWH EC , who were not on therapy. The main dysregulation seemed to occur in complexes I, III, and IV of the OXPHOS pathway. FBA identified altered flux in several glycolytic intermediates like pyruvate, αKG, glutamate, and fructose-6-phosphate in PLWH ART compared with PLWH EC and HC. Long-term cART also affected the redox homeostasis in T lymphocytes. The in vitro pharmacological inhibition of the OXPHOS complexes in the latent lymphocytic cell model suggested a role of the complex IV in latency reversal, complex I, III, and IV in apoptosis, and complex IV in immunosenescence. Altered glutaminolysis (i.e., glutamine lysed to glutamate) and increased plasma glutamate have been observed in several cohorts from both high income (Gelpi et al, 2021) and low-and middleincome countries (Gelpi et al, 2021) and are required for optimal HIV-1 infection of CD4 + T cells (Clerc et al, 2019) . Glutaminolysis is the primary pathway fueling the TCA cycle and OXPHOS in naïve and memory T cell subsets which are critical factors for immune recovery in successfully treated PLWH (Rosado-Sánchez et al, 2019). HIV-1 infection is more common in T cells with elevated glycolysis and OXPHOS and inhibition of these metabolic activities can block HIV-1 replication and reservoir transactivation (Valle-Casuso et al, 2019). Impairment of the metabolic steps preceding OXPHOS can also result in lipid accumulation in macrophages (Castellano et al, 2019) . Enhanced glycolysis and OXPHOS are characteristics of CD8 + T cell exhaustion (Rahman et al, 2021) . However, long-term molecular immune pathogenic consequences of successful cART have not yet been evaluated. In our study, we identified system-level up-regulation of OXPHOS as the main characteristic of PLWH on long-term cART. When comparing with PLWH EC , an up-regulation of OXPHOS, and to certain extent glycolysis, was observed in PLWH ART . Given that HIV-1 preferentially selects cells that have elevated cellular OXPHOS and glycolysis for infection and replication, reservoir seeding (Hegedus et al, 2014; Palmer et al, 2014; Valle-Casuso et al, 2019) and cell-to-cell spread of HIV-1, this metabolic environment permit ongoing replication during cART (Sigal et al, 2011) . Therefore, we hypothesize that up-regulation of OXPHOS in PLWH ART was the reason behind the relatively larger HIV-1 reservoir in long-term successfully treated infection compared with PLWH EC with natural control of viral replication. This metabolic modulation could potentially be a barrier to the post-treatment control of viral replication.

A recent seminal study showed that a higher HIV-1 viral set point in untreated patients during acute HIV-1 infection correlated positively with OXPHOS and that in vitro pharmacological inhibition of complex I (by rotenone or metformin) and complex III (by antimycin A) suppressed viral replication and immunometabolism through an NLRX1 and FASTKD5-dependent mechanism (Guo et al, (H) Quantification of H2A.X (S139). The graph shows fold change (Fc) of protein expression in relation to respective control after normalization to β-Actin. Experiments were performed in three biological replicates. Significance was determined using two-tailed t test (P < 0.05 with * < 0.05, ** < 0.033, *** < 0.002) and represented with mean and SD. Significance for each drug is compared with respective control. See also Figs S4 and S5. Source data are available online for this figure.

GSMM for HIV-1 control Ambikan et al.

https://doi.org/10.26508/lsa.202201405 vol 5 | no 9 | e202201405 2021). Furthermore, we recently observed that blocking glycolysis with 2-deoxyglucose (2-DG) increased cell death in lymphocytic and pre-monocytic HIV-1 latent cell models (Gelpi et al, 2021) , in line with other studies (Valle-Casuso et al, 2019; Guo et al, 2021) . These studies indicate a critical role for glycolysis and OXPHOS in HIV-1 immuno-pathogenesis. In the present study, we observed that latently infected cells treated with antimycin resulted in increased markers of apoptosis in latent J-Lat 10.6 cells compared with the parental Jurkat cells. This indicates an increased preferential cell death of the latently infected cells without latency reversal. Only inhibition of complex IV by arsenic trioxide showed a small degree of latency reversal. Both the transcriptomic data and the in vitro assays indicated the role of complex I, III, and IV as essential components of the electron transport chain for generation of ATP and cellular energy requirements. Complexes I and III have a role in ROS production and are essential in inflammatory macrophages and T helper 17 (T H 17) cells while also playing a vital role in lymphocyte activation, proliferation, and differentiation (Yin & O'Neill, 2021) . Recently, it has been shown that complex III is crucial for the suppressive function of Tregs (Weinberg et al, 2019) . Our DCQ identified increased frequency of Tregs in PLWH ART compared with PLWH EC which is in line with recent findings (Caetano et al, 2020) ; however, PLWH EC can present more activated Tregs (Gaardbo et al, 2014; Caetano et al, 2020) . Finally, an earlier study reported the Cox-II enzyme leads to reduced T-cell apoptosis in HIV-1 infected cells (Tripathy & Mitra, 2010) . In contrast, our study indicated pharmacological inhibition of the complex IV with arsenic trioxide increased apoptosis (as measured by the annexin V) both in latent J-Lat 10.6 cells and non-latent Jurkat. Interestingly, inhibition of complex IV in J-Lat 10.6 cells also showed latency reactivation which could potentially be linked to apoptosis. In our FBA, we identified altered flux in pyruvate, glutamate, and αKG in the PLWH ART compared with PLWH EC and HC. Recently, we identified a higher level of glutamate in PLWH ART in several cohorts compared with HC (Gelpi et al, 2021) . The level was even higher in PLWH ART with metabolic syndrome (Gelpi et al, 2021 ). Blood glutamate level has been reported to be higher in PLWH with dementia (Ferrarese et al, 2001) . Reducing the blood glutamate concentrations with blood glutamate scavengers like pyruvate facilitates the efflux of glutamate from the brain to the blood. This can limit the neurotoxic effect of glutamate (Boyko et al, 2012) and has been reported to effectively improve neurological recovery in traumatic brain injury (Gottlieb et al, 2003; Zlotnik et al, 2007; Boyko et al, 2012) . The coordination between glutamate and pyruvate and its neuroprotective role in chronic HIV-1 infected patients on therapy needs further studies to understand neurological complications in HIV infection after successful treatment.

Although immune cell senescence decreases the overall cellular activity, it is associated with a high metabolic need, usually by increasing aerobic glycolysis. In the case of our lymphocytic cell culture model, we detected an enrichment of the senescent marker CD57 compared with the parental cell line, indicative of increased chronic activation of latently infected cells. Furthermore, we detected increased levels of DNA damage (H2A.X [S139]) (Mah et al, 2010) , decreased proliferation (Ki-67) (Lawless et al, 2010) , and DNA replication (PCNA) (González-Magaña & Blanco, 2020) after OXPHOS inhibition. Earlier studies have shown how OXPHOS inhibition in human fibroblasts induced senescence (Stöckl et al, 2006) . Therefore, high plasticity of metabolic reprogramming could induce an increase in glycolysis during OXPHOS inhibition which could potentially be coupled to induction of senescence in the HIV-1 latent cells during the suppressive therapy.

Our study also showed that ROS was increased in patients on long-term (median 19 yr) compared with short-term (median 7 yr) of suppressive therapy. This could be linked to the use of the older nucleoside reverse transcriptase inhibitors (NRTIs) like zidovudine (AZT), stavudine (d4T), or didanosine (ddI) as a part of the initial treatment regimen. The cell's epigenetic state is closely associated with ROS-induced oxidative stress due to mitochondrial damage and altered OXPHOS (Guillaumet-Adkins et al, 2017) . It is known that antiretrovirals such as AZT, d4T, and ddI can cause mitochondrial damage, ultimately altering OXPHOS (Pinti et al, 2006) . Recent molecular studies have reported that PLWH on treatment has epigenetic age acceleration (Gross et al, 2016) compared with the non-infected individuals that can partially be reversed with cART initiation (Esteban-Cantos et al, 2021) . Therefore, understanding the biological mechanism of potential accentuated aging in PLWH on long-term successful therapy who were exposed to earlier generation treatment regimen and dysregulated metabolic profile could potentially provide a clinical intervention strategy to improve the quality of life of PLWH ART .

In conclusion, our study indicated a system-level up-regulation of OXPHOS and, to a certain extent, glycolysis in PLWH ART compared with the PLWH EC . Furthermore, we show how this up-regulation could play a role in latent reservoir dynamics and immunosenescence in HIV-1-infected individuals with long-term successful therapy. Pharmacological inhibition of the OXPHOS complexes could have a role in latency reversal, apoptotic properties, and immunosenescence in latently infected cells. Further studies are warranted to elucidate the molecular mechanisms underlying the observed shift in OXPHOS in PLWH ART and how its coordination with glutaminolysis can lead to immune dysregulation during successful therapy.

The study population includes three groups of PLWH, with two groups as suppressed viremia (PLWH ART and PLWH EC , n = 19 each), and one group with viremia (PLWH VP herein, n = 19). In addition, we enrolled 19 HC. The study was approved by the regional ethics committees of Stockholm (2013 Stockholm ( /1944 Stockholm ( -31/4 and 2009 Stockholm ( /1485 and amendment (2019-05585 and 2019-05584, respectively) and performed in accordance with the Declaration of Helsinki. All participants gave informed consent. The patient's identity was anonymized and delinked before analysis.

PBMCs were used for RNA-sequencing (RNA-Seq) using Illumina HiSeq2500 or NovaSeq6000 as described by us . Differential gene expression analysis was performed using the R/Bioconductor package DESeq2 v1.26.0 (DOI: 10.18129/ B9.bioc.DESeq2). Gene list enrichment analysis for cART-specific genes was performed using enrichr module of python package GSEAPY v0.9.16 (Subramanian et al, 2005; Chen et al, 2013) and MSigDB hallmark gene set v7.4. GSEA between PLWH ART and PLWH EC was performed using GSEA v4.1.0 software (Subramanian et al, 2005) and MSigDB hallmark gene set v7.4. Metabolomics data were generated using the Metabolon HD4 (Metabolon Inc.) .

DCQ by measuring the proportion of different cells in each sample was performed using the deconvolution algorithm adapted from Estimating the Proportions of Immune and Cancer cells (Chen et al, 2013) . The reference gene expression profile consists of gene-level expression data of 18 blood cell types and it is based on Human Protein Atlas version 20.1 and Ensembl version 92.38. Signature genes for the 18 blood cell types in the reference profile were downloaded from CellMarker (Zhang et al, 2019) and Pan-glaoDB (Franzén et al, 2019) . The transcript per million (TPM) transformed gene expression data of all genes from the samples were used in the procedure along with reference profile and signature gene list to estimate the cell proportion.

Significantly regulated genes (adjusted P < 0.05) in all the pair-wise comparisons among the four cohorts were used to derive the cARTspecific genes. The list of significant genes in each of the comparisons was considered as individual sets and various set operation procedures were used for the derivation. The set operations performed are represented below. ART = {z | z 2 X 1 or z 2 X 2 or z 2 X 3 } NULL = {z | z 2 Y 1 or z 2 Y 2 or z 2 Y 3 } ART-specific genes = {z | z 2 ART and z Ï NULL} where, X 1 = {z | z is gene regulated in HC versus PLWH ART } X 2 = {z | z is gene regulated in PLWH EC versus PLWH ART } X 3 = {z | z is gene regulated in PLWH VP versus PLWH ART } Y 1 = {z | z is gene regulated in HC versus PLWH VP } Y 2 = {z | z is gene regulated in PLWH EC versus PLWH VP } Y 3 = {z | z is gene regulated in HC versus PLWH EC } GSMM, FBA, and essentiality analysis Group-specific human GSMMs were reconstructed by integrating transcriptomics data on human reference GSMM obtained from Metabolic Atlas (Robinson et al, 2020 ). The metabolic model reconstruction was performed using task-driven Integrative Network Inference for Tissues (tINIT) algorithm (Agren et al, 2012 (Agren et al, , 2014 Robinson et al, 2020) . The algorithm creates a context-specific model by selecting only reaction that can carry flux based on the provided transcript expression table (transcript per million). The reconstructed models were then checked for biological feasibility by analyzing their capacity to carry out 56 essential metabolic tasks. FBA was performed using MatLab function solveLP from RAVEN toolbox v2.4.0 (Wang et al, 2018) and ATP hydrolysis as objective function. Plasma metabolomics data were used as a reference to constrain the exchange reactions in the model assuming that exchange reaction fluxes were relatively influenced by availability of extracellular metabolites. We calculated log 2scaled changes of exchange metabolites against the control cohort, and it was used proportionally to compute the reaction bounds.

Network topology analysis was performed on the metabolic networks generated for the cohorts. The metabolic networks were created by drawing edges between reactants, products, and enzymatic genes of each of the reactions, which showed significant (>10 −5 ) and varying flux values among the cohorts. The networks were then analyzed using igraph toolkit. The absolute value of the flux scaled between 0 and 1 was used as edge weight. Leiden algorithm (Traag et al, 2019 ) was used to identify communities and the betweenness centrality of all the nodes was computed. Nodes were ranked based on their centrality measurement. Nodes with high centrality were considered as most influential for the existence and functioning of the network.

R package ggplot2 v3.3.2 (Wickham, 2016) was used to create all bubble plots, scatter plots, and boxplots. R/Bioconductor package ComplexHeatmap v2.2.0 (Gu et al, 2016) was used to create all the heat maps. Network diagrams were drawn in Cytoscape ver 3.6.1 (Shannon et al, 2003) . Venn diagrams were generated using the online tool InteractiVenn (Chen & Boutros, 2011) .

Total DNA was extracted from PBMCs using QIAamp DNA mini kit (QIAGEN) according to manufacturers' instructions. HIV-1 DNA quantification was performed using Internally Controlled qPCR (IC-qPCR) as described by Vicenti et al (2018) . In brief, total HIV-1 DNA was quantified in PLWH ART (n = 17) and PLWH EC (n = 14) using 500 ng of DNA in duplicates. Quantification was performed using Takara Universal Mastermix (Takara) on an ABI 7500F using primers (β globin F; AGGGCCTCACCACCAACTT, β globin R; GCACCTGACTCCTGAGGAGAA, HXB2 F; GCCTCAATAAAGCTTGCCTTGA, HXB2 R; GGCGCCACTGCTAGA-GATTTT) and probes (β globin; HEX-ATCCACGTTCACCTTGCCCCACA-TAM, HXB2; FAM-AAGTAGTGTGTGCCCGTCTG-MGBEQ) targeting β globin and HIV-1 (HXB2) and normalized to β globin levels.

The latency cell model J-Lat clone 10.6 (NIH HIV reagent program) was used together with its parental cell line Jurkat. Cells were cultured in StableCell RPMI 1640 (Sigma-Aldrich) supplemented with 10% fetal bovine serum (Gibco) and 20 U/ml penicillin and 20 μg/ml streptomycin (Gibco) at 37°C and 5% CO 2 .

Cytotoxicity assays were performed for metformin (Sigma-Aldrich), arsenic trioxide (Sigma-Aldrich), oligomycin (Sigma-Aldrich), antimycin (Sigma-Aldrich), and aTOS (Sigma-Aldrich) (Fig S3A) . Experimental concentrations with low cytotoxicity were chosen and assayed for 24 h. All assays were performed in biological triplicates and analyzed for viability using flow cytometry, as described below.

PBMCs were subjected to flow cytometry analysis. Samples were thawed in 37°C water bath and washed with flow cytometry buffer (PBS + 2% FBS + 2 mM EDTA). Total cellular ROS levels were measured using the CellROX Deep Red Flow Cytometry Assay Kit (Invitrogen) according to the manufacturer's instructions. Briefly, 750 nM of CellROX deep red reagent was added to PBMCs and incubated for 1 h at 37°C, protected from light. The cells were then stained with Live/Dead fixable near-IR dye (Invitrogen), and cell surface markers were detected by incubating cells with anti-CD3 (clone OKT3, BD Bioscience), anti-CD4 (clone SK3; BD Bioscience), anti-CD8 (clone HIT8a; BioLegend), anti-CD14 (clone M5E2; Bio-Legend), and anti-CD16 (clone 3G8; BD Bioscience) for 20 min on ice in flow cytometry buffer. Positive and negative controls for ROS measurement were performed by incubating cells with either tertbutyl hydroperoxide (200 µM) or N-acetyl cysteine (5 mM) for 45 min at 37°C before the addition of CellROX deep red reagent. All cells were fixed with 2% paraformaldehyde before acquiring on BD FACS Symphony flow cytometer (BD Bioscience). Compensation setup was performed using single-stained controls prepared with antibody-capture beads: anti-mouse Ig, κ/negative control compensation particles set (BD Biosciences) for mouse antibodies and ArC amine-reactive compensation bead kit (Invitrogen) for use with LIVE/DEAD fixable dead cell stain kits.

Flow cytometry for cell lines was conducted by extracellular staining using anti-CD57 (clone HNK-1; BioLegend) and LIVE/DEAD Near-IR viability stain (Invitrogen) followed by fixation using ki-67 fixation/permeabilization kit (eBioscience). Intracellular staining was performed using anti-Ki-67 (clone Ki-67; BioLegend) and anti-PCNA (clone PC10; BioLegend). Analysis of apoptosis was performed using Annexin-V Alexa647 conjugate (Thermo Fisher Scientific) staining in combination with LIVE/DEAD Near-IR viability stain (Invitrogen) prior fixation using 4% PFA. Samples were acquired on BD FACS Fortessa (BD Bioscience). Flow cytometry data were analyzed and compensated with FlowJo 10.6.2 (TreeStar Inc.) and statistical analysis was performed using Mann-Whitney U test or two-tailed t test in Prism 9.3.0 (GraphPad Software Inc.).

Cells were lysed in RIPA buffer (Sigma-Aldrich) supplemented with 1× PhosSTOP (Sigma-Aldrich) and 2× cOmplete protease inhibitor cocktail (Roche) on ice for 30 min. Protein estimation was performed using DC protein assay (Bio-Rad Laboratories) and 37.5-48 μL of protein run in each well on NuPage 4-12% BisTris 20 well, 1 mm precast gels (Thermo Fisher Scientific) and transferred using the iBlot transfer system (Invitrogen) with iBlot PVDF Transfer stack (Invitrogen). Membranes were incubated with Phospho-Histone H2A.X (Ser139) (Cell Signaling Technology) and β-Actin (Sigma-Aldrich). Secondary antibody incubation was performed using Dako Immunoglobulins/HRP (Aglient Technologies) and membranes developed using ECL (Amersham) on ChemiDoc (Bio-Rad Laboratories). Relative protein quantification was performed using ImageLab 6.0.1 (Bio-Rad Laboratories) and statistical significance using a two-tailed t test in Prism 9.3.0 (GraphPad Software Inc.).

The raw RNA sequencing (RNAseq) data have been deposited in the NCBI/SRA with PRJNA420459. The metabolomics data are available from dx.doi.org/10.6084/m9.figshare.19747582. All the codes are available at github: github.com/neogilab/LongART.

Supplementary Information is available at https://doi.org/10.26508/lsa. 202201405.

Reconstruction of genome-scale active metabolic networks for 69 human cell types and 16 cancer types using INIT

Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling

Integrative proteotranscriptomic and immunophenotyping signatures of HIV-1 Elite control phenotype: A cross-talk between glycolysis and HIF signaling

Multiomics personalized network analyses highlight progressive immune disruption of central metabolism associated with COVID-19 severity

Plasma metabolic signature and abnormalities in HIV-infected individuals on long-term successful antiretroviral therapy

The effect of blood glutamate scavengers oxaloacetate and pyruvate on neurological outcome in a rat model of subarachnoid hemorrhage

HIV-1 elite controllers present a high frequency of activated regulatory T and Th17 cells

Plasma metabolomics identifies lipid abnormalities linked to markers of inflammation, microbial translocation, and hepatic function in HIV patients receiving protease inhibitors

HIV infection and latency induce a unique metabolic signature in human macrophages

Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool

VennDiagram: A package for the generation of highly-customizable Venn and Euler diagrams in R

Entry of glucose-and glutamine-derived carbons into the citric acid cycle supports early steps of HIV-1 infection in CD4 T cells

Epigenetic age acceleration changes 2 years after antiretroviral therapy initiation in adults with HIV: A substudy of the NEAT001/ ANRS143 randomised trial

Increased glutamate in CSF and plasma of patients with HIV dementia

PanglaoDB: A web server for exploration of mouse and human single-cell RNA sequencing data

Immunoregulatory T cells may be involved in preserving CD4 T cell counts in HIV-infected long-term nonprogressors and controllers

The central role of the glutamate metabolism in long-term antiretroviral treated HIV-infected individuals with metabolic syndrome

Human PCNA structure, function and interactions

Blood-mediated scavenging of cerebrospinal fluid glutamate

Methylome-wide analysis of chronic HIV infection reveals five-year increase in biological age and epigenetic targeting of HLA

Complex heatmaps reveal patterns and correlations in multidimensional genomic data

Epigenetics and oxidative stress in aging

Multi-omics analyses reveal that HIV-1 alters CD4 + T cell immunometabolism to fuel virus replication

HIV-1 pathogenicity and virion production are dependent on the metabolic phenotype of activated CD4+ T cells

Senolysis by glutaminolysis inhibition ameliorates various age-associated disorders

Recommendations for analytical antiretroviral treatment interruptions in HIV research trials-report of a consensus meeting

Quantitative assessment of markers for cell senescence

GammaH2AX as a molecular marker of aging and disease

Unbiased metabolomics links fatty acid pathways to psychiatric symptoms in people living with HIV

Lipid profiles and APOE4 allele impact midlife cognitive decline in HIV-infected men on antiretroviral therapy

Increased glucose metabolic activity is associated with CD4+ T-cell activation and depletion during chronic HIV infection

Persistent metabolic changes in HIV-infected patients during the first year of combination antiretroviral therapy

Anti-HIV drugs and the mitochondria

EPIC: A tool to estimate the proportions of different cell types from bulk gene expression data

Elevated glycolysis imparts functional ability to CD8 + T cells in HIV infection

An atlas of human metabolism

Glutaminolysis and lipoproteins are key factors in late immune recovery in successfully treated HIV-infected patients

Cytoscape: A software environment for integrated models of biomolecular interaction networks

Glycolysis downregulation is a hallmark of HIV-1 latency and sensitizes infected cells to oxidative stress

Cell-tocell spread of HIV permits ongoing replication despite antiretroviral therapy

Beyond the polymerase-γ theory: Production of ROS as a mode of NRTI-induced mitochondrial toxicity

Sustained inhibition of oxidative phosphorylation impairs cell proliferation and induces premature senescence in human fibroblasts

Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles

Differential modulation of mitochondrial OXPHOS system during HIV-1 induced T-cell apoptosis: Up regulation of complex-IV subunit COX-II and its possible implications

Cellular metabolism is a major determinant of HIV-1 reservoir seeding in CD4+ T cells and offers an opportunity to tackle infection

Development of an internally controlled quantitative PCR to measure total cell-associated HIV-1 DNA in blood

Integrative lipidomics and metabolomics for system-level understanding of the metabolic syndrome in long-term treated HIV-infected individuals

RAVEN 2.0: A versatile toolbox for metabolic network reconstruction and a case study on Streptomyces coelicolor

Mitochondrial complex III is essential for suppressive function of regulatory T cells

Elegant Graphics for Data Analysis

A network-based approach reveals the dysregulated transcriptional regulation in non-alcoholic fatty liver disease

The role of the electron transport chain in immunity

Combined metabolic activators therapy ameliorates liver fat in nonalcoholic fatty liver disease patients

Transcriptomics and targeted proteomics analysis to gain insights into the immune-control mechanisms of HIV-1 infected elite controllers

Quantitative humoral profiling of the HIV-1 proteome in elite controllers and patients with very long-term efficient antiretroviral therapy

CellMarker: A manually curated resource of cell markers in human and mouse

Brain neuroprotection by scavenging blood glutamate

License: This article is available under a Creative Commons License

The authors acknowledge support from the National Genomics Infrastructure in Genomics Production Stockholm funded by Science for Life Laboratory, the Knut and Alice Wallenberg Foundation and the Swedish Research Council, and SNIC/Uppsala Multidisciplinary Center for Advanced Computational Science for assistance with massively parallel sequencing and access to the UPPMAX computational infrastructure. The study is funded by the Swedish Research Council grants 2017-01330, 2018-06156, and 2021-01756 to U Neogi and 2017-05848 and 2020-02129 to A Sönnerborg.

AT Ambikan: data curation, formal analysis, investigation, visualization, methodology, and writing-original draft. S Svensson-Akusjärvi: data curation, formal analysis, visualization, methodology, and writing-original draft. S Krishnan: visualization, methodology, and writing-review and editing. M Sperk: formal analysis, methodology, and writing-review and editing. P Nowak: data curation, investigation, and writing-review and editing. J Vesterbacka: data curation, investigation, and writing-review and editing. A Sönnerborg: conceptualization, resources, funding acquisition, investigation, project administration, and writing-review and editing. R Benfeitas: supervision, visualization, methodology, project administration, and writing-review and editing. U Neogi: conceptualization, resources, supervision, funding acquisition, visualization, methodology, writing-original draft and project administration.

The authors declare that they have no conflict of interest.