Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) provides an invaluable resource of high-throughput gene expression data that can be integrated and analyzed. Datasets from GEO that were considered eligible in our analysis met the following criteria: 1) Datasets with T-cell transcriptomic data from patients who developed aGvHD later (pre-aGvHD hereafter) and those who did not develop aGvHD (non-aGvHD hereafter); 2) Datasets with group information for each sample (pre-aGvHD vs. non-aGvHD); 3) Datasets with information about the technology and platform used for studies; and 4) Datasets with more than 10,000 probes. Based on the above criteria, two aGvHD datasets were downloaded from the GEO database. Details of each dataset are shown in Table S1 .

Bias of high-throughput experiments were commonly derived from heterogeneity and variables. In this study, we recruited datasets from different platforms and samples and handled by individual researchers. To avoid the possible unreliable results, samples from the two datasets were integrated to increase the number of samples (37 pre-aGvHD samples vs. 48 non-aGvHD samples), followed by batch normalization using R package “sva” (11). Next, gene differential analysis (|LogFC| >1, P < 0.05) was performed by comparing T cells isolated from pre-aGvHD and non-aGvHD blood samples using R package “limma”. Afterward, volcano map and heatmap were depicted by “limma” package and “pheatmap” package, respectively, in the R computing environment (12).

Gene Set Variation Analysis (GSVA) was applied to calculate the activities of metabolic and signaling pathways for each sample (13). Pathway signatures were obtained from MSigDB of the Broad Institute (https://www.gsea-msigdb.org/gsea/index.jsp) and Reactome annotation (https://reactome.org) (14, 15). For comparison of GSVA scores, the T-cell expression data were multiplied by 1 for aGvHD-dependent lines or by −1 for aGvHD-independent lines to reflect the direction of aGvHD dependency (positive for aGvHD dependency and negative for aGvHD independency). The data were then standardized to z-scores across samples for comparison and the creation of correlation matrix heatmaps.

The correlation of gene expression level and pathway score derived from GVSA was evaluated by Pearson correlation. |Correlation coefficient| >0.4 and P < 0.05 were considered statistically significant.

We used the CIBERSORT tools (16) to identify immune cell profile from the RNA expression data.

TFs potentially driving the expression of glycolytic genes and metabolic reprogramming were predicted by the module “UCSC_TFBS” under the “Protein_Interactions” function of DAVID (https://david.ncifcrf.gov/home.jsp) and Pearson correlation analysis. TFs regulating the expression of individual gene (such as PFKP, ENO3, or GAPDH) were predicted through R package “TFBSTools” (17) and JASPAR database (http://jaspar.genereg.net/). The target genes of E2F6 were predicted using hematopoietic cell-specific chromatin immunoprecipitation followed by sequencing (ChIP-seq) data collected in Cistrome data browser (18).

Glycolysis-associated genes, PFKP, ENO3, and GAPDH, were applied to develop an aGvHD risk model. Nomogram was constructed to predict aGvHD risk with multivariate logistic regression analysis. To assess the calibration of the risk model, calibration curves were plotted. To quantify the discrimination performance of the risk model, Harrell’s C-index was measured, and receiver operating characteristic (ROC) curve was drawn. Furthermore, the risk model was subjected to bootstrapping validation (1,000 bootstrap resamples) to calculate a relatively corrected C-index. To determine the clinical usefulness of the risk model, decision curve was plotted through calculating net benefits across different threshold probabilities.

CD3⁺ T cells were collected from 5 non-aGvHD patients and 9 pre-aGvHD patients by EasySep™ Human T Cell Enrichment Kit (Stemcell Technology, Cat # 19051). For 5′ single-cell RNA sequencing (RNA-seq) data, raw reads obtained from the 10× Genomics single-cell RNA-seq platform were demultiplexed and mapped to the human reference genome GRCh38 using the CellRanger software (version 3.0.2) (https://support.10xgenomics.com/single-cell-gene-expression/software) with default parameters. In this study, CD3⁺ cells were retained for downstream analysis. Finally, our study included 48,718 CD3⁺ T cells. Normalized gene expression values were plotted for each cell as violin plot in R.

For the analysis of cell surface molecules, single-cell suspensions were prepared. CD4⁺ memory T cells (CD4⁺CD45RO⁺) were detected by flow cytometry.

The study was approved by the Ruijin Hospital Ethics Committee, and all processes were consistent with Helsinki Declaration standards.

The differentially expressed genes (DEGs) and GSVA results were displayed with P values, fold changes, and ranks. An unpaired two-tailed Student’s t-test (for two group comparisons) or a one-way ANOVA was performed, and the Wilcoxon rank-sum test was performed using R package ggplot2. The results of multivariate logistic regression were displayed with odds ratio (OR) and P values. The strength of the Pearson correlation was displayed with correlation coefficient as the following guide: 0.00–0.19, “very weak”; 0.20–0.39, “weak”; 0.40–0.59, “moderate”; 0.60–0.79, “strong”; 0.80–1.0, “very strong.” P value <0.05 was statistically significant. All statistical tests and graphing were performed by R. In the figures, statistical significance was shown as follows: *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.

To explore transcriptional metabolic reprogramming of T cells before full onset of aGvHD, we integrated the T-cell transcriptomic data across different studies (19, 20). The raw data and platform information of GSE4624 and GSE73809 were downloaded from GEO database. Both studies collected T cells prior to aGvHD onset (pre-aGvHD group) and on matched time in control patients who were without aGvHD (non-aGvHD group) and concluded that gene expression signature of alloreactive T cells had a dominant effect on the development of aGvHD. After annotation, 8,455 genes in GSE4624 (GPL3639) and 30,905 genes in GSE73809 (GPL17586) were obtained ( Table S1 ). To reduce the possibility of false positive results, the two datasets were batch normalized and integrated ( Figure S1 ). The expression patterns of metabolic pathway genes have been determined to reflect actual metabolic activities in patients well (14). To analyze the metabolic activities of each sample, we examined the gene sets of seven super-pathways based on Reactome annotation (14, 15) including amino acid metabolism (348 genes), carbohydrate metabolism (286 genes), energy metabolism (110 genes), lipid metabolism (766 genes), nucleotide metabolism (90 genes), TCA cycle (148 genes), and vitamin and cofactor metabolism (168 genes) with GSVA approach ( Table S2 ) (13). Next, we compared the enrichment scores of seven metabolic signatures between T cells from pre-aGvHD and non-aGvHD samples. The results showed that carbohydrate, nucleotide, TCA cycle, and amino acid metabolism were significantly enriched in pre-aGvHD samples ( Figure 1A , Table S3 ). Of note, carbohydrate metabolism was the most significantly enriched (P = 0.0003). To identify the metabolic pathways dysregulated in pre-aGvHD samples, we examined sub-metabolic gene sets of the four enriched super-pathways based on Reactome and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation. The sub-pathways with high activity in pre-aGvHD samples were mainly glycolysis-related, such as glycolysis, pentose phosphate pathway, fructose and mannose metabolism, galactose metabolism, and respiratory electron transport ( Figure 1B ), of which glycolysis had the most extensive correlation with other enriched sub-pathways ( Figure 1C ). These results suggested that glycolysis might be the crucial metabolic pathway. The extensively activated pathways might be the results of glycolysis activation before aGvHD full onset.

To further identify the genes involved in the activation of glycolysis employed by T cells in the pre-aGvHD state, we performed differential expression analyses of genes derived from REACTOME_Glycolysis and KEGG_Glycolysis_gluconeogenesis gene sets. Compared with T cells isolated from non-aGvHD samples, 18 glycolytic genes were upregulated and six genes were downregulated in pre-aGvHD group (|Fold Change| >1 and P < 0.05; Figure 2A , Table S4 ). The top 5 upregulated genes were PFKP, SOD1, ENO3, GAPDH, and STMN1. Among them, PFKP, ENO3, and GAPDH have well-known functions in glycolysis. Phosphofructokinase (encoded by PFKP) catalyzes the phosphorylation of fructose 6-phosphate (F6P) to fructose 1,6-bisphosphate (F16BP) by ATP, which is the most important control step of glycolysis (21). Enolase (encoded by ENO3), also known as phosphopyruvate dehydratase, catalyzes the transformation of 2-phosphoglycerate (2-PG) to phosphoenolpyruvate (PEP) (22). Glyceraldehyde-3-phosphate dehydrogenase (encoded by GAPDH) is a glycolytic enzyme that catalyzes the critical step by converting glyceraldehyde 3-phosphate (G3P) into 1,3-bisphosphoglycerate (1,3-BPG). Inhibition of glyceraldehyde-3-phosphate dehydrogenase (GAPDH) downregulates glycolysis in both myeloid and lymphoid cells, preventing immune activation (23–25). Participation of upregulated genes in glycolysis and their interaction with other metabolic pathways are visualized in Figure 2B .

In the pathogenesis of aGvHD, donor-derived alloreactive T cells are activated by recognition of host antigens, then adopt a pathogenic effector phenotype and migrate to target organs. Activated T cells require more energy to proliferate and mediate effector functions. However, the underlying relationship between metabolic reprogramming and alloreactive T-cell activation is still poorly understood (10). To assess the biological relevance of metabolic reprogramming, we evaluated T-cell subsets and T-cell functional signatures by GSVA. We first scored seven T-cell subsets for their relative abundance by using CIBERSORT (16) and evaluated the correlation of the T-cell subsets with metabolic reprogramming and differentially expressed glycolytic genes. Activated CD4 memory T cell was identified to be the top subset that positively correlated with glycolysis-related pathways and glycolytic genes ( Figures 3A, B ). The expressions of TXN, ALDOB, GAPDH, STMN1, PSMC4, ENO3, PRKACB, and PFKP had strong correlations with activated CD4 memory T subset ( Figure 3B ). Furthermore, we validated the positive correlation between CD4 memory T subset with aGvHD development and glycolysis activation in our patient cohort. We found that the proportion of CD4 memory T cells was significantly higher in pre-aGvHD samples than their counterparts in non-aGvHD group (13.16% ± 6.585% vs. 1.378% ± 0.8093%, P < 0.05) ( Figure S2A ). Single-cell RNA sequencing data revealed significantly higher expression of glycolytic genes PFKP and GAPDH in CD4 memory T cells of pre-aGvHD patients than non-aGvHD patients ( Figure S2B ), supporting the important function of glycolytic activation in CD4 memory T cells during the aGvHD development. Next, we examined enrichment of T-cell functional signatures derived from gene ontology (GO) for validation. T-cell proliferation and migration signatures were positively associated with glycolysis-related pathways and glycolytic genes ( Figures S3A, B ). The expressions of ALDOB, ENO3, PRKACB, STMN1, TXN, PFKP, GAPDH, ALDOA, PFKP, AGRN, PSMC4, and SOD1 had significant correlations with T-cell proliferation and migration ( Figure S3B ). Overall, these results suggested that metabolic activity was intrinsically coupled with T-cell activation pathways.

Metabolic reprogramming is largely determined by signaling pathways and TFs. TCR, the phosphatidylinositol 3-kinase (PI3K)/protein kinase B (AKT), mTOR, MYC signaling, and TFs such as IRF4, SREBP, PGC1α, HIF1α, ATF4, and E2F have been widely involved in anabolic or catabolic metabolism, including glycolysis and REDOX balance (14, 26–28). Thus, we compared the enrichment scores of hallmark signaling and TCR signaling pathways between T cells from pre-aGvHD and non-aGvHD samples. The results showed that MYC, mTORC1, TCR, Hedgehog, and Wnt/β-catenin signaling pathways were significantly enriched in pre-aGvHD groups ( Figure S4 , Table S5 ). To identify specific driver signaling pathways, we performed correlation analyses of metabolic reprogramming and differentially expressed glycolytic genes with hallmark signaling and TCR signaling pathways. TCR, mTORC1, and MYC signaling pathways were identified as drivers with top positive correlations ( Figures 4A, B ). The expressions of PFKP, STMN1, PSMC4, GAPDH, TXN, ENO3, and PRKACB had stronger correlations (Correlation coefficients >0.5 and P < 0.05) with TCR, mTORC1, and MYC signaling activity ( Figure 4B ).

Gene expression is generally regulated by TFs directly, which consequently influences metabolic activity. To elucidate TFs that potentially drive the expression of glycolytic genes and metabolic reprogramming, DEGs between pre-aGvHD and non-aGvHD groups ( Figures S5, S6 ) were utilized for the prediction of TFs by the module “UCSC_TFBS” under the “Protein_Interactions” function of DAVID (https://david.ncifcrf.gov/home.jsp). TFs predicted with adjusted P < 0.05 were considered to be significantly enriched ( Table S6 ). The enriched TFs and other well-known metabolism-regulating TFs (26, 27, 29) whose expression data are available to access in this study (E2F1, E2F4, E2F6, MYC, ATF4, ATF6, NFKB1, STAT3, LPIN1, HIF1A, RPS6KB2, YY1, SREBF2) were submitted for the following analyses. We performed correlation analyses of metabolic reprogramming and differentially expressed glycolytic genes with the abovementioned TFs. E2F6 was identified as the top positively correlated TF ( Figures S7A, B ), which was consistent with the positive correlation of E2F target signature with metabolism reprogramming and glycolytic gene expression ( Figures 4A, B ). Furthermore, expressions of MYC, MTOR, and E2F6 were proven to be higher in pre-aGvHD T cells than those in non-aGvHD T cells by our own patient cohort ( Figure S8A ). These results suggest that TCR, mTORC1, and MYC signaling pathways and E2F6 are potential “master regulators” of metabolic reprogramming and glycolytic activity.

Regarding the central function of glycolysis in alloreactive T-cell activation in pre-aGvHD samples, we used PFKP, ENO3, and GAPDH, which were among the most significantly upregulated glycolytic genes and their expressions were proven to be higher in pre-aGvHD T cells than those in non-aGvHD T cells ( Figure S8B ), to construct an aGvHD predictive model. The predictive model that incorporated the three predictors was developed and presented in the nomogram ( Figure 5A ). The area under the ROC curve was 0.806, showing a high predictive value ( Figure 5B ). The C-index of the nomogram was 0.806 (95% CI: 0.713–0.899) for the cohort and was confirmed to be 0.786 through bootstrapping validation, suggesting that the risk model had good discriminative ability. The calibration curve demonstrated that the nomogram had good concordance to predict aGvHD risk in this cohort ( Figure 5C ). The decision curve showed that the aGvHD predictive model provided superior net benefit when clinical decision thresholds were between 14% and 91% ( Figure 5D ). Furthermore, the predictive value of the aGvHD risk model was validated by our own patient cohort ( Figures S8C, D ).

Using TFBSTools and JASPAR database, we predicted E2F1, E2F4, and E2F6 to be possible TFs promoting the expression of PFKP, ENO3, and GAPDH. Moreover, hematopoietic cell-specific ChIP-seq data further supported E2F6 as a TF regulating the expression of PFKP, ENO3, and GAPDH ( Figure S9 ), which was in accordance with the positive correlation between E2F6 and their expressions ( Figure S7B ) ( 18). Furthermore, we confirmed that T cells with higher E2F6 expression harbored increased expression of PFKP, ENO3, and GAPDH in our patient cohort ( Figure S10 ). Taken together, these results suggest that TCR, mTORC1, and MYC signaling pathways might promote the expression of PFKP, ENO3, and GAPDH through TF E2F6 to activate glycolysis and its related pathways in alloreactive T cells during aGvHD development ( Figure 6 ).