Single-cell transcriptome reveals age-related variations in splenic endothelial cell clusters

To elucidate the transcriptional changes associated with aging in mouse splenic endothelial cells, we isolated ECs from mouse spleen tissues from 4 young (2-month-old, female) and 4 old (2-year-old, female) C57/B16 mice. We barcoded and sequenced the ECs using 10x Genomics-based single-tube protocol and excluded pericytes (Pdgfrb) or immune cells (Ptprc), smooth muscle cells (Acta2), fibroblasts (Col1a1), and erythrocytes (Hba-a1, Hba-a2, Hbb-bs) as per previous analyses [22] and EC markers (Pecam1) have been selected. A total of 10,467 ECs (Young: 5,053, Old: 5,414) were prepared for downstream analysis after strict quality control, doublets checking (Fig. S4) and gene expression filtering (Fig. 1A). Following normalization, unsupervised graph-based clustering partitioned the cells into groups as visualized by the uniform manifold approximation and projection (UMAP) plots, including the clusters UMAP plots (Fig. 1B) and the young and old group UMAP plots (Fig. S1 A).

Fig. 1

Construction of single cell sequencing atlas of splenic endothelial cells. (A) Flow chart of scRNA-seq and bioinformatics analysis of the young and old splenic endothelial cells. (Young, n = 5053; Old, n = 5414). (B) UMAP plot showing different celltypes in mouse splenic endothelial cells. Young and old groups are sharing the same clusters: Capillary arterial, Capillary venous, Capillary1, Capillary2, Artery, Immunology1, Immunology2, Immunology3, Immunology4, Proliferating. (C) Feature plots display the expression profiles of celltype-specific marker genes for different clusters in mouse splenic ECs. (D) Heatmap showing the top 30 marker genes of specific clusters in the old group and the enrichment function annotations of each are on the right. (E) Heatmap showing the particular regulon in different clusters of splenic ECs.(F-I) Representative micrographs of old mouse splenic sections, stained for an EC marker (Cd31) and Cd14(Monocytes, F), Cd79a (Immunology3, G), Cd52 (Immunology, H), Fcerla (Mast cell, I) and counterstained with DAPI

Apart from the conventional clusters of splenic ECs identified in previous study [22], such as Capillary arterial (Gpihbp1), Capillary venous (Sema3d), Capillary1 (Glul), Capillary2 (Mal), and Artery (Stmn2) (Fig. 1C), we identified 5 additional clusters in both young and old ECs. One of these, the Proliferating EC, which highly expressed the marker Mki67, has been associated with vascular EC regeneration in previous studies [22, 23]. Notably, our cluster marker investigation revealed four distinct clusters of ECs with immunological relevance, as many of their markers are related to immune genes. We conducted an overlap analysis of the top 100 markers of these immunology-related ECs and classic immune cells (Fig. 3B). The results revealed that each subtype of immunology-related ECs displayed a strong association with specific types of immune cells: Immunology1 with dendritic cells (DC), Immunology2 with Mast cell, Immunology3 with B cells, and Immunology4 with Macrophages. They were subsequently named Immunology1 (H2-Eb1), Immunology2 (Fcerla), Immunology3 (Cd79a), and Immunology4 (Lyz2) (Fig. 1C).

We identified a set of Top50 markers for each cell type (Table S1). The cell proportions of different clusters in the young and old splenic groups (Fig. S1 A) indicated that the number of Capillary Venous ECs increase with aging, which shares the same situation in the previous study [23]. Other clusters such as Capillary1, Capillary2 and Capillary Arterial are declined with aging.

Gene Ontology (GO) analysis of the Top 30 marker genes unveiled the functional characteristics linked to each specific cell type (Fig. 1D). The function “branching involved in blood vessel morphogenesis” was enriched for Capillary1, indicating its involvement in the development of vessel branches. For Capiilary2, the function “regulation of transcription from RNA polymerase II promoter in response to stress” was enriched. In addition, the Capillary Arterial, Capillary venous and Artery clusters were highly associate with vascular development. The Proliferating ECs showed a strong ability to regulate the mitotic cell cycle, confirming its role in vascular proliferation.

