Machine learning combined with omics-based approaches reveals T-lymphocyte cellular fate imbalance in abdominal aortic aneurysm | BMC Biology

Single-cell transcriptome profiles of human AAA

As shown in the analysis process (Fig. 1A), data from 4 aneurysmal abdominal aorta samples in the database were included in this research. Following data integration and quality filtering, 7343 cells were projected onto a UMAP plot based on integrative unsupervised clustering analysis (Fig. 1B), of which 16 populations were obtained and 9 main cell types were identified and assigned to each population based on CellMarker database and SingleR algorithm (Fig. 1C). The levels of cell type-specific canonical markers were presented in Fig. 1D.

The highest cellular population was T cells followed by monocytes/macrophages, whereas smooth muscle cell (SMC) and endothelial cell (EC) populations were present in small numbers (Fig. 1E), in agreement with previous findings of cellular responses of AAA cells at the tissue and single-cell level in mice [11, 12]. We applied the cell-chat analysis to infer intercellular signal transduction during AAA progression. T lymphocytes, with the heaviest weights, played central roles in the intersected cellular communication network (Fig. 1F, Additional file 1: Fig. S1D). Cell cycle phase scores were assigned based on the gene expression signatures for the G2, S, and G2/M phases. Immune cells presented high proliferation scores (Fig. 1G), and the highest proliferating cells were observed in the T subpopulation (Additional file 1: Fig. S1F). In contrast, SMCs and ECs showed reduced proliferation activities, evidenced by the lowest proliferation scores (Additional file 1: Fig. S1F).

The top markers (sorted by average Log (fold change)) for a cell type were identified relative to all other types. The identified neural cells featured the high expression of WFDC2, FXYD2, and KLK1 (Additional file 1: Fig. S1A). We demonstrated that the immune cell lineages are comprised: (i) B cell (clusters 2, 5); (ii) CD4 + T (clusters 0, 1, 3); (iii) CD8 + T (clusters 7, 10); (iv) mast cell (cluster 15); (v) NK (cluster 11); (vi) monocyte/macrophage (clusters 4, 6, 8). Moreover, non-immune cell lineages were observed, namely endothelial cell (EC) (cluster 12), smooth muscle cell (SMC) (cluster 13), and neural cell (Schwann cell) (cluster 9) (Additional file 1: Fig. S1B) expressing NGFR associated protein [14]. Furthermore, the significant correlation of the 2000 most variable genes was identified to delineate the relationship among cell clusters in AAA. Histiocytes were clustered together, and the close relationships between B and macrophages as well as T and NK were demonstrated (Additional file 1: Fig. S1C). The functional analysis suggested proliferation was suppressed in histiocytes, and the immune state was activated in immune cell lineages for AAA patients (Additional file 1: Fig. S1E). The sustained immune response inhibited SMC proliferation and upregulated cellular impairment, promoting AAA progression [15]. Moreover, SMC showed high expression of collagen and proteoglycan genes (such as FN1, COL1A2, DCN) (Additional file 1: Fig. S1A). SMC and EC both exhibited cellular damage and collagen-containing extracellular matrix degeneration hallmarks, suggesting they underwent phenotypic transformation in AAA (Additional file 1: Fig. S1E). After hypothetical ligand-receptor pairing among each subpopulation was sorted, a large number of TIMP1-CD63/TIGB1-collagen genes fibrosis axis [16] between T-cell and SMC were revealed (Additional file 1: Fig. S1D), indicating the aggravated microenvironment within the vessel wall. An additional supplementary table (Additional file 2: Table S1) catalogs the identified cell populations and specifies their details.

T lymphocyte transcriptomes revealed phenotypic and functional heterogeneity in AAA

