Early-to-mid stage PD shows a more-cytotoxic & late-differentiated immune profile
In the discovery analysis, we systematically scrutinized various immune subsets and their functional states in 28 PD patients (25 iPD aged 60-70 years and three genetic PD patients with mutations in GBA or PINK1) and 24 matched healthy controls (HC) (refer to “cohort design” in “Methods”, Supplementary Tables 1 and 2; for simplicity, ‘PD patients’ are hereafter shortened as PD). This was realized by investigating 37 different innate and adaptive immune subsets and more than 700 combinatorial T-cell features, using a 35-marker mass cytometry (also known as Cytometry by Time of Flight or CyTOF) panel and five panels of multiple-color flow-cytometry (FCM) composed of 33 lineage and functional T-cell markers (Fig. 1a, Supplementary Tables 3 and 4). We recruited the participants from the ongoing nation-wide Luxembourg Parkinson’s study with more than 800 PD and 800 HC22 (https://www.parkinson.lu/research-participation) and controlled for several major confounding factors, cytomegalovirus (CMV) serostatus, medications and comorbidities, known to affect the immune system, to ensure that our observations are PD-specific (Fig. 1a and Supplementary Table 1). Furthermore, we narrowed the patients to those with early-to-mid stage disease [Hoehn and Yahr (H&Y) staging scale: mean = 2.3, ranging from 1.5 to 3.0; most of them ≤2.5, except for five participants with a scale of 3]. Most of the selected patients had a disease duration within 10 years from clinical diagnosis while three of them had a disease duration of 12, 13 and 19 years, respectively.
a Graphical representation of the cohort and experimental setup. b PCA plot showing no distinct immunological fingerprint of PD based on the entire peripheral immune system. c Volcano plot highlighting the most significantly (p < 0.05, fold change >1.3) decreased and increased subsets in PD vs HC. The horizontal dashed line corresponds to –log10(0.05), while the two vertical dashed lines correspond to log2 value of −0.3785 or 0.3785. The marker definition for some of the highlighted subsets is provided here or other panels: CD4 naive (CD45RA+CCR7+), CD4 Tfh (CD45RA–CXCR5+) and naive B (IgD+CD27–) (Supplementary Fig. 1). Scatter dot plots showing the frequency of total CD8 T cells (d), CD8 CD45RA+CCR7– (the simplified gating for TEMRA) (e), CD8 TCM (CD45RA−CCR7+, central memory) (f) and CD57 MSI in CD8 TEMRA (g). Scatter dot plots showing the frequency of CD8+ NKT (CD8+CD4–CD3+CD19–CD56+TCRgd−, h), CD4+ NKT (CD8–CD4+CD3+CD19–CD56+TCRgd−, i) and CD4–CD8– (DN) NKT (CD8−CD4−CD3+CD19−CD56+TCRgd−, j) and CD57 MSI in total NKT (k). l, Frequency of CD56hiCD57– immature NK cells. m Representative viSNE plot from either HC or PD highlighting the expression levels of CD45RA, CCR7, CD27 and CD57 in total CD8 T cells. The arrow indicates the area of CD8 CD45RA+CCR7–CD27− (TEMRA). Scatter dot plots showing the frequency of neutrophils (CD66b+CD45midCD16+CD294−, n), eosinophils (CD66b+CD45midCD16–CD294+, o) and innate lymphoid cell type 2 (ILC2, p). q Representative cytometry plots showing the expression of CD294 and CD117. The number showing the frequency among living CD45+ singlets. The results in (c–l, n–p) were analyzed using unpaired two-tailed Student’s t test without adjustments for multiple comparisons. Data are presented as mean ± standard deviation (s.d.). Each symbol represents the measurement from one individual (d–l, n–p). ns or unlabeled, not significant; all significant P-values are indicated. CyTOF mass cytometry, CMV cytomegalovirus, NKT Natural killer T cells, PBMC peripheral blood mononuclear cells, HC healthy controls, n = 24, PD patients with Parkinson’s disease, n = 28; MSI Median signal intensity. Panel a was created with BioRender.com. Source data are provided as a Source Data file.
As a strong cryopreservation effect has been observed on the FCM readouts of T-cell markers23, we performed the cytometric analysis on fresh blood (Supplementary Fig. 1). Our CyTOF analysis did not show a different general immunological fingerprint of PD vs HC, based on the entire peripheral immune system (Fig. 1b). Nevertheless, several immune cell types, especially the T-cell compartments, were altered (Fig. 1c). Total classical αβ T cells were modestly reduced in PD (Supplementary Fig. 2A), reflected by a decrease of total CD4 T cells (Supplementary Fig. 2B), whereas the γδ T cells were unchanged among total living CD45+ cells (Supplementary Fig. 2C). The decreased frequency of total CD4 T cells was mainly due to a reduction in CD4+CXCR5+ T follicular helper cells (Tfh) (Supplementary Fig. 2D), CD45RA+CCR7+ naive (Supplementary Fig. 2E) and CD45RA–CCR7+ central memory (TCM) CD25− conventional T cells (Tconv) (Supplementary Fig. 2F), but not CD45RA−CCR7− effector memory (TEM) Tconv (Supplementary Fig. 2G). Although total CD8 T cells showed no difference in PD vs HC (Fig. 1c, d), the CD8 naive/memory subset composition displayed alterations (Fig. 1c). Our unbiased volcano plot analysis (Fig. 1c) showed that the frequency of cytotoxic terminally-differentiated effector memory T cells (CD45RA+CCR7−, TEMRA)24,25 was substantially increased among total CD8 T cells (Fig. 1c, e), whereas the frequency of CD8 TCM was reduced in PD (Fig. 1f). The proportion of CD8 naive and TEM showed similarity between PD and HC (Supplementary Fig. 2H, I). Furthermore, the expression levels (median signal intensity, MSI) of CD57, a marker for terminal differentiation, among CD8 TEMRA showed a trend to be increased (p = 0.0663) in PD (Fig. 1g). A trend to be increased was also observed for the frequency of CD57+ cells among CD8 TEMRA (p = 0.0568) (Supplementary Fig. 2J), further indicating a late-differentiated CD8 T-cell profile. Moreover, another cytotoxic cell type, natural killer T cells (NKT) also exhibited a late-differentiated state, as reflected by an increased frequency of CD8+ NKT (Fig. 1h), among total NKT26, especially in women. Meanwhile, the frequency of less-differentiated CD4+ (Fig. 1i) and CD4−CD8− (also known as double negative, DN, Fig. 1j) NKT was either decreased or intact, respectively. The female-biased relative increase in PD might be explained by a lower fraction of CD8+ NKT in healthy females than males in their sixties (Supplementary Fig. 2K). CD8+ NKT also expressed higher levels of CD57 (Fig. 1k). The frequency of CD56highCD57– immature NK was also decreased (Fig. 1l). As CD8 T-cell composition was considerably changed, while total CD8 T cells were intact, we performed an unsupervised analysis on gated CD8 T cells. Indeed, our unsupervised viSNE analysis confirmed an enhanced frequency of CD8 TEMRA (CD45RA+CCR7−CD27−) (Fig. 1m). In brief, we observed a more-cytotoxic and late-differentiated immune profile in early-to-mid stage PD, as reflected by several relevant subsets, such as CD8 TEMRA, NKT and NK.
With the 35-marker CyTOF analysis in whole blood, we could assess many other immune subsets, such as granulocytes (neutrophils, eosinophils and basophils), monocytes (classical, intermediate and non-classical), dendritic cells (myeloid DC and plasmacytoid DC, known as mDC and pDC respectively), NK (immature and late), B cells (naive, memory, plasma cells) and innate lymphoid cells (ILCs: ILC1, ILC2 and ILC3) (Supplementary Fig. 1A–M and Supplementary Table 5). Most of them did not show any significant change in PD vs HC regarding the frequency among total CD45+ cells or among the relevant parent gates (Supplementary Table 5). Consistent with a more-cytotoxic profile, an increased frequency of neutrophils was observed in PD, especially for males (Fig. 1n, Supplementary Fig. 2L). The heightened frequency of neutrophils was accompanied with a reduction in eosinophils, especially for males (Fig. 1o, Supplementary Fig. 2M). Unlike the scenario of CD8+ NKT, the gender-biased changes in these two types of granulocytes in PD cannot be simply explained by the pre-existing difference between male and female HC (Supplementary Fig. 2L, M). In the meantime, basophils were unchanged in PD vs HC (Supplementary Fig. 2N). The reduced eosinophils and increased neutrophils are conforming to two recent works27,28, respectively, where although mainly routine whole blood counts have been analyzed.
Unexpectedly, the frequency of ILC229 in PD was decreased to almost half of that in HC (Fig. 1p, q), while ILC1 and ILC3 showed no difference (Supplementary Fig. 2O, P). The reduction was even more pronounced in males, which could be explained by a higher number of ILC2 already in male vs female HC in their sixties (Supplementary Fig. 2Q). The reduction in ILC2 is consistent with the observed decrease in eosinophils, as ILC2 control eosinophil homeostasis at least in typical Type II immunopathology (i.e., allergy)30. Finally, we found a slight but significant decrease in the frequency of IgD+CD27− naive B cells (Supplementary Figs. 1 and 2R). Collectively, our unbiased comprehensive immunophenotyping analysis reveals a more-cytotoxic and late-differentiated immune profile, while the frequency of immature NK, eosinophils and ILC2 was decreased in early-to-mid stage PD.
Early-to-mid stage PD exhibits a heightened effector profile in CD8 T cells
Based on our T-cell-dominant observations by CyTOF, we next analyzed T cells in depth using five FCM panels with a total of 33 T-cell-relevant markers (Supplementary Table 4), the combinations of which gave rise to ~700 features. The PCA determined a distinct immunological fingerprint in iPD vs HC (with the exception of one PD and one HC labeled as “PD9” and “HC17” in Fig. 2a). The three genetic PD were not identified as outliers compared to iPD in the PCA plot based on the comprehensive T-cell analysis. Alike to the CyTOF data, no difference was observed in the frequency of total CD8 T cells among fresh peripheral blood mononuclear cells (PBMC) in iPD vs HC (Supplementary Fig. 3A). Different from the CyTOF data, we did not observe any significant difference in the frequency of total T cells and CD4 T cells between the two groups (Supplementary Fig. 3A). This might be due to the exclusion of granulocytes (accounting for the majority of immune cells in whole blood) in PBMC. Although the overall frequency of major T-cell populations remained unchanged in the FCM analysis (Supplementary Fig. 3A), the PCA results indicated the existence of changes in specific T-cell subsets (Fig. 2a).
a PCA plot showing a distinct immunological fingerprint of PD vs HC based on T-cell combinatorial features using the FCM analysis. b Volcano plot showing the most significantly (p < 0.05, fold change >1.4) decreased and increased subpopulations in PD vs HC. The horizontal dashed line corresponds to the value of 1.3 (p = 0.05), while the two vertical dashed lines correspond to the log2 value of −0.485 or 0.485 (fold change of 1.4). c Scatter dot plots (left) and representative FCM plots (right) showing the increase in CD8 TEMRA (CD45RO−CD45RA+CCR7−, the simplified gating strategy for TEMRA) for all (left) or females (middle) or males (right). The TEMRA gate was highlighted in blue dashed rectangle in FCM plots. The combination of markers used to define TEMRA was described in the y-axis title. d Scatter dot plots showing the ratios between CD8 TEMRA and CD8 TCM (CD45RO+CCR7+CD27+). e-i Scatter dot plots (upper) and representative FCM plots (lower) showing the frequency of T-bethigh (e), CD45RO−T-bet+ (f), CD45RO−Ki67+ (g), CD45RO–Helios+ (h) and CD45RO–CD57+ (i) among CD8 T cells. In e, blue and magenta dashed rectangle highlights T-bethigh and T-betmid cells, respectively, while gold dashed rectangle shows Eomes+ cells. The percentage of the interest was enlarged. j–l Serological levels of GZMA (j) for all (left) or females (middle) or males (right), GZMB (k) and perforin (l). The results in b–i were analyzed using unpaired two-tailed Student’s t test while the results in j–l were analyzed using unpaired two-tailed Mann–Whitney nonparametric test without adjustments for multiple comparisons. Data are presented as mean ± standard deviation (s.d.). Each symbol represents the measurement from one individual (c–l). ns or unlabeled, not significant; all significant P-values are indicated. HC healthy controls, n = 24; PD, patients with Parkinson’s disease, n = 28; FCM flow cytometry; GZMA/GZMB Granzyme A/B, TCM central memory. For c, d, j, female HC, n = 10; female PD, n = 8; male HC, n = 14; male PD, n = 19. Source data are provided as a Source Data file.
Among the most significantly changed immune subpopulations, the FCM analysis again found a strong increase in the frequency of TEMRA among total CD8 T cells (Fig. 2b, c), confirming our CyTOF results. CD8 TEMRA cells re-express CD45RA, but lose the expression of CD45RO, CCR7 and CD27. Following a simplified gating strategy, we were able to pinpoint that CD45RA+CD45RO−CCR7− cells (CD8 TEMRA) were increased among total CD8 T cells in iPD (Fig. 2c). To more strictly identify CD8 TEMRA, we next included CD27 in another staining panel. Similar to that from the simplified gating strategy (Fig. 2c), the frequency of CD45RO−CCR7−CD27− effector CD8 T cells was also increased (Supplementary Fig. 3B). As the fraction of CD45RA/CD45RO double-negative (mean: ~5%) or double-positive (mean: ~5%) cells was tiny among CD8 T cells (Supplementary Fig. 3C), most of the CD45RO–CCR7–CD27– effector CD8 T cells (Supplementary Fig. 3B) should be CD45RA+CCR7−CD27− CD8 TEMRA, and thus also referred as CD8 TEMRA in certain analyses. Meanwhile, TCM (CD45RO+CCR7+CD27+) and transitional memory (TM) (CD45RO+CCR7–CD27+) CD8 T cells were reduced (Supplementary Fig. 3D, E), while naive (CD45RO−CCR7+CD27+) CD8 T cells showed no difference between PD and HC (Supplementary Fig. 3F). The reduction in CD8 TCM and TM cells was in line with a lower frequency of long-lived memory (KLRG1−CD127+)31 CD8 T cells (Supplementary Fig. 3G) and consistent with the CyTOF results (Fig. 1f). Although total CD8 cells were not changed, the mean ratios of CD8 TEMRA over TCM were substantially enhanced from 12.05 in HC to 42.90 in iPD (Fig. 2d). Jointly, we firmly demonstrate that the portion of CD8 TEMRA is increased while the fraction of CD8 TCM and TM is contracted in the peripheral blood of early-to-mid stage iPD, indicative of an imbalance in the differentiation process within CD8 compartments.
We next investigated the functional states of various T-cell subsets. T-bet is an essential marker for effector CD8 lymphocyte functions32,33. Consistently, CD8 T cells from iPD vs HC also displayed a higher frequency of T-bethigh and CD45RO−T-bet+ CD8 T cells (Fig. 2e, f). T-bet+ cells were mainly CD45RO−, but not CD45RO+ cells, indicating that those cells were mostly CD45RA+ terminally differentiated cells, based on the largely mutually exclusive relationship between CD45RA and CD45RO expression (Supplementary Fig. 3C). The increased CD8 effector function was further consolidated by an augmented proliferation and activation levels among CD45RO− CD8 T cells, as quantified by the expression of Ki67 (Fig. 2g) and HELIOS34 (Fig. 2h), respectively. The increase of CD45RO−CD57+ (Fig. 2i) and CD57+ (Supplementary Fig. 3H) portions among CD8 T cells further supports an advanced differentiation state35 of CD8 T cells in iPD.
We further analyzed whether CD8 T cells from iPD show any signs of exhaustion36, senescence37 or other major dysfunctionality, as our participants were mainly over 60. We did not observe any sign of exhaustion in the expression of several key T-cell exhaustion markers such as PD-1, CTLA4 and LAG3 (Supplementary Fig. 3I–K). Moreover, the CD8 T-cell senescence marker KLRG1 was also unchanged between iPD and HC (Supplementary Fig. 3L). The activation marker ICOS was significantly decreased in CD8 T cells (Supplementary Fig. 3M). However, this decline only reflected the decrease in the frequency of CD8 TCM cells (Supplementary Fig. 3D), as only the frequency of ICOS+CD45RO+, but not of ICOS+CD45RO− cells was lessened (Supplementary Fig. 3N). We also observed a decreased expression of the amino acid transporter CD98 (Supplementary Fig. 3O), nevertheless being independent of CD45RO expression levels (Supplementary Fig. 3P). Together, these data support that CD8 T cells in early-to-mid stage iPD exhibit a terminally-differentiated but non-exhausted and non-senescent state.
Effector CD8 T cells tend to migrate to the tissues where an active immune response occurs38. Therefore, we also assessed the expression of several lymphocyte-trafficking relevant chemokine receptors, such as CXCR3, CCR4 and CCR6. In the blood of iPD, we observed a lower frequency of CXCR3+ and CCR4+ cells and their corresponding expression levels (geometric mean, MFI) among total CD8 T cells, whereas CCR6 showed no difference (Supplementary Fig. 4A, B). The decrease in CXCR3 and CCR4 expression might be attributable to the decline in the frequency of TCM cells among CD8 T cells (Supplementary Fig. 3D), since among all the four subsets, CD8 TCM cells displayed the highest expression levels of these two chemokine receptors (Supplementary Fig. 4C). The expression of the three individual receptors on CD8 TEMRA was low relative to that on other CD8 memory subsets, but without showing clear difference in iPD vs HC (Supplementary Fig. 4C). The frequency of CCR4 and CCR6 co-expressing cells was also decreased in iPD (Supplementary Fig. 4D), although the cells expressing both receptors was sparse among CD8 T cells (<1% in average). Accordingly, the frequency of cells lacking all the tested three chemokine receptors (CXCR3, CCR4 and CCR6) was increased in iPD vs HC (Supplementary Fig. 4D). We further asked whether CD8 TEMRA express one of the major homing receptors allowing them to properly migrate to the central nervous system (CNS). As integrin alpha 4 (also known as CD49d) is the major brain homing factor of peripheral CD8 T cells, controlling trafficking of CD8 T cells into the CNS39,40, we analyzed CD49d expression. However, we did not find any significant difference in CD49d expression or the frequency of CD49d+ cells between iPD and HC among various subsets of CD8 T cells, including CD8 TEMRA (Supplementary Fig. 4E, F). In brief, our data suggest that CD8 TEMRA have an uncompromised potential to migrate into the CNS.
Since CD8 T cells and NKT are known to have cytotoxic potential and also secrete cytotoxic effector molecules, we asked whether a more-cytotoxic state was already reflected by serological parameters. To this end, we measured two of the most abundant extracellular cytotoxic effector molecules, granzyme A (GZMA) and granyzme B (GZMB), together with the membrane-pore forming molecule perforin. Highly encouragingly, circulating GZMA was significantly elevated in iPD, especially in females (Fig. 2j), while the levels of both GZMB and perforin were comparable between iPD and HC (Fig. 2k, l). In short, our serological data also suggest a more-cytotoxic immune milieu in iPD.
Although naive CD4 cells displayed no difference, the frequency of CD4 TCM cells was slightly, but significantly decreased (Supplementary Fig. 5A, B). Moreover, the portion of intermediate CD4 (CCR7−CD27+CD45RO−)41 T cells was increased in iPD (Supplementary Fig. 5C). Unlike CD8 T cells, the frequency of effector and TEMRA among total CD4 T cells remained unchanged (Supplementary Fig. 5D, E). Together, CD8, but not CD4 T cells, favor a terminally-differentiated effector program over the generation of a long-lived central memory T-cell profile in iPD.
The combination of the chemokine receptors, such as CXCR3, CCR4 and CCR6 can also be used to distinguish the CD4 T helper lineages Th1, Th2 and Th17 (Supplementary Fig. 5F). By analyzing the combinations of these chemokine receptors and also assessing the expression of the CD4 Th master transcription factors T-bet, GATA3 and RORγT, we observed neither a significant change in the frequency of Th subsets, nor in the ratios between those subsets in iPD vs HC (Supplementary Fig. 5G–J). No change in any Th subset was also in agreement with our CyTOF results (Supplementary Table 5). We also did not observe any significant difference in the Th1 and Th2 cytokines, such as GM-CSF, IFN-γ, IL-5 and IL-13 measured in the sera of our cohort (Supplementary Fig. 5K). Furthermore, since IL-33 plays a crucial role in regulating the ILC2 response and serves as an alarmin in CNS42, we quantified the circulating IL-33 level and found it was similar between PD and HC (Supplementary Fig. 5K). We further measured several CNS-homing relevant chemokines, such as IL-8, IP10 and MCP-1. None of them showed any significant difference (Supplementary Fig. 5K). Together, our data from various viewpoints firmly suggest that the CD4 composition is less affected in early-to-mid stage iPD.
Early-to-mid stage PD shows a reduced frequency of CD8 Treg but not CD4 Treg
CD4 regulatory T cells (CD4 Treg) play an important role in suppressing effector T-cell responses to avoid an overshooting immune reaction that could harm surrounding self-tissues43,44,45. Previous reports have shown a reduced frequency of CD4 Treg and/or an impaired suppressive capability of Treg in PD and other neurodegenerative diseases7,46,47. Through our FCM analysis, however we did not observe any significant change in the frequency of CD4+FOXP3+ T cells (CD4 Treg) among CD4 T cells between PD and HC (Fig. 3a). Unexpectedly, we found that the frequency of CD8+FOXP3+ Treg (CD8 Treg) was reduced in iPD, especially for females (Fig. 3b, c). CD8 Treg are able to suppress the responses of effector cells, including CD8 cells [reviewed elsewhere48,49]. In line with our observation of reduced CD8 Treg, the expression of other markers relevant to CD8 Treg, such as CD25 (IL2RA) and CD122 (IL2RB)50, was also reduced in CD8 T cells of iPD vs HC (Fig. 3b). The ratio between the frequency of CD8 TEMRA and of CD8 Treg was increased, with a mean of 70.62 in PD vs 30.37 in HC, further highlighting an effector-dominant, and consequently more ‘active’ CD8 T-cell compartment, especially for females (Fig. 3d). Although the frequency of CD4 Treg was unchanged in our iPD, this data alone could not exclude the possibility of a compromised suppressive function of CD4 Treg in iPD. Therefore, we also analyzed several functional markers, such as CD45RO and phospho-S6 [pS6, reflecting mTORC1 activity51]. The levels of both CD45RO and pS6 were decreased among CD4 Treg in iPD (Fig. 3e). Despite the reduction in those markers, the expression of FOXP3 and CTLA4, which are decisive for maintaining suppressor function45, remained unchanged among CD4 Treg (Fig. 3f). These data indicate that it is mainly CD8 Treg cellularity that was impaired in early-to-mid stage, especially female iPD, whereas the CD4 Treg frequency and suppressive capacity were unlikely affected.
a Scatter dot plot showing the frequency of FOXP3+ cells among CD4 T cells. b Scatter dot plots showing the frequency of FOXP3+ cells for all or females or males, CD25+ cells and CD122+ cells among CD8 T cells. c Representative FCM plots showing the reduced frequency of FOXP3+ CD8 T cells (CD8 Treg) in PD vs HC. d Scatter dot plots showing the ratio between CD8 TEMRA (CD45RO–CCR7−CD27−) and CD8 Treg for all (left) or females (middle) or males (right). e Scatter dot plots showing the frequency of CD45RO+ and CD45RO+pS6+ among CD4 Treg. f Scatter dot plots showing MFI of FOXP3 and CTLA4 among total CD4+FOXP3+ Treg. The results (a, b, d–f) were analyzed using unpaired two-tailed Student’s t test. Data are presented as mean ± standard deviation (s.d.). Each symbol represents the measurement from one individual (a, b, d–f). ns or unlabeled, not significant; all significant P-values are indicated. FCM flow cytometry, MFI geometric mean (geomean), Treg regulatory T cells, HC healthy controls, n = 24, PD patients with Parkinson’s disease, n = 28. For b and d, female HC, n = 10; female PD, n = 8; male HC, n = 14; male PD, n = 19. Source data are provided as a Source Data file.
TEMRA & the ratios between TEMRA and specific subsets as early diagnosis markers
Considering the substantial difference in CD8 TEMRA between iPD and HC, we evaluated the potential of using CD8 TEMRA as a peripheral diagnostic biomarker. We first analyzed the possible correlation between the frequency of CD8 TEMRA and various available quantitative clinical parameters, such as age, disease duration, H&Y staging scale, UPDRS-III (Unified Parkinson’s Disease Rating Scale III), LEDD (Levodopa equivalent daily dose) and MOCA (Montreal Cognitive Assessment). Although no significant correlation was observed between CD8 TEMRA and most of those clinical data, the frequency of CD8 TEMRA in PD showed a significant negative correlation with disease duration from initial symptoms (Fig. 4a) or clinical diagnosis (Fig. 4b). This suggests that CD8 TEMRA populations might be more involved in the early rather than the later stage of iPD. Samples from early-to-mid stage iPD could already be well distinguished from that of HC with an area under the curve (AUC) of 0.7662, based on CD8 TEMRA frequency alone (Fig. 4c). As the frequency of CD8 TEMRA negatively correlated with disease duration, we applied another receiver operating characteristic (ROC) analysis by focusing on patients diagnosed only within 5 years. An excellent AUC value of 0.8580 [the category of AUC ranges is defined elsewhere52] was obtained (Fig. 4d).
Correlation between the frequency of CD8 TEMRA (CD45RO–CCR7−CD27−) and the disease duration from initial symptoms (a) or clinical diagnosis (b). Coefficient and P-value based on Spearman correlation (n = 28). ROC analysis based on the frequency of CD45RA+CD45RO−CCR7− for all PD (n = 28) (c) or patients diagnosed within 5 years (n = 11) (d) vs all HC (n = 24). e The frequency of CD8 TEMRA (CD45RA+CD45RO−CCR7−) vs that of CD8 Treg. f The frequency of CD8 TEMRA (CD45RA+CD45RO−CCR7−) quantified by FCM vs the frequency of ILC2 quantified by CyTOF. g The frequency of CD8 (CD45RA+CD45RO−CCR7−) vs the ratios between CD8 TEMRA (CD45RO–CCR7−CD27−) and CD8 Treg. The dashed red circle highlights the PD-dominant area. h ROC analysis based on the frequency of CD8 TEMRA (CD45RA+CD45RO–CCR7–) only from female PD (n = 8) or HC (n = 10). i, j ROC analysis based on the ratios between the frequency of CD8 TEMRA (CD45RO–CCR7−CD27−) and of CD8 Treg, for all iPD (n = 28) vs HC (n = 24) (i) or with only female iPD (n = 8) and HC (n = 10) (j). k–n ROC analysis based on the ratios between the frequency of TEMRA (CD45RA+CD45RO−CCR7−) and of TCM (CD45RO+CCR7+CD27+), for all iPD (n = 28) (k) or patients diagnosed within 5 years (n = 11) (l) vs all HC (n = 24) or only from female (m) iPD (n = 8) and HC (n = 10) or only from male (n) iPD (n = 19) and HC (n = 14). The scatter dot plot was displayed in the right panels (l–n). P-values were based on unpaired two-tailed Student’s t test. Data are presented as mean ± standard deviation (s.d.). Each symbol represents one individual. All the PD (n = 28, where one female PD sample was excluded as the same patient visited twice within a short period) and HC (n = 24) were used unless specified. P-values displayed under each ROC analysis test the null hypothesis of the AUC = 0.50 with a two-tailed test. AUC area under curve, ILC2 type 2 innate lymphoid cells, ROC receiver operating characteristic, Treg regulatory T cells, TCM central memory. Source data are provided as a Source Data file.
Since the frequency of CD8 Treg and ILC2 also showed a difference between iPD and HC, we assessed whether those subsets could provide additional diagnostic power (Fig. 4e, f). As the ratio between CD8 TEMRA and CD8 Treg might represent the CD8 effector potential in the given individual, we also integrated the ratios and the frequency of CD8 TEMRA (Fig. 4g). Neither the frequency of CD8 Treg, the frequency of ILC2, nor the ratios between CD8 TEMRA and CD8 Treg substantially increased the potential of implementing CD8 TEMRA as a biomarker, always leaving seven HC samples mixed in the ‘PD area’ (Fig. 4g). A detailed analysis into the association between CD8 TEMRA to Treg ratios and the frequency of CD8 TEMRA revealed that six out of the seven ‘indistinguishable’ HC samples were males’. Considering this and other gender-biased results, we performed another ROC analysis with only female samples. Encouragingly, an outstanding AUC value of 0.9125 was achieved (Fig. 4h). The ROC analysis using the ratio between CD8 TEMRA and CD8 Treg also showed an excellent diagnostic potential when only females were analyzed (Fig. 4i, j). As the ratio between CD8 TEMRA and TCM was substantially increased in iPD, we also performed the ROC analysis using this ratio. Even with all the patients, an excellent AUC value of 0.8503 was obtained (Fig. 4k). With those diagnosed within five years or only females, an outstanding AUC value of 0.9470 or 0.9250 was accomplished, respectively (Fig. 4l, m). Even if only males were analyzed, an AUC value of 0.8120 was obtained (Fig. 4n). In brief, our data firmly support that CD8 TEMRA alone as well as the ratios between CD8 TEMRA and CD8 TCM or Treg are reliable, but yet unrecognized, easily-accessible cellular markers for early diagnosis of iPD, especially females.
Validation in another subcohort & discovery of enhanced cytotoxic pathways in CD8 T cells
Given the high clinical diagnostic potential, we further validated our key findings regarding CD8 T cells in another subcohort (12 HC vs 11 iPD) of the Luxembourg Parkinson’s study22 using the FCM analysis (Supplementary Tables 2 and 6). The total CD8 T cells were unchanged among total CD3 or living lymphocyte singlets (Fig. 5a). Although all the samples of the validation subcohort were cryopreserved, the frequency of CD8 TEMRA was still clearly increased while the frequency of CD8 TCM was decreased in early-to-mid stage iPD vs HC (Fig. 5b, c). The ratio between CD8 TEMRA and TCM was still much higher in early-to-mid stage iPD than in HC (Fig. 5d). However, the ratio was generally lower in cryopreserved materials relative to fresh blood samples, possibly due to the cryopreservation effects23,53.
a Scatter dot plots showing the frequency of total CD8 T cells among CD3 or total living lymphocyte singlets using the FCM analysis of cryopreserved PBMC from another subcohort of the Luxembourg Parkinson’s study (11 iPD vs 12 HC, all CMV+ except for two seronegative HC). Scatter dot plots showing the frequency of CD8 TEMRA (b) or TCM (c) among total CD8 T cells using different marker combinations. Of note, one sample was excluded from PD as it was identified as an outlier by the default setting of the ROUT method of Graphpad. d Scatter dot plots showing the ratios between CD8 TEMRA and CD8 TCM using different marker definitions. Scatter dot plots showing the frequency of GZMA+ (e), GZMB+ (f), perforin+/PRF1+ (g) or GZMK+ cells (h) among various CD8 T subsets. i Scatter dot plots showing the frequency of GZMA+GZMB+PRF1+GZMK– cells among total CD8 T cells or CD8 subsets. j Gating strategy was used to identify different CD8 subsets by FCM (BD FACSymphonyTM S6). k Representative FCM plots showing the expression of GZMA, GZMB, PRF1 and GZMK among CD8 TEMRA. l Scatter dot plots showing MFI of perforin among total CD8 cells or CD8 subsets. The results in a–i, l were analyzed using ordinary one-way ANOVA test with two-stage linear step-up procedure of Benjamini/Krieger/Yekutieli correction. q-values (FDR) were displayed. Data are presented as mean ± standard deviation (s.d.). Each symbol represents the measurement from one individual. all q-values are indicated. CMV Cytomegalovirus, FCM flow cytometry, GZMA/B/K Granzyme A/B/K, HC healthy controls, n = 12 including 5 males; PD patients with Parkinson’s disease, n = 11 including 8 males. Source data are provided as a Source Data file.
To further investigate whether CD8 T cells harbor a higher cytotoxic potential in iPD, we analyzed the expression of several cytotoxicity-related markers within CD8 subsets using the FCM analysis. The frequency of GZMA−, GZMB− and perforin (PRF1)-expressing cells gradually increased along the differentiation trajectory of CD8 T cells from naive (Tn) to TCM, TEM and TEMRA (Fig. 5e–g). The frequency of GZMB-expressing cells was selectively enhanced only in TEM and TEMRA of iPD vs HC. In contrast, the frequency of GZMK-expressing cells was reduced in TEM and TEMRA (Fig. 5h). As the frequency of circulating exhausted-like GZMK+CD8+ increases during ageing54, a reduction in GZMK expression might point towards less-exhausted and more-active ‘younger’ CD8 T cells in iPD vs HC. We next quantified the combinations defined by the expression of several cytotoxic markers in CD8 T cells, as the synergistic effects might be required between relevant cytotoxic molecules to effectively perform a cytotoxic function55. Encouragingly, the frequency of GZMA+GZMB+PRF1+GZMK– cells was increased among total CD8 T cells, especially in TEMRA and TEM in iPD vs HC (Fig. 5i). Indeed, a largely mutually exclusive expression pattern between GZMK and other examined cytotoxic molecules (GZMA, GZMB and PRF1) was observed within CD8 TEMRA (Fig. 5j, k). Although the frequency of PRF1-expressing cells was not significantly changed in CD8 subsets, the MFI of PRF1 still showed a trend to be augmented in iPD vs HC among total CD8 T cells, TEM and TEMRA (Fig. 5l). Collectively, these cytometry data reveal that the cytotoxic function of CD8 T cells, especially the synergy among cytotoxic molecules of CD8 TEMRA and TEM, is enhanced in early-to-mid stage iPD.
To more comprehensively characterize the transcriptional program in CD8 T cells and gain mechanistic insights, we performed single-cell RNA-sequencing (scRNA-seq) based on the four sorted subsets of CD8 T cells. The reason to first sort the four subsets prior to scRNA-seq is that the decisive markers (CD45RA/CD45RO isoforms) of our target cell types are encoded by the same gene, which otherwise cannot be easily distinguished by the widely-employed shotgun sequencing methods. Considering our female-biased results in CD8 T cells, we selected only females for scRNA-seq (Supplementary Table 2). After various quality control steps, a total of 24,832 fluorescence-activated cell sorting (FACS)-isolated individual CD8 T cells from cryopreserved PBMC of iPD (n = 5) and HC (n = 4) were kept for further analysis (Fig. 6a, Supplementary Fig. 6A). Mapping the four sorted subsets with the two-dimensional uniform manifold approximation and projection (UMAP) showed a clear separation between CD8 Tn, TCM and the remaining CD8 subsets while TEM and TEMRA were partially overlapped (Fig. 6b). As expected, a higher expression of CCR7, CD27, CD28, IL7R and SELL/CD62L was observed in CD8 Tn and TCM compared with the later-differentiated ones (i.e., CD8 TEM and TEMRA) (Fig. 6c). The same expression pattern was detected for TCF7 (also known as TCF1) and TOX between terminally-differentiated exhausted CD8 T cells56 and CD8 TEMRA compared to less-differentiated subsets (Fig. 6c). In line with our FCM results, the expression of several cytotoxic effector molecules, such as GZMA, GZMB, GZMH, GZMN, PRF1, NKG7, KLRD1, KLRG1, and especially GNLY, was enhanced in CD8 TEM and TEMRA of iPD vs HC (Fig. 6d). The portions of the cells co-expressing two or more various cytotoxic or effector molecules, such as GZMA, GZMB, PRF1 and IFNG, were also increased in iPD (Fig. 6e). This pattern mainly occurred in CD8 TEM and TEMRA without pre-stimulation (comparing Fig. 6b with Fig. 6f–i, Supplementary Fig. 6B–D). The enrichment analysis revealed that the pathways involved in graft-versus-host disease, antigen processing and presentation (e.g., CD74 and TAPBP), leukocyte trans-endothelial migration and cell adhesion (e.g., ITGB2, ITGB7, ITGAL, CD2 and CD226) and NK-mediated cytotoxicity (e.g., GZMB, PRF1, KLRD1, ITGAL, ITGB2 and PLCG2) were most altered among the upregulated genes in CD8 TEMRA of iPD vs HC (Fig. 6j–m). The gene THEMIS is a positive regulator of TCR signaling strength of peripheral CD8 T cells in response to low-affinity self-peptide MHC and cytokines57. THEMIS was among the most upregulated genes in CD8 TEMRA of iPD (Fig. 6j). Concurrently, several genes, critically involving in the mRNA de-adenylation-dependent degradation (enrichment P = 7.63E-8 among downregulated genes, Fisher’s exact test, GO process), such as CNOT6L58, ZFP3659 and ZFP36L2, were substantially downregulated in CD8 TEMRA of iPD vs HC (Fig. 6J, Supplementary Fig. 6E). The downregulation of this decay process potentially increases the mRNA stability in CD8 TEMRA to be better equipped for the fast-track translational steps.
a Schematic for CD8 scRNA-seq. The other gates are in Supplementary Fig. 6a. b UMAP showing distribution of n = ~25,000 cells. Violin plots of selected genes defining CD8 naive/memory subsets (c), signifying cytotoxicity (d) and involving in top-ranked pathways enriched among upregulated genes in CD8 TEMRA (m). e Balloon plot showing the percentages of cells co-expressing the indicated markers. f, h UMAP showing joint density of GZMA+GZMB+PRF1+ (f) and GZMA+GZMB+PRF1+IFNG+ (h) in all CD8 cells. UMAP showing the individual GZMA+GZMB+PRF1+ (g) or GZMA+GZMB+PRF1+IFNG+ (i) cells among CD8 subsets. For visual comparability, random downsampling was employed to display the same number of cells for different conditions and subsets. j Volcano plots showing the expression changes in CD8 TEMRA. The selected top up- or downregulated DEGs are labeled in red or blue. Vertical dashed line indicates the log2FC value of 0.5 or −0.5, while the horizontal red line indicates –log10(0.05). k Top-ranked enriched KEGG pathways among upregulated DEGs in CD8 TEMRA of iPD vs HC. l Top-ranked GSEA pathways in upregulated DEGs in CD8 TEMRA. The lower part showing the rank distribution of the genes involved in the indicated pathway. The list on the right showing the leading-edge genes. n UMAP showing the unsupervised clustering analysis for n = ~7200 CD8 TEMRA cells. The clustering results were split for iPD or HC cells. Cells were randomly down-sampled to the same number for iPD and HC. o Violin plots of selected markers distinguishing clusters. p Heatmap of selected top DEGs in the clusters of CD8 TEMRA. P-values in j and p were analyzed using two-tailed nonparametric Wilcoxon Rank Sum test adjusted based on Bonferroni correction. In p, only the genes with adjusted P-value ≤0.05 were considered. P-values in k and l were analyzed using the one-tailed Fisher’s exact test and the one-tailed empirical phenotype-based permutation test, respectively. DEG differentially expressed genes, FACS fluorescence-activated cell sorting, FC fold change, GZMA/B/K Granzyme A/B/K, IFNG interferon gamma, UMAP uniform manifold approximation and projection. Panel a was created with BioRender.com.
Our clustering analysis within CD8 TEMRA identified two clusters (C0 and C1) (Fig. 6n). C0 expressed genes encoding proteins with effector or TCR signal activating functions [CAMK4, GZMK, COTL160 and THEMIS] (Fig. 6o). Although GZMK can be considered an exhaustion marker either during ageing54 or among CCR7+ TCM61, it is also used to identify TEM-like populations in CCR7− cells. In line with the latter, the two key exhaustion markers HAVCR2/TIM3 and CTLA4 were barely expressed in CD8 TEMRA (Supplementary Fig. 6F). The expression of the other exhaustion markers TIGIT or LAG3 was either much lower in C0 or comparable between patients and controls, respectively (Supplementary Fig. 6F). Meanwhile, PDCD1/PD-1 mRNA was undetectable in our samples. Thus, our data support C0 as a more TEM-like and non-exhausted subset. C1 was characterized by a high expression of several NK functional markers [KLRC2, KLRC3, FCGR3A/CD16, IKZF2/HELIOS and TYROBP/DAP1262] (Fig. 6o). The analysis of differentially-expressed genes (DEG) in iPD vs HC further demonstrated that the genes involved in NK-mediated cytotoxicity, such as GZMB, PRF1 and FCGR3A, were among the most upregulated ones in C1 (Fig. 6p). Several MHC type II molecules (e.g., HLA-DRB1, HLA-DQA1 and HLA-DQB1), known CD8 T-cell activation markers63 as well as another NK cytotoxic effector molecule (KLRD1) and leukocyte trans-endothelial migration and adhesion molecules (ITGB2 and ITGAL) were upregulated in C0 and, to a lesser extent in C1, of CD8 TEMRA (Fig. 6p). Together, these scRNA-seq data support that the cytotoxic pathways are enhanced in CD8 TEMRA from early-to-mid stage iPD vs HC.
Similar pathways enriched for the upregulated genes in iPD vs HC were observed in CD8 TEM (Supplementary Fig. 6G–I). Our volcano plot (Supplementary Fig. 6G) showed that, AOAH (acyloxyacyl hydrolase), a lysosomal enzyme known to detoxify LPS64 and several established cytotoxic molecules (GZMH, NKG7 and KLRG1), were substantially upregulated in CD8 TEM of iPD vs HC. Our clustering analysis within CD8 TEM identified three clusters (C0, C1 and C2) (Supplementary Fig. 6J). Both C0 and C2 of CD8 TEM produced IFNG mRNA without pre-stimulation, while only C1 expressed the markers of CD161+/KLRB1+ tissue-homing CD8 cells, including IL23R, CCR6, SLC4A10, ZBTB16, COLQ and CEBPD65,66 (Supplementary Fig. 6K). C0 expressed high levels of both cytotoxic molecules (GZMB, GZMH and GNLY) and the terminal-differentiation-related transcription factor ZEB267, indicative of a more terminally-differentiated status (Supplementary Fig. 6K). C2 showed an expression profile reminiscent of NK-like cells (including KLRC2, KLRC3, TYROBP/DAP12, LYN and IKZF2/HELIOS). Our DEG analysis shows that many cytotoxic molecules such as GZMA, GZMB, GZMH, GZMM, NKG7, GNLY and KLRD1 as well as several MHC-II genes were upregulated in at least one of the clusters of CD8 TEM in iPD vs HC (Supplementary Fig. 6L). TGFB1 and its downstream signaling molecules, such as SMAD7, SKI and SKIL, which inhibit the expression of CD8 key effector molecules68, were all substantially downregulated in all the clusters of CD8 TEM (Supplementary Fig. 6M). This again indicates a more-active effector status. Concomitantly, pro-survival or anti-apoptotic genes such as LMNA, GADD45B and BCL2 were downregulated, indicative of more short-lived CD8 TEM in iPD (Supplementary Fig. 6M). In total, the cytotoxic pathways in CD8 TEM were also upregulated in early-to-mid stage iPD vs HC.
Reprogrammed CD8 Tn & TCM to favor terminal differentiation in early-to-mid stage iPD
Since the ratios between CD8 TEMRA and TCM were substantially enhanced in iPD, while total CD8 T cells were kept intact, we reasoned that the CD8 differentiation process might be affected. To unbiasedly quantify the CD8 differentiation process, we performed pseudotime trajectory analysis in our CD8 scRNA-seq datasets using the tool Slingshot69. Among the most significantly-differentially-expressed genes along the CD8 differentiation trajectory, the marker genes known to change over the CD8 differentiation stages, such as IL7R, LEF1 and GZMA, were confirmed by our unsupervised analysis (Fig. 7a). Encouragingly, the mean pseudotime increased from pre-sorted CD8 Tn to TCM, TEM and TEMRA (Fig. 7b). The relative fractions of cells for several mid-to-later pseudotime windows (>60) were much higher in CD8 T cells of iPD vs HC, while those in earlier windows (20-50) were lower (Fig. 7c, d). Hence, our unsupervised pseudotime analysis also demonstrates an accelerated differentiation process within CD8 T compartments of early-to-mid stage iPD.
a Heatmap showing the top 20 most-increased and -decreased genes along the predicted pseudotime trajectory ranked using log2FC (if p-value < 0.01) (n = ~25,000 cells). Each column represents one single cell. b Two-dimensional representation showing the CD8 memory T-cell differentiation trajectory using UMAP (upper). Lower, the ridge plots of pseudotime distribution among different CD8 subsets along differentiation trajectory. c Percentages of each binned psedutotime window along the pseudotime trajectory in iPD (n = ~12,000 cells) or HC (n = ~13,000 cells). d Stacked barplot showing the relative fraction of cells at the given predicted differentiation stage (during the given pseudotime window) between iPD and HC. e Unsupervised clustering analysis of CD8 TCM. Percentages of the two clusters and the ratios between two clusters were displayed. For comparability in visualization, cell numbers were randomly down-sampled to the same number for different conditions and subsets. f Violin plots of selected marker genes distinguishing different clusters within CD8 TCM. g Heatmap of selected top up- or downregulated genes in CD8 TCM from iPD vs HC. h Balloon plot showing the percentage of cells co-expressing the indicated markers in the given subset. If a cell shows the read count equal to or higher than 1 for each of the markers in the indicated combination, it is regarded as the cell co-expressing the set of markers. i UMAP plot showing joint density of CCL5 and GZMK (left) and of GZMA and GZMK (right) in all the individual CD8 T cells. Arrows refer to the area of TCM co-expressing the indicated markers. For a, the one-tailed ANOVA method was used to test whether any of the spline coefficients (for pseudotime fitting) is non-zero. In g, differential expression was analyzed using two-tailed nonparametric Wilcoxon Rank Sum test based on Bonferroni correction (only those with adjusted P-value ≤0.05 were considered). Each dot represents one single cell in (b, e, i). FC fold change, HC healthy controls, PD patients with Parkinson’s disease, UMAP uniform manifold approximation and projection.
As pseudotime analysis showed an effect already at an earlier differentiation stage, we next assessed transcriptome in CD8 TCM to characterize potential mechanisms. An unsupervised analysis within CD8 TCM identified two clusters (C0 and C1). The more-active cluster C0 was characterized by high levels of cytotoxic or proinflammatory effector molecules (e.g., GZMA and CCL5), while the resting cluster C1 was devoid of the expression of those effector molecules. The ratios between the more-active cluster C0 and the relatively-resting cluster C1 almost doubled in iPD than HC (Fig. 7e, f). The integrin ITGA4, known to have a higher expression in CD8 TEM compared to TCM, was highly expressed in C0 (Fig. 7f). The less-active C1 also expressed the TCR inhibitory docking adapter (GAB2)70 and the transcription factor KLF7 that is highly expressed in CD8 Tn71 (Fig. 7f). Moreover, C1 over-expressed the antigen-specific memory T-cell marker CD40LG in iPD vs HC (Fig. 7f), again indicative of C1 as a quiescent TCM subset. Meanwhile, MAP3K1/MEKK1, the negative regulator of virus-specific CD8 T cells72, was mainly expressed in C1, especially for iPD (Fig. 7f). In summary, our clustering analysis alone already indicates a more-active status of CD8 TCM in iPD.
To identify the potential transcriptional change in iPD vs HC, we further performed DEG analysis within CD8 TCM. First, a decrease in the fraction of cells expressing GAB2 was observed in C1, further indicating a more-active status of CD8 TCM (Fig. 7f). Meanwhile, THEMIS was upregulated in both clusters of CD8 TCM of iPD (Supplementary Fig. 7A). Several critical transcription factors (ID3, LEF1 and ETS1) that favor more the memory T-cell development71,73,74 were also upregulated (Fig. 7f, g, Supplementary Fig. 7A). Furthermore, our DEG analysis showed that three out of five somatic linker histone H1 family genes (HIST1H1C/H1.2, HIST1H1D/H1.3 and HIST1H1E/H1.4)75 and several other histone genes showed a substantial upregulation in CD8 TCM of iPD (Fig. 7g). Specific subtypes of the liner histone genes (H1.2 and H1.4) have been recently shown to facilitate the differentiation of another terminally-differentiated immune cell type (i.e., neutrophils)76. This indicates a reprogramming of the chromatin-access status within CD8 TCM as the epigenetic landscape is known to substantially change during the CD8 memory differentiation process77. Although CCR7+GZMK+ cells have been identified as exhausted-like CD8 memory T-cell progenitors61, none of the analyzed exhaustion markers (TIGIT, TIM3/HAVCR2, LAG3 and CTLA4) showed an expression in GZMK+CCL5+GZMA+ C0 in our TCM dataset (Supplementary Fig. 7A). We further examined whether GZMA, CCL5 and GZMK are co-expressed within the same cells of CD8 TCM. The scRNA-seq confirmed that CCL5 and GZMK were co-expressed in the CD8 TCM cluster area (Fig. 7h, i). GZMK and GZMA were also co-expressed, although only in a small fraction of CD8 TCM (Fig. 7h, i). Consistent with our FCM data, the expression between GZMK and GZMB was largely mutually exclusive (Supplementary Fig. 7B). In short, our unbiased scRNA-seq from different angles demonstrate that CD8 TCM were already more ready, which might promote the terminal differentiation in early-to-mid stage iPD.
To evaluate whether the accelerated differentiation originates even earlier, we next systematically analyzed scRNA-seq data within CD8 Tn. The volcano plot shows that the early-differentiation transcription factor LEF174 was among the most highly-upregulated genes in total CD8 Tn of iPD vs HC (Supplementary Fig. 7C). The pathway enrichment analysis revealed that RUNX1-regulated transcription related to myeloid differentiation was ranked top (odd ratio of 84.68, P-value = 4.71E-4) (Supplementary Fig. 7D). More precisely, RUNX1 and RUNX2, reported to favor CD8 terminal differentiation78, were upregulated in CD8 Tn of iPD vs HC (Supplementary Fig. 7E, F). Concurrently, RUNX3, known to suppress the expression of both RUNX1 and RUNX2 genes and terminal differentiation79, was slightly reduced in CD8 Tn (Supplementary Fig. 7F). TCR signaling was also among the top-ranked pathways enriched for upregulated genes in CD8 Tn of iPD vs HC (Supplementary Fig. 7D). Our clustering analysis identified two clusters in CD8 Tn (Supplementary Fig. 7G). C0, the main cluster (>95% of CD8 Tn), expressed RUNX2, while C1 expressed the effector molecule CCL5 and more abundantly the other two markers (e.g., FYN and ITGB1/CD29) that are highly expressed in memory CD8 T cells relative to CD8 Tn80,81. This indicates C1 as already a more memory-like CD8 subset (Supplementary Fig. 7H). The DEG analysis within different clusters also showed that LEF1 and THEMIS were upregulated in CD8 Tn of iPD vs HC. Moreover, similar to the scenario in other CD8 subsets, several members of the linker histone H1 family genes showed a substantial increase in both clusters of CD8 Tn, while CNOT6L was still among the most downregulated genes (Supplementary Fig. 7H). Together, our unbiased scRNA-seq analysis indicates that early-to-mid stage iPD already initiates and reprograms the transcriptional machinery that facilitates accelerated differentiation from as early as CD8 Tn.