We further investigated the molecular mechanisms underlying EC phenotypic differentiation using single-cell regulatory network inference and clustering (SCENIC) [19]. We analyzed the regulons exerting the most significant impact on each type of EC (Fig. 1E), elucidating the regulation strength of each regulon for different cell types. Double immunostaining of IF and FISH for an EC marker (CD31) and the markers of these specialized immunology EC phenotypes validated the scRNA-seq data (Fig. 1F–I).

Age-related cellular and molecular characteristics of murine splenic ECs

To further elucidate the impact of aging on the molecular function of the splenic endothelium, we analyzed the overall differentially expressed genes (DEGs) in young and old splenic ECs (Fig. 2A and Table S2). We discovered that functions such as “response to interferon-gamma,” “cytokine-mediated signaling pathway,” and “response to tumor necrosis factor” were enriched for upregulated DEGs, indicating an increased expression of inflammatory signatures and a shift towards a pro-inflammatory phenotype (Fig. 2B). These findings align with previous research [24, 25]. Functions such as “endothelium development,” “response to growth factor,” and “blood circulation” were enriched in downregulated DEGs, suggesting vascular dysfunction in senescent endothelial cells [24].

Fig. 2

Cellular and molecular characteristics of aged splenic ECs. (A) Volcano plot showing aging-associated up- and down- regulated differentially expressed genes (DEGs) in all celltypes (adjusted P-value < 0.05,|LogFC| >0.25), high expression of up-regulated gene are labeled. (B) Barplot showing GO terms of the overall DEGs between young and old group in all celltypes. (C) Circos plot showing aging-associated up- and down- regulated DEGs, each of the connecting curve showing the gene is whether up- or down- regulated between two clusters. (D) Dot plots showing the top five celltype-specific DEGs of different clusters. Only those with annotations are showed. Up-regulated genes are colored in red while the down-regulated ones are in blue

To further investigate cell type-specific alterations in gene expression, we identified key cell types and molecular mechanisms affected by splenic senescence (Table S2). We observed the highest numbers of both up- and down-regulated DEGs in Immunology1 (up:198, down:109), Immunology4 (up:128, down:174), and Proliferating (up:169, down:187) clusters (Fig. S1 C), suggesting these cell types may be most affected by aging. The chord plot (Fig. 2C) indicated more overlapping DEGs in classic clusters (such as Capillary Arterial, Capillary Venous, Capillary1, Capillary2, and Artery) than in unique clusters like Immunology1, Immunology2, Immunology3, Immunology4, and Proliferating. The top 5 DEGs in different clusters revealed an increased immune signature, as illustrated by the overexpression of H2-Aa, H2-Ab1, H2-Eb1, Cd74 in clusters such as Capillary Arterial, Capillary Venous, Capillary1, and Immunology1 (Fig. 2D, Fig. S1 C). This is consistent with existing literature [26,27,28,29]. Cxcl12 was found to be downregulated in almost every cluster of splenic ECs (Fig. S1 D), which could affect the ability of immune cell homing for the splenic ECs [10, 30].

Cellular and molecular characteristics of immunology-related ECs

We have delineated four distinct subpopulations of ECs with immunological relevance, each demonstrating a strong association with a particular immune cell type (Fig. 3B). Immunology1 cells express dendritic cell markers such as Fabp4 and H2-Eb1, suggesting functional similarities to DC, while Immunology2 cells express Fcerla alongside other mast cell markers, indicating mast cell-like characteristics. Immunology3 cells are characterized by the expression of B cell markers including Cd19, Cd79a, and Cd79b, and Immunology4 cells are marked by macrophage-associated genes such as Lyz2 and Cd68 (Fig. 3B).

To further investigate their relationship with corresponding immune cells and the age-related changes in their cellular and molecular characteristics, we integrated splenic immune cell data from the Tabula Muris datasets [28, 31] (B cells, Mast cells, Macrophages) and splenic DCs from a prior study [32]. The feature plot (Fig. 3A) illustrates that the positions of these immunology-related ECs are distinct from those of classical ECs, suggesting functional divergence. We then assessed the correlation between classical ECs, immunology-related ECs, and immune cells using a defined methodology. Classical ECs (Capillary arterial, Capillary venous, Capillary1, Capillary2, Artery, and Proliferating) showed strong inter-correlations with coefficients exceeding 0.8 (Fig. S1 E). In the young cohort, the correlation between classical ECs and immunology-related ECs ranged from 0.5 to 0.8, lower than within the classical group, whereas in the older cohort, the correlation coefficients were higher. The correlation analysis between immunology-related ECs and immune cells (Fig. 3C) revealed that Immunology1 and DC had the strongest association, particularly in the older group. In the young cohort, B cell and Immunology3, as well as Mast cells with Immunology2 and Macrophages with Immunology4, shared stronger relationships.