All T lymphocyte populations were isolated and combined (Fig. 2A). Initially, singleR identified 5 clusters of T cells, and unsupervised re-clustering further revealed 8 phenotypes (Fig. 2B). Each phenotype exhibited varying gene expression patterns, and top phenotype-specific markers were presented in Fig. 2C. The activation status propensities were further assessed, revealing 5 CD4⁺ clusters (T_1, T_2, T_4, T_5, T_7) and 3 CD8⁺ clusters (T_3, T_6, T_8) (Additional file 1: Fig. S2A). T_1, T_2, and T_7, with strong correlation (Fig. 2D), were characterized by the upregulated expression of chemokine (C–C motif) receptor 7 (CCR7), lymphoid-enhancer-binding factor 1 (LEF1), recombinant human CD62L/SELL protein (SELL), and transcription factor 7 (TCF7), indicating the stem-cell and central-memory-like phenotypes. T_1 and T_2, with the largest proportion (Additional file 1: Fig. S2B), were related to functions associated with V(D)J recombination and lymphocyte proliferation (Fig. 2E), displaying a naïve state. T_7 mediated important cellular functions such as cell recruitment and differentiation (Fig. 2E). This cellular expressed tumor necrosis factor (TNF) and Th1-related cytokine (Additional file 1: Fig. S2A), representing the Th1-like central-memory subcluster. Based on the expression of specific classical markers forkhead box P3 (FOXP3), interleukin 2 receptor subunit alpha (IL2RA), T cell immunoreceptor with Ig and ITIM domains (TIGIT), and cytotoxic T-lymphocyte associated protein 4 (CTLA4), T_4 was labeled as Treg, whose protective role has been demonstrated in animal models of AAA [17]. Effector-memory Th1-like subcluster was observed in T_5 that expressed C-X-C motif chemokine receptor 3 (CXCR3), C-X-C motif chemokine receptor 4 (CXCR4), C–C motif chemokine receptor 6 (CCR6), and TNF, interferon gamma (IFNG) (Additional file 1: Fig. S2A). T_5 had the most extensive and strong interactions among all T subclusters (Fig. 2D). Because of its high levels of granzyme K (GZMK), granzyme A (GZMA), killer cell lectin like receptor K1 (KLRK1), and CXCR4, T_3 was identified as effector-memory CD8⁺ T cells, featuring cell killing and response to INFG and TNF. Terminally differentiated cytotoxic CD8⁺ T profile was present in T_6 that expressed granzyme B (GZMB), natural killer cell granule protein 7 (NKG7), Fc gamma receptor IIIa (FCGR3A), and fibroblast growth factor binding protein 2 (FGFBP2). T_8, a small pool of cells preferentially releasing noktochor (NKT), has acquired cytotoxicity and killing ability after activation (Fig. 2E, Additional file 1: Fig. S2A).

Furthermore, T_3 and T_5 exhibited a high proliferating state, indicating that they played a critical role in AAA pathogenesis (Fig. 2F). The intracellular communication network was established to investigate the complex T-cell phenotype interaction, which revealed the strong cross-talk between T_5 and the remaining populations (Fig. 2G). Pathway-integrated analyses demonstrated that T_1 and T_2 were mainly involved in Wnt-signaling and proliferation-associated pathways. NF-ĸB, transforming growth factor beta (TGFB), and mitogen activated kinase-like protein (MAPK) signalings and oxidative-phosphorylation-hallmark were upregulated in T_5. T_3 was implicated in TNF signaling and T_6 showed selective enrichment in fluid shear stress-associated pathway (Additional file 1: Fig. S2C, D). Upregulated TNFRSF4 co-expression in T_4 was detected (Additional file 1: Fig. S2A), which has been reported to be a key negative regulator of Treg [18, 19]. Furthermore, T_4 was significantly involved in P53 signaling and apoptosis pathways (Additional file 1: Fig. S2C). Collectively, these results suggested a potential correlation between the human AAA state and the impairment of immune suppressor ability of Treg.

Trajectory reconstruction recapitulated T-cell fate decision and expression kinetics

To fully understand the dynamics of T-cell infiltration over the disease progression, we performed a pseudotime analysis based on the transcriptional changes during the cell continuous development process using Monocle2 algorithm, with the tendency of transition among different phenotypes revealed (Fig. 3A, Additional file 1: Fig. S3A). Differentiation trajectory revealed a main bifurcation event that led to 2 cell infiltration fates (Fig. 3B). Fate_2 comprised of T_3 and T_6, while fate_1 comprised of T_4 in the pathogenesis of AAA (Fig. 3C). The details of infiltration fates were shown in the supplementary document (result and discussion, Additional file 1: Fig. S3B–D). On the basis of cell development and differentiation, all T-cells can also be categorized into 8 subsets by monocle dimensionality reduction algorithm, supporting our findings that the presence of 8 major different cell phenotypes in heterogeneous T-cell populations during human AAA progression (Additional file 1: Fig. S3B). In AAA, T cells bifurcated into vastly 3 infiltration states at various time points (Additional file 1: Fig. S3C). T_1 and T_2, endowed with high plasticity, primarily dominate the root of chronological trajectory (starting state 3), which was assigned as the least mature cellular state in the pseudotime and precursors of other phenotypes. The remaining phenotypes were predominantly in terminal infiltration states (Additional file 1: Fig. S3D).

