Multiomic profiling of chronically activated CD4+ T cells identifies drivers of exhaustion and metabolic reprogramming

CD4+ memory T cells exhibit progressive exhaustion phenotype during chronic stimulation in vitro

While a complex phenotype in the context of a tumor or infection, T-cell exhaustion can be replicated in culture models by repeatedly exposing primary T cells to antigen or CD3/CD28 activation [55], and has been used for discovery of mediators of the exhausted state [56–59]. We adapted this protocol to profile the molecular kinetics of human CD4+ T-cell exhaustion, focusing on memory CD4+ T cells to minimize cellular heterogeneity, given their abundance and critical functions in the tumor microenvironment. We chose the memory subset not only for their functionality in anti-tumor responses, but also as a starting cell population that would allow us a balance of specificity/homogeneity (versus whole CD4+ T cell population) and practicality (allowing us to isolate a large enough population to culture to perform our experiments). These cells also expand well in response to antigen stimulation, a benefit for our downstream analysis material needs.

To define the molecular basis of dynamic cellular state transitions of the initiation and progression of T-cell exhaustion, our pipeline employed multi-pronged quantitative proteomic, phosphoproteomic, and metabolomic workflows to unbiasedly analyze changes in global protein expression, phosphorylation, and metabolite abundance. CD4+ CD45RO+ T cells were isolated from 2 healthy donors (donor #1 primarily shown, but donor #2 expression and protein agreement can be found in S1 Fig. Complete data set in S1 Data). Donor #2 samples were collected at slightly differing time points from donor #1 to capture different expression dynamics during chronic stimulation, but the key time points of 0, 6, and 12 days were kept the same for comparison between data sets. T cells were expanded prior to start of the chronic stimulation protocol, whereby cells were next subjected to chronic CD3 and CD28 stimulation using magnetic bead-bound antibodies in the presence of low levels of IL-2 (to maintain cell viability, S2 Fig) over the course of 12 days (Fig 1A). Media and antigen (CD3/CD28 Dynabeads) were replenished every 2 days, and cell aliquots were collected for phenotyping, cytokine profiling, and multiomic analysis at select intervals (Fig 1B). We profiled cell cultures harvested at representative time points (days 6, 8, 10, and 12), along with cells sustained in a “resting” condition (i.e., stimulation for 2 days before antigen removal followed by culturing without antigen until day 10 in cultures parallel to those undergoing repeated stimulation) (Fig 1A).

Fig 1. In vitro chronic CD4+ T-cell activation model replicates exhaustion phenotype while providing material for multiomics workflow.

(A) Diagram depicting the culture model used to generate time course samples from healthy human CD4+ memory T cells. Cells were stimulated with CD3/CD28 Dynabeads for 12 days, with fresh beads plus full media change every other day (fresh media alone added on alternate days). Samples were collected on days highlighted in bold. (B) Multiomics workflow (see Methods for details). Experiments were performed in biological triplicate, with proteomics, phosphoproteomics, and metabolomics done on same samples. (C) ELISA-based quantification of secreted cytokines (IL-2, IFN-γ, TNF-α) on days 0, 6, 10, and 12, 3 replicates per time point. (D) Cell counts were taken each time point, and average cell division compared to previous time point calculated. Cell viability remained stable throughout at ~75%–85%. (E) Flow cytometric analysis of inhibitory receptors PD-1, TIM-3, and TIGIT at resting, day 6, and day 12 time points. Data normalized to mode for accurate population comparisons to account for differences in number of cells analyzed. The data underlying Fig 1E can be found in the supplemental flow cytometry files uploaded to flowrepository.org. Figure created with BioRender.com.

https://doi.org/10.1371/journal.pbio.3002943.g001

Conditioned media was recovered at each time point for 3 prominently secreted cytokines (IL-2, IFN-γ, TNF-α) 8 to 10 h after re-stimulation to confirm decreased secretion as typically seen in exhaustion. Consistent with previous studies [55,56], we observed a marked decrease in cytokine secretion by 10 days with almost complete loss of these functions by 12 days (Figs 1C and S3), along with decreased cell proliferation (Fig 1D), and increased expression of well-known inhibitory receptors (PD-1, TIM-3, LAG-3, TIGIT) (Figs 1E and S1A). IL-2 secretion was down-regulated earliest, reaching almost undetectable levels by day 6, whereas IFN-γ and TNF-α showed maximally reduced levels by day 12 (Fig 1C). Cell proliferation followed a similar trajectory, with high cell division until day 10, followed by decreased cell division rates until proliferation became stagnant by day 12 (Fig 1D). Some cell death was observed throughout the protocol, especially in the early stages of stimulation (day 2 through day 6 or 8), leveling out later in the protocol (day 8 to 10), and slightly increasing by day 12 as the cells enter later stages of chronic stimulation. As expected, flow cytometry measurements revealed pronounced up-regulation of cell surface levels of PD-1, TIM-3, TIGIT, LAG-3, CD160, and CD137 (Figs 1E and S4), with some level of heterogeneity of marker co-expression (S4 Fig).

Taken together, these phenotypic results indicated that the CD4+ T cells were functionally exhausted by day 12. With confidence in our model, we went forward with our systematic characterization of samples collected over a time course, allowing us to capture how CD4+ memory T cells respond to chronic antigen simulation.

Quantitative proteomic analysis of chronically stimulated CD4+ memory cells illuminates temporal signaling events associated with exhaustion

Utilizing quantitative mass spectrometry using isobaric stable-isotope-based multiplex labeling (TMT), we measured dynamic changes in protein expression during the exhaustion time course. Post-filtering and normalization (S5A Fig), we detected and quantified with high confidence changes in the relative abundance a total of 8,401 proteins across all time points in the first donor sample (S5C Fig). Two comparisons of most physiological relevance were comparison of day 12 (termed “late exhaustion”) versus either control (termed “resting”) or day 6 cultures (termed “pre-exhaustion”). For these comparisons, 6,094 and 6,122 proteins, respectively, were deemed to be significantly differentially expressed (adjusted p-value Fig 2A). Among the proteins identified, many were exhaustion-associated (PDCD1, LAG3, ENTPD1, TOX, etc.), important CD4+ T-cell factors (CD40L, EOMES, GNG4, etc.), or T-cell signaling modulatory proteins (SLAM proteins, CD200, ICOS, CD27, CD28, etc.).

Fig 2. Global proteomic and phosphoproteomic profiling reveals signaling dynamics while recapitulating exhaustion markers.

(A) Differential protein expression in late exhaustion (day 12) versus resting, and late exhaustion (day 12) versus early exhaustion (day 6) populations. (B) Expression patterns of select T-cell exhaustion-associated proteins. Rows z-score normalized. (C) Expression patterns of putative exhaustion signature markers defined by Beltra and colleagues [60]. Rows z-score normalized. (D) Enrichment analysis of proteomics data using an exhaustion-associated geneset and a CD8+ T-cell exhaustion versus memory geneset derived from differential expression of RNA sequencing data (Wherry and colleagues [61]). (E) Global proteomic expression grouped according to expression profile by hierarchical clustering. Genes from each cluster were analyzed by Enrichr for pathway enrichment results. Rows z-score normalized. (F) Differential patterns of phosphosite abundance in late exhaustion (day 12) versus resting cells, and late exhaustion (day 12) versus early exhaustion (day 6). (G) Kinase-substrate enrichment analysis of differential phosphosites between late (day 12) and early (day 6) exhaustion. (H) Phosphorylation motif analysis showing predicted kinase classes implicated for observed changes in phosphorylation between late (day 12) and early (day 6) exhaustion. The data underlying this figure can be found in either S1 Data under proteomics or under phosphoproteomics, and D and G can be directly found in S3 Data. Figure created with BioRender.com.

https://doi.org/10.1371/journal.pbio.3002943.g002

As an initial assessment of reliability, we examined our differential expression data for known markers and drivers of exhaustion (Fig 2B). Specifically, we detected significant changes in the expression levels of the highly characterized inhibitory receptors PD-1, LAG-3, TIM-3, and TIGIT, that are known to be concurrently expressed on exhausted T cells [62], as well as multiple other previously reported markers. GNG4, recently identified to demarcate exhausted CD4+ T-cell population [63], was also significantly up-regulated in the time course (Fig 2B).

Expression profiles of the 4 major exhaustion-associated inhibitory receptors were significantly up-regulated across the time course (FDR S1 Data and S4 Data), thereby validating the model. Likewise, virtually all additional markers such as CD244 (2B4/SLAMF4), ENTPD1 (CD39), CD38, SLAMF7, SIGLEC10, and BTLA were significantly up-regulated by late exhaustion, while CTLA-4 and CD44, however, were not significantly changed. While also associated with T-cell activation, the persistent elevated co-expression of factors such as PD-1, ENTPD1, and CD38, is strongly linked to immune exhaustion [62].

Likewise, we observed significant changes in the abundance of previously described transcriptional regulators of T-cell exhaustion, including TOX, BLIMP-1 (PRDM1), EOMES, BATF, NR4A1/3, IRF4, ID2, and NFAT-family proteins. Notably, the transcription factor TOX, which is required for entry of CD8+ T cells into the exhausted state [23,64,65], showed progressive up-regulation in each time point relative to resting cell cultures. Additional evidence of exhaustion arises from our observation that Interleukin-7 Receptor Subunit Alpha (IL7R) protein expression is down-regulated in our chronically stimulated cultures. Specifically, this receptor is highly expressed in memory cells, as we observed in our resting condition, and is down-regulated when cells become exhausted [24]. We also observed an increase in both granzyme B and perforin at day 12, showing similar slight gain of function as seen in terminally exhausted CD4+ [51,52] and CD8+ T cells [53].

We compared our data set to exhaustion transcriptomic (mRNA) profiles reported in other studies. While CD4+ T-cell exhaustion data sets are largely lacking in the public domain, Beltra and colleagues used flow cytometry and RNA-seq to identify 4 subsets of CD8+ T-cell exhaustion states termed “terminal,” “intermediate,” “progenitor 1,” and “progenitor 2” [60]. Comparing our protein expression data to these gene clusters, we observed that day 12 chronically stimulated CD4+ cells closely resemble the “terminal” exhausted CD8+ cells (Fig 2C). We also observe some agreement with the “intermediate” phenotype, which is likely explained by heterogeneity innate in bulk cell population analysis.

Next, we employed Geneset Enrichment Analysis (GSEA) to further compare our proteomic results to previously reported exhaustion gene profiles. Using 2 genesets, one containing curated exhaustion proteins from the literature and the other from RNA-seq of differential gene expression measured in exhausted versus memory CD8+ T cells in mice [61]. We observe that over time during chronic stimulation, the normalized enrichment score (NES) of these signatures within our proteomic data set increases and reaches its peak by the final day 12 time point (Fig 2D).

Recapitulation of the dynamics of known exhaustion-associated proteins provided confidence in the model system. To verify the consistency of our findings, we nevertheless evaluated the experimental results obtained from the replicate second donor, again evaluating the surface expression of well-defined exhaustion markers, cytokine phenotyping, and unbiased quantitative proteomics for global protein expression. Notably, we observed high overall agreement in the patterns of inhibitory receptor expression based on flow cytometry and proteome expression (S1 Fig), pointing to the overarching robustness of our main findings.

Up-regulation of PRC2-mediated gene repression, prostanoid metabolism, and other select other pathways in exhausted cells

We next examined proteome dynamics at a global level to look for trends across pathways and biological functions. Hierarchical clustering of the global set of proteomic profiles pointed to 4 major clusters of protein expression throughout the progression to exhaustion (Fig 2E). These patterns were defined by: (1) “Activation/Intermediate Chronic,” wherein protein expression was elevated in early chronic stimulation time points (days 6 to 8) but down by late chronic stimulation/exhaustion as well as resting cultures; (2) “Chronic Activation/Exhaustion,” which corresponds to elevated expression in late chronic stimulation time points (days 10 to 12) relative to other time points, corresponding to more terminal/late exhaustion; (3) “Resting,” where protein expression is highest in resting cells, and then decreases by day 6, remaining low through day 12; and (4) “Resting/Chronic Activation Shared,” which shows an inverse pattern of the first cluster, with decreased protein expression at intermediate time points relative to resting and late chronic stimulation/exhaustion.

The “Activation/Intermediate Chronic” cluster was predictably enriched for biological processes and pathways linked to T-cell activation (Fig 2E), including cell cycle control, pyrimidine metabolism, translation, and mRNA transport that activated T cells deploy when undergoing proliferation following initial antigen exposure. The mTORC pathway was also active during this stage, promoting cell proliferation, survival, growth, and migration. mTORC, along with Myc, also controls cellular metabolism, prompting T cells to ramp up energy production via glycolysis, a critical early step during T-cell activation. Interestingly, we observed proteosome-dependent degradation, which has been linked to regulating T-cell metabolism and fate specification in CD8+ T cells [66].

In the “Chronic Activation/Exhaustion” cluster, we observed enrichment of cell cycle factors, in this case reflecting the marked down-regulation of proliferation in this cell population (for a comparison of the different proteins related to cell cycle that appear in “Chronic Activation/Exhaustion” and “Activation/Intermediate Chronic” clusters, see S1 Data). We also observed enrichment of genes involved in RNA degradation, which could be explained by RNA turnover being used as a means for turning off gene expression in cells undergoing exhaustion [67]. We also found supporting evidence that NFAT-mediated target gene expression is occurring, in agreement with the up-regulation of NFAT-related proteins in the proteome (NFATC1, NFATC2, NFATC3, NFAT5). Interestingly, chromatin remodeling PRC2 Complex, which mediates long-term gene silencing through modification of histone tails, was elevated in day 12 cultures, suggesting a regulatory role in controlling gene suppression during chronic stimulation. We also observed pathways related to cell death, which may drive a portion of the cells to undergo apoptosis under exhaustion-inducing conditions. Lastly, we found evidence of metabolic change, most intriguingly pointing to C20 Prostanoid Biosynthesis. Specifically, the C20 class of prostanoids are comprised of prostaglandins, thromboxanes, and prostacyclins, which have been implicated in the regulation of the immune response in T cells [68,69].

The “Resting” cluster was enriched for common T-cell processes. Some were, however, related to T cell activation, which may be caused by the fact these proteins can be down-regulated by the chronic time points to negatively control activation. These pathways include IL-2 signaling, cell–cell adhesion pathways, modulators of T-cell activation and TCR signaling, and some metabolic functions such as oxidative phosphorylation, a resting T cells primary energy source. Example proteins from these pathways that are up in resting and go down with chronic activation include IL2RB, CD3E and CD3G, LCK, CD4, and many OXPHOS proteins such as those in the ATP synthase complex like ATP5B/E/I. Notably, this observation suggests however that exhausted T cells are not using oxidative phosphorylation to the same extent that resting cells are, which is consistent with the findings that impaired mitochondrial function and OXPHOS is associated with T-cell exhaustion [70].

The last cluster, “Resting/Chronic Activation Shared,” represents proteins that are up-regulated in both resting and late chronic stimulation relative to the earlier stimulation time points and has the largest number of components. Proteins in this grouping are involved in cell maintenance and suppression of activation, proliferation, migration/cell–cell interactions, and other biological processes that are shared by the 2 phenotypes. Furthermore, we also found enrichment of signaling factors linked to CD4+ TCR pathway, calcineurin-dependent NFAT signaling, epigenetic chromatin organization, and HDAC-mediated signaling, further pointing to the importance of coordination of metabolic changes with epigenetic and transcriptional control through signaling cascades. We should note that evidence of calcineurin-dependent NFAT signaling pathway was found in both “Resting/Chronic Activation Shared” (CREBBP, MAP3K1, YWHAB, PRKCB, BAD, PRKCE, NFATC3, NFATC2, PRKCA, NFATC1, AKAP5, CABIN1, FKBP1A, PPP3CA, MAPK9, PPP3CB, CAMK4, CHP1, EP300, MEF2D, MAPK3) and “Chronic Activation/Exhaustion” (ITCH, ERG3, MAF, IFNG, IRF4, CASP3, IL2RA, SLC3A2, CBLB, PTGS2, TNF, JUNB). While the proteins contributing to the signals differed in each cluster, a majority of these proteins (from both clusters) are up by day 12, suggesting that day 12 is where the strongest signal comes from.

Posttranslational regulation of exhaustion factors and signaling pathways

To directly assess changes in signaling pathways arising from chronic stimulation, we performed phosphoproteomic analysis of the same samples to capture global protein phosphorylation changes during the progression of exhaustion. After stringent filtering, we quantified 11,728 unique high-confidence phosphosites on 2,685 phosphoproteins. A total of 7,499 phosphorylation events (on 2,100 proteins) were found to be significantly differentially regulated (FDR Fig 2F).

Reflecting the dramatic impact of TCR-mediated signaling events, most changes occurred between resting and day 6 cells, including phosphorylation of many TCR pathway components, cell cycles, and cytoskeletal proteins. We observed less significant phosphorylation changes between day 6 and 8, but increased phosphorylation of exhaustion- and CD4+ memory cell-associated transcription factor EOMES, as well as decreased phosphorylation of MAP kinases. Significant phosphorylation changes were observed between days 8 and 10 where increased EOMES phosphorylation was sustained. Changes in phosphorylation were minimal by day 12 compared to day 10, with only tens of phosphosites significantly changing.

Other notable T-cell regulatory proteins include NFATs, known to exert key transcriptional roles in T cells, which were highly phosphorylated in our data set, indicating that many sites significantly changed over time. For example, we detected and quantified 17 phosphosites on NFAT1 (NFATC2), some of which were significantly up-regulated in late exhaustion (day 12) (S99, S110, S243, S330, S757/759, T224) or pre-exhaustion (day 6) (S230, S236, S243, S757/759, S814, and T224). We also picked up 9 phosphosites on NFAT2 (NFATC1), some that have been functionally annotated and were significantly differentially regulated. Serine 245, significantly hyperphosphorylated in late exhaustion compared to pre-exhaustion, is a known substrate of kinases PKACA and PIM1, impacting intracellular localization, transcriptional activity, and cell motility. We also quantified 9 unique phosphosites on NFAT4 (NFATC3), of which S163, S165, and S207 were detected and have annotated function. However, these phosphosites were not significantly different between time points, except for S163, which is dephosphorylated in day 8 and day 10 versus resting, and S207, which was dephosphorylated on day 6 and 8 compared to resting. S163 is known to be phosphorylated and inhibited by JNK2 [71], affecting localization and transcription. S207 is known to be phosphorylated by CK1A, which also controls subcellular localization. Such observed phosphorylation patterns suggest that NFAT4 is cytoplasmic during the resting cell state but enters the nucleus during the early stages of chronic activation (day 6/8) before returning to the cytoplasm during the terminal exhaustion state (day 10/12).

We observed increased phosphorylation on proteins such as EOMES, PTPN14 (PEZ), NAB2, and AKAP5 in late exhaustion (day 12) compared to resting (Fig 2F). EOMES was hyperphosphorylated at many sites, especially in the T-box associated domain between serine 630 and serine 655. AKAP5, which interacts with important T-cell signaling proteins such as PKA, PKC, and calcineurin and inhibits IL-2 transcription by disrupting calcineurin-dependent dephosphorylation of NFAT [72,73], was observed to have 2 sites that were highly phosphorylated over time with chronic stimulation. Chronic antigen exposure induced pronounced phosphorylation of the phosphatase PEZ (Fig 2F), encoded by PTPN14, and a major regulator of adheren junctions and cell–cell adhesion. NAB2 is also seen highly phosphorylated, which known to enhance T cell function. As increased phosphorylation is observed in parallel with loss of effector function, it is possible this protein is inhibited by phosphorylation of these sites.

A subset of phosphosites with increased fold change and significance appears slightly separated from the main cluster (Fig 2F). Notably, these phosphosites tend to follow the same phosphorylation dynamics of up-regulation in all stimulation time points compared to resting but staying consistent from day 6 until 12. Corresponding phosphosites were on proteins enriched for pathways such as cell proliferation (G2/M phase control, mitotic spindle formation) and T-cell maturation (Enrichr, q value

Global phosphosite trends reveal pathway dynamics during chronic stimulation

To glean insights into signaling kinetics, we clustered each phosphosite profile, which revealed 4 basic expression patterns. The first cluster showed a phosphorylation pattern that was low in resting cells but stably hyperphosphorylated from pre- to late-exhaustion. This cluster corresponded to regulators of the cell cycle, transcriptional repression via RanBP2, control of gene expression by vitamin D receptor (VDR), histone methylation (epigenetic control), and YY1 activity, an important regulator of exhaustion [56] (S6D Fig). VDR, a nuclear receptor that is known to play an important role in T-cell function, regulates hundreds of genes [74]. Interestingly, we observed up-regulation of VDR upon T-cell activation by day 6, but only a slight increase by day 12 of chronic activation. Upon binding the active form of its ligand D3, VDR translocates into the nucleus of T cells and binds to promotor regions of PDCD1 (PD-1), TIM3, and TIGIT, inhibiting their expression, thereby decreasing the exhaustion phenotype [75]. Hence, the up-regulation we observed could act as an “exhaustion modifier,” supported by studies showing reduced disease severity and T-cell exhaustion in SARS-CoV-2 patients with higher levels of vitamin D [76,77].

In the next cluster phosphorylation peaks by day 6 to 8 and becomes dephosphorylated again by late exhaustion. This cluster corresponded to various biological responses expected to emerge from activated T cells such as cell cycle control, mRNA processing, mTOR signaling, pyrimidine biosynthesis and salvage, PLK1 signaling, and Aurora B signaling.

The subsequent cluster represents proteins heavily phosphorylated in resting cells, which were dephosphorylated by day 6, and generally stayed at low levels through to late exhaustion. These factors corresponded to IL-2 signaling (FYN, GAB2, CBL, etc.), the TCR pathway (LAT, FYN, MAPK1, SOS1, RAF1), PI3K/mTOR/AKT signaling (AKT1S1, CAMK4, MAPK1, MYC, etc.), protein secretion (GBF1, GNAS, GOLGA4, etc.), and Myc targets (SSB, TRA2B, TRIM28, U2AF1, etc.), all of which are known to be down-regulated during chronic stimulation and exhaustion [10].

Similar to the cluster in our proteomics data set, the last cluster generally followed a pattern of dephosphorylation by day 6, gradually restoring phosphorylation to resting levels by day 12. These components were enriched for chromatin remodeling enzymes, as well as components of T-cell activation and mRNA processing.

Differential kinase activity in chronically stimulated CD4+ T cells

To determine which protein kinases exerted the largest impact in driving these signaling events and, by inference, exhaustion progression, we leveraged knowledge of annotated substrates and phosphorylation motif information to infer specific kinase or kinase-family activity in our model.

First, we performed kinase-substrate enrichment analysis on our phosphosite data (Fig 2G). Relative to pre-exhaustion (day 6), we observed an increase in ATM, MAPKAPK2 (MK2), PRKG1, AKT1, and PRKDC kinase activity by late exhaustion (day 12), and a sharp decrease in CDK2, CDK1, MAP2K1, CDK6, CDK4, CDK9, HIPK2, PLK1, GSK3A, CDK5, and CSNK2A1 kinase activity (Fig 2G). Inhibition of the top hit, ATM, has already been shown to restore exhausted CD8+ T-cell function in an HCV-infection context [78]. Together with p38, MAPKAPK2 (MK2) signals in response to cell stress or TCR/CD28 engagement [79]; chronically stimulated T-cells experience both. AKT1 has many roles in T-cell function and development, but in late-stage differentiated T cells, inhibition of AKT1 has been shown to enhance tumor elimination [80,81]. Up-regulation of PRKDC, linked to DNA damage repair, may reflect increased cell stress experience during T-cell exhaustion. Conversely, kinases prominently down-regulated in terminally exhausted (day 12) cells were cell cycle-dependent kinases, consistent with concomitant decrease in cell proliferation.

Additionally, we used phosphorylation motifs that were differentially phosphorylated to predict kinase-class participation of these changes (Fig 2H). When we compared late exhaustion to pre-exhaustion, we observed a large percentage of phosphorylation events were caused by kinases GSK3, MAPK, WW domain CK1, CDK, and PKA, among others, providing further insight into kinase activity that takes places between pre- and late-exhaustion.

Expression trajectories implicate CD276, FLT-1, and SLC proteins in CD4+ T-cell exhaustion

To hone in on individual proteins (versus pathway level) that are likely contributing to the exhausted phenotype, we devised a partitioning strategy for sifting through the data with the aim of extracting physiologically relevant groupings. Since proteins with similar expression patterns often play roles in similar processes [82], we examined profile similarity to observe coherent patterns during progression of exhaustion. We first divided proteins by functional class implicated in T-cell exhaustion, specifically cell surface proteins and transcription factors. Next, we applied k-means clustering to group each type according to their expression patterns to identify clusters enriched for exhaustion-associated proteins to implicate other potential regulators exhibiting similar profiles. Notably, the resulting tight groupings revealed sets of highly correlated cell surface proteins and transcription factors that show low expression in resting cells and increasingly elevated expression throughout the progression of exhausted phenotypes reaching peak expression at day 12 (Fig 3A and 3B).

Fig 3. Hierarchical clustering of cell surface and transcriptional regulator protein expression reveals clusters enriched for known exhaustion-associated proteins.

(A) Cell surface protein cluster (left) enriched for exhaustion-associated proteins, with known (green) and novel (yellow) protein profiles highlighted (right). (B) Transcriptional regulator protein clusters (left) enriched for exhaustion-associated proteins, with known (green) and novel (yellow) protein profiles displayed (right). The data underlying this figure can be found in S1 Data. Figure created with BioRender.com.

https://doi.org/10.1371/journal.pbio.3002943.g003

Strikingly, clusters that were characterized by low expression in resting cells and increasing expression throughout progression of the exhaustion phenotype (Fig 3A and 3B) were significantly enriched for many known exhaustion-associated factors including inhibitory receptors PD-1, TIM-3, LAG-3, TIGIT, and ENTPD1, and exhaustion-associated transcription factors TOX, BLIMP-1, BATF, and NFATC1. Such observations implicate candidate novel associations not previously linked to T-cell exhaustion that may play a role as drivers or markers of CD4+ T-cell exhaustion based on the coherence of their expression profiles. We also observed a cluster that follows a similar profile in our kinase clustering, which contains a variety of kinases from different classes, including MAP kinases, CDKs, NEKs, and STKs, among others (S4 Data). When we performed an integrative analysis (using eXpression2Kinases [83]) to investigate whether these kinases or transcription factors may be regulating each other, we found that 2 transcription factors from the clusters shown in Fig 3B were enriched for targets in our kinase cluster (cluster 3 in S4 Data) that follows the same expression pattern; NFYB targeted 11 out of 21 kinases from this list, and TCF3 targeted 5 out of 21 kinases (S5 Data), suggesting these transcription factors may be playing a role in the regulation of kinases from this cluster.

Besides the most widely known receptors (PD-1, LAG-3, TIM-3, TIGIT), the cluster of cell surface proteins (84 genes total) also contained a few recently characterized inhibitory receptors (CD244/SLAMF4, SLAMF7, SIGLEC10), as well as many novel proteins that have not previously been associated with chronic stimulation/exhaustion (S6 Data). CD244, also known as SLAMF4, is a member of the signaling lymphocytic activation molecule (SLAM) family and has been shown to have both stimulatory and inhibitory activities [84,85]. However, in the context of chronic infections, it has been suggested that CD244 limits CD8+ effector responses [86]. Furthermore, this gene has also been previously associated with exhaustion in patients with chronic lymphocytic leukemia as a marker of T-cell exhaustion. However, cells displaying CD244 still produced some cytokine while still losing cytotoxicity and proliferation capacity [87]. The self-ligand SLAMF7 was recently implicated as a player in T-cell exhaustion, and its signaling was shown to phosphorylate STAT1/3, expression of multiple inhibitory receptors, and expression of transcription factors associated with T-cell exhaustion [88]. SLAMF7 signaling activity is mainly controlled by the presence/absence of cytoplasmic adaptor proteins SH2D1A/B, which follows the same expression pattern in our data as the receptor itself (SHD1A detected, significantly up-regulated by day 12, FDR 89,90] and was identified as being up-regulated on activated CD4+ T cells [91], but has not yet been directly linked to exhaustion. In our data set, this gene is not only up-regulated on pre-exhausted cells (day 6), but even more highly expressed by late exhaustion (day 12), suggesting that SIGLEC10 is involved in the exhaustion phenotype of these CD4+ T cells.

Among the proteins in the cell surface cluster were a few receptors/cell surface proteins that have been linked to dampening of T-cell responses in various aspects but have not yet been linked to exhaustion. CD276 (B7-H3) is a member of the B7 superfamily, which is a group of receptors that play modulatory roles during T-cell responses. This molecule, originally suggested to have both costimulatory and inhibitory effects in murine studies, was shown to be inhibitory in human T cells [92]. Although shown to down-regulate T-cell activation, proliferation, and cytokine production, CD276 has not yet been linked to chronic T-cell stimulation or exhaustion. We observed significant up-regulation over time, with resting cells expressing the least amount of CD276, pre-exhausted cells (day 6) expressing an intermediate amount, and exhausted cells (day 12) expressing the most, suggesting that this molecule may be up-regulated by CD4+ T cells in response to repeated antigen stimulation to dampen T-cell responses.

Another protein of interest in this cluster is Fms Related Receptor Tyrosine Kinase 1 (FLT-1). This receptor is a part of the vascular endothelial growth factor receptor (VEGFR) family, specifically bound and activated by VEGF-A, which is often secreted by tumors into the microenvironment and drives TOX-dependent T-cell exhaustion [93]. Since we observed up-regulation over time in exhausted CD4+ T cells following the same expression pattern as other inhibitory receptors associated with exhaustion, up-regulation of FLT-1 may contribute to this cell state.

Independent marker validation of in vivo phenotypes using human clinical samples

To investigate if the exhaustion patterns and markers seen in vitro also occur in vivo, we first analyzed the surface expression of CD276 and FLT-1 on T cells present in peripheral blood. Five HIV+ and 5 HIV- donor peripheral blood mononuclear cell (PBMC) samples were analyzed by spectral flow cytometry using validated antibodies specific for these factors (Figs 4 and S7). While expression of these 2 markers on CD4+ T cells in both HIV–positive and HIV–negative donors was low overall, we found that CD4+ T cells positive for either marker (CD276 or FLT-1) showed concomitantly significantly higher PD-1 expression, establishing both that these markers do indeed appear on cells in vivo, and implying that they may mark functionally exhausted T cells. We found this trend to be true regardless of HIV status, and we did not see a significant difference between HIV–positive and HIV–negative donors. However, this may be due to the effectiveness of modern HIV treatments, as these patients showed similar T cell subset frequencies and protein expression as the HIV–negative donors.

Fig 4. Flow cytometric analysis of HIV+ and HIV- primary T cells reveals increased PD-1 expression in CD276+ and FLT-1+ CD4+ T cells.

Mean fluorescence intensity of PD-1 in CD4+ T cells from both HIV+ and HIV- patients (grouped due to observing same effect in both) for CD4+ cells that are CD276+ or FLT-1+. The data underlying this figure can be found in the supplemental flow cytometry files uploaded to flowrepository.org. Figure created with BioRender.com.

https://doi.org/10.1371/journal.pbio.3002943.g004

Lastly, we observed an enrichment of solute carrier (SLC) family proteins (S6 Data), which mediate import and export of diverse small molecules (e.g., ions, peptides, nutrients) need for metabolism. As T cells are heavily regulated by metabolic changes [11,12,94], these transporters may contribute to exhaustion via modulation of metabolism, again hinting at the major metabolic shifts immune cells undergo.

Transcriptional regulators altered in chronically stimulated CD4+ memory T cells

Significant and long-lasting metabolic changes are often strengthened and controlled via transcription. Transcriptional control of exhaustion is also known to be critical and therefor the same clustering strategy as described for the cell surface proteins was next applied to transcriptional regulators to identify potential novel transcription factors or co-factors. We observed 4 major protein expression patterns (S7 Data). We again saw a promising cluster that displayed the pattern of low in resting and increasing during chronic stimulation (Fig 3B). Similar to observations of cell surface proteins clusters, this cluster was enriched for exhaustion-associated transcription factors (TOX, EOMES, NFACT1 and 3, BATF, NR4A1, EGR2). We also observed other exhaustion-associated factors in a cluster that follows a somewhat similar pattern, but expression comes up earlier and stays high, and the actual expression patterns of the factors found here show expression profiles that might fit into the previous cluster. The exhaustion-associated factors found here were BLIMP-1 (PRDM1), ID2, IRF4, and NFIL3 (Fig 3B). These results suggest that there may be novel proteins following similar expression profiles identified here.

One protein identified in this cluster was MAF BZIP Transcription Factor (c-Maf, MAF), a leucine zipper-containing transcription factor. While not as well studied as factors like TOX or BLIMP-1, the role of MAF has been investigated in CD8+ T cells, where its overexpression leads to repression of cytokines IL-2 and IFN-γ production, and an increased expression of T cell exhaustion genes, such as PD-1 [62,95]. In another study, MAF was suggested to even be a driver of exhaustion, where CD8+ T cells also had increased anti-tumor capacity when MAF was knocked out [96]. Lastly, MAF induced exhaustion-associated genes in mice, such as PD-1, TIGIT, TIM-3, and LAG-3, in cooperation with the known exhaustion transcription factor, BLIMP-1 (PRDM1) [97]. However, MAF is not up-regulated by day 12 in our second donor, conflicting with this evidence.

Another interesting transcription factor found in this cluster is transcription factor zinc finger E-box binding homeobox 2 (ZEB2), which is a member of the Zfh1 family of 2-handed zinc finger/homeodomain proteins. While its significance in CD4+ T cells has not been closely studied, ZEB2 is a transcriptional repressor and is involved in terminal differentiation of CD8+ T cells (effector and memory) [98,99]. It has been previously reported that ZEB2 is up-regulated in murine CD4+ T cells by day 7 during an LCMV infection [100]. While ZEB2 appears to be an indicator of terminal effector-like cells, ZEB2 may also be playing a role in the CD4+ T cell chronic infection dysfunctional phenotype, directly repressing both IL7R and IL-2, which are 2 indicators of exhausted CD4+ T cell phenotype. ZEB2 has also been shown to correlate with the expression of ID2 and T-bet, suggesting their involvement in a larger transcriptional network [100]. Further, ZEB2 has also been found in a gene signature differentially implicated in exhaustion versus memory, which was also grouped alongside TOX and EOMES [101].

Besides individual factors of interest, we also observed enrichment of factors suggesting signaling pathways and complex up-regulation, for example, downstream components of non-canonical NF-κB signaling. The downstream transcriptional program of non-canonical NF-κB signaling is driven by a heterodimer of RELB and p52 (NFKB2), both of which are highly significant in the same cluster as transcription factors such as TOX, EOMES, and NFATC1. However, we did not find the same significant up-regulation of canonical NF-κB signaling factors p50 (NFKB1), REL, and RELA. In non-canonical NF-κB signaling, p100 is a pre-cursor to p52, and it is designated for processing into p52 through phosphorylation. We observed up-regulated phosphorylation of serine 707 of this protein at day 12, but decreased phosphorylation during days 6 to 8, which has been annotated as being phosphorylated by GSK3b and promotes it degradation by FBXW7α. Such an observation perhaps suggests that this pathway is active during early exhaustion and may be primed but inhibited through GSK3 activity in late exhaustion. Such characteristics could also be part of the mechanism of inhibition by GSK3b and its ability to somewhat improve functionality of exhausted cells when inhibited [102,103]. Target genes of non-canonical NF-κB pathway have been shown to include BLIMP-1 (via RelB, albeit in tumor cells) [104], suggesting this gene may be an early regulator of the exhaustion transcriptional network.

Another observation in this cluster was the enrichment of subunits of the polycomb repressive complex 2 (PRC2), providing further evidence from our pathway enrichment results where PRC2 complex repressive activity was first identified. We detected 8 out of 9 members of this complex (JARID2, AEBP2, EED, EZH1, EZH2, SUZ12, RBBP7, RBBP4), all of which were significantly up-regulated by late exhaustion (day 12) when compared to resting (S7 Data). The PRC2 complex is an epigenetic modifier, which functions by methylating H3K27 histones on target genes, usually suppressing expression. EZH2 and the PRC2 complex have been well studied and are known to play a role in development, proliferation, and function of T cells, as well as potential immunomodulators [105]. Recently, it was found that PD-L1 expression on exhausted CD4+ T cells was regulated by PRC2 and SWI/SNF [106]. While the opposite was observed in tumor cells, with PRC2 up-regulation leading to down-regulation of PD-L1, it is interesting to see that anti-CLTA-4 immune checkpoint inhibitor therapy is actually improved when given in conjunction with EZH2 therapies [107].

Network model of dysregulation in exhausted CD4+ memory T cells identifies p300 as a key node driving phenotype

To identify key drivers of the chronically stimulated dysfunctional phenotype we defined molecularly, we incorporated both our proteomic and phosphoproteomic data sets by using a causal network analysis across time that allows us to predict proteins in our model that are responsible for large changes of expression or phosphorylation.

Using known and annotated interactions such as protein–protein interactions or phosphorylation events, we examined the networks of differential phospho/proteins in each time point which factored in the previous comparison starting from resting cells. We solved a Steiner tree problem, connecting the largest number of differential features (leaves) to upstream regulators (connectors). In our final network (Fig 5A), which represents late exhaustion (day 12) of chronic stimulation, we observed several clustered subnetworks, with the largest revolving around p300 (EP300) in the top left corner (Fig 5B) with connections to a few important players in the T-cell exhaustion phenotype, including transcription factors EOMES, BATF, TBX21, and the cytokine IFN-γ. Furthermore, p300 is connected to PTPN2, a phosphatase known to inhibit T-cell activation by dephosphorylating key signaling proteins such as LCK, further suggesting p300 may play a role in T-cell function. In our proteomic data set, p300 is also significantly increased over time, peaking at day 12 chronic stimulation (Fig 5B).

Fig 5. Causal network analysis predicts p300 as regulator of exhaustion phenotype and p300 inhibition leads to down-regulation of inhibitory receptors PD-1 and TIM-3.

(A) Causal network analysis leveraging annotated interactions (transcription factor-target, kinase-substrate, phosphatase-substrate) as a framework to evaluate protein expression changes (late exhaustion (day 12) versus resting, taking into account previous time point comparisons) to predict key nodes (regulators) impacting most changes seen across the data set. (B) Pruned version of a large subnetwork with p300 (black outline) at main hub. Orange circles mark proteins of interest mentioned in main text. (C) Depiction of experimental design of p300 inhibition studies. (D) Flow cytometric analysis of PD-1 and TIM-3 showing shifts in fluorescent intensity of cell populations in treated versus untreated (cell number normalized to mode). Cell numbers shown for each plot are 20,312 for p300 treated and 30,438 for DMSO, which corresponds to all CD4+ T cells from the sample (excluding doublets and debris, see Methods for gating), corresponding to >95% of sample. (E) Bar chart of mean fluorescent intensity showing significant changes (T test) in PD-1 and TIM-3 during p300 inhibition versus controls (antigen removed for 2 days, or treatment with T-cell activation inhibitor Dasatinib); *, p-value p-value S8 Data and S9 Data (A, B), supplemental flow cytometry files uploaded to flowrepositroy.org (D), and the S3 Data (E, F). Figure created with BioRender.com.

https://doi.org/10.1371/journal.pbio.3002943.g005

To test if p300 is in fact playing a role in the phenotype, we inhibited p300 using A-485, a small molecule inhibitor. Specifically, cells were chronically stimulated for 10 days, and then treated with drug or vehicle for 2 days alongside the final 2 days of stimulation (Fig 5C), and then analyzed cells for basic exhaustion-associated phenotype using flow cytometry. Cells treated with p300 inhibitor showed significantly reduced expression levels of inhibitory receptor and exhaustion markers PD-1 and TIM-3 (Fig 5D and 5E). While cells still had expression of these receptors, cells treated with p300 inhibitor overall had on average significantly lower expression, and increased number of cells with low expression (Fig 5D). To better quantify this change, we looked at mean fluorescent intensity (MFI) values of drug treated versus DMSO treatment and observed significant decrease in receptor expression (Fig 5E). As 2 controls, we treated cells with Dasatinib to block T-cell activation or removed beads from day 10 until day 12 (rested for 2 days). Since p300 is known to play a role in cell cycle, we wanted to ascertain if the effects seen with p300 inhibition were not solely a consequence of blocking T-cell activation. If this was the case, we would expect PD-1/TIM-3 expression of p300-inhibitor treated cells to be equivalent to either resting or treating with Dasatinib, which was not the case. Further evidence that p300 inhibition was not blocking T-cell activation was that activation marker CD137 was unchanged in p300 inhibitor treated cells versus DMSO but was decreased in rested and Dasatinib conditions (S8 Fig).

P300 is a histone acetyltransferase cofactor, which imparts its effects via histone remodeling, usually activating bound genes. In our p300 subnetwork, we also observed that this gene has connections to HDAC2, which is a class I HDAC, along with HDAC1 and 3, suggesting further epigenetic regulation. T-cell fate and function is also known to be strongly controlled by epigenetic regulation, prompting us to investigate which genes p300 and other similar epigenetic-modifying cofactors are regulating.

Control of transcriptional changes influenced by epigenetic co-factor binding dynamics

Epigenetic control is mediated by complexes of regulatory cofactors and transcription factors. Our analyses have identified roles for epigenetic cofactors (e.g., PRC2, p300, HDACs) and TFs, therefore, we sought to determine whether TF-cofactor interactions might further indicate key regulatory complexes driving T-cell exhaustion. To assess recruitment of epigenetic modifiers by sequence-specific DNA-binding factors, we performed biochemical assays using the CASCADE (Comprehensive ASsessment of Complex Assembly at DNA Elements) microarray platform [108]. In CASCADE, recruitment of cofactors by TFs is assayed by profiling the recruitment of cofactors to double-stranded DNA microarrays comprising hundreds of annotated TF-binding site sequences. We applied nuclear extracts from acutely stimulated (2 days) and 12-day chronically stimulated cells to the microarrays and measured cofactor recruitment using fluorescently labeled antibodies. Acutely stimulated cells were used here as a control here to allow us to see how complexes form during healthy T-cell challenge compared to that of exhausted cell challenge. We determined multiple cofactor recruitment motifs for each cofactor, which we used to infer the TF family to which a given cofactor was recruited by searching against TF motif databases (see Methods). Alongside p300, we opted to probe both KAT2B, a histone acetyltransferase known to associate with p300 that has been implicated as a marker in T-cell exhaustion [109,110], and a mixture of class I HDACs (HDAC1/2/3), which have also been linked to T-cell exhaustion [111–113] and were connected in our p300-centered subnetwork (Fig 5B).

Generally, we observed changes in recruitment for all 3 profiled cofactors between acute and chronically stimulated T cells, with p300 exhibiting the most pronounced changes in recruitment by TF families (Fig 5F and S10 Data): P300 (14 families), HDACs (10 families), KAT2B (10 families). Consistent with loss of TCR signaling in the exhausted state, we observed loss of cofactor recruitment to IRF/STAT and NF-kB family (REL, NFKB) TFs. Notably, we also observed coherent enhanced recruitment of all 3 cofactors to the TBX family, which contains the major exhaustion regulator EOMES. Strikingly, however, the strongest recruitment signal for p300 in exhausted cells is the NFAT family.

Furthermore, we examined changes in TF-cofactor interactions related to inhibitor receptor expression. We observed increased recruitment of HDACs to PKNOX/TGIF family, both known regulators of PD-1 [114,115]. We also observed p300 recruitment to multiple nuclear receptors families (including RXRA) that are implicated as regulators of TIM-3 [116]. We likewise observed increased binding of KAT2B to another nuclear receptor group, including family member NR4A1 reported to regulate TIM-3, and possibly other inhibitory receptors [117–119]. Increased recruitment of HDACs to the MEIS/PBX family was also seen in the exhausted cells, where MEIS2 is predicted to control TIM-3 [117].

Lastly, we observed enhanced recruitment of p300 to EBOX family in chronically stimulated cells. Of the 3 TFs in this motif family (XBP1, CREB3, CREB3L1), our proteomics data had shown significant changes in XBP1 protein abundance in exhausted cells. Interestingly, XBP1 has been shown to suppress antitumor activity of T cells [120], promote PD-1 expression, and contribute to translating metabolic remodeling to changes in gene expression [12]. Specifically, XBP1 regulates mitochondrial activity, and can be activated by a decrease in glycolysis, as seen in exhaustion.

Together, these results demonstrate changes in TF-cofactor complexes as T cells progress to the exhausted state, supporting our proposed role for p300 in concert with that of known regulators (e.g., EOMES) that likely dictate the metabolic shifts exhibited by exhausted CD4+ T cells.

Chronic stimulation of CD4+ memory T cells causes cell-intrinsic metabolic remodeling

Our proteomic profiling and network modeling revealed changes in metabolic pathways, as discussed in our pathway enrichment and network analyses results. Due to the large role metabolism plays in T-cell function and differentiation, we performed a systematic metabolic enrichment analysis. Such measurements quantify potential metabolic alterations that associate with chronic stimulation under nutrient rich culture conditions, allowing us to characterize intrinsic cellular responses rather than secondary changes due to changing nutrient availability. Applying the network inference tool MOMENTA [121] to our proteomic data, we observed significant up- and down-regulation of key metabolic pathways when we compared late exhausted cells to resting cells (Fig 6A). We observed changes associated with the electron transport chain and glycolysis in late exhaustion (day 12), providing further evidence of down-regulation of these energy sources in a nutrient-independent manner.

Fig 6. Metabolic pathway changes in late exhaustion (day 12) chronically stimulated CD4+ T cells.

(A) Metabolic pathway enrichment predictions of proteomic profiles obtained from late exhaustion (day 12) stimulated cells versus resting cells. (B) Heatmaps of expression levels of proteins annotated to individual predicted pathways. Rows z-score normalized. (C) Heatmap of top 1% differential metabolite features measured between late exhaustion (day 12) stimulated and resting cells. Rows z-score normalized. The data underlying this figure can be found in S3 Data (A, B) and S1 Data (C). Figure created with BioRender.com.

https://doi.org/10.1371/journal.pbio.3002943.g006

We also observed several pathways pointing to increased, possibly dysfunctional, fatty acid metabolism. Mitochondrial fatty acid oxidation itself is not significantly up-regulated; however, peroxisomal FAO is, along with carnitine shuttling (a key step in importing fatty acids into the mitochondria) and metabolism of omega3 and omega6 fatty acids. Peroxisomal FAO serves a different biological function that mitochondrial FAO; peroxisomal FAO is required to break down very long chain fatty acids, prostaglandins, leukotrienes, and other more complex fatty acids. While carnitine does not play a role in peroxisomal FAO, it is necessary to transport end products into the mitochondria for use in the citric acid cycle and OXPHOS. Another difference is where spare electrons are fed during oxidation, which normally goes into the ETC, but peroxisomes donate theirs directly to oxygen (O₂) in the cell, which creates reactive oxygen species (ROS). It is known that exhausted T cells often have increased oxidative stress, usually associated with mitochondrial dysfunction; however, this new mechanism may also be contributing. One common way cells adapt to this oxidative stress is through production of antioxidating pathways, such as selenoamino acid metabolism (Fig 6B), which is seen up-regulated and known to impact T-cell function via its antioxidant properties [122].

Other top pathways include phosphatidylinositol phosphate and hyaluronan metabolism. Phosphatidylinositol phosphate metabolism, known for generating second messengers with important roles in T-cell function and development [123], plays an important role in the PI3K pathway driving T-cell activation. The down-regulation we observed suggests a possible mechanism underlying exhausted T-cell hyporesponsiveness. Hyaluronan synthesis, also down-regulated, is required for IL-2-mediated cell proliferation [124]. IL-2 secretion by activated CD4+ T cells is rapidly down-regulated by chronic stimulation in our culture model (Fig 1C). The down-regulation of hyaluronan would make chronically stimulated cells less responsive to proliferation-inducing cytokines like IL-2.

Consistent with FAO modulation, the most significant hit up-regulated in late exhaustion (day 12) relative to pre-exhaustion (day 6) is lysine metabolism, which along with methionine, is needed for carnitine biosynthesis, as well as 3-oxo-10R-octadecatrienoate beta-oxidation. These results again point to exhausted cells utilizing FAO as a preferential energy source, as previously suggested [11].

We also observed ROS detoxification, which is expected and shares some overlap in functionality with selenoamino acid metabolism. We also observed new pathways emerge, including prostaglandin formation, further encouraging our observation of “C20 Prostanoid Biosynthesis” from our cluster enrichment analysis, and possibly related to peroxisomal FAO, where prostanoids are broken down for energy. Specifically, we see prostaglandin formation from linoleic acid (linoleate), and we also observe up-regulated linoleate metabolism. Inversely, we observed significant pathways that were down in day 12/up in day 6, including glycolysis, pyruvate metabolism, and metabolic pathways of nucleic and amino acids, likely as day 6 cells were still activated and proliferating.

Direct measurement of metabolite dynamics verifies metabolic remodeling

The putative large-scale alterations in metabolic pathway activity that were inferred by our functional enrichment analysis prompted us to verify these predicted changes using an independent experimental approach. To this end, we performed global metabolic profiling using untargeted mass spectrometry to directly assess changes in the levels of key metabolites in these same cell populations (see Methods). Large changes in the steady-state levels of endogenous cellular small molecules were observed over the exhaustion time course (Fig 6C), reflecting profound changes in cellular metabolism.

Consistent with our proteomic predictions, metabolite level functional enrichment confirmed our MOMENTA predictions, as well as implicated other metabolic alterations. Comparing our differential metabolite patterns to the human MFN metabolomic model (curated combination of KEGG, BiGG, and the Edinburgh model) [125], we observed significant perturbations (p-value

Hence, our orthogonal approach strongly supports and further refines our findings that chronically stimulated cells intrinsically down-regulate glycolysis and ETC, while breaking down longer/more complex fatty acids through peroxisomal FAO, and up-regulating pathways countering the increased oxidative stress from this altered metabolic state.

Metabolic pathway regulation in inhibitory-receptor positive HIV-donor T cells reflect metabolic shifts seen in vitro

To validate that the changes observed in the in vitro culture system mirror what occurs in vivo, we obtained T cells from 4 HIV+ donors, which were pooled to obtain enough material for downstream analyses. Cells were FACS-sorted based on inhibitory receptor (IR) expression (PD-1+ and TIGIT+) and CD4/CD8, and bottom-up proteomics was performed. Due to low cell count post-sort, proteomic coverage was shallow, again highlighting the need for the cell culture model. Yet while most exhaustion-associated proteins were not detected, a few (NFATC2, RUNX3, TBX21) were found to be up-regulated in CD8+ IR+ T cells, while others such as DMNT1, were prominently elevated in CD4+ IR+ T cells.

Although certain exhaustion-associated proteins were not directly detected, pathway enrichment analysis using MOMENTA showed high agreement of top metabolic pathways up-regulated in HIV+ IR+ T cells with those altered in day 6 and day 12 (versus resting) in our in vitro model (Table 1). These include elevation in pathways related to adenine and adenosine salvage pathways and retinoate biosynthesis.

We also investigated other pathways enrichment encompassing the top 100 most significant factors up-regulated in HIV+ IR+ T cells and found that components related to prostaglandin synthesis, fatty acid metabolism, and p53 pathways were all significantly enriched in exhausted T cells in vivo (Table 2), again matching what was seen in our in vitro exhaustion model.

Up-regulation of exhaustion and identified genes in breast cancer tumor lymphocytes in vivo

To ensure the exhaustion protein expression we observed in our in vitro model also occurs more generally in CD4+ T cells in vivo, we accessed a publicly available single-cell RNA sequencing dataset and plotted marker expression (using Single Cell Portal [126]) in tumor infiltrated CD4+ and CD8+ T lymphocytes displaying an exhausted phenotype (Fig 7A). Strikingly, in this breast cancer data set [127], we observed up-regulation of a significant number of exhaustion-associated factors detected in our model (Fig 7B), most notably in a putative exhausted CD4+ subpopulation demarcated by GNG4 expression and increased inhibitory receptor expression (Fig 7A, left red box on UMAPs). Inhibitory receptors showed similar levels to that of exhausted CD8+ T cells (Fig 7A, right red box on UMAPs). Receptors such as PD-1, TIGIT, and LAG-3, as well as exhaustion-associated transcriptional factors TOX, BLIMP-1 (PRDM1), and NFATC1, among others, were up-regulated in vivo in agreement with our model.

Next, we looked to see if transcripts corresponding to proteins identified in our analysis were likewise expressed in the same CD4+ cell subpopulation. Indeed, we found that despite low general global correlation between mRNA and protein levels that many of the novel markers we identified in our in vitro model were also up-regulated in the tumor resident T cells, including p300, ZEB2, VDR, and many SLC transporters (Fig 7B).

Overall, this data suggests that many of the proteins we identified in the chronic T-cell stimulation model are likewise elevated in vivo in exhausted CD4+ (and often simultaneously in exhausted CD8+) T cells in a tumor context.

Study limitations

While our culture model allowed for in-depth multiomic profiling and nutrient-independent observations, it does have some limitations. The in vitro conditions (nutrient availability, cell–cell interactions, etc.) and method of continuous TCR activation (via CD3/CD28 Dynabeads) are artificial and do not fully replicate what immune cells experience in vivo. However, the up-regulation of many known exhaustion markers and regulators strongly supports the utility of the model. In order to more closely replicate physiological conditions, human serum was used in place of fetal bovine serum. We also recognize these results were largely reproduced in a small number of human donors, and so may not generalize more broadly due to genetic and environmental variability.

Another confounder is cellular heterogeneity, given that bulk omics methods were used here, as T-cell populations are often heterogeneous. In particular, we observed evidence of variability by flow cytometry based on co-expression levels of inhibitory receptors. Although our quantitative mass spectrometry-based approach measures the average expression across all the subpopulations, we believe our results captured trends that are representative across the entire population. With advances in proteomics, phosphoproteomics, and metabolomics methods, the ability to further hone in on specific subpopulations will one day be feasible and make for interesting follow-up studies.

We also acknowledge that our phosphoproteomics method detected primarily phosphorylation sites that have not yet been functionally annotated. While lack of curated functional significance and incomplete kinase-substrate knowledge are current bottlenecks in this field. However, the data it is still insightful as signaling changes were predicted at the pathway level, and observed enriched phosphorylation of proteins in the same system serves as an indicator of differentially regulated activity.

In closing, this work provides a massive resource of multiple data types temporally profiling the effects of CD4+ T-cell exhaustion/chronic stimulation for the broader research community. Our data allows comparison to existing results and ongoing in vivo studies, points to potential new targets and avenues of research, including the identification of regulatory factors that can potentially be modulated to either reduce exhaustion or potentially even to induce the exhausted state in hyperactive contexts such as autoimmune disease. With the major findings of our study and the limitations in mind, this resource opens the door for follow-up studies. Further research into the mechanisms of control of p300 during chronic T-cell activation should permit dissection of the relationships to other important upstream and downstream proteins. Such analyses could identify new druggable targets with a more specific mode of action as targeting p300 itself is detrimental to normal cell function due to cell-cycle interference. Likewise, as metabolism dictates T-cell function, deploying tool compounds to investigate the impact of the metabolic changes we detected, most notably preferential alterations in peroxisomal fatty acid oxidation, and the associated regulators driving these alterations, could reveal actionable targets.

Human primary PBMC isolation

Blood was obtained through New York Biologics, Inc. It is collected from de-identified healthy donors and shipped refrigerated directly to our lab. The blood is tested for potentially dangerous viral or bacterial infections. PBMCs were isolated from blood packs using density centrifugation with Histopaque-1077 (Sigma-Aldrich, 10771). PBMCs are washed multiple times with PBS before being frozen in 10% DMSO + 90% HI-FBS (HyClone) via Mr. Frosty container placed at −80°C overnight, and then transferred to long-term liquid nitrogen storage. All data produced from these samples were analyzed anonymously.

T-cell isolation

Frozen PBMCs were thawed quickly at 37°C in a water bath until just an ice core remained. Warm media was added and then returned to a 15 ml Falcon tube until PBMCs were thawed, and then spun down. After aspiration, cells were resuspended in PBS with 1% BSA (Thermo, 15260037) and 50 mM EDTA, and T cells were isolated using untouched memory CD4+ T-cell isolation kit (Miltenyi, 130-091-893). T cells were then plated on appropriately sized tissue culture plates in AIM-V (Thermo, 12055091) with 5% heat-inactivated human AB serum (Valley Biomedical, HP1022HI) and incubated at 37°C (5% CO₂) for 2 h to remove any contaminating monocytes. Cells in suspension were then collected and used for expansion or experiment.

T-cell expansion

To obtain enough cells/material for mass spectrometry-based approaches, T cells were expanded before use. After T-cell isolation, T cells were plated in AIM-V with 5% heat-inactivated human AB serum (hereby referred to as T-cell media) at a density of 1 × 10⁶ cells per ml. For CD4+ T cells, 30 U/ml of IL-2 was added, and 25 μl of ImmunoCult (CD3/CD28, Stemcell Technologies) was added per million cells. On day 3, cells were diluted approximately 4- to 8-fold for a concentration of 1–2.5 × 10⁵ cells per ml. On day 5, cells were again diluted to a concentration of 1–2.5 × 10⁵ cells per ml. On day 7, cells were diluted to as close to 4-fold as possible to achieve between 1–6 × 10⁵ cells per ml. Cell counts and viability were then monitored each day until day 9 or 10 and were harvested when proliferation rate and viability began to decline. These cells were then used for experimental procedures.

T-cell in vitro exhaustion

Once isolated and expanded, T cells were plated at 1 × 10⁶ cells per ml, with IL-2 added where specified. CD3/CD28 Dynabeads (Thermo, 11132) were added at a 1:1 cell to bead ratio. On day 2, cells were collected, spun down, resuspended in 1 ml of T-cell media and beads were removed using a magnet. Beads were then washed 2 more times, each with more agitation than the last to dissociate any cells stuck to beads. For “resting” condition samples, cells were then plated in fresh medium at 1 × 10⁶ cells per ml, and for “chronic” stimulation samples, cells were re-plated at the same density, but had fresh Dynabeads added at 1:1. Media and fresh beads were then added every 2 days to the chronic stimulation samples until end of time course (day 12).

Small molecule inhibition cultures

Cells were cultured through the in vitro exhaustion protocol described above until day 10. After the day 10 media and antigen refresh, cells were treated with 1 to 10 μm of A-485 p300 inhibitor, or with Dasatinib to block T-cell activation. “Rested” condition cells were also plated without Dynabeads to allow for 2 days of rest before analysis. All conditions were then cultured for 2 more days before analysis via flow cytometry.

Cytokine release assay

On time points where cytokine release was tested, 5 × 10⁵ cells were plated in triplicate per replicate in 1 ml of media without IL-2 in a 12-well plate. CD3/CD28 Dynabeads were added at 1:1. Cells were incubated for 6 to 10 h before media was collected and frozen at −80°C. Once all time points were acquired, media samples were thawed and analyzed via ELISA (Thermo, BMS221INST, KHC3011, EHIFNG) or Proteome Profiler Human Cytokine Array Kit (R&D Systems, ARY005B).

Mass spectrometry sample preparation

Samples were taken from −80°C storage and immediately lysed using guanidine hydrochloride lysis buffer (6 M guanidine hydrochloride, 100 mM (pH 8.5) Tris, 40 mM chloroacetamide, 10 mM TCEP, PhosStop). Bradford was performed on each sample to obtain protein concentrations. Protein mixtures were then enzymatically digested using trypsin protease at a 1:50 to 1:100 μg of trypsin to protein ratio. Digestion was performed overnight (16 h) at 37°C with gentle shaking. Samples were then acidified using final concentration of 1% formic acid to stop digestion. Each sample was then desalted using Sep-Pak desalting columns and dried down using a speed-vac before peptide concentration was measured using Pierce peptide quantification assay (Thermo, 23275). Approximately 100 μg of each sample was then labeled using TMTpro 16-plex kit (Thermo, A44520). Samples were then pooled after labeling and dried down in a speed-vac. Next, the pooled samples were resuspended in 2% acetonitrile, 0.1% ammonium hydroxide and loaded onto the HPLC (Agilent 1100). Reversed-phase chromatography was performed with a 53-min gradient, where 48 fractions were collected starting after 5 minutes of loading. Every 12th fraction was then pooled, and dried down, for a final total of 12 fractionated samples; 5% of these were taken and used for proteomic analysis, and the rest was used for phosphopeptide enrichment and subsequent phosphoproteomic analysis.

Phosphopeptide enrichment

Fractionated and dried samples were thawed after storage at −80°C on the bench for 5 min. Each fraction was then resuspended in 800 μl of sample buffer (80% acetonitrile, 0.5% trifluoroacetic acid) and spun at 1,600 rpm on a temperature-controlled shaker at 20°C for 15 min. Resuspended fractions were then loaded onto 96-well plates and placed in the KingFisher (Thermo Fisher) sample purification system. Six other plates were prepared and placed into the KingFisher alongside the sample plate: (1) Plate with comb tip to cover magnetic rods. (2) Bead plate, containing washed iron beads (10 μl beads/100 μg of peptide, Cube Biotech, 31505-Fe) in wash buffer (80% acetonitrile, 0.1% trifluoroacetic acid). (3) First wash plate containing 800 μl of wash buffer. (4) Second wash plate containing 800 μl of wash buffer. (5) First elution plate containing 200 μl of elution buffer (50% acetonitrile, 7.5% ammonium hydroxide). (6) Second elution plate containing 200 μl of elution buffer. The KingFisher was programmed to pick up the iron beads using the magnetic rods and incubate them with sample for 20 min with agitation. Beads (now with bound phosphopeptides) were then transferred to the first wash plate and incubated for 5 min with agitation, and then again for the second wash plate. Finally, beads were transferred to the elution plates and incubated for 5 min with agitation each. The 2 eluates for each fraction were then collected and pooled before drying down in a SpeedVac.

Metabolite extraction and clean-up from cellular supernatant

If samples were intended for metabolomics in parallel with proteomics, samples were not immediately lysed with guanidine lysis buffer as described above. Instead, samples were snap frozen after washing with PBS in a chemically resistant microcentrifuge tubes (e.g., Eppendorf). A mixture of pre-chilled (−20°C) methanol/acetonitrile/water (40/40/20, v/v) was prepared for liquid–liquid extraction (LLE). Next, we added buffer at 3 to 5× volume of sample. Samples were then vortexed for 30 s and snap frozen by placing in liquid nitrogen until frozen (1 min). Samples are thawed on ice and sonicated for 3 cycles of 5-min pulses using a probe sonicator (40 kHz, Ultra Autosonic) at 10% power. The samples were then incubated at −20°C for 1 h. Samples were next centrifuged for 15 min at max speed (14,000 rpm) at 4°C. The supernatant containing metabolites were transferred to a new tube, while protein (precipitate) was kept at −80°C until ready, at which point guanidine lysis buffer was added and brought through the proteomic workflow. Metabolite supernatant was dried in a speed-vac at 30°C prior to solid phase micro-extraction (SPME) clean up.

Dried metabolite extracts were resuspended in PBS and transferred to a 96-well plate. The 96-blade SPME system was next washed in ethanol:water (70:30, v/v) for 30 min and preconditioned for 30 min in methanol:water (50:50, v/v). Samples were then extracted for 1 h, and the blades with bound metabolites were washed for 20 s in water (MS grade, Thermo Fisher) and metabolites desorbed in acetonitrile:water (50:50,v/v) for 1 h. The desorption solvent was then evaporated in a speed-vac until dry and stored at −80°C until ready for analysis prior to LC/MS, at which point it was resuspended in 2% acetonitrile and separated on the nano-liquid chromatography mass spectrometer (nLC-MS/MS).

LC-MS analysis of peptides and phosphopeptides

Fractionated, multiplexed samples were resuspended in 2% acetonitrile and 0.1% formic acid (mobile phase A) and transferred into a 96-well plate to be loaded into our Easy nanoLC1200 HPLC system and analyzed by our Exploris 480 with FAIMS Pro (Thermo Fisher). The peptide solutions were first loaded onto a reversed-phase nanotrap column (75 mm i.d. 3 2 cm, Acclaim PepMap100 C18 3 mm, 100 A°, Thermo Fisher) and separated by an EASY-Spray column (ES803A, Thermo Fisher) using a gradient (6% to 17% over 77 min, then 17% to 36% over 45 min for peptides, 6% to 19% over 58 min, then 19% to 36% over 34 min for phosphopeptides) of mobile phase B (0.1% formic acid, 80% acetonitrile) at a flow rate of 250 nL/minute. Positive ion mode was used throughout, with a capillary temperature of 275°C and a spray voltage of 2,100 V. Our DDA method chose the 12 most abundant ions for fragmentation (Exploris 480 NCE 33%), with FAIMS cycling through −50, −57, and −64 volts for each full scan. Precursor scans were acquired at a 120,000 FWHM resolution with a maximum injection time of 120 ms in the Orbitrap analyzer. The following 0.8 s were dedicated to fragmenting the most abundant ions at the same FAIMS compensation voltage, with charge states between 2 and 5, via HCD (NCE 33%) before analysis at a resolution of 45,000 FWHM with a maximum injection time of 60 ms. Phosphopeptides were analyzed the same besides the shortened gradient described and increased injection time of 150 ms.

LC-MS analysis of metabolites

LC-MS was performed using an Orbitrap Q-Exactive Exploris 480 mass spectrometer (Thermo Scientific). Metabolites (2 μls) were loaded onto a pre-column (75 mm i.d. × 2 cm, 3 μm) and then separated on a C18 capillary column (75 mm i.d. × 2 cm, 2 μm, 100 Å, Thermo Fisher Scientific) over a 45-min gradient. Mobile phase A was 2% acetonitrile and mobile phase B was 80% of acetonitrile. The gradient consisted of 2% to 60% mobile phase B for 20 min, was increased to 95% mobile phase B over 10 min, and maintained at 95% mobile phase B for 15 min with flow rate of 300 nL/min. The MS was operated in switching mode to acquire data on both positive and negative mode in the same run. The spray voltage was set to 2.1 kV and −1.8 kV, respectively. Full scan covers m/z 60 to 1,000 at a resolution of 60,000. The AGC target was set to 300% and the maximum ion injection time was set to 25 ms. MS2 scans were performed at 15,000 resolution with a maximum injection time of 64 ms using stepped collision energies of 10, 20, and 40. Dynamic exclusion was enabled using a time window of 10 s.

Prize-collecting Steiner forest

The Steiner forest problem is a method of network optimization in a weighted graph G(V, E) that is composed of a node set V and edge set E. Furthermore, we consider a function p(v) that assigns a prize to each node v ∈ V and a function w(e) that assigns a weight to each edge e ∈ E. Formally, we seek a forest of trees F(V_F,E_F) that minimizes the objective function

where c(e) = 1 –w(e), β is a scaling factor that affects the number of prized nodes included in the optimal forest, μ is a parameter that penalizes hub nodes (nodes with high degree), κ is the number of trees in the forest, and ω is a parameter that controls the number of trees. Prior to computing the optimal forest, a dummy node is attached to a subset of the nodes in G. Once the optimization is complete, the dummy node and all of its artificial edges are removed to reveal a forest with each tree in the forest rooted at a node that was connected to the dummy node. We used the Omics Integrator package to solve the Steiner forest problem [138]. As the Steiner forest algorithm requires an input network, we assembled a directed molecular network capturing transcription factor–target, kinase–substrate, phosphatase–substrate as well as other interactions from Reactome [139], Kegg [140], Signor [141], and ligand–receptor interactions [142].

For each interaction e in the input network, the Steiner forest problem requires the assignment of a weight or probability w(e), suggesting that edges with low cost c(e) = 1 –w(e) are more likely to be selected in the optimal forest. Given a directed edge e = (x, y), where the node x is the tail or the source of the interaction and node y is the head, we defined that the weighting function is the reciprocal of the outdegree of x as . Consequently, edges that involve tail genes with high outdegree will be more likely to be removed during the Steiner forest optimization. Such a step allowed us to more effectively penalize hub nodes, indicating higher confidence that edges selected in the optimal forest are more specific to T-cell exhaustion instead.

As for a set of priced genes, we considered genes that were significantly differentially expressed or phosphorylated in each time step after stimulation if their corresponding FDR lg₂ of the corresponding fold changes in a given time point t compared to the unstimulated case, .

In addition, Omics Integrator requires the user to assign values to certain parameters that affect the topology of the optimal forest. The parameters of interest are μ, β, ω, and D. We assigned a value of 0 to parameter μ because we found μ to be too punitive, and the resulting optimal forest would either be empty or leave out all but a couple of nodes with non-zero prizes, even for small values μ~10⁻⁴. In addition, our edge weight-assignment in (1) already punishes hub nodes that might be overrepresented in the literature. Furthermore, we set β = 1 and ω = 1. Lastly, we assigned a value of 5 to D. The values chosen for μ, β, ω, and D were constant across time points because any differences in network size and structure among the optimal forests could then be attributed to the experimental data and the input networks and not to different parameter values.

Nuclear extraction

The nuclear extract protocols are as previously described [108,143,144], with modifications as detailed below. Nuclear extractions were performed on day 2 of the exhaustion protocol for acutely stimulated cells and day 12 of the exhaustion protocol for exhausted cells. To harvest nuclear extracts from the stimulated CD4+ memory T cells, cells were collected with media in 50 ml conical tubes and centrifuged at 500×g for 5 min at 4°C. After centrifugation, the supernatant was aspirated off leaving behind the cell pellet. Cell pellets were resuspended in ice-cold 1× PBS (Cytvia, Catalogue #SH3028.02) and 0.1 mM Protease Inhibitor (Sigma-Aldrich, Catalogue #P8340). If there were multiple cell pellets from the same condition, they were combined into 1 tube at this step. Cells were centrifuged at 500×g at 4°C and the supernatant was aspirated off. To lyse the plasma membrane, the cells were resuspended in Buffer A and incubated for 10 min on ice (10 mM HEPES, pH 7.9, 1.5 mM MgCl, 10 mM KCl, 0.1 mM Protease Inhibitor, Phosphatase Inhibitor (Santa-Cruz Biotechnology, Catalogue #sc-45044), 0.5 mM DTT (Sigma-Aldrich, Catalogue #4315)). After the 10 min incubation, Igepal detergent (final concentration of 0.1%) was added to the cell and Buffer A mixture and vortexed for 10 s. To separate the cytosolic fraction from the nuclei, the sample was centrifuged at 500×g for 5 min at 4°C to pellet the nuclei. The cytosolic fraction was collected into a separate microcentrifuge tube. The pelleted nuclei were washed gently with Buffer A, which was added and removed without disrupting the pellet. The pelleted nuclei were then resuspended in Buffer C (20 mM HEPES, pH 7.9, 25% glycerol, 1.5 mM MgCl, 0.2 mM EDTA, 0.1 mM Protease Inhibitor, Phosphatase Inhibitor, 0.5 mM DTT, and 420 mM NaCl) and vortexed for 30 s. To extract the nuclear proteins (i.e., the nuclear extract), the nuclei were incubated in Buffer C for 1 h while mixing at 4°C on Life Technologies HulaMixer (Orbital: OFF, Reciprocal: 90°/30 s, Vibro: 5°/5 s). To separate the nuclear extract from the nuclear debris, the mixture was centrifuged at 21,000×g for 20 min at 4°C. The nuclear extract was collected in a separate microcentrifuge tube and flash frozen using liquid nitrogen. Nuclear extracts were stored at −80°C.

Protein-binding microarray (PBM)

Microarray DNA double stranding and PBM protocols are as previously described [108,143–145]. Any changes to the previously published protocols are detailed. Double-stranded microarrays were pre-wetted in HBS (20 mM HEPES, 150 mM NaCl) containing 0.01% Triton X-100 for 5 min and then de-wetted in an HBS bath. Microarrays were incubated with 2% non-fat dried milk solution (Lab Scientific, Cat # M0841) in HBS for 1 h in the dark. The array was then rinsed in an HBS bath containing 0.1% Tween-20 and subsequently de-wetted in an HBS bath. Next, the array was incubated with nuclear extract (420 μg of nuclear extract for 4 × 180 K array design) for 1 h in the dark in a binding reaction buffer (20 mM HEPES, pH 7.9, 100 mM NaCl, 1 mM DTT, 0.2 mg/ml BSA, 0.02% Triton X-100, 0.4 mg/ml salmon testes DNA (Sigma-Aldrich, Catalogue #D7656)). The array was then again rinsed in an HBS bath containing 0.1% Tween-20 and de-wetted in an HBS bath. After the protein incubation, the array was incubated for 20 min in the dark with 20 μg/ml primary antibody for the COF of interest (S1 Table). The primary antibody was diluted in 2% non-fat dried milk in HBS. After the primary antibody incubation, the array was rinsed in an HBS bath containing 0.1% Tween-20 and de-wetted in an HBS bath. Microarrays were then incubated with 20 μg/ml of either alexa488 or alexa647 conjugated secondary antibody (S1 Table) for 20 min in the dark. The secondary antibody was diluted in 2% milk in HBS. The array was rinsed in bath containing HBS and 0.1% Tween-20. The array was then place in a Coplin Jars twice for 3 min in 0.05% Tween-20 in HBS and once for 2 min in HBS in Coplin jars to remove excess antibody. After the washes, the array was de-wetted in an HBS bath. Microarrays were scanned with a GenePix 4400A scanner and fluorescence was quantified using GenePix Pro 7.2. Exported fluorescence data were normalized with MicroArray LINEar Regression.

Replicate PBM experiments for each COF in both acute and exhausted T-cell conditions were as follows: KAT2B (4 replicates each condition); HDAC mix (HDAC 1, II, III) (4 replicates each condition); P300 (2 replicates each condition). Reported TF-COF interactions in each conditions had to be seen in 2 replicate experiments.

PBM microarray design and analysis

The microarray design contained a set of 346 nonredundant TF-binding models from the JASPAR 2018 core vertebrate motif set and was obtained using the JASPAR 2018 R bioconductor package. TF-binding models refer to the consensus sequence for each motif and all single-nucleotide variants (SVs) of the consensus, which allows us to directly construct motifs from our binding data. Each transcription factor binding site (TFBS) (i.e., consensus or SV sequences) was embedded within a 34-base pair (bp) long target region, within the context of a 60-bp long DNA probe of the form: GC cap (2 bp) + target region (34 bp) + primer sequence (24 bp). Within the target region, we embedded each consensus sequence with a variable number of padding bases to account for different length motifs, using the format: 5′ pad (2 bp) + consensus sequence (L bp) + 3′ variable pad (32–L bp). The pad sequences (both 5′ and 3′) were generated randomly with the constraint that sequential positions contain non-repeating nucleotides. Background target DNA probes were also included in the design to estimate background fluorescence intensities in the experiments. These regions were selected at 34-base genomic segments from the human genome (hg38).

Cofactor recruitment motifs modeled using single variant (SV) probes was performed as previously described [108,143]. DNA probe log fluorescence values were first z-transformed using the fluorescence values from the background probes. DNA motifs were then determined for all TF-binding models (i.e., sets of consensus and associated SV probe sequences). For each TF-binding model, the base preference at each position was evaluated using Δz-scores calculated as the difference relative to the median z-score obtained across all possible nucleotides at that position (i.e., the consensus and the 3 SV probes). This z-score-based binding motif captures the differential impact of base variants on the log fluorescence binding, which is analogous to a change in binding energy. The z-score matrices can be transformed into position weight matrices (PWMs) to facilitate comparison against motif databases using the approach previously described [108,143]. Motifs were only considered significant if they passed the following criteria. (1) z-score of the consensus probe for each TF-binding models >0.5; (2) the average information content (IC) of at least one 5-bp window within the motif must > 1; (3) the adjusted p-value of the best motif match against JASPAR motif database >0.05. Motif matching was performed using TomTom (Meme suite) [146]. Motif strength is the mean of Δz-score for the top 15% of probes for each TF-binding model. Individual TF-binding motifs were clustered into TF family groups based on highly similar and generally indistinguishable motifs. Matched motifs were called into their TF family cluster. The motif strength for each TF family group is the mean of all the significant motifs within that TF family cluster.

Antibodies used for these assays can be found in S1 Table.

Fluorescence-activated cell sorting (FACS)

Briefly, ~2 × 10⁸ frozen PBMCs from 4 people living with HIV (PWH) and an uninfected person were thawed, washed, and stained with Zombie NIR viability dye (1:400, Biolegend) in phosphate-buffered saline (DPBS) for 15 min at room temperature (RT). PBMCs from the 4 PWH were combined into 1 tube for staining and sorting while PBMCs from the uninfected individual were used for necessary fluorescence experimental controls. Viability staining was quenched with FACS buffer (DPBS + 0.5% bovine serum albumin + 2 mM ethylenediaminetetraacetic acid), PBMCs were washed, and then pre-stained with a Fc receptor blocking solution (human TruStain FcX, 1:50, Biolegend) for 10 min at RT. Next, PBMCs were stained using an 8-parameter antibody cocktail (CD8 (BUV805, SK1, BD Biosciences), PD-1 (BV421, EH12.2H7, Biolegend), CD3 (BV510, OKT3, Biolegend), TIGIT (PE-Cy7, A15153G, Biolegend), CD4 (Alexa Fluor 700, RPA-T4, Biolegend), CD14 (APC/Fire750, 63D3, Biolegend), and CD19 (APC/Fire750, SJ25C1, Biolegend)) combined with monocyte blocker (1:20, Biolegend) and brilliant stain buffer plus (1:10, BD Biosciences) for 30 min on ice. After cell surface staining, PBMCs were washed twice with FACS buffer, resuspended in FACS buffer, and kept on ice until acquisition. All staining was performed in the dark and cells were sorted using a FACSAriaTM cell sorter (BD Biosciences). Full stained PBMCs from PWH were sorted into 4 populations of cells based on IR expression on T cells: CD4+ PD-1+ +/- TIGIT+ (denoted IR+) cells, CD4+ PD-1- TIGIT- (denoted IR-), CD8+/- CD4- (non-CD4 CD3+) PD-1+ +/- TIGIT+ (denoted IR+), and CD8+/- CD4- (non-CD4 CD3+) PD-1- TIGIT- (denoted IR-); ~1–2.5 × 106 cells per population were collected for downstream proteomic analysis. Additional information for all antibodies used for FACS can be found in S2 Table.

Additional demographic information regarding subjects used for sorting

The 4 people living with HIV infection ranged from 26 to 56 years old, included 2 males and 2 females, 2 identified as black/African American, and 2 identified as Hispanic; all 4 individuals were part of the Boston University HIV/Aging cohort with an enrollment requirement of effective viral suppression with undetectable HIV-1 RNA for a minimum of 6 month (

Data reporting

All relevant data and supplemental data are within the paper and its Supporting Information files excluding raw flow cytometry and mass spectrometry files, which can be found in the following locations: Flow cytometry files are freely accessible at http://flowrepository.org using the identifier “FR-FCM-Z8G6,” and the mass spectrometry data files have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the data set identifier PXD057703. For more information, please see reference below.

Perez-Riverol Y, Bai J, Bandla C, Hewapathirana S, García-Seisdedos D, Kamatchinathan S, Kundu D, Prakash A, Frericks-Zipper A, Eisenacher M, Walzer M, Wang S, Brazma A, Vizcaíno JA (2022). The PRIDE database resources in 2022: A Hub for mass spectrometry-based proteomics evidences. Nucleic Acids Res 50(D1):D543-D552 (PubMed ID: 34723319).