Fig. 3

Cellular and molecular characteristics of four immunology-related ECs. (A) Left: Feature plots showing the positions of four subtypes of immunology-related ECs. Right: Bar chart illustrating the most enriched Gene Ontology (GO) term within these subtypes. The X-axis represents the number of genes associated with the functions, and the bars are color-coded based on Log10P values. (B) Heatmap showing the shared markers between immunology-related ECs and Immune cells (DCs, Mast cells, B cells, Macrophages and monocytes). (C) Correlation matrix displaying the relationships between immunology-related ECs and Immune cells. The upper matrix, shaded in grey, illustrates the correlations between Old immunology-related ECs and Immune cells, while the lower matrix depicts the correlations for young ones. Each correlation coefficient is labeled within the corresponding circle in the matrix. (D) Upper: Boxplot showing the Immunology-related gene (IRG) score of Young and Old regular ECs subtypes and immunology ECs. ****p < 0.001. Lower: Boxplot showing the immune cell homing (ICH) genes score of Young and Old regular ECs subtypes and immunology ECs. ****p < 0.001. (E) Boxplot showing the Antigen Presenting gene score of young and old immunology-related ECs. ****p < 0.001. (F) Boxplot showing the Allergy-related gene score of young and old immunology-related ECs. ****p < 0.001. (G) Boxplot showing the B cell Activation gene score of young and old immunology-related ECs. ****p < 0.001. (H) Boxplot showing the Phagocytosis gene score of young and old immunology-related ECs. ****p < 0.001

Then we further analyzed the differential expression genes between immunology-related ECs and their relative immune cells (Fig. S5 A-D). We calculated the DEGs between them and the top5 DEGs are shown. The presence of endothelial cell markers Gm42418 among all the upregulated DEGs in four immunology-related EC clusters indicates their characteristics. Besides, immune cell related genes also presented in the up-regulated DEGs in each immunology-related EC clusters such as Ly6a in Immunology1(Fig. S5 A), Mcpt8 in Immunology2 (Fig. S5 B), Ighm and Igkc in Immunology3(Fig. S5 C) and Adgre1 in Immunology4(Fig. S5 D), illustrating that immunology-related ECs possesses attributes of both endothelial cells and immune cells, representing an intermediate state of endothelial cells while retaining the conservatism of endothelial cells.

To gain a deeper understanding of immunology-related ECs, we compared them with classical ECs using the Immunology-Related Gene Set Score [14] (Fig. 3D, Methods). It is clear that the scores of immunology-related ECs exceed those of classical ECs in both young and old groups, reinforcing our hypothesis that immunology-related ECs possess enhanced immune regulatory functions. When comparing within classical and immunology-related ECs, the scores for the young group are higher than those for the old group, aligning with previous studies that suggest a decline in the immune regulatory functions of ECs with age [23].

Previous research has highlighted the role of ECs in the recruitment and homing of immune cells [33, 34].By interacting with circulating innate and adaptive immune cells and regulating their extravasation from the bloodstream into the tissue parenchyma, ECs could play a crucial role in controlling tissue and lymph node inflammation. Therefore, we calculated the Immune Cell Homing (ICH) score [15] in classical ECs and immunology-related ECs (Fig. 3D). In both young and old groups, classical ECs exhibit higher scores than immunology-related ECs, suggesting that immunology-related ECs have transitioned away from their endothelial characteristics and have adopted properties more akin to immune cells.

Given that ECs are considered semi-professional APCs [10], we evaluated the antigen presentation score between classical ECs and immunology-related ECs to discern any differences (Fig. S1 F). The higher scores of both young and old immunology-related ECs compared to classical groups support the notion that immunology-related ECs may possess stronger immunomodulatory functions than classical ECs.