BEAM demonstrated dynamic expression of genes governing diversification of T-cell fates and construction of the Lineage trajectories. Unsupervised clustering of these cell fate genes revealed 3 clusters with different kinetics along the pseudotime (Fig. 3D). Genes in cluster_2 possessed decisive roles in fate_1 specification while cluster_1 determined Fate_2 as infiltration tendency. Assessment of the canonical markers based on BEAM revealed a shift in T cells toward immune-potentiating or immunosuppressive fate with progression (Fig. 3E, F). The proinflammatory interferon response and cytotoxic makers such as C–C motif chemokine Ligand 3 (CCL3), interferon gamma (FNG), interferon regulatory factor 1 (IRF1), perforin 1 (PRF1), GZMA, NKG7, and AHNAK nucleoprotein (AHNAK) were expressed mainly in Fate_2. Inhibitory marker TIGHT and cytotoxic T-lymphocyte associated protein 4 (CTLA4) were expressed only in Fate_1 with a reduced acceleration compared to the Fate_2’s markers. The modest expression of CCR7 and selectin L (SELL) straddled both pre-branch and Fate_1 branches, suggesting a relatively immature state and insignificant activation in T cells at this regulatory cell fate.

AAA immune microenvironment characterization

We used a large number of independent RNA-seq datasets, including all genome enrichment analysis (GSEA) results, to verify scRNA-seq-derived derived hypotheses in subsequent analyses (Figs. 4, 5, and 6, Additional file 1: Figs. S4 and S5). GSEA-calculated statistics revealed significant activation of specific pathways in AAA patients, including “antigen processing and presentation,” “autoimmune thyroid disease,” “T-cell receptor (TCR) signaling,” and “Th1/2 cell differentiation” (Fig. 4A). This result suggested that the specific Ag–driven T-cell autoimmunity is responsible for the pathogenesis. Cardiac muscle contraction, fatty acid elongation, and oxidative phosphorylation were suppressed (Fig. 4A). Disease ontology analysis also unveiled a strong link between AAA and few disease types, including vasculitis, collagen disease, and autoimmune disease (Additional file 1: Fig. S4A). Hallmark signaling computed via GSVA indicated that the immunity/inflammation-related pathways were upregulated while metabolism-related pathways were inhibited. Interferon (IFN), KRAS proto-oncogene (KRAS), TNF, NF-κB, IL6-JAK-STAT3, and IL2–STAT5 were the highly expressed cytokines and signaling molecules (Fig. 4B). Based on the PCA analysis, the immune infiltration landscape revealed a significant difference in the immune status (Additional file 1: Fig. S4B). The outcomes indicated that AAA was associated with elevated infiltration of immunocytes, especially activated T cells. T-helper cells, especially Th1, had a higher and more significant infiltration profile than other immune cells. Infiltration of activated B cell, mast cell, and eosinophil was also found in AAA samples (Fig. 4C, Additional file 1: Fig. S4C). Abdominal aortic wall tissues contained high levels of immune functional factors and enhanced immune response (Fig. 4D, Additional file 1: Fig. S4C). Activated T cells were most likely to be regulated by IFN-gamma-response, TNF-signaling-via–NF-κB, and IL6-JAK-STAT3-signaling (P < 0.001) (Fig. 4E). Genes of major categories, including MHC molecules, immunomodulators/checkpoints (CP), effector cells, and immunosuppressive cells, were assessed to identify immunophenoscore (IPS)—an effective predictor of immune response and therapy (Fig. 4F, Additional file 1: Fig. S4D).

Aneurysmal tissue showed a higher IPS than that of control samples (P < 0.01), indicating increased potency of immunogenicity in AAA (Fig. 4G). AAA showed a higher score for antigen presentation and effector cells and lower score for immunosuppressive cells than those of controls, suggesting a disorder of immunity regulation in AAA, especially immunosuppressive cell dysfunction in disease progression (Fig. 4H). This hypothesis agreed with the results from previous analysis that enhanced immunity with deregulated T regulatory function was identified in aneurysm wall microenvironment [20].

Biological implications underlying TIRS

To bring further insight into the role of genes as cell fate determinants, following the significant cellular state-related differential gene set (968 genes, P < 0.0001), BEAM cell fate leader gene set (437 genes), and differentiation-trajectory ordering gene set (375 genes) extracted out, a total of 374 T-cell infiltration state-related genes were obtained which were differentially expressed across pseudotime through the intersection of the above three gene sets. Further, eBayes-moderated ANOVA t-test was conducted on the basis of abdominal aorta wall bulk data. Eventually, the most dysregulated gene set (FDR < 0.05) in aneurysm tissues shared 50 genes overlapping with 374 T-cell infiltration state-related genes (Additional file 1: Fig. S5A), which were integrated as final T-cell infiltration regulatory gene signature (TIRS) (Additional file 1: Fig. S5B).

A coregulatory network was established to map TIRS to bulk-cell level. A high density of interactions was determined among TIRS (P < 0.001), showing the regulatory mechanisms underlying the pathogenesis and T-cell fate based on this unique gene set. The vast majority of TIRS components were co-expressed and upregulated in the aneurysmal region during the progression of infiltrated T-cell cellular differentiation. TIRS that governs the different cell fates was further clustered into 3 modules. Most components (red) were identified as the determinant module responsible for the transition of T-cell into Fate_2 (immune-potentiating). Fate_2 module genes (red) showed more statistically significant expression dysregulation than Fate_1 module genes (orange) (Fig. 5A).

To investigate the functional significance of TIRS, we performed the over-representation analysis (ORA) by applying all MSigDB sets. TIRS was involved in modules related to T-cell datasets, cellular activation and differentiation, autoimmunity, and inflammatory responses (Fig. 5B). Figure 5C illustrates the chromosomal locations and average expression of TIRS in aortic bulk tissue. As an unsupervised machine learning technique, PCA was performed to summarize TIRS into a smaller number of high-order components that outlined significant differences between AAA samples and controls (Fig. 5D). To better understand the correlation between TIRS and patient status, we characterized further the TIRS by extracting significant components to generate a continuous variable PCA-T-score (PTS) for each sample. Integrative analysis proved the significant difference in PTS between healthy controls and disease status despite the different aneurysmal size (P < 0.001), where lower PTS could be an indicator of AAA status (Fig. 5E). Bar plot based on the non-parametric tests further demonstrated the percent weight of different arterial dilation degree between low-PTS group (P < 0.00001) and high-PTS group (P < 0.05). Samples with greater disease severity were primarily distributed in low-PTS group compared to high-PTS group (large size: 44% vs 36%, small size: 42% vs 36%, health sample: 0% vs 22%, P < 0.01), suggesting a potential association between TIRS and individual status (Fig. 5F). Similar findings were also drawn in another aneurysm sample dataset (Additional file 1: Fig. S5C).

A significant connection between PTS and infiltrating degree of multiple immunocyte subpopulations, especially activated CD4⁺T and EM_CD8⁺T in abdominal aorta tissue, was identified (Additional file 1: Fig. S5D). Global TIRS expression levels between high- and low-immune infiltration tissues also showed significant differences (Additional file 1: Fig. S5E, Fig. 5G). The scatter plot further demonstrated that PTS decreased with elevated level of immune infiltration. A noticeable trend for an inverse correlation between PTS and IPS were also identified (Fig. 5H). IPS showed statistically significant differences between the high- and low-score groups (Fig. 5I), highlighting a close association between TIRS and the status of aneurysm immune microenvironment such as antigen recognition/presentation and cellular activation. T-cell activation (the process of proliferation, differentiation, and transformation of T cells stimulated by antigens) increases cell surface PD-1 expression [21]. PD-1 neutralizing antibodies and inhibitors suppressed aortic tissue inflammation and decreased vascular SMC apoptosis and vessel wall calcification in murine AAA, thereby alleviating AAA progression [22]. Here we also showed that PD-1 was expressed significantly higher in the low-PTS group than that in the high-PTS group (P < 0.0001) (Additional file 1: Fig. S5F). Taken together, TIRS may correlate with the therapeutic response to immunotherapy and serve as a promising pathological biomarker in AAA.

Detection and evaluation of key biomarkers