The findings discussed above led us to consider whether the four identified immunology-related ECs might align with the concept of immune cell-like ECs (EndICLT) [11]. We proceeded to evaluate the unique immunoregulatory functions of these four immunology-related ECs in relation to their associated immune cells. For instance, DCs are renowned for their role as efficient APCs [35], prompting us to score this function among the four immunology-related ECs. The highest antigen presentation score was observed in Immunology1 ECs in both young and old groups (Fig. 3E), suggesting a strong functional resemblance to DCs. This similarity extends beyond shared markers to encompass functional attributes, with aging appearing to enhance antigen presentation capabilities.

For Immunology2, we calculated the allergy gene score, given the central role of Mast cells in initiating allergic immune responses (Fig. 3F) [36]. Immunology1, Immunology2, and Immunology4 scored similarly, reflecting the complex interplay of various immune cells and molecules involved in allergic responses. In terms of B cell activation, Immunology3 demonstrated the strongest performance, underscoring its close association with B cells, and this function was preserved with aging. For Immunology4, we evaluated the phagocytosis score, and as anticipated, it scored the highest among the immunology-related EC groups. The older group achieved lower scores compared to the younger group, aligning with previous studies that suggest the phagocytic function of Macrophages is impaired with age [37].

In order to explore if the age-induced changes in immunology-related ECs are consistent with the immune cell in peripheral blood. We used the peripheral blood scRNA-seq data in the young and aged mouse from article “Single-cell transcriptomics of peripheral blood in the aging mouse”(GSE120505) [38]. After quality control, we reclustered and identified 10 clusters of cells in human peripheral blood cells (Fig. S6 A). The subcluster related to immunology-related EC (DC, Basophil, B cell and Macrophage) was extracted and aged-related DEGs of each cluster were analyzed (Fig. S6 B-E).

In peripheral blood, age-related changes in DCs are linked to blood coagulation, while in DC-resembling ECs (Immunology1), they are associated with antigen processing and extracellular matrix organization. Three intersecting genes were found among the upregulated DEGs in peripheral blood DCs and Immunology1(Fig. S6 B). One of them is Ecm1, upregulated in the aging heart, contributes to cardiac fibroblast stimulation and fibrosis in aging and myocardial infarction [39]. In peripheral basophils, upregulated DEGs are linked to leukocyte chemotaxis, while in mast cell-resembling ECs (Immunology2), they are related to vasculature development. There were no overlapping genes between the age-upregulated DEGs in basophils and Immunology2 (Fig. S6 C), indicating heterogeneous changes during aging despite their shared origin. For peripheral macrophages, the function of regulating response to external stimuli was upregulated. In macrophage-resembling ECs (Immunology4), age-related DEGs were associated with blood vessel endothelial cell migration. Psen1, relevant to autosomal-dominant Alzheimer’s disease, was among the overlapping genes (Fig. S6 E) [40].

These results provide evidence that peripheral blood immune cells and immunology-related endothelial cells share some similar age-induced changes in both functional and gene phenotypes such as immune response and chemotaxis. However, the aging-related changes in clusters of immunology-related endothelial cells are more focused on their regulatory function in vessel development.

Overall, Immunology-related ECs are not only express markers similar to immune cells but also retain their immunomodulatory functions. This further substantiates the notion that they are undergoing EndICLT.

Age-related molecular changes along differentiation trajectories

Previous research suggests that cells undergoing EndICLT do not fully differentiate into immune cells but rather occupy an intermediate state between ECs and immune cells [11]. To investigate whether the four immunology-related ECs fall along differentiation trajectories, we conducted pseudotime analysis using Monocle 2 [41, 42]. We first analyzed all splenic ECs and immune cells in both young and old groups (Fig. 4A). The analysis delineated three major differentiated cell groups: (1) ECs; (2) B cells; (3) Macrophages and DCs. It appeared that the four immunology-related EC groups were positioned along these trajectory paths.

Fig. 4

The pseudotime trajectory reveals relationship between splenic immunology-related ECs and immune cells. (A) Pseudotime analysis conducted on all ten ECs cell types (Capillary arterial, Capillaryvenous, Capillary1, Capillary2, Artery, Immunology1, Immunology2, Immunology3, Immunology4, Prolifearting) in splenic ECs (young and old integrated) and 4 groups of Immune cells (DC, Mast cells, B cells and Macrophages). The points are colored by celltypes. (B) Pseudotime trajectory analysis of Immunology4 EC related subtypes in young group: Capillary1, Capillary2, Proliferating, Immunology4 and Macrophages. (C) Pseudotime trajectory analysis of Immunology4 EC related subtypes in old group: Capillary1, Capillary2, Proliferating, Immunology4 and Macrophages. (D) Heatmap showing the expression profiles along the pseudotime of Time differential genes (q value < 1e-4) in the trajectory of Immunolog4 EC related subtypes in old group including Capillary1, Capillary2, Proliferating, Immunology4, and Macrophage, which were divided into three clusters with the expression pattern and the gene in cluster 3 represented on the right. (E) The cluster dendrogram of co-expression in the trajectory of Immunology4 EC related subtypes in old group. (F) Dotplot showing the expression of the co-expression modules in Immunology4 EC related subtypes in old group. (G) The kME plot of the co-expression modules in Immunology 4 EC related subtypes. Top 10 hub genes of each co-expression module of Immunology4 EC related subtypes in old group are visualized on the right. (H) The pseudotime trajectory for co-expression module. The Y-axis represents the kME score of each hub genes. (I) Left, pie plot showing overlapped genes between cluster 3 and hub genes of co-expression module 2 of the Immunology 4 EC related subtypes. Right, bar plot showing enriched GO terms and KEGG of the overlapped genes listed in the middle

To gain a clearer understanding of the status of these four groups, we divided the analysis into four subsets (Fig. S2 A, B). The first subset included Capillary1, Capillary2, Proliferating, Immunology1, and DCs, with Proliferating ECs as the starting point due to their regenerative potential (Fig. S2 A-i). Some Immunology1 cells were located on the branch between ECs and DCs, in State 3, suggesting a potential transitional status (Fig. S2 B-I). The second subset comprised Capillary1, Capillary2, Proliferating, Immunology2, and Mast cells. In the older group, Proliferating and Immunology2 cells branched towards the Mast cells trajectory, suggesting that aging may disrupt EC differentiation (Fig. S2 A-ii, Fig. S2 B-ii). Immunology2 cells also occupied a middle status between ECs and Mast cells.

The third subset included Capillary1, Capillary2, Proliferating, Immunology3, and B cells (Fig. S2 A-iii, Fig. S2 B-iii), with all Immunology3 cells positioned in an intermediate state between ECs and immune cells. The fourth subset contained Capillary1, Capillary2, Proliferating, Immunology4, and Macrophages. Similar to the third subset, all Immunology4 ECs were situated in an intermediate state (Fig. S2 A-iv, Fig. S2 B-iv). Notably, the proximity of Immunology4 ECs to the branch point of the trajectory prompted us to consider potential molecular changes in Immunology4 that could facilitate the transition from classical ECs to Macrophages.

To further explore the age-related changes, we separated the Immunology4-related ECs into young and old groups and conducted individual trajectory analyses. The trajectory analysis revealed distinct patterns between the young and old groups (Fig. 4B, C). We used the cellAlign algorithm [21] to analyze the distinctions between the trajectories of old and young. First, we analyzed the global alignment to quantify overall similarity in expression throughout the trajectory (Fig. S7 A, B). Then we used the young dataset as the reference to map the difference between the two datasets(Fig. S7 C).The results indicated that during pseudotime, most of the genes in young and old groups are conserved while some of the genes are unique. Secondly, we analyzed the local alignment to identifies regions of the trajectory that are closer to each other (Fig. S7 D). The regions are remained conservative within the 0-100 range while in the later stages of the trajectory, they begun to differentiate. This clearly demonstrates the differences in trajectories between young and old group, highlighting specific gene changes in the differentiation trajectory associated with aging.

Within the older group, Immunology4 ECs were found to occupy an intermediate state between ECs and immune cells, a pattern not observed in the young group (Fig. 4C). This observation suggests that aging may push Immunology4 into a transitional state.

To substantiate our hypothesis, we initially examined the gene expression profiles at the branch point in the aged group of Immunology4 ECs related subsets (Fig. 4D). The genes were classified into three groups according to their expression dynamics. It was clear that the majority of genes in Cluster 1 displayed high expression levels in Cell Fate 2, while those in Cluster 2 were predominantly expressed in Cell Fate 2 as well. In contrast, Cluster 3 exhibited high expression levels before the branch point.

Subsequently, we identified co-expression modules using the R package hdWCGNA [17] to verify whether gene expression alters or remains consistent throughout the transition process from ECs to immune cells in the aged group. We constructed a co-expression network and visualized it using a dendrogram (Fig. 4E), and three co-expression modules were calculated. Genes within Module 1 showed elevated expression levels in Proliferating ECs and Capillary1 and Capillary2 ECs, while those in Module 2 and 3 displayed heightened expression in Macrophages (Fig. 4F).

We then computed the Module Eigengenes in each module, a metric frequently used to summarize the gene expression profile of an entire co-expression module. The module eigengenes were computed by performing principal component analysis (PCA) on the subset of the gene expression matrix comprising each module [17]. The top 10 hub genes of Module 1 (Fig. 4G-upper) are primarily involved in the regulation of endothelial cell proliferation and regulation of vascular permeability. For Module 2, the hub genes mainly focus on the regulation of lymphocyte proliferation and positive regulation of leukocyte-mediated cytotoxicity (Fig. 4F-middle). Module 3 hub genes are related to the positive regulation of the mitotic cell cycle phase transition. We analyzed how the module eigengenes change throughout the pseudotime trajectories for each co-expression module using the hdWGCNA function PlotModuleTrajectory.

We conducted pseudotime trajectory analysis with Monocle2 and studied module dynamics throughout the cellular transitions from Proliferating ECs to Macrophages. The trajectories indicated that Module 1 was turning off their expression programs throughout the transition from ECs to Macrophages, while Module 2 was turning on in the process. These results paralleled the observation that its endothelial function declined during the transition while its immune phenotype was upregulated. The Module 3 genes remained static during the transition (Fig. 4H). This finding suggested that the genes in Module 2 may contribute to the regulation of the transition.

To identify potential transitional genes, we overlapped genes in Cluster 3 and Module 4’s top 50 hub genes, and 18 genes were detected (Fig. 4I). The GO enrichment showed the function of overlapped genes mainly about regulation of inflammatory response (Gpx1, Alox5ap, Ctss), extracellular matrix binding (Ctss, Lgals3), antigen processing and presentation (Ifi30, Ctss), and phagocytic and endocytic vesicle (Ctss, Pld4, Ftl1). Previous research has shown that Ctss is involved in inducing the release of inflammatory cytokines and leading to endothelial dysfunction in hyperglycemic conditions, and another study demonstrated that decreasing Ctss can ameliorate age-related dry eye. These findings suggested that aging may activate Ctss and initiate the EndICLT process (Fig. S2 C).

Besides, we also analyzed the transcription factors in the old Immunology4 subsets (including Capillary1, Capillary2, Proliferating, Immunology4 and Macrophage) used pySCENIC [19]. After constructing the co-expression modules, we calculated the cell-type specific regulators(Fig. S8 A, B).Then to explore the TFs changes in trajectory, we integrate transcription factor activity into the pseudotime matrix(the pseudotime was normalized). According to our pervious trajectory (Fig. 4C), the transition occurred in the time interval between 0.35 and 0.45 (from 0 to 1). We specifically investigated immunology4-specific transcription factors and plotted the transcription factor activity over pseudotime to see whether any specific changes near transition time(Fig. S8 C).We found that Ltf showed a specific upregulation around 0.4(Fig. S8 D). It is well-known that LTF could activate the NF-κB signaling pathway, promoting macrophage activation [43]. For example, LTF binding to CD14 receptor competes with the bacterial LPS (product of dying bacteria) [44] and can attenuate NF-κB-induced transcription of genes for various inflammatory mediators [45]. Also, it is reported that high LTF expression might contribute to meniscal aging and degeneration through the NF-κB signaling pathway [46].Therefore, we could speculate that the endothelial-immunology transition is regulated by Ltf.

Profiling the impact of aging on spleen ECs cellular communication and molecular signaling pathways

To explore cell-to-cell interactions between the young and old groups, we utilized the R package CellChat for analysis [20]. In line with observations from other studies [23, 47], the aged group exhibited stronger signal communication in terms of both quantity and intensity (Fig. 5A, B). Compared to the young group, inflammation pathways such as Ccl, Cxcl, and Mif were upregulated (Fig. 5A-lower, C). Pathways in the old Immunology2 ECs were the most upregulated (Fig.