TIRS was subjected to the machine learning-based process for feature selection. Taking the intersection of LASSO and SVM-RFE, 4 key biomarkers were identified, including FOSB, JUNB, cystatin F (CST7), and TBC1 domain family member 4 (TBC1D4) (Fig. 6A). The expression dysregulations of these key biomarkers were verified in the abdominal aortic wall, perivascular adipose tissue from AAA patients. In Abdominal Aortic Wall Dataset 1, the expressions of FOSB, JUNB, TBC1D4, and CST7 were significantly higher in the abdominal aortic wall compared to the control aortic wall (Fig. 6B). The biomarker discrimination efficacy was evaluated using ROC curves, with AUCs for FOSB, JUNB, CST7, and TBC1D4 recorded at 0.911, 0.917, 0.926, and 0.955, respectively (Fig. 6C). The clinical impact curve (CIC) suggested that AAA patients could benefit from FOSB, JUNB, and a two-gene combination at high-risk thresholds from 0 to 1 (Fig. 6D). In Abdominal Aortic Wall Dataset 2, similar results were observed, with higher expressions of FOSB, JUNB, TBC1D4, and CST7 compared to the control aortic wall. The AUCs were 0.982 for FOSB, 0.911 for JUNB, 0.893 for CST7, and 0.857 for TBC1D4. The CIC also indicated potential benefits for AAA patients from FOSB, JUNB, and the two-gene combination at high-risk thresholds (Fig. 6E–G). In Perivascular Adipose Tissue Dataset 3, FOSB and JUNB expressions were elevated in dilated perivascular adipose tissue compared to non-dilated tissue, while CST7 and TBC1D4 showed no significant differences (Fig. 6H). AUCs for FOSB, JUNB, CST7, and TBC1D4 were 0.956, 0.887, 0.506, and 0.514, respectively (Fig. 6I). The CIC indicated that AAA patients could also benefit from FOSB, JUNB, and the two-gene combination at high-risk thresholds in this dataset (Fig. 6J).

We also confirmed the dysregulation of key biomarkers in clinical samples from AAA patients, analyzing both abdominal aortic wall and peripheral blood (Fig. 7A). In the Abdominal Aortic Wall Inhouse Dataset 1, FOSB and JUNB levels were significantly higher in the AAA aortic wall compared to the control, while CST7 and TBC1D4 showed no significant difference (Fig. 7B). The AUCs for FOSB, JUNB, CST7, and TBC1D4 were 1.000, 1.000, 0.600, and 0.800, respectively (Fig. 7C). The clinical impact curve (CIC) indicated that AAA patients could benefit from FOSB, JUNB, and the two-gene combination at high-risk thresholds from 0 to 1 (Fig. 7D). In the Peripheral Blood Inhouse Dataset 2, FOSB, JUNB, CST7, and TBC1D4 expressions were elevated in AAA patients compared to controls (Fig. 7E). The AUCs for these biomarkers were 0.989 for FOSB, 0.906 for JUNB, 0.933 for CST7, and 0.736 for TBC1D4 (Fig. 7F). The CIC also suggested potential benefits for AAA patients from FOSB, JUNB, and the two-gene combination at high-risk thresholds in this dataset (Fig. 7G). Overall, individuals with elevated FOSB and JUNB levels may possess a high risk of AAA onset or progression and require special attention and timely interventional therapy. We tested this hypothesis in a mouse model of AAA. We produced mouse AAA model using the BAPN + angII method. The gross anatomical images showed obvious production of auxiliary aortic aneurysm. We confirmed the existence of AAA using the hematoxylin and eosin (H&E) or Verhoeff’s dye liquor and Van Gieson’s dye liquor (EVG) staining. H&E staining results showed that the AAA model group mice had abdominal aorta whose three-layer structure of the tube wall is disordered, severely damaged, with inflammatory cell infiltration and endometrial hyperplasia. H&E and EVG stainings displayed elastic fibers in the abdominal aortic aneurysm wall of AAA model group mice obvious degradation and appear fracture (Fig. 7H, I). Using immunohistochemistry to measure FOSB and JUNB signals, we found that the expression of FOSB and JUNB in AAA mice increased compared with those in control mice (Fig. 7J–L). ELISA results also showed that serum FOSB and JUNB were higher in AAA than the control group (Additional file 1: Fig. S6A